[Edu-sig] <FUN>more math with Python</FUN>

[Kirby Urner]

>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