Since Swift Numerics will be heavily based on numerical computing and math, I want to make this proposal for an add-on module on numerical analysis. Example use cases include (but are not limited to) the following:

• Finding the roots (or zeroes) of continuous functions
• Interpolating polynomials
• Approximating the definite integrals of continuous functions
• Solving linear systems
• Solving initial value problems

NOTE: Complex numbers are required for implementing Muller’s method.

P.S.: As someone with a background in applied mathematics, thank you @stephentyrone for introducing this wonderful Swift package!

Take a look at e.g. the Tensorflow Github Repo's "topics" :)

For the reasons discussed in this forums thread it would be slightly advantageous to make the Complex type and module have different names. This requires renaming one of them.

This would allow people who want to work with a single concrete complex type to do something like:

typealias Complex = ComplexType<Float> // if we rename the type
typealias Complex = ComplexModule.Complex<Float> // If we rename the module

They would then still be able to get at the generic type as either ComplexType or ComplexModule.Complex.

As you can see, renaming the type makes this specific use case a tiny bit neater, but makes the normal generic usage a tiny bit messier (while renaming the module only changes the import statement in routine use, which is less invasive).

Should we include some linting tool to enforce Swift syntax? E.g. SwiftLint?

In order to write generic mathematical functionality against the ElementaryFunctions protocol, it is necessary to be able to use basic arithmetic. This includes the infix operators (+, -, *, /) and the prefix operators (+, -).

Therefore, the ElementaryFunctions protocol should either refine SignedNumeric, or if it is intended that non-numeric types be able to conform, the arithmetic operators should be included in the protocol.

print(Foundation.cbrt(1000.0)) // 10.0
print(Double.pow(1000.0, 1.0/3.0)) // 9.999999999999998

Calculation through "pow" gives a different result.

Currently, Complex.init?(length: 0, phase: .infinity) produces nil.

However, when the length of a Complex number is identically 0, the phase is irrelevant and the value is always exactly (0, 0). This is even true for non-finite phase, because the length is still 0.

Thus, the polar initializer for Complex should produce (0, 0) when the input length is 0, regardless of the phase.

The current implementation of Complex has both real and imaginary as get-only. These properties should have setters.

pull-request: #44

Since Swift supports Unicode in identifiers, it is possible to improve the appearance of the code by adding aliases for operators, such as: × ⋅ ÷ ∶ (RATIO) - (MINUS SIGN), etc.

What do you think of it?

Add a flexibly-sized bigint type that grows as required, conforming to BinaryInteger.

I made this library called NDArray a while ago and got it to a POC state. I know ShapedArray is on the radar so I wanted to know if I could contribute the code to a longer term project.

