Imperative Programming Part 1.
Use a
Double base so that the reciprocal is not truncated by integer division. Repeated multiplication gives the requested linear-work approximation to Euler's number.def eaLinear(n: Int): Double = {
require(n > 0)
val base = 1.0 + 1.0 / n.toDouble
var i = 0
var result = 1.0
while (i < n) {
result *= base
i += 1
}
result
}In exact arithmetic over the real numbers, the loop invariant is and . It holds initially, one multiplication establishes the next exponent, and at exit gives the specified power. The nonnegative loop variant is . There are exactly loop multiplications, hence in asymptotic analysis, and one division to form the base. With floating-point arithmetic, the invariant describes the intended exact result; the actual output has the corresponding rounding errors.
def eaFast(n: Int): Double = {
require(n > 0)
var remaining = n
var base = 1.0 + 1.0 / n.toDouble
var result = 1.0
while (remaining > 0) {
if (remaining % 2 == 1) result *= base
base *= base
remaining /= 2
}
result
}Writing , the exact-arithmetic loop invariant is and . For an even integer exponent , replace the base by its square and the exponent by . For an odd integer exponent , first multiply the accumulator by the base and then make the same replacement. Both preserve the product. The loop variant is
remaining, which strictly decreases while positive. At exit the remaining exponent is zero, so the returned value approximates .There are iterations, squarings and one accumulator multiplication for each set bit of : between and multiplications. Thus the bound is as , or including the small case . All power variables are
Double; no unintended Int squaring occurs.Maintain the reciprocal factorial term rather than calculating large factorials separately:
def eb(n: Int): Double = {
require(n > 0)
var i = 0
var term = 1.0
var sum = 1.0
while (i < n) {
i += 1
term /= i.toDouble
sum += term
}
sum
}The exact-arithmetic loop invariant is , , and . It holds for because . Dividing by the new index produces the next reciprocal factorial; adding it extends the partial sum. The loop variant is . At exit the sum includes all terms from through , using exactly divisions and no multiplications in the loop. In particular, the initial
1.0 accounts for the zero-index term. The work bound is , with constant auxiliary space.Avoid forming the enormous powers in the original ratios. Put , and ; algebra givesUse simultaneous exponentiation by squaring to compute the two powers with one loop and a shared remaining exponent:
def ec(n: Int): Double = {
require(n > 1)
val dn = n.toDouble
val a = 1.0 + 1.0 / dn
val b = 1.0 + 1.0 / (dn - 1.0)
var remaining = n - 1
var u = a
var v = b
var pa = 1.0
var pb = 1.0
while (remaining > 0) {
if (remaining % 2 == 1) {
pa *= u
pb *= v
}
u *= u
v *= v
remaining /= 2
}
((dn + 1.0) * a) * pa - dn * pb
}The exact-arithmetic loop invariant is , and . The even integer and odd integer cases preserve both equalities in parallel. The loop variant is
remaining. On termination, the displayed expression is exactly the desired approximation in the ideal arithmetic model. For , the loop uses at most multiplications, followed by three more: multiplications and one loop.Converting
n to Double before adding one also avoids overflowing n + 1 at Int.MaxValue. The algebra avoids power overflow but does not remove catastrophic cancellation: the final subtraction combines quantities of order to obtain a result of order . For sufficiently large , floating-point rounding can outweigh the mathematical truncation error. For example, at n == 1000000 this Double computation returns approximately 2.7180521293, although the exact formula is approximately 2.71828182845916. Arbitrarily many requested decimal digits require increased arithmetic precision, not simply a larger Int argument.Use the reciprocal-factorial sum for increasing accuracy. Let be the requested number of decimal places and use absolute error tolerance proportional to . Expanding gives and , so these are genuine first nonzero error orders. The first approximation needs of order . Even its fast implementation therefore uses multiplications, while the linear version uses exponentially many in . The third approximation needs of order , but its fast powering still uses multiplications, with two parallel accumulators and a cancellation-prone subtraction.
For the Taylor series sum, Stirling's formula gives , so a sufficient term count grows as . Its divisions therefore grow more slowly than the multiplication counts of either fast-power alternative. More concretely, its positive tail satisfiesThe first omitted term also gives a lower bound . This is the reciprocal-factorial series remainder bound, supplying an explicit stopping certificate rather than relying only on Big O notation. It follows by bounding every successive tail denominator factor below by .
The sum also avoids subtraction of nearly equal large values and has small positive terms, making it numerically preferable. For a particular small accuracy request, constants and actual term counts may be compared directly; an asymptotic error bound alone does not prove a universal finite-input operation ranking. With ordinary
Double, sufficiently many tiny positive terms stop changing the rounded sum, so beyond machine precision a higher-precision version is necessary regardless of the formula.
Articles by others on the same topic
There are currently no matching articles.