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

[Kirby Urner]

At 07:14 PM 10/1/2001 -0400, Tim Peters wrote:
>[Tim]
> >> 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.
>
>[Kirby Urner]
> > Not following this.
> >
> >    >>> def altdrs(n):
> >             return round(fact(n)/e)
> >
> > returns the same answer as drs(n) right down to n=1.
>
>Then 1 certainly isn't the smallest n for which equality fails <wink>.


OK, I get it now:

   >>> def notsame():
           n = 1
           while altdrs(n)==drs(n):
              n = n+1
           return n

   >>> notsame()
   19
   >>> drs(19)
   44750731559645106
   >>> altdrs(19)
   44750731559645112.0

I.e. at n=18, we're still in the money:

   >>> drs(18)
   2355301661033953
   >>> altdrs(18)
   2355301661033953.0


> > ...
> >            while float(v) <> phi:
>
>Note that loops like this may never terminate:


Very true.  Would be better to have some loop that detects
when there's no difference between successive iterations
and cuts it off there -- assumes the series is convergent
in a nice way.

Or you could do some abs(target-current)<epsilon, but have
to choose epsilon with some care.

>It's a pragmatic ("programming") issue, not a semantic issue:  the purpose
>of special-casing 1 and -1 within F.__mul__ would not be to deliver a
>different answer, but to return the correct answer *faster*.


OK, I understand.

My code doesn't do any case checking for 1, -1, but does a lot
of instance checking, as my fraction objects may have matrices
or polynomials to deal with (which doesn't mean I've implemented
rational expressions -- just means you can multiply a polynomial
by (1/2) and have the coefficients be fractions:

    >>> p
    x**6 + 2*x**5 + 3*x**4 + 4*x**3 + 4*x**2 + 2
    >>> p*F(1,2)
    (1/2)*x**6 + x**5 + (3/2)*x**4 + 2*x**3 + 2*x**2 - 0*x + 1
                                                       ^^^

Oops, I see a bug I need to deal with).

Kirby