In fact, it's such a good approximation that
drs(n) == round(n!/e)
where round(x) rounds x to the nearest integer, for all n >= 1.
A very excellent point.
If the students are advanced enough, it's an interesting exercise to find the smallest n for which the limitations of floating-point arithmetic cause round(n!/e) to return a different result than your exact computation using rationals.
Not following this.
def altdrs(n): return round(fact(n)/e)
returns the same answer as drs(n) right down to n=1. A related question (that I understand) would be how far do we need to push the series in brackets to get 1/e within the limitations of floating point arithmetic:
from math import e def e_neg1(): f = 0 sign = 1 t = 0 while float(f)<>(1/e): f += sign * F(1, fact(t)) sign = -sign t += 1 return (t-1,f)
e_neg1() (18, (138547156531409/376610217984000))
I get 18 terms, with the accompanying fraction. >>> 138547156531409/376610217984000 # (import future division) 0.36787944117144233 >>> 1/e 0.36787944117144233 This is related to the thread wherein we get phi to the limits of floating point, using the continued fraction: phi= 1 + 1 -------------- 1 + 1 -------- 1 + 1 ---- 1 + ... It's related because evaluating a continued fraction involves recursive use of the Fraction object F:
def evalcf(qs): """ Evaluate continued fraction with partial quotients qs """ if len(qs)==1: return F(qs[0],1) return F(qs[0],1) + F(1,evalcf(qs[1:]))
evalcf([1,1,1,1,1]) (8/5) evalcf([1,1,1,1,1,1,1,1,1]) (55/34) float(evalcf([1,1,1,1,1,1,1,1,1])) 1.6176470588235294
from math import sqrt phi = (1 + sqrt(5))/2 def getphi(): qs = [] v = 0 while float(v) <> phi: qs.append(1) v = evalcf(qs) return (len(qs),v)
getphi() (40, (165580141/102334155))
And later you flipped in a different way, by testing the parity of t. A common idiom in numeric programs is to multiply by 1 or -1 instead, like so:
Slaps forehead. Yes, of course! This is what you see most often.
f = 0 sign = 1 for t in range(n+1): f += sign * F(1, fact(t)) sign = -sign return fact(n) * f
As a *programming* issue, it's thus interesting to consider whether F.__mul__ should or shouldn't try to recognize the special cases of multiplication by 1 and -1. There's isn't a "one size fits all" answer, so it's not for students with low tolerance for uncertainty <wink>.
Not following again. Like, of course it should, right?
f = F(1,2) f (1/2) f*1 (1/2) f*-1 (-1/2)
What other behavior would I want? Kirby