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/scripts/karatsuba.py`) that
  29 will 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.
  66
  67 Its complexity is `O(log(n)*n^2)` as it requires one division per iteration, and
  68 it doubles the amount of correct digits per iteration.
  69
  70 ### Sine and Cosine (`bc` Math Library Only)
  71
  72 This `bc` uses the series
  73
  74 ```
  75 x - x^3/3! + x^5/5! - x^7/7! + ...
  76 ```
  77
  78 to calculate `sin(x)` and `cos(x)`. It also uses the relation
  79
  80 ```
  81 cos(x) = sin(x + pi/2)
  82 ```
  83
  84 to calculate `cos(x)`. It has a complexity of `O(n^3)`.
  85
  86 **Note**: this series has a tendency to *occasionally* produce an error of 1
  87 [ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't
  88 any way around it; [this article][7] explains why calculating sine and cosine,
  89 and the other transcendental functions below, within less than 1 ULP is nearly
  90 impossible and unnecessary.) Therefore, I recommend that users do their
  91 calculations with the precision (`scale`) set to at least 1 greater than is
  92 needed.
  93
  94 ### Exponentiation (`bc` Math Library Only)
  95
  96 This `bc` uses the series
  97
  98 ```
  99 1 + x + x^2/2! + x^3/3! + ...
 100 ```
 101
 102 to calculate `e^x`. Since this only works when `x` is small, it uses
 103
 104 ```
 105 e^x = (e^(x/2))^2
 106 ```
 107
 108 to reduce `x`.
 109
 110 It has a complexity of `O(n^3)`.
 111
 112 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 113 their calculations with the precision (`scale`) set to at least 1 greater than
 114 is needed.
 115
 116 ### Natural Logarithm (`bc` Math Library Only)
 117
 118 This `bc` uses the series
 119
 120 ```
 121 a + a^3/3 + a^5/5 + ...
 122 ```
 123
 124 (where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small
 125 and uses the relation
 126
 127 ```
 128 ln(x^2) = 2 * ln(x)
 129 ```
 130
 131 to sufficiently reduce `x`.
 132
 133 It has a complexity of `O(n^3)`.
 134
 135 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 136 their calculations with the precision (`scale`) set to at least 1 greater than
 137 is needed.
 138
 139 ### Arctangent (`bc` Math Library Only)
 140
 141 This `bc` uses the series
 142
 143 ```
 144 x - x^3/3 + x^5/5 - x^7/7 + ...
 145 ```
 146
 147 to calculate `atan(x)` for small `x` and the relation
 148
 149 ```
 150 atan(x) = atan(c) + atan((x - c)/(1 + x * c))
 151 ```
 152
 153 to reduce `x` to small enough. It has a complexity of `O(n^3)`.
 154
 155 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 156 their calculations with the precision (`scale`) set to at least 1 greater than
 157 is needed.
 158
 159 ### Bessel (`bc` Math Library Only)
 160
 161 This `bc` uses the series
 162
 163 ```
 164 x^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ...
 165 ```
 166
 167 to calculate the bessel function (integer order only).
 168
 169 It also uses the relation
 170
 171 ```
 172 j(-n,x) = (-1)^n * j(n,x)
 173 ```
 174
 175 to calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`.
 176
 177 **Note**: this series can also produce errors of 1 ULP, so I recommend users do
 178 their calculations with the precision (`scale`) set to at least 1 greater than
 179 is needed.
 180
 181 ### Modular Exponentiation (`dc` Only)
 182
 183 This `dc` uses the [Memory-efficient method][8] to compute modular
 184 exponentiation. The complexity is `O(e*n^2)`, which may initially seem
 185 inefficient, but `n` is kept small by maintaining small numbers. In practice, it
 186 is extremely fast.
 187
 188 ### Non-Integer Exponentiation (`bc` Math Library 2 Only)
 189
 190 This is implemented in the function `p(x,y)`.
 191
 192 The algorithm used is to use the formula `e(y*l(x))`.
 193
 194 It has a complexity of `O(n^3)` because both `e()` and `l()` do.
 195
 196 ### Rounding (`bc` Math Library 2 Only)
 197
 198 This is implemented in the function `r(x,p)`.
 199
 200 The algorithm is a simple method to check if rounding away from zero is
 201 necessary, and if so, adds `1e10^p`.
 202
 203 It has a complexity of `O(n)` because of add.
 204
 205 ### Ceiling (`bc` Math Library 2 Only)
 206
 207 This is implemented in the function `ceil(x,p)`.
 208
 209 The algorithm is a simple add of one less decimal place than `p`.
 210
 211 It has a complexity of `O(n)` because of add.
 212
 213 ### Factorial (`bc` Math Library 2 Only)
 214
 215 This is implemented in the function `f(n)`.
 216
 217 The algorithm is a simple multiplication loop.
 218
 219 It has a complexity of `O(n^3)` because of linear amount of `O(n^2)`
 220 multiplications.
 221
 222 ### Permutations (`bc` Math Library 2 Only)
 223
 224 This is implemented in the function `perm(n,k)`.
 225
 226 The algorithm is to use the formula `n!/(n-k)!`.
 227
 228 It has a complexity of `O(n^3)` because of the division and factorials.
 229
 230 ### Combinations (`bc` Math Library 2 Only)
 231
 232 This is implemented in the function `comb(n,r)`.
 233
 234 The algorithm is to use the formula `n!/r!*(n-r)!`.
 235
 236 It has a complexity of `O(n^3)` because of the division and factorials.
 237
 238 ### Logarithm of Any Base (`bc` Math Library 2 Only)
 239
 240 This is implemented in the function `log(x,b)`.
 241
 242 The algorithm is to use the formula `l(x)/l(b)` with double the `scale` because
 243 there is no good way of knowing how many digits of precision are needed when
 244 switching bases.
 245
 246 It has a complexity of `O(n^3)` because of the division and `l()`.
 247
 248 ### Logarithm of Base 2 (`bc` Math Library 2 Only)
 249
 250 This is implemented in the function `l2(x)`.
 251
 252 This is a convenience wrapper around `log(x,2)`.
 253
 254 ### Logarithm of Base 10 (`bc` Math Library 2 Only)
 255
 256 This is implemented in the function `l10(x)`.
 257
 258 This is a convenience wrapper around `log(x,10)`.
 259
 260 ### Root (`bc` Math Library 2 Only)
 261
 262 This is implemented in the function `root(x,n)`.
 263
 264 The algorithm is [Newton's method][9]. The initial guess is calculated as
 265 `10^ceil(length(x)/n)`.
 266
 267 Like square root, its complexity is `O(log(n)*n^2)` as it requires one division
 268 per iteration, and it doubles the amount of correct digits per iteration.
 269
 270 ### Cube Root (`bc` Math Library 2 Only)
 271
 272 This is implemented in the function `cbrt(x)`.
 273
 274 This is a convenience wrapper around `root(x,3)`.
 275
 276 ### Greatest Common Divisor (`bc` Math Library 2 Only)
 277
 278 This is implemented in the function `gcd(a,b)`.
 279
 280 The algorithm is an iterative version of the [Euclidean Algorithm][10].
 281
 282 It has a complexity of `O(n^4)` because it has a linear number of divisions.
 283
 284 This function ensures that `a` is always bigger than `b` before starting the
 285 algorithm.
 286
 287 ### Least Common Multiple (`bc` Math Library 2 Only)
 288
 289 This is implemented in the function `lcm(a,b)`.
 290
 291 The algorithm uses the formula `a*b/gcd(a,b)`.
 292
 293 It has a complexity of `O(n^4)` because of `gcd()`.
 294
 295 ### Pi (`bc` Math Library 2 Only)
 296
 297 This is implemented in the function `pi(s)`.
 298
 299 The algorithm uses the formula `4*a(1)`.
 300
 301 It has a complexity of `O(n^3)` because of arctangent.
 302
 303 ### Tangent (`bc` Math Library 2 Only)
 304
 305 This is implemented in the function `t(x)`.
 306
 307 The algorithm uses the formula `s(x)/c(x)`.
 308
 309 It has a complexity of `O(n^3)` because of sine, cosine, and division.
 310
 311 ### Atan2 (`bc` Math Library 2 Only)
 312
 313 This is implemented in the function `a2(y,x)`.
 314
 315 The algorithm uses the [standard formulas][11].
 316
 317 It has a complexity of `O(n^3)` because of arctangent.
 318
 319 [1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm
 320 [2]: https://en.wikipedia.org/wiki/Long_division
 321 [3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring
 322 [4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number
 323 [5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
 324 [6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place
 325 [7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
 326 [8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method
 327 [9]: https://en.wikipedia.org/wiki/Root-finding_algorithms#Newton's_method_(and_similar_derivative-based_methods)
 328 [10]: https://en.wikipedia.org/wiki/Euclidean_algorithm
 329 [11]: https://en.wikipedia.org/wiki/Atan2#Definition_and_computation