Initial work on docc setup.

Author: stephentyrone

Sources/ComplexModule/Documentation.docc/Complex.md

@@ -0,0 +1,51 @@
+# ``Complex``
+
+## Topics
+
+### Real and imaginary parts
+
+A `Complex` value is represented with two `RealType` values, corresponding to
+the real and imaginary parts of the number:
+```swift
+let z = Complex(1,-1) //  1 - i
+let re = z.real       //  1
+let im = z.imaginary  // -1
+```
+All `Complex` numbers with a non-finite component is treated as a single
+"point at infinity," with infinite magnitude and indeterminant phase. Thus,
+the real and imaginary parts of an infinity are nan.
+```swift
+let w = Complex<Double>.infinity
+w == -w               // true
+let re = w.real       // .nan
+let im = w.imag       // .nan
+```
+See <doc:Infinity> for more details.
+
+- ``init(_:_:)``
+- ``init(_:)-5aesj``
+- ``init(imaginary:)``
+- ``real``
+- ``imaginary``
+
+### Magnitude and norms
+
+See the article <doc:Magnitude> for more details.
+
+- ``magnitude``
+- ``length``
+- ``lengthSquared``
+- ``normalized``
+
+### Polar representations
+
+- ``init(length:phase:)``
+- ``phase``
+- ``length``
+- ``polar``
+
+### Conversions from other types
+
+- ``init(_:)-4csd3``
+- ``init(_:)-80jml``
+- ``init(exactly:)-767k9``

Sources/ComplexModule/Documentation.docc/ComplexModule.md

@@ -0,0 +1,63 @@
+# ``ComplexModule``
+
+Types and operations for working with complex numbers.
+
+## Representation
+
+The `Complex` type is generic over an associated `RealType`; complex numbers
+are represented as two `RealType` values, the real and imaginary parts of the
+number.
+```
+let z = Complex<Double>(1, 2)
+let re = z.real
+let im = z.imaginary
+```
+
+### Memory layout
+
+A `Complex` value is stored as two `RealType` values arranged consecutively
+in memory. Thus it has the same memory layout as:
+- A Fortran complex value built on the corresponding real type (as used
+by BLAS and LAPACK).
+- A C struct with real and imaginary parts and nothing else (as used by
+computational libraries predating C99).
+- A C99 `_Complex` value built on the corresponding real type.
+- A C++ `std::complex` value built on the corresponding real type.
+Functions taking complex arguments in these other languages are not
+automatically converted on import, but you can safely write shims that
+map them into Swift types by converting pointers.
+
+## Real-Complex arithmetic
+
+Because the real numbers are a subset of the complex numbers, many
+languages support arithmetic with mixed real and complex operands.
+For example, C allows the following:
+```c
+#include <complex.h>
+double r = 1;
+double complex z = CMPLX(0, 2); // 2i
+double complex w = r + z;       // 1 + 2i
+```
+The `Complex` type does not provide such mixed operators. There are two
+reasons for this choice. First, Swift generally avoids mixed-type
+arithmetic, period. Second, mixed-type arithmetic operators lead to
+undesirable behavior in common expressions when combined with literal
+type inference. Consider the following example:
+```swift
+let a: Double = 1
+let b = 2*a
+```
+If we had a heterogeneous `*` operator defined, then if there's no prevailing
+type context (i.e. we aren't in an extension on some type), the expression
+`2*a` is ambiguous; `2` could be either a `Double` or `Complex<Double>`. In
+a `Complex` context, the situation is even worse: `2*a` is inferred to have
+type `Complex`.
+
+Therefore, the `Complex` type does not have these operators. In order to write
+the example from C above, you would use an explicit conversion:
+```swift
+import ComplexModule
+let r = 1.0
+let z = Complex<Double>(0, 2)
+let w = Complex(r) + z
+```

Sources/ComplexModule/Documentation.docc/Infinity.md

@@ -0,0 +1,48 @@
+# Zero and infinity
+
+Semantics of `Complex` zero and infinity values, and important considerations
+when porting code from other languages.
+
+Unlike C and C++'s complex types, `Complex` does not attempt to make a 
+semantic distinction between different infinity and NaN values. Any `Complex`
+datum with a non-finite component is treated as the "point at infinity" on 
+the Riemann sphere--a value with infinite magnitude and unspecified phase.
+
+As a consequence, all values with either component infinite or NaN compare
+equal, and hash the same. Similarly, all zero values compare equal and hash
+the same.
+
+## Rationale
+
+This decision has some drawbacks,¹ but also some significant advantages.
+In particular, complex multiplication is one of the most common operations with
+a complex type, and one would like to be able to use the usual naive arithmetic
+implementation, consisting of four real multiplications and two real additions:
+```
+(a + bi) * (c + di) = (ac - bd) + (ad + bc)i
+```
+`Complex` can use this implementation, because we do not differentiate between
+infinities and NaN values. C and C++, by contrast, cannot use this
+implementation by default, because, for example:
+```
+(1 + ∞i) * (0 - 2i) = (1*0 - ∞*(-2)) + (1*(-2) + ∞*0)i
+                    = (0 - ∞) + (-2 + nan)i
+                    = -∞ + nan i
+```
+`Complex` treats this as "infinity", which is the correct result. C and C++
+treat it as a nan value, however, which is incorrect; infinity multiplied
+by a non-zero number should be infinity. Thus, C and C++ (by default) must
+detect these special cases and fix them up, which makes multiplication a
+more computationally expensive operation.²
+
+### Footnotes:
+¹ W. Kahan, Branch Cuts for Complex Elementary Functions, or Much Ado
+About Nothing's Sign Bit. In A. Iserles and M.J.D. Powell, editors,
+_Proceedings The State of Art in Numerical Analysis_, pages 165–211, 1987.
+
+² This can be addressed in C programs by use of the `STDC CX_LIMITED_RANGE`
+pragma, which instructs the compiler to simply not care about these cases.
+Unfortunately, this pragma is not often used in real C or C++ programs
+(though it does see some use in _libraries_). Programmers tend to specify
+`-ffast-math` or maybe `-ffinite-math-only` instead, which has other
+undesirable consequences.

Sources/ComplexModule/Documentation.docc/Magnitude.md

@@ -0,0 +1,137 @@
+# Magnitude and norms
+
+Introduction to the concept of norm and discussion of the `Complex` type's
+`.magnitude` property.
+
+In mathematics, a *norm* is a function that gives each element of a vector
+space a non-negative length.¹
+
+Many different norms can be defined on the complex numbers, viewed as a
+vector space over the reals. All of these norms have some basic
+properties. If we use *‖z‖* to represent any norm of *z*, it must:
+- be *subadditive* (a.k.a. satisfy the triangle inequalty):
+  ‖z + w‖ ≤ ‖z‖ + ‖w‖ for any two complex numbers z and w.
+- be *homogeneous*
+  ‖az‖ = |a|‖z‖ for any real number a and complex number z.
+- and be *positive definite*
+  ‖z‖ is zero if and only if z is zero.
+
+The three most commonly-used norms are:
+- 1-norm ("taxicab norm"):
+  ‖x + iy‖₁ = |x| + |y|
+- 2-norm ("Euclidean norm"):
+  ‖x + iy‖₂ = √(x² + y²)
+- ∞-norm ("maximum norm" or "Чебышёв [Chebyshev] norm")²:
+  ‖x + iy‖ = max(|x|,|y|)
+
+> Exercise:
+> 1. Check that these properties hold for one of the three norms
+     that we just defined. (Hint: write z = a+bi and w = c+di,
+     and use the fact that the absolute value is a norm on the
+     real numbers, and therefore has the same property).
+
+The `Complex` type gives special names to two of these norms; `length`
+for the 2-norm, and `magnitude` for the ∞-norm.
+
+## Magnitude:
+
+The `Numeric` protocol requires us to choose a norm to call `magnitude`,
+but does not give guidance as to which one we should pick. The easiest choice
+might have been the Euclidean norm; it's the one with which people are most
+likely to be familiar.
+
+However, there are good reasons to make a different choice:
+- Computing the Euclidean norm requires special care to avoid spurious
+  overflow/underflow (see implementation notes for `length` below). The
+  naive expressions for the taxicab and maximum norm always give the best
+  answer possible.
+- Even when special care is used, the Euclidean and taxicab norms are
+  not necessarily representable. Both can be infinite even for finite
+  numbers.
+  ```swift
+  let big = Double.greatestFiniteMagnitude
+  let z = Complex(big, big)
+  ```
+  The taxicab norm of `z` would be `big + big`, which overflows; the
+  Euclidean norm would be `sqrt(2) * big`, which also overflows.
+  
+  But the maximum norm is always equal to the magnitude of either `real`
+  or `imaginary`, so it is necessarily representable if `z` is finite.
+- The ∞-norm is the choice of established computational libraries, like
+  BLAS and LAPACK.
+
+For these reasons, the `magnitude` property of `Complex` binds the
+maximum norm:
+```swift
+Complex(2, 3).magnitude    // 3
+Complex(-1, 0.5).magnitude // 1
+```
+
+## Length:
+
+The `length` property of a `Complex` value is its Euclidean norm.
+
+```swift
+Complex(2, 3).length    // 3.605551275463989
+Complex(-1, 0.5).length // 1.118033988749895
+```
+
+Aside from familiarity, the Euclidean norm has one important property
+that the maximum norm lacks:
+- it is *multiplicative*:
+  ‖zw‖₂ = ‖z‖₂‖w‖₂ for any two complex numbers z and w.
+
+> Exercises:
+> 1. Why isn't the maximum norm multiplicative?
+  (Hint: Let `z = Complex(1,1)`, and consider `z*z`.)
+> 2. Is the 1-norm multiplicative?
+
+### Implementation notes:
+
+The `length` property takes special care to produce an accurate answer,
+even when the value is poorly-scaled. The naive expression for `length`
+would be `sqrt(x*x + y*y)`, but this can overflow or underflow even when
+the final result should be a finite number.
+```swift
+// Suppose that length were implemented like this:
+extension Complex {
+  var naiveLength: RealType {
+    .sqrt(real*real + imaginary*imaginary)
+  }
+}
+
+// Then taking the length of even a modestly large number:
+let z = Complex<Float>(1e20, 1e20)
+// or small number:
+let w = Complex<Float>(1e-24, 1e-24)
+// would overflow:
+z.naiveLength // Inf
+// or underflow:
+w.naiveLength // 0
+```
+Instead, `length` is implemented using a two-step algorithm. First we
+compute `lengthSquared`, which is `x*x + y*y`. If this is a normal
+number (meaning that no overflow or underflow has occured), we can safely
+return its square root. Otherwise, we redo the computation with a more
+careful computation, which avoids spurious under- or overflow:
+```swift
+let z = Complex<Float>(1e20, 1e20)
+let w = Complex<Float>(1e-24, 1e-24)
+z.length // 1.41421358E+20
+w.length // 1.41421362E-24
+```
+
+### Footnotes:
+
+¹ Throughout this documentation, "norm" refers to a
+  [vector norm](https://en.wikipedia.org/wiki/Norm_(mathematics)).
+  To confuse the matter, there are several similar things also called
+  "norm" in mathematics. The other one you are most likely to run into
+  is the [field norm](https://en.wikipedia.org/wiki/Field_norm).
+  
+  Field norms are much less common than vector norms, but the C++
+  `std::norm` operation implements a field norm. To get the (Euclidean)
+  vector norm in C++, use `std::abs`.
+
+² There's no subscript-∞ in unicode, so I write the infinity norm
+  without the usual subscript.

Sources/RealModule/Documentation.docc/Augmented.md

@@ -0,0 +1,20 @@
+# ``Augmented``
+
+## Overview
+
+Consider multiplying two Doubles. A Double has 53 significand bits, so their
+product could be up to 106 bits wide before it is rounded to a Double result.
+So up to 53 of those 106 bits will be "lost" in that process:
+```swift
+let a = 1.0 + .ulpOfOne // 1 + 2⁻⁵²
+let b = 1.0 - .ulpOfOne // 1 - 2⁻⁵²
+let c = a * b           // 1 - 2⁻¹⁰⁴ before rounding, rounds to 1.0
+```
+Sometimes it is necessary to preserve some or all of those low-order bits;
+maybe a subsequent subtraction cancels most of the high-order bits, and so
+the low-order part of the product suddenly becomes significant:
+```swift
+let result = 1 - c      // exactly zero, but "should be" 2⁻¹⁰⁴
+```
+Augmented arithmetic is a building-block that library writers can use to
+handle cases like this more carefully.

Sources/RealModule/Documentation.docc/RealModule.md

@@ -0,0 +1,13 @@
+# ``RealModule``
+
+<!--@START_MENU_TOKEN@-->Summary<!--@END_MENU_TOKEN@-->
+
+## Overview
+
+<!--@START_MENU_TOKEN@-->Text<!--@END_MENU_TOKEN@-->
+
+## Topics
+
+### <!--@START_MENU_TOKEN@-->Group<!--@END_MENU_TOKEN@-->
+
+- <!--@START_MENU_TOKEN@-->``Symbol``<!--@END_MENU_TOKEN@-->

Sources/RealModule/Documentation.docc/Relaxed.md

@@ -0,0 +1,46 @@
+# ``Relaxed``
+
+## Overview
+
+Because of rounding, and the arithmetic rules for infinity and NaN values,
+floating-point addition and multiplication are not associative:
+```swift
+let ε = Double.leastNormalMagnitude
+let sumLeft  = (-1 + 1) + ε  //  0 + ε = ε
+let sumRight =  -1 + (1 + ε) // -1 + 1 = 0
+
+let ∞ = Double.infinity
+let productLeft  = (ε * ε) * ∞  // 0 * ∞ = .nan
+let productRight =  ε * (ε * ∞) // ε * ∞ = ∞
+```
+For some algorithms, the distinction between these results is incidental; for
+some others it is critical to their correct function. Because of this,
+compilers cannot freely change the order of reductions, which prevents some
+important optimizations: extraction of instruction-level parallelism and
+vectorization.
+
+If you know that you are in a case where the order of elements being summed
+or multiplied is incidental, the Relaxed operations give you a mechanism
+to communicate that to the compiler and unlock these optimizations. For
+example, consider the following two functions:
+```swift
+func sum(array: [Float]) -> Float {
+  array.reduce(0, +)
+}
+
+func relaxedSum(array: [Float]) -> Float {
+  array.reduce(0, Relaxed.sum)
+}
+```
+when called on an array with 1000 elements in a Release build, `relaxedSum`
+is about 8x faster than `sum` on Apple M2, with a similar speedup on Intel
+processors, without the need for any unsafe code or flags.
+
+### multiplyAdd
+
+In addition to `Relaxed.sum` and `Relaxed.product`, `Relaxed` provides the
+``multiplyAdd(_:_:_:)`` operation, which communciates to the compiler that
+it is allowed to replace separate multiply and add operations with a single
+_fused multiply-add_ instruction if its cost model indicates that it would
+be advantageous to do so. When targeting processors that support this
+instruction, this may be a significant performance advantage.