Implement IEEE 754 Decimal128 (and possibly Decimal64 as well, though it's less important). This has a dependency on fixed-width integer types because we want to be able to leverage UInt128 arithmetic for the implementation.

Currently, dividing Complex(0, 0) by itself returns Complex(0, 0), but it should return .infinity to represent the undefined result.

The encoding would change from {"x":0,"y":0} to [0,0] in JSON.

SE-0239 was originally pitched with CodingKeys, but in a discussion of the implementation, Tony Parker wrote:

This is one of those cases where I think we should just use two unkeyed values. ClosedRange is well defined and I don't think it's likely for it to grow the need for more keys. Same goes for the other Range types.

This would probably also make the encoded range types more parsable by other languages when encoded to JSON.

... I feel very strongly that we should prefer a compact representation for fundamental and mathematically-defined types like this.

Instead of the synthesized conformance, we could use the following:

extension Complex: Decodable where RealType: Decodable {
  public init(from decoder: Decoder) throws {
    var unkeyedContainer = try decoder.unkeyedContainer()
    let x = try unkeyedContainer.decode(RealType.self)
    let y = try unkeyedContainer.decode(RealType.self)
    self.init(x, y)
  }
}

extension Complex: Encodable where RealType: Encodable {
  public func encode(to encoder: Encoder) throws {
    var unkeyedContainer = encoder.unkeyedContainer()
    try unkeyedContainer.encode(x)
    try unkeyedContainer.encode(y)
  }
}

Rational Number Data Type

This is an example proposal for a rational number module. A good example use of a rational number
would for when a user wishes to simplify symbolic expressions without needing to worry about
floating point errors. I believe this type would be most useful as an internal currency for symbolic
calculations, where performance is sacrificed for accuracy. This type would also allow for the
division of integers without explicit calculations needing to be performed. Some issues will arise
such as the reduction of fractions to simplest form. Creating a performant algorithm for finding
common factors will be a must. Another important issue will be the conversion of both terminal and
repeating decimal numbers into simplified fractions.

The Rational type would conform to the Real protocol and could be used for both
the real and imaginary components of a Complex value.

Simple Example

struct Rational {
    
    let numerator: Int
    let denominator: Int
    
    // todo - fix this so that infinities are handled properly
    var decimal: Double {
        denominator != 0 ? Double(numerator)/Double(denominator) : .infinity
    }
    
    public init(_ numerator: Int, _ denominator: Int) {
        self.numerator = numerator
        self.denominator = denominator
    }
    
    /// Need to test and refactor this to be more robust.
    var isNegative: Bool {
        if numerator >= 0 && denominator > 0 {
            return false
        } else if numerator <= 0 && denominator < 0 {
            return false
        } else {
            return true
        }
    }
}


// MARK: - Equatable Conformance
extension Rational: Equatable {
    // todo: come back and adjust this because using the decimal is a bit of a naive
    // approach to checking equivalence.
    static func == (lhs: Rational, rhs: Rational) -> Bool {
        return lhs.decimal == rhs.decimal
    }
}


// MARK: - Comparable Conformance
extension Rational: Comparable {
    // todo: Come back and fix for same reason as Equatable.
    static func < (lhs: Rational, rhs: Rational) -> Bool {
        return lhs.decimal < rhs.decimal
    }
}


// MARK: - Int and Rational Math

/// Compute the value of a rational number to the power of an integer.
/// The three main cases are handled
/// * power greater than 0
/// * power less than 0
/// * power equal to 0
/// - important- Relies on a for loop to multiply the numerator and denominator by its self (power-1) times. This could be very slow for large powers.
func pow(base: Rational, power: Int) -> Rational {
    if power == 0 {
        return 1
    } else {
        var num = base.numerator
        var denom = base.denominator
        for _ in 1...abs(power)-1 {
            num *= base.numerator
            denom *= base.denominator
        }
        
        return power > 0 ? RationalNumber(num, denom) : RationalNumber(denom, num)
    }
}

// MARK: AdditiveArithmetic Conformance
extension Rational: AdditiveArithmetic {
    static var zero: Rational = 0
    
    static func += (lhs: inout Rational, rhs: Rational) {
        lhs = lhs+rhs
    }
    
    static func -= (lhs: inout Rational, rhs: Rational) {
        lhs = lhs-rhs
    }
    
}

// MARK: Numeric Conformance 
extension Rational: Numeric {
    var magnitude: Rational {
        return RationalNumber(abs(numerator), abs(denominator))
    }
    
    init?<T>(exactly source: T) where T : BinaryInteger {
        self.numerator = Int(source)
        self.denominator = 1
    }
    
    static func *= (lhs: inout Rational, rhs: Rational) {
        lhs = lhs*rhs
    }
}




// MARK: - (Rational, Rational) Math
extension Rational {
    static func *(lhs: Rational, rhs: Rational) -> Rational {
        return Rational(lhs.numerator*rhs.numerator,
                              lhs.denominator*rhs.denominator)
    }
    
    static func /(lhs: Rational, rhs: Rational) -> Rational {
        return Rational(lhs.numerator*rhs.denominator,
                              lhs.denominator*rhs.numerator)
    }
    
    static func +(lhs: Rational, rhs: Rational) -> Rational {
        return Rational(lhs.numerator*rhs.denominator + lhs.denominator*rhs.numerator,
                              lhs.denominator*rhs.denominator)
    }
    
    static func -(lhs: Rational, rhs: Rational) -> Rational {
        return Rational(lhs.numerator*rhs.denominator - lhs.denominator*rhs.numerator,
                              lhs.denominator*rhs.denominator)
    }
    
    static func ^(lhs: Rational, rhs: Int) -> Rational {
        return pow(base: lhs, power: rhs)
    }
    
    static prefix func -(input: Rational) -> Rational {
        return Rational(-input.numerator, input.denominator)
    }
}


// MARK: - ExpressibleByIntegerLiteral Conformance
extension Rational: ExpressibleByIntegerLiteral {
    init(integerLiteral value: Int) {
        numerator = value
        denominator = 1
    }
}

// MARK: - Sequence Extension
extension Sequence where Element: Numeric {
    func product() -> Element {
        return reduce(1, *)
    }
}

Tests

final class RationalTests: XCTestCase {
    
    // MARK: - Math operations
    
    func testScalarMultiplication() {
        let scalar: Rational = 3
        let rational: Rational = Rational(2, 3)
        
        let standard: Rational = Rational(6, 3)
        
        XCTAssertEqual((scalar*rational).numerator, standard.numerator)
        XCTAssertEqual((scalar*rational).denominator, standard.denominator)
    }
    
    func testMultiplication() {
        let left = Rational(5, 2)
        let right = Rational(3, 4)
        
        let standard = Rational(15, 8)
        
        XCTAssertEqual((left*right).numerator, standard.numerator)
        XCTAssertEqual((left*right).denominator, standard.denominator)
    }
    
    func testDivision() {
        let left = Rational(5, 2)
        let right = Rational(3, 4)
        
        let standard = Rational(20, 6)
        
        XCTAssertEqual((left/right).numerator, standard.numerator)
        XCTAssertEqual((left/right).denominator, standard.denominator)
    }
    
    func testAddition() {
        let left = Rational(5, 2)
        let right = Rational(3, 4)
        
        let standard = Rational(26, 8)
        
        XCTAssertEqual((left+right).numerator, standard.numerator)
        XCTAssertEqual((left+right).denominator, standard.denominator)
    }
    
    func testSubtraction() {
        let left = Rational(5, 2)
        let right = Rational(3, 4)
        
        let standard = Rational(14, 8)
        
        XCTAssertEqual((left-right).numerator, standard.numerator)
        XCTAssertEqual((left-right).denominator, standard.denominator)
    }
    
    func testExponentiation() {
        let left = Rational(5, 2)
        let power: Int = 3
        
        let standard = Rational(125, 8)
        
        let output = pow(base: left, power: power)
        
        XCTAssertEqual(output.numerator, standard.numerator)
        XCTAssertEqual(output.denominator, standard.denominator)
    }
    
    func testNegativeExponentiation() {
        let left = Rational(5, 2)
        let power: Int = -3
        
        let standard = Rational(8, 125)
        
        let output = pow(base: left, power: power)
        
        XCTAssertEqual(output.numerator, standard.numerator)
        XCTAssertEqual(output.denominator, standard.denominator)
    }
    
    func testZeroExponentiation() {
        let left = Rational(5, 2)
        let power: Int = 0
        
        let standard = Rational(1, 1)
        
        let output = pow(base: left, power: power)
        
        XCTAssertEqual(output.numerator, standard.numerator)
        XCTAssertEqual(output.denominator, standard.denominator)
    }
    
    func testExponentiationOperator() {
        let left = Rational(5, 2)
        let power: Int = 3
        
        let standard = Rational(125, 8)
        
        let output = left^power
        
        XCTAssertEqual(output.numerator, standard.numerator)
        XCTAssertEqual(output.denominator, standard.denominator)
    }
    
    // MARK: - Equivalence and Comparison
    
    func testEquivalence() {
        let left = Rational(4, 2)
        let right = Rational(8, 4)
        
        XCTAssertEqual(left, right)
    }
    
    // MARK: - AdditiveArithmetic
    
    func testInOutAddition() {
        var left = Rational(4, 2)
        let right = Rational(3, 1)
        
        let standard = Rational(10, 2)
        
        left += right
        
        XCTAssertEqual(left, standard)
    }
    
    func testReduce() {
        let numbers = [Rational(2, 2), Rational(3, 3), Rational(4, 4)]
        
        let standard = Rational(24, 24)
        
        XCTAssertEqual(numbers.product(), standard)
    }
}
  

I never had heterogeneous real-complex addition or subtraction, and I've removed the multiply and division operators because of type-inference issues discussed in the README. This is a minor inconvenience for some use cases, and we should investigate ways in which we could make them available--especially multiplication and division, which

• are used more frequently because of the natural interpretation of C as vector space over R
• would allow us to elide a lot of spurious arithmetic operations an get some perf wins

Hey! Very excited about this package (as you might have seen from my three PRs #33 #34 #36 )!

Earlier this summer/fall I started playing around with Abstract Algebra in Swift, I repackage it as a SPM package this hour, just to be able to show it to you fellow math lovers. I just called it AbstractAlgebra - for lack of imagination. It is quite crude and probably ridiculously slow but shows some ideas.

In general, what fits into this package? I fully understand if e.g. Abstract algebra does not fit...

Thanks!

Add a ShapedArray type and/or protocol suitable for use as the common currency type for vectors / linear algebra / image processing / tensors, like ndarray is for numpy.

This should hew pretty closely to what S4TF already does with their ShapedArray type; it doesn't need to be identical, but we need to have a good reason for any changes.

Provide a family of [U]Int128, 256, ... up to a reasonable threshold where flexibly-sized bignum becomes competitive (most likely 1024 bits). This should be done as generically as possible so that support for additional sizes can be added as needed with minimal complexity.

These types should conform to FixedWidthInteger.

In the standard library, Swift uses “unsafe” to denote operations that are memory-unsafe. This sets an expectation which I think we should respect and follow.

Currently, Complex has a member named unsafeLengthSquared. This operation is memory-safe, so in accordance with Swift’s established precedent it should not have “unsafe” in its name.

I propose renaming it to lengthSquared, and documenting that the calculation could overflow (or underflow).

This also helps with performance tuning.

A new module providing additional functionality on the standard library integer types (and more generally, any fixed-width integer type), including:

• - bitwise rotation
• - full add / sub (with carryin/carryout)
• - saturating arithmetic (functions, not operators)
• - division with rounding control
• - shift with rounding control
• - "true" mod function (Euclidean division)
• - GCD
• - division by constant (precompute magic multiply and shift a la libdivide)
• - modulus by constant (use Lemire's method)

I wonder if we might rather define our trig functions in terms of an Angle type:

public struct Angle<T: Real> {
  public var radians: T
  public init(radians: T) { self.radians = radians }
  public static func radians(_ val: T) -> Angle<T> { .init(radians: val) }
  
  public var degrees: T { radians * 180 / .pi }
  public init(degrees: T) { self.init(radians: degrees * .pi / 180) }
  public static func degrees(_ val: T) -> Angle<T> { .init(degrees: val) }
}

Browsing SE-0246, this alternative design doesn't seem to have been considered, but I think it would help to make our functions more useful and self-documenting. Depending on what you're doing, it can be more natural to think in terms of degrees rather than radians. In C, people tend to bring these utility functions/macros by themselves, whereas other languages like Java include them as part of the standard library.

The design of the type given above means that creating and consuming an Angle in radians in essentially free.

Clearly, we cannot make this change to the ElementaryFunctions protocol requirements, as we need the conforming type to provide those implementations. However, we can do this for the free functions - either in addition to, or in place of the regular functions which take a <T:ElementaryFunctions> argument.

Since the Numerics modules cover complex numbers and operations, and you mentioned here (#6) that ShapedArrays are within the scope of the project, would quaternion math and possibly other algebras fall under the scope of Swift Numerics as well?

Use cases include everywhere 3D rotation is used, from games and visualization to simulation and control systems. A well-tested and well-documented implementation would be a big boon for a new math library. Some quick googling shows this is not in the C++ STL or a standard part of NumPy, but it is a standard part of MATLAB and part of Accelerate (maybe better suited as part of an update/rewrite of Accelerate?).

It would be great to be able to work with (brand new shiny big) binary integers and using these to define bounded number types. Number bounded by range, e.g. a BigInt bound from -273 to infinity, representing Kelvin.

I've given an attempt to do this in my SPM Package "DelTal".

/// `BUNCInt` is short _B_ound _U_nsigned _N_amed _C_ategorized _Int_eger
public typealias Kelvin = BUNCInt<NoBound<BigUInt>, KelvinName, Temperatur>

public struct KelvinName: IntegerName {
    public static let nameOfInteger = "Kelvin"
}

public struct Temperatur: IntegerCategory {
    public static let nameOfCategory = "Temperature"
}

public extension Kelvin {
    enum Water {}
}

/* Verbose (needed) syntax for: `extension Kelvin.Water` */
public extension BUNCInt.Water where Name == KelvinName, Category == Temperatur {
    static var meltsAt: Kelvin { .init(magnitude: 273) }
    static var boilsAt: Kelvin { meltsAt + 100 }
}

It allows for expressive and safe wrapping of primitives according to Object-Calisthenics rule #3, and fits perfectly with Swift philosophy of safety.

It would really be amazing if developers were able to get compile-time errors inline, just like, defining this overflowing UInt8 yields a compile-time

let uint8: UInt8 = 1337

Integer literal '1337' overflows when stored into 'UInt8'

Being able to get that for custom defined integer types (or typealiases) would be... awesome!

Currently .magnitude uses the two-norm (sqrt(real*real + imag*imag)), because it's the most obvious choice, but we don't need to use it, and there are some good reasons why we should consider using either the 1-norm (|real| + |imag|) or ∞-norm (max(|real|, |imag|)) instead.

• The two-norm requires special care to avoid spurious overflow/underflow, but the naive expressions for the 1-norm or ∞-norm are always correct.
• Even when care is used, near the overflow boundary the two-norm may not be representable. The ∞-norm in particular is always representable (because of this, the ∞-norm would seem to be the obvious alternative; it gives the nicest API surface).
• Both the 1-norm and ∞-norm are much easier to compute for operators than the 2-norm (O(n) vs. "no closed form expression, but usually O(n^3) iterative methods"), so it would be nice to establish a precedent of .magnitude binding one of these cheaper-to-compute norms.

Because Numeric doesn't specify the semantics, there's no problem changing what we do, but this is a change that would be good to make soon before too many people start using it. We should still provide the two-norm (because it's often useful), but we would simply do it under a different name (length, norm, euclideanNorm, etc).

All of the usual ElementaryFunctions operations should be defined for Complex types. Implementation details may be found in Kahan's "Branch Cuts for Complex Elementary Functions; or Much Ado About Nothing's Sign Bit".

We can simplify the implementation somewhat for Swift complex numbers, because:
(a) we do not need to worry about fenv flags
(b) we do not need to bother with sign of zero/inf, because we're not trying to be precise with those cases (they do not pull their weight, complexity-wise).

I just took a look at this project, and have some suggestions for slight tweaks.

• magnitude should return RealType.magnitude instead of RealType.

• IntegerLiteralType should alias to RealType’s version.

• The Codable conditional conformance could be split into separate Encodable and Decodable conditional conformances.

• The debugDescription could be unconditional.

    extension Complex: CustomDebugConvertible {
        public var debugDescription: String {
            String(describing: Self.self) + “(“ + String(reflecting: x) + “, “ + String(reflecting: y) + “)”
        }
    }
  

A plotting module would be a great component of Swift Numerics. And since this is Swift, a plotting module could utilize Metal for high performance visualization and animation. The Python community provides many examples of plotting libraries such as Matplotlib, Plotly, HoloViews, and many more.

There are a lot of aspects of trigonometry that are easier to express using τ instead of using 2π (see also: https://tauday.com).

Any thoughts on including tau as a constant in this library?

Update SE-0259 for inclusion as a module in swift-numerics, and use it for our tests.

Fields are additive groups, which means they have negation. Therefore the AlgebraicField protocol should refine SignedNumeric.

Currently, Complex.magnitude is calculated using the abs() free function instead of the .magnitude property on its real and imaginary parts. I don’t know if there is any difference at runtime (I suspect not, as the conditional checks in Swift.abs() should get optimized away), but from a code-reading standpoint it seems reasonable to prefer implementing the Complex.magnitude property in terms of the RealType.magnitude property.

At present, the implementations of these operations convert their Int arguments to the floating-point type. This works fine for small integers (the common case), but isn't quite right for large integers (which will round on conversion). This is especially problematic because it can convert odd integers to even integers, which will change the sign of results when the floating-point argument is negative.

Hey, I already asked in #38 about the scope of swift-numerics, but since I have not got a reply yet I keep on pitching some math stuff.

Last summer I wrote EquationKit(just gave SPM support to it) which offers multivariate polynomial with support for partial differentiation:

let polynomial = (3*x + 5*y - 17) * (7*x - 9*y + 23)
print(polynomial) // 21x² + 8xy - 50x - 45y² + 268y - 391)
let number = polynomial.evaluate() {[ x <- 4, y <- 1 ]}
print(number) // 0

let y＇ = polynomial.differentiateWithRespectTo(x)
print(y＇) // 42x + 8y - 50
y＇.evaluate() {[ x <- 1, y <- 1 ]} // 0

let x＇ = polynomial.differentiateWithRespectTo(y)
print(x＇) // 8x - 90y + 268
 x＇.evaluate() {[ x <- 11.5,  y <- 4 ]} // 0

Would that fit in here? I saw #42 and would I would mention this.

Currently, the Complex type uses the word phase in a variety of places.

The technical term for the angle of a complex number is its “argument”, however that word already has a quite different meaning in Swift and programming languages in general, so it would be confusing to use here.

The commonplace, readily-understood term for the angle of a complex number, is, well, its “angle”. There is little to no chance of confusion if we use “angle”.

In contrast, “phase” is far more esoteric, perhaps to the point of jargon. While a mathematician, an electrical engineer, and a quantum physicist (hmm, that reminds me of a joke…) will know what it means, there is a high chance that many programmers in the wider audience will be unfamiliar with the term.

I think it would be most clear to use “angle” instead of “phase” for the Complex type.

Implement the IEEE-754 Float16 type, conforming to BinaryFloatingPoint. Use x86 and arm half<->float conversions, and armv8.2 half-precision arithmetic when available.

So that we do not have to run tests with -c release, and to separate the two more cleanly.

One of Swift's main philosophies is "Clarity is more important than brevity", thus ought we not to spell the math functions out?

In fact the whole standard library is filled with verbose method with clear naming. E.g.

Swift Double's func truncatingRemainder(dividingBy other: Double) -> Double has a verbose and clear naming, compared to Python's not so clear divmod

Thus we maybe ought to rename

func erf

to

func errorFunction

cosh -> hyperbolicCosine or cosineHyperbolic etc.

Offer shorthand syntax as optional package (product)

We could also keep the shorthand syntax for those who want to use that. Maybe that idea is controversial though. We could even let the name of the implementing functions be verbose and create "shorthand syntax" versions in a separate SPM product. And users who think shorthand syntax is neeater can import NumericsShorthandSyntax or import NumbericsBrief or similar.

I think a module that performs symbolic calculations with expressions could be very useful. This would manipulate algebraic expressions by simplifying, factoring, expanding, etc., providing exact answers. Ideally with something like this, you would be able to define your own variables and functions along with rules to simplify them. I also think that this is something that should be including because it is something that is far from trivial to implement.

With Swift's custom operators and strong typing, this will be extremely powerful.

There are plenty of use case for this. Firstly, as previously, this would allow you to solve equations with exact value rather than estimations. Moreover, this would allow for the definition of types whose properties themselves are expressions, thus allowing you to perform calculations and easily derive formulas. For example, one could define a point type with expressions as their coordinates, enabling to easily get the formula for the distance between two points.

struct Point {
    var x, y: Expression
}

let p1 = Point(x: "x1", y: "y1")
let p2 = Point(x: "x2", y: "y2")
let lengthFormula = sqrt((p1.x - p2.x) ^ 2 + (p1.y - p2.y) ^ 2)

print(lengthFormula)
// Prints "√((x1 - x2)² + (y1 - y2)²)"

Expressions would also be able to do things like isolate for variables:

var expression: Expression = "x ^ 2 + 1 = 3y + 5"
expression.isolate(for: "x")

print(expression)
// Prints "x = ±√(3y + 4)"

For mathematical computations, something like this has the potential to really make Swift stand out, opening it to numerous applications.

We should have a protocol that floating-point types as well as Complex conform to, which includes multiplication and division. This will enable people to write more useful generic algorithms for those types.

One option is a protocol (perhaps named Field) that refines Numeric and includes division. Then algorithms could be written against either Field or Field & ElementaryFunctions.

ElementaryFunctions was added to apple/swift with conformances for SIMD types. It seems those conformances were dropped when ElementaryFunctions was moved to apple/swift-numerics.

Is it reasonable to re-add SIMD{n}: ElementaryFunctions conformances to apple/swift-numerics?

Context: tensorflow/swift#404

release plz

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

/// # 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<R: Real, C: Complex>(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
        }
    }
}

Currently, the polar initializer Complex.init?(length:phase:) is failable. However, Complex already has a representation for invalid results, namely (.nan, .nan), which is treated as the point at infinity.

Thus, it would be reasonable for the polar initializer to be non-failable, and always produce (length * .cos(phase), length * .sin(phase)). In the case that either length or phase is non-finite, the result would be Complex.infinity.

This is the same result one would get by first calculating the sine and cosine of the phase as a RealType then constructing the Complex result from its components, so it seems reasonable that it’s what users will expect.

Hi, my framework uses the Numerics package. We developed a model using the SwiftRT API and a version using TensorFlow. On the Mac we get the exact same results. The results were recorded and put into the unit tests to assure a constant point of reference.

On Ubuntu and gLinux, TensorFlow produces the exact same results. SwiftRT using the Numerics package almost produces the same results up to 0.001 difference. I am wondering if the Numerics package is relying on system library functions which are inconsistent across platforms?

If so, this is a big problem. Any insights?

Thanks, Ed

It is common to need a random sample from either a real or complex normal distribution.

The standard library provides static random(in:) methods which sample uniformly, in a constrained extension of BinaryFloatingPoint. (The constraint is RawSignificand: FixedWidthInteger.)

I propose that we add similar functionality for sampling from a normal distribution on Real and Complex with the corresponding constraints.

Hi, I can successfully build on Linux but I get the following warnings from NumericsShims.h:
It would be great to have a clean build :)

Thanks, Ed

$ swift build
In file included from /usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/NumericsShims.c:19:
/usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/include/NumericsShims.h:137:10: warning: implicit declaration of function 'lgammaf_r' is invalid in C99 [-Wimplicit-function-declaratio
  return lgammaf_r(x, signp);
         ^
 warning: implicit declaration of function 'lgamma_r' is invalid in C99 [-Wimplicit-function-declaration]
  return lgamma_r(x, signp);

2 warnings generated.
<module-includes>:1:10: note: in file included from <module-includes>:1:
#include "/usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/include/NumericsShims.h"
         ^
/usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/include/NumericsShims.h:137:10: warning: implicit declaration of function 'lgammaf_r' is invalid in C99
  return lgammaf_r(x, signp);
         ^
<module-includes>:1:10: note: in file included from <module-includes>:1:
#include "/usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/include/NumericsShims.h"
         ^
/usr/local/google/home/ewconnell/swiftrt/.build/checkouts/swift-numerics/Sources/NumericsShims/include/NumericsShims.h:259:10: warning: implicit declaration of function 'lgamma_r' is invalid in C99
  return lgamma_r(x, signp);
         ^

The natural place to write many generic algorithms is the protocol composition AlgebraicField & ElementaryFunctions, as it provides all the basic arithmetic operations including negation and division, as well as powers, roots, and transcendental functions.

To facilitate and encourage this practice, we should provide a convenient name for that composition.