contrib/bc/manuals/algorithms.md

   1 # Algorithms
   2
   3 This `bc` uses the math algorithms below:
   4
   5 ### Addition
   6
   7 This `bc` uses brute force addition, which is linear (`O(n)`) in the number of
   8 digits.
   9
  10 ### Subtraction
  11
  12 This `bc` uses brute force subtraction, which is linear (`O(n)`) in the number
  13 of digits.
  14
  15 ### Multiplication
  16
  17 This `bc` uses two algorithms: [Karatsuba][1] and brute force.
  18
  19 Karatsuba is used for "large" numbers. ("Large" numbers are defined as any
  20 number with `BC_NUM_KARATSUBA_LEN` digits or larger. `BC_NUM_KARATSUBA_LEN` has
  21 a sane default, but may be configured by the user.) Karatsuba, as implemented in
  22 this `bc`, is superlinear but subpolynomial (bounded by `O(n^log_2(3))`).
  23
  24 Brute force multiplication is used below `BC_NUM_KARATSUBA_LEN` digits. It is
  25 polynomial (`O(n^2)`), but since Karatsuba requires both more intermediate
  26 values (which translate to memory allocations) and a few more additions, there
  27 is a "break even" point in the number of digits where brute force multiplication
  28 is faster than Karatsuba. There is a script (`$ROOT/karatsuba.py`) that will
  29 find the break even point on a particular machine.
  30
  31 ***WARNING: The Karatsuba script requires Python 3.***
  32
  33 ### Division
  34
  35 This `bc` uses Algorithm D ([long division][2]). Long division is polynomial
  36 (`O(n^2)`), but unlike Karatsuba, any division "divide and conquer" algorithm
  37 reaches its "break even" point with significantly larger numbers. "Fast"
  38 algorithms become less attractive with division as this operation typically
  39 reduces the problem size.
  40
  41 While the implementation of long division may appear to use the subtractive
  42 chunking method, it only uses subtraction to find a quotient digit. It avoids
  43 unnecessary work by aligning digits prior to performing subtraction and finding
  44 a starting guess for the quotient.
  45
  46 Subtraction was used instead of multiplication for two reasons:
  47
  48 1.      Division and subtraction can share code (one of the less important goals of
  49         this `bc` is small code).
  50 2.      It minimizes algorithmic complexity.
  51
  52 Using multiplication would make division have the even worse algorithmic
  53 complexity of `O(n^(2*log_2(3)))` (best case) and `O(n^3)` (worst case).
  54
  55 ### Power
  56
  57 This `bc` implements [Exponentiation by Squaring][3], which (via Karatsuba) has
  58 a complexity of `O((n*log(n))^log_2(3))` which is favorable to the
  59 `O((n*log(n))^2)` without Karatsuba.
  60
  61 ### Square Root
  62
  63 This `bc` implements the fast algorithm [Newton's Method][4] (also known as the
  64 Newton-Raphson Method, or the [Babylonian Method][5]) to perform the square root
  65 operation. Its complexity is `O(log(n)*n^2)` as it requires one division per
  66 iteration.
  67
  68 ### Sine and Cosine (`bc` Only)
  69
  70 This `bc` uses the series
  71
  72 ```
  73 x - x^3/3! + x^5/5! - x^7/7! + ...
  74 ```
  75
  76 to calculate `sin(x)` and `cos(x)`. It also uses the relation
  77
  78 ```
  79 cos(x) = sin(x + pi/2)
  80 ```
  81
  82 to calculate `cos(x)`. It has a complexity of `O(n^3)`.
  83
  84 **Note**: this series has a tendency to *occasionally* produce an error of 1
  85 [ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't
  86 any way around it; [this article][7] explains why calculating sine and cosine,
  87 and the other transcendental functions below, within less than 1 ULP is nearly
  88 impossible and unnecessary.) Therefore, I recommend that users do their
  89 calculations with the precision (`scale`) set to at least 1 greater than is
  90 needed.
  91
  92 ### Exponentiation (`bc` Only)
  93
  94 This `bc` uses the series
  95
  96 ```
  97 1 + x + x^2/2! + x^3/3! + ...
  98 ```
  99
 100 to calculate `e^x`. Since this only works when `x` is small, it uses
 101
 102 ```
 103 e^x = (e^(x/2))^2
 104 ```
 105
 106 to reduce `x`. It has a complexity of `O(n^3)`.
 107
 108 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 109 their calculations with the precision (`scale`) set to at least 1 greater than
 110 is needed.
 111
 112 ### Natural Logarithm (`bc` Only)
 113
 114 This `bc` uses the series
 115
 116 ```
 117 a + a^3/3 + a^5/5 + ...
 118 ```
 119
 120 (where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small
 121 and uses the relation
 122
 123 ```
 124 ln(x^2) = 2 * ln(x)
 125 ```
 126
 127 to sufficiently reduce `x`. It has a complexity of `O(n^3)`.
 128
 129 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 130 their calculations with the precision (`scale`) set to at least 1 greater than
 131 is needed.
 132
 133 ### Arctangent (`bc` Only)
 134
 135 This `bc` uses the series
 136
 137 ```
 138 x - x^3/3 + x^5/5 - x^7/7 + ...
 139 ```
 140
 141 to calculate `atan(x)` for small `x` and the relation
 142
 143 ```
 144 atan(x) = atan(c) + atan((x - c)/(1 + x * c))
 145 ```
 146
 147 to reduce `x` to small enough. It has a complexity of `O(n^3)`.
 148
 149 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 150 their calculations with the precision (`scale`) set to at least 1 greater than
 151 is needed.
 152
 153 ### Bessel (`bc` Only)
 154
 155 This `bc` uses the series
 156
 157 ```
 158 x^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ...
 159 ```
 160
 161 to calculate the bessel function (integer order only).
 162
 163 It also uses the relation
 164
 165 ```
 166 j(-n,x) = (-1)^n * j(n,x)
 167 ```
 168
 169 to calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`.
 170
 171 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 172 their calculations with the precision (`scale`) set to at least 1 greater than
 173 is needed.
 174
 175 ### Modular Exponentiation (`dc` Only)
 176
 177 This `dc` uses the [Memory-efficient method][8] to compute modular
 178 exponentiation. The complexity is `O(e*n^2)`, which may initially seem
 179 inefficient, but `n` is kept small by maintaining small numbers. In practice, it
 180 is extremely fast.
 181
 182 [1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm
 183 [2]: https://en.wikipedia.org/wiki/Long_division
 184 [3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring
 185 [4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number
 186 [5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
 187 [6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place
 188 [7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
 189 [8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method