Equation Solver Module

I think a set of equation solvers from linear to quartic should be implemented and optimized for getting exact solutions as well as approximations. I have made a simple example for a `cubicSolve` function

``` Swift

/// # Cubic Equation Solver
///
/// Solves the cubic equation of the form `ax^3 + bx^2 + cx + d = 0` - **eq: 1**
/// Should follow some variation of the following algorithm
/// - important: a, b, c, d should all be real numbers.
/// - reference: [Cubic Equation](https://en.wikipedia.org/wiki/Cubic_equation)
///
/// Procedure:
/// 1. Divide both sides of equation by a.
/// i. In the case that a == 0 fallback onto the `quadraticSolve`.
/// ii. other wise move on to step 2.
/// The equation should now look as follows x^3 + (b/a)x^2 + (c/a)x + d/a = 0.
/// 2. Calculate the discriminant D
/// i. Let r_1 = b/a , r_2 = c/a , and r_3 = d/a.
/// ii. Let `G = (3*r_2 - r_1^2)/9` and `F = (9*r_1*r_2 - 27*r_3 - 2*r_1^3)/54`
/// iii. Then `D = G^3 + F^2`
/// a. If `D > 0` , one real root with 2 complex conjugate roots.
/// b. if `D = 0` , all roots are real and atleast 2 are repeated.
/// c. If `D < 0`, all roots are real and unequal.
/// A. `D > 0`, where `cbrt` and `sqrt` are the cuberoot and squareroot respectively.
/// 1. The only real solution: ` x_1 = cbrt(F + sqrt(D)) + cbrt(F - sqrt(D)) - (1/3)*r_1`
/// 2. First Complex solution : `x_2 = Complex(-(1/2)(cbrt(F + sqrt(D)) + cbrt(F - sqrt(D))) - (1/3)*r_1 , (cbrt(3)/2)*cbrt(F + sqrt(D)) - cbrt(F - sqrt(D)))`
/// 3. Second Complex solution: `x_3 = Complex(-(1/2)(cbrt(F + sqrt(D)) + cbrt(F - sqrt(D))) - (1/3)*r_1 , -(cbrt(3)/2)*cbrt(F + sqrt(D)) - cbrt(F - sqrt(D)))`
/// B. For `D = 0` and `D < 0 ` Reduce Equation to Depressed Cubic
/// 1. Make the substitution of `x = t - b/(3a)` in **eq: 1** resulting in:
/// `t^3 + p*t + q = 0` - **eq: 2** where:
/// `p = (3*a*c - b^2)/(3*a^2)` and `q = (2b^3 - 9*a*b*c + 27*a^2*d)/(27*a^3)`
/// i. If `p = q = 0` then all roots `t_1 = t_2 = t_3 = 0`
/// ii. If `4*p^3 + 27*q^2 = 0` and `p != 0` then `t_1 = 3*q/p` and ` t_2 = t_3 = -3*q/(2*p)`
/// iii. If `D < 0` then `t(k) = 2*sqrt(-p/3)*cos((1/3)*acos((3*q/(2*p)*sqrt(-3/p))) - 2*pi*k/3)` for `k = 1,2,3`
/// 2. Convert back to x using `x(k) = t(k) - b/(3*a)`
///
/// - note: General Optimizations can be performed for trigonometric identities or for making approximations based upon values of a ,b , c, and d
///
///
///
///
func cubicSolve(a: R, b: R, c: R, d: R) -> [C] {
if a == 0 { return quadraticSolve(a: b, b: c, c: d) }

let r_1 = b/a
let r_2 = c/a
let r_3 = d/a

let g = (3.0*r_2 - r_1^2.0)/9.0 // G
let f = (9*r_1*r_2 - 27*r_3 - 2*r_1^3)/54 // F

let d = g^3 + f^2 // discriminant
if d > 0 {
let x_1 = cbrt(f + sqrt(d)) + cbrt(f - sqrt(d)) - (1/3)*r_1
let x_2 = Complex(-(1/2)(cbrt(f + sqrt(d)) + cbrt(f - sqrt(d))) - (1/3)*r_1 , (cbrt(3)/2)*cbrt(f + sqrt(d)) - cbrt(f - sqrt(d)))
let x_3 = Complex(-(1/2)(cbrt(f + sqrt(d)) + cbrt(f - sqrt(d))) - (1/3)*r_1 , -(cbrt(3)/2)*cbrt(f + sqrt(d)) - cbrt(f - sqrt(d)))
return [x_1, x_2, x_3]

} else if d = 0 {
let p = (3*a*c - b^2)/(3*a^2)
let q = (2b^3 - 9*a*b*c + 27*a^2*d)/(27*a^3)
let conversionFactor = -(-b/(3*a))

if p == 0 && q == 0 {
return [0, 0, 0]
} else if p != 0 && 4*p^3 + 27*q^2 == 0 {

let t_1 = 3*q/p
let t_2 = -3*q/(2*p)
let t_3 = -3*q/(2*p)
return [t_1 + conversionFactor, t_2 + conversionFactor, t_3 + conversionFactor]
}

} else if d < 0 {
return (1...3).map {
2*sqrt(-p/3)*cos((1/3)*acos((3*q/(2*p)*sqrt(-3/p))) - 2*pi*$0/3) + conversionFactor
}
}
}
```

Start with the proposed `cubicSolve` signature and algorithm in issue #42. Clarify the intended linear-to-quartic scope, exact versus approximate behavior, numeric type requirements, and expected handling of degenerate cases before implementation. Done should mean that the agreed solver scope and public API are implemented and verified.

Written by the indexing model from the issue text.