Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 06:32 PM 1/7/2007, Terry Carroll wrote:
>I may add this algorithm to the cookbook.

You should.

Dick




___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Sun, 7 Jan 2007, Terry Carroll wrote:

> ...Say you want to get an approximation of 0.096.  Well, a good first
> approximation is 1/10.  But 1/10 is 0.010

Um, that should be "...is 0.100"


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Sat, 6 Jan 2007, Dick Moores wrote:

> Well, I have to admit I don't understand your code at all. But I see it
> works.

A continuing fraction is an expression that looks like this:

1
A +  --
   1
  B + -
 1
 C + --
1
  D + -
 E + etc.

i.e., A + (1 / (B + (1 / (C + (1 / (D + (1 / E + ... )))

That is, the denominator of the fraction itself has a term that is a 
fraction; and that fraction itself has a term that is a fraction, etc., 
either to infinity or to some finite point.

The basic idea of using a continued fraction to get a rational 
approximation is something like this.   Say you want to get an 
approximation of 0.096.  Well, a good first approximation is 1/10.  But 
1/10 is 0.010.  That's too big, or, to put another way, the denominator is 
too small.  It should really be 1/(10+something), where "something" is 
some fraction 1/x (if it were larger than 1, then we'd be starting from 
11+ something).

The next best approximation is to figure out that "something" is.   It 
turns out, it's 1/2, i.e. .096 = about 1/(10+(1/2)), or about 0.9624.  
That's still too big, so that 1/2 we added should really be 1/(2+ 
something).

The next "something" we calculate is again 1/2, so it's 
1/(10+(1/(2+(1/2)), or .096154.  Getting closer, but the same thing.

Again, the next something is 1/2, so we get 1/(10+(1/(2+(1/(2+1/2))), 
which is 0.096, exactly.

So, in that butt-ugly ascii I have above, A=10, B = 2, C = 2 and D = 2,
and it turns out not to have been necessary to go any further.

(I admit I cheated above and picked 0.096 as an example, because some 
other examples, like .093 and .098 actually have a "something" that's 
equal to 1).

Now, I actually don't understand all of the math in the article at 
http://mathworld.wolfram.com/ContinuedFraction.html , which I based this 
on; but enough of it is clear that I could basically just take certain 
equations and translate them into Python:


> a, r, p, q = [], [], [], []

This is a list of terms; "a" is the successive denominator, i.e. compared
to my description above a[0] is A; a[1] is B; a[2] is C, etc.  r[n] is
used to calculate a[n]; it actually could have been dispensed with.  p[n]
and q[n] are the value of the continued fraction after n iterations.

> a.append(None)
> r.append(None)
> p.append(None)
> q.append(None)

This is just so I can assign to a[n], r[n], etc.; it makes the code more 
readable.

> if n == 0: r[n] = x
> else: r[n] = 1/(r[n-1]-a[n-1])

These statements are equations 8 & 9 from the Wolfram page cited above.

> a[n] = int(r[n])

This is equation 10.

> if n == 0:
> p[n] = a[0]
> q[n] = 1

These are effectively equations 25 & 26 for the special case of n == 0, 
i.e. equation 24.

> elif n ==1:
> p[n] = a[n]*p[n-1] + 1
> q[n] = a[n]

These are equations 25 & 26 for the special case of n == 1, taking into 
account equations 22 to get values for what would have otherwise looking 
for p[-1], i.e. an entry before the first element in the list.

> else:
> p[n] = a[n]*p[n-1] + p[n-2]
> q[n] = a[n]*q[n-1] + q[n-2]

This is the general case for equations 25 & 26.

The rest is just cranking through the iterations.

> (BTW your pi is a bit off but I used yours, instead of math.pi, which 
> is 3.1415926535897931 . 

Oops.  I transposed the "79" and "89"; that's what I get for going from 
memory.

I may add this algorithm to the cookbook.

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 10:40 AM 1/6/2007, Dick Moores wrote:
>At 11:21 PM 1/4/2007, Terry Carroll wrote:
> >On Wed, 3 Jan 2007, Dick Moores wrote:
> >
> > > Be that as it may, farey() is an amazing program.
> >
> >Not to beat this subject to death, but the comment at the bottom of
> >http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317 about
> >continued fractions piqued my interest.  I'm no mathematician, but I
> >encountered continued fractions a long time ago and was fascinated by
> >them.  So I read the URL pointed to,
> >http://mathworld.wolfram.com/ContinuedFraction.html , and came up with the
> >following:
> >
> >#
> >
> >def cf(x, tol=0.0001, Trace=False):
> > """
> > Calculate rational approximation of x to within tolerance of tol;
> > returns a tuple consisting of numerator and denominator p/q
> > Trace=True causes iterated results to be shown
> > """
> > a, r, p, q = [], [], [], []
> > Done = False
> > n = 0
> > if Trace: print "x:%f tol:%f" % (x, tol)
> > while not Done:
> > a.append(None)
> > r.append(None)
> > p.append(None)
> > q.append(None)
> > if n == 0: r[n] = x
> > else: r[n] = 1/(r[n-1]-a[n-1])
> > a[n] = int(r[n])
> > if n == 0:
> > p[n] = a[0]
> > q[n] = 1
> > elif n ==1:
> > p[n] = a[n]*p[n-1] + 1
> > q[n] = a[n]
> > else:
> > p[n] = a[n]*p[n-1] + p[n-2]
> > q[n] = a[n]*q[n-1] + q[n-2]
> > if Trace:
> > print "n:%d a:%d p:%d q:%d approx:%f" % \
> >   (n, a[n], p[n], q[n], float(p[n])/q[n])
> > if abs(float(p[n])/q[n] - x) < tol:
> > Done = True
> > num = p[n]; denom = q[n]
> > n += 1
> > return (num, denom)
> >
> >#
> >
> >Here's a result for pi:
> >
> > >>> print cf(3.14159265357989,0.001, Trace=True)
> >x:3.141593 tol:0.00
> >n:0 a:3 p:3 q:1 approx:3.00
> >n:1 a:7 p:22 q:7 approx:3.142857
> >n:2 a:15 p:333 q:106 approx:3.141509
> >n:3 a:1 p:355 q:113 approx:3.141593
> >n:4 a:292 p:103993 q:33102 approx:3.141593
> >(103993, 33102)
> >
> >i.e., the first 5 approximations it came up with were 3/1, 22/7, 333/106,
> >355/113 and a whopping 103993/33102.
> >
> >For the 0.36 example you used earlier:
> >
> > >>> print cf(0.36, .01, Trace= True)
> >x:0.36 tol:0.01
> >n:0 a:0 p:0 q:1 approx:0.00
> >n:1 a:2 p:1 q:2 approx:0.50
> >n:2 a:1 p:1 q:3 approx:0.33
> >n:3 a:3 p:4 q:11 approx:0.363636
> >(4, 11)
> > >>>
> >
> >it went right from 1/3 to 4/11 (0.363636), skipping the 3/8 (0.375) option
> >from the farey series.
> >
> >But this continued fraction algorithm is ill-suited to answer the question
> >"what's the closest fraction with a denominator < N", because it doesn't
> >try to find that, it jumps further ahead with each iteration.
> >
> >Anyway, I thought you might find it interesting based on our discussion.
>
>Terry,
>
>Well, I have to admit I don't understand your code at all. But I see it works.
>
>I modified one of my functions of frac.py, and came up with
>
>===
>from __future__ import division
>import time, psyco
>
>psyco.full()
>
>def d(number):
>  import decimal
>  decimal.getcontext().prec = 16
>  return decimal.Decimal(str(number))
>
>def bestFracForMinimumError(decimal, minimumError):
>  denom = 0
>  smallestError = 10
>  count = 0
>  while True:
>  denom += 1
>  num = int(round(d(decimal) * d(denom)))
>  error = absd(num)) / d(denom)) - d(decimal)) /
>d(decimal)) * d(100)
>  if d(error) <= d(smallestError):
>  count += 1
>  smallestError = d(error)
>  q = d(num)/d(denom)
>  print "%d/%d = %s has error of %s per cent" % (num,
>denom, q, smallestError)
>  if d(smallestError) <= d(minimumError):
>  print "count is", count
>  break
>
>=

I just realized that I had used "<=" in my code where I should have 
used "<".  So after make these changes, the results also changed, of 
course. See them at 
http://www.rcblue.com/Python/PartOfReplyToTerryOnTutorList-2.txt

Dick

>You can see the results of both
>bestFracForMinimumError(3.14159265357989, 0.0002)
>
>(BTW your pi is a bit off but I used yours, instead of math.pi, which
>is 3.1415926535897931 . Also, I needed 0.0002 in order to produce
>your 103993/33102)
>
>and
>
>bestFracForMinimumError(.36, .01)
>
>at 



___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 11:21 PM 1/4/2007, Terry Carroll wrote:
>On Wed, 3 Jan 2007, Dick Moores wrote:
>
> > Be that as it may, farey() is an amazing program.
>
>Not to beat this subject to death, but the comment at the bottom of
>http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317 about
>continued fractions piqued my interest.  I'm no mathematician, but I
>encountered continued fractions a long time ago and was fascinated by
>them.  So I read the URL pointed to,
>http://mathworld.wolfram.com/ContinuedFraction.html , and came up with the
>following:
>
>#
>
>def cf(x, tol=0.0001, Trace=False):
> """
> Calculate rational approximation of x to within tolerance of tol;
> returns a tuple consisting of numerator and denominator p/q
> Trace=True causes iterated results to be shown
> """
> a, r, p, q = [], [], [], []
> Done = False
> n = 0
> if Trace: print "x:%f tol:%f" % (x, tol)
> while not Done:
> a.append(None)
> r.append(None)
> p.append(None)
> q.append(None)
> if n == 0: r[n] = x
> else: r[n] = 1/(r[n-1]-a[n-1])
> a[n] = int(r[n])
> if n == 0:
> p[n] = a[0]
> q[n] = 1
> elif n ==1:
> p[n] = a[n]*p[n-1] + 1
> q[n] = a[n]
> else:
> p[n] = a[n]*p[n-1] + p[n-2]
> q[n] = a[n]*q[n-1] + q[n-2]
> if Trace:
> print "n:%d a:%d p:%d q:%d approx:%f" % \
>   (n, a[n], p[n], q[n], float(p[n])/q[n])
> if abs(float(p[n])/q[n] - x) < tol:
> Done = True
> num = p[n]; denom = q[n]
> n += 1
> return (num, denom)
>
>#
>
>Here's a result for pi:
>
> >>> print cf(3.14159265357989,0.001, Trace=True)
>x:3.141593 tol:0.00
>n:0 a:3 p:3 q:1 approx:3.00
>n:1 a:7 p:22 q:7 approx:3.142857
>n:2 a:15 p:333 q:106 approx:3.141509
>n:3 a:1 p:355 q:113 approx:3.141593
>n:4 a:292 p:103993 q:33102 approx:3.141593
>(103993, 33102)
>
>i.e., the first 5 approximations it came up with were 3/1, 22/7, 333/106,
>355/113 and a whopping 103993/33102.
>
>For the 0.36 example you used earlier:
>
> >>> print cf(0.36, .01, Trace= True)
>x:0.36 tol:0.01
>n:0 a:0 p:0 q:1 approx:0.00
>n:1 a:2 p:1 q:2 approx:0.50
>n:2 a:1 p:1 q:3 approx:0.33
>n:3 a:3 p:4 q:11 approx:0.363636
>(4, 11)
> >>>
>
>it went right from 1/3 to 4/11 (0.363636), skipping the 3/8 (0.375) option
>from the farey series.
>
>But this continued fraction algorithm is ill-suited to answer the question
>"what's the closest fraction with a denominator < N", because it doesn't
>try to find that, it jumps further ahead with each iteration.
>
>Anyway, I thought you might find it interesting based on our discussion.

Terry,

Well, I have to admit I don't understand your code at all. But I see it works.

I modified one of my functions of frac.py, and came up with

===
from __future__ import division
import time, psyco

psyco.full()

def d(number):
 import decimal
 decimal.getcontext().prec = 16
 return decimal.Decimal(str(number))

def bestFracForMinimumError(decimal, minimumError):
 denom = 0
 smallestError = 10
 count = 0
 while True:
 denom += 1
 num = int(round(d(decimal) * d(denom)))
 error = absd(num)) / d(denom)) - d(decimal)) / 
d(decimal)) * d(100)
 if d(error) <= d(smallestError):
 count += 1
 smallestError = d(error)
 q = d(num)/d(denom)
 print "%d/%d = %s has error of %s per cent" % (num, 
denom, q, smallestError)
 if d(smallestError) <= d(minimumError):
 print "count is", count
 break

=

You can see the results of both
bestFracForMinimumError(3.14159265357989, 0.0002)

(BTW your pi is a bit off but I used yours, instead of math.pi, which 
is 3.1415926535897931 . Also, I needed 0.0002 in order to produce 
your 103993/33102)

and

bestFracForMinimumError(.36, .01)

at 

Dick


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Wed, 3 Jan 2007, Dick Moores wrote:

> Be that as it may, farey() is an amazing program. 

Not to beat this subject to death, but the comment at the bottom of 
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317 about 
continued fractions piqued my interest.  I'm no mathematician, but I 
encountered continued fractions a long time ago and was fascinated by 
them.  So I read the URL pointed to, 
http://mathworld.wolfram.com/ContinuedFraction.html , and came up with the 
following:

#

def cf(x, tol=0.0001, Trace=False):
"""
Calculate rational approximation of x to within tolerance of tol;
returns a tuple consisting of numerator and denominator p/q
Trace=True causes iterated results to be shown
"""
a, r, p, q = [], [], [], []
Done = False
n = 0
if Trace: print "x:%f tol:%f" % (x, tol)
while not Done:
a.append(None)
r.append(None)
p.append(None)
q.append(None)
if n == 0: r[n] = x
else: r[n] = 1/(r[n-1]-a[n-1])
a[n] = int(r[n])
if n == 0:
p[n] = a[0]
q[n] = 1
elif n ==1:
p[n] = a[n]*p[n-1] + 1
q[n] = a[n]
else:
p[n] = a[n]*p[n-1] + p[n-2]
q[n] = a[n]*q[n-1] + q[n-2]
if Trace:
print "n:%d a:%d p:%d q:%d approx:%f" % \
  (n, a[n], p[n], q[n], float(p[n])/q[n])
if abs(float(p[n])/q[n] - x) < tol:
Done = True
num = p[n]; denom = q[n]
n += 1
return (num, denom)

#

Here's a result for pi:

>>> print cf(3.14159265357989,0.001, Trace=True)
x:3.141593 tol:0.00
n:0 a:3 p:3 q:1 approx:3.00
n:1 a:7 p:22 q:7 approx:3.142857
n:2 a:15 p:333 q:106 approx:3.141509
n:3 a:1 p:355 q:113 approx:3.141593
n:4 a:292 p:103993 q:33102 approx:3.141593
(103993, 33102)

i.e., the first 5 approximations it came up with were 3/1, 22/7, 333/106, 
355/113 and a whopping 103993/33102.

For the 0.36 example you used earlier:

>>> print cf(0.36, .01, Trace= True)
x:0.36 tol:0.01
n:0 a:0 p:0 q:1 approx:0.00
n:1 a:2 p:1 q:2 approx:0.50
n:2 a:1 p:1 q:3 approx:0.33
n:3 a:3 p:4 q:11 approx:0.363636
(4, 11)
>>>

it went right from 1/3 to 4/11 (0.363636), skipping the 3/8 (0.375) option 
from the farey series.  

But this continued fraction algorithm is ill-suited to answer the question
"what's the closest fraction with a denominator < N", because it doesn't
try to find that, it jumps further ahead with each iteration.

Anyway, I thought you might find it interesting based on our discussion.

(Yeah, I know, I didn't choose the formats well on those print 
statements.)

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 08:20 AM 1/4/2007, Terry Carroll wrote:
>On Wed, 3 Jan 2007, Dick Moores wrote:
>
> > Terry, I just noticed that farey(0.36, 10) returns (1, 3), a pretty
> > big miss, IMO. The correct fraction with smallest error and maximum
> > denominator of 10 is 3/8, which I'm proud to say my klunky frac.py
> > () produces.
>
>Dick, good catch.  Consider adding a comment to the Cookbook page pointing
>out that the method sometimes does not produce the incorrect answer,
>including for that input.

Done.

Dick



___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Wed, 3 Jan 2007, Dick Moores wrote:

> Terry, I just noticed that farey(0.36, 10) returns (1, 3), a pretty 
> big miss, IMO. The correct fraction with smallest error and maximum 
> denominator of 10 is 3/8, which I'm proud to say my klunky frac.py 
> () produces.

Dick, good catch.  Consider adding a comment to the Cookbook page pointing
out that the method sometimes does not produce the incorrect answer, 
including for that input.

I note that the web page for the recipe says only that its produces "a 
rational," not necessarily the closest rational, but the description in 
the book says it's the closest:

  You have a number v (of almost any type) and need to find a rational
  number (in reduced form) that is as close to v as possible but with a 
  denominator no larger than a prescribed value. (Recipe 18.13)


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 01:17 PM 1/2/2007, Terry Carroll wrote:
>On Mon, 1 Jan 2007, Dick Moores wrote:
>
> > bestFracForMinimumError() is only a small part of a program I wrote
> > long ago, called frac.py
>
>Dick, if your goal is to have a routine to get the fraction with the least
>possible error (as opposed to learing how to use Decimal), have a look at
>this recipe from the Python Cookbook:
>
>  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317

Terry, I just noticed that farey(0.36, 10) returns (1, 3), a pretty 
big miss, IMO. The correct fraction with smallest error and maximum 
denominator of 10 is 3/8, which I'm proud to say my klunky frac.py 
() produces.

Be that as it may, farey() is an amazing program. I appreciate the 
fast replies from you and Danny to my plea for explication. I'm still 
working through them.

Dick





___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Danny Yoo]

>> Dick, if your goal is to have a routine to get the fraction with the least
>> possible error (as opposed to learing how to use Decimal), have a look at
>> this recipe from the Python Cookbook:
>>
>>  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317
>
> Terry, that is truly ingenious. Is there an explication anywhere of
> exactly how it works?


Hi Dick,

On a first glance, it looks like it's doing a binary search on the Farey 
series on some order.  (Actually, it's not, but we'll get to that later in 
this post.) Farey sequences have an entry on Wikipedia:

 http://en.wikipedia.org/wiki/Farey_number

Look at the terms of F_8, for example.  If you take two points with some 
other point between them, say 1/5 and 2/7, notice that:

   1 3   2
   -  <  --  <   -
   5 12  7

If we just do the funny thing by adding the numerators and denominators 
--- by taking the "mediant" --- we end up with another term in the Farey 
sequence.  This is very cute.


Wait.  Actually, what the algorithm is doing doesn't appear to be 
searching through a particular Farey sequence of order n.  Instead, what 
it's doing is much simpler: it appears to be just taking advanatage of the 
in-betweenness property of "mediants".

 http://en.wikipedia.org/wiki/Mediant_%28mathematics%29

Oh well, that still works.  The algorithm seems to be misnamed, though: I 
think it should really be described as "inexact to rational via mediant 
approximation."


The trick that the algorithm is really a binary-search, using the 
definition of "mediant" to find "midpoints".  Whenever we see something 
like:


while we haven't found the answer:
 We know that the answer's somewhere between the lower and upper
 bounds.  (precondition)

 midpoint = some calculation combining the lower and upper bounds

 if the midpoint is too big:
 move the upper bound
 elif the midpoint is too small:
 move the lower bound
 else
 we've found the answer

 At the end of this, we guarantee our answer's still between
 the lower and upper bounds.  (postcondition)


then we should suspect a binary search.  In the case of the algorithm in 
the Cookbook, we can see that it follows this structure very closely.

Concretely, when they say:

 if v * mediant[1] > mediant[0]: ...

we can do simple equational reasoning to see that this is really 
saying:

 if v > (mediant[0] / mediant[1]): ...

aka: "if the value we're looking for is bigger than the mediant, move the 
lower bound up."  The mediant is being represented by a 2-tuple 
(numerator, denominator), so with that, you should be able to read the 
case analysis off more easily.


Best of wishes!
___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Wed, 3 Jan 2007, Dick Moores wrote:

> At 01:17 PM 1/2/2007, Terry Carroll wrote:
>
> >  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317
> 
> Terry, that is truly ingenious. Is there an explication anywhere of 
> exactly how it works?

There is in the printed copy of the Python Cookbook.  I happened to have 
read it last week, which is the only reason I knew about it.  I don't have 
the book handy, but here's what I remember of it.

It's based on the concept of the Farey sequence of fractions.  I may 
explain it imperfectly, but it's basically the sequence of fractions with 
integer denominators and numerators, each reduced to the lowest common 
denominator, and then arranged in numerically ascending order.  The 
sequence has an order N, where N is the highest integer used in the 
denominator.

It's easier to show an example.  Here's the Farey sequence for order 4; 
I'll include a decimal approximation to show the basis for the ordering:

 0/1,   1/4,   1/3,   1/2,   2/3,   3/4,   1/1

 0.00.25   0.33  0.500.67   0.75   1.00


It's been observed (by a guy named Farey, hence the name) that, for any
subsequence of three consecutive fractions in the sequence, if expressed
as A/B, C/D, E/F; then C/D = (A+E)/(B+F).

For example, in the above seqence, take the three-fraction subsequence 
{1/2, 2/3, 3/4}.  In this,  A=1, B=2, C=2, D=3, E=3, F=4, so:

  (A+E)/(B+F) = (1+3)/(2+4) = 4/6 = 2/3 = C/D.

The inner fraction (C/D) is called the mediant.

If I understand the Python algorithm correctly, it takes advantage of this
process by essentially constructing a Farey sequence of order 1 (which is
essentially {0/1, 1/1}, or {0, 1}), and then calculating the mediant
between those two points, thereby constructing a subsequence of order 2;
then another mediant between that mediant and one of its neighbors (which
neighbor is chose by considering whether the decimal fraction you seek to
approximate is greater than or less than the mediant), and then doing this
continuously, calculating subsquences of orders 3, 4, etc, until it
reaches the desired precision.

I just checked, and (big surprise) Wikipedia has an article on the Farey
sequence:

http://en.wikipedia.org/wiki/Farey_number

There's apparently more to it than I described above.  They give the same 
formula I use, but where I use A, B, C, D, E, F, they use A, B, P, Q, C, 
D, respectively.

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 01:17 PM 1/2/2007, Terry Carroll wrote:
>On Mon, 1 Jan 2007, Dick Moores wrote:
>
> > bestFracForMinimumError() is only a small part of a program I wrote
> > long ago, called frac.py
>
>Dick, if your goal is to have a routine to get the fraction with the least
>possible error (as opposed to learing how to use Decimal), have a look at
>this recipe from the Python Cookbook:
>
>  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317

Terry, that is truly ingenious. Is there an explication anywhere of 
exactly how it works?

Dick

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 10:39 AM 1/2/2007, Kent Johnson wrote:
>Dick Moores wrote:
>>from decimal import Decimal as D
>>def bestFracForMinimumError(decimal, minimumError):
>>  denom = 0
>>  while True:
>>  denom += 1
>>  num = round(D(str(decimal)) * D(str(denom)))
>>  error = abs(str((str(D(num) / D(str(denom))) -
>This looks backwards^^

Thanks for catching that.


>Don't you need D(str(num)) ? Then converting it back to a str before 
>you call abs will not work.
>
>Your old approach of
>def D(num):
>   return Decimal(str(num))
>
>would probably make for more readable code and fewer errors.

Yes, I went back to it.

>I'm not sure this approach will work, though, if you are trying to 
>get higher precision, I would think you would have to do all the 
>calculations in Decimal, not going to floats for num. I admit I 
>haven't thought it through, though. I think you can use 
>Decimal.quantize() instead of round().

Thanks very much, Kent. You got me back on track.

Dick



___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Tue, 2 Jan 2007, Terry Carroll wrote:

> Dick, if your goal is to have a routine to get the fraction with the least 
> possible error (as opposed to learing how to use Decimal), have a look at 
> this recipe from the Python Cookbook:
> 
>  http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317

BTW, to get the higher precision in Decimal:

import decimal

def farey(v, lim):
  '''Named after James Farey, an English surveyor.
  No error checking on args -- lim = max denominator,
  results are (numerator, denominator), (1,0) is infinity
  '''
  if v < 0:
n,d = farey(-v, lim)
return int(-n),int(d)
  z = lim-lim   # get 0 of right type for denominator
  lower, upper = (z,z+1), (z+1,z)
  while 1:
mediant = (lower[0] + upper[0]), (lower[1]+upper[1])
if v * mediant[1] > mediant[0]:
if lim < mediant[1]: return (int(upper[0]), int(upper[1]))
lower = mediant
elif v * mediant[1] == mediant[0]:
if lim >= mediant[1]: return (int(mediant[0]), int(mediant[1]))
if lower[1] < upper[1]: return (int(lower[0]), int(lower[1]))
return (int(upper[0]), int(upper[1]))
else:
if lim < mediant[1]: return lower
upper = mediant

dec = decimal.Decimal("0.345765988765560057657654654")
me = decimal.Decimal("0.01")
max_denom = 1/me

(num, denom) = farey(dec,max_denom)
error = abs(decimal.Decimal(str(float(num)/denom)) - dec) *100
print "%s = %d/%d with %s per cent error" % (dec, num, denom, error)


Which gives:

0.345765988765560057657654654 = 878844001/2541730620 with 
4.3994234234534600E-11 per cent error

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

On Mon, 1 Jan 2007, Dick Moores wrote:

> bestFracForMinimumError() is only a small part of a program I wrote 
> long ago, called frac.py 

Dick, if your goal is to have a routine to get the fraction with the least 
possible error (as opposed to learing how to use Decimal), have a look at 
this recipe from the Python Cookbook:

 http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317

> dec = .345765988765560057657654654
> me = .01
> num, denom, error = bestFracForMinimumError(dec, me)
> print "%f = %d/%d with %f per cent error" % (dec, num, denom, error)

The parameters are different from yours, but as I understand it, your 
equivalent is:

>>> dec = .345765988765560057657654654
>>> me = .01
>>> max_denom = 1/me
>>> (num, denom) = farey(dec,max_denom)
>>> error = abs((num/denom) - dec) *100
>>> print "%f = %d/%d with %f per cent error" % (dec, num, denom, error)
0.345766 = 40258524/116432863 with 0.00 per cent error
>>>

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Kent Johnson]

Dick Moores wrote:
> from decimal import Decimal as D
> 
> def bestFracForMinimumError(decimal, minimumError):
>  denom = 0
>  while True:
>  denom += 1
>  num = round(D(str(decimal)) * D(str(denom)))
>  error = abs(str((str(D(num) / D(str(denom))) - 
This looks backwards^^

Don't you need D(str(num)) ? Then converting it back to a str before you 
call abs will not work.

Your old approach of
def D(num):
   return Decimal(str(num))

would probably make for more readable code and fewer errors.

I'm not sure this approach will work, though, if you are trying to get 
higher precision, I would think you would have to do all the 
calculations in Decimal, not going to floats for num. I admit I haven't 
thought it through, though. I think you can use Decimal.quantize() 
instead of round().

HTH
Kent

> D(str(decimal))) / str(D(str(decimal)) * d("100"
>  if error <= D(minimumError):
>  break
>  return int(num), D(denom), error
> 
> dec = D(".34576598876876867756765765")
> 
> me = D(".0001")
> 
> print bestFracForMinimumError(dec, me)
> 
> 
> Traceback (most recent call last):
>File "fracSimple2-c.py", line 17, in 
>  print bestFracForMinimumError(dec, me)
>File "fracSimple2-c.py", line 8, in bestFracForMinimumError
>  error = abs(str((str(D(num) / D(str(denom))) - D(str(decimal))) 
> / str(D(str(
> decimal)) * d("100"
>File "E:\Python25\lib\decimal.py", line 578, in __new__
>  "First convert the float to a string")
> TypeError: Cannot convert float to Decimal.  First convert the float 
> to a string
> 
> I don't understand this TypeError. Seems to me that I've converted 
> EVERYTHING in that line 8 to a string.
> 
> Dick
> 
> ___
> Tutor maillist  -  Tutor@python.org
> http://mail.python.org/mailman/listinfo/tutor
> 
> 


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

from decimal import Decimal as D

def bestFracForMinimumError(decimal, minimumError):
 denom = 0
 while True:
 denom += 1
 num = round(D(str(decimal)) * D(str(denom)))
 error = abs(str((str(D(num) / D(str(denom))) - 
D(str(decimal))) / str(D(str(decimal)) * d("100"
 if error <= D(minimumError):
 break
 return int(num), D(denom), error

dec = D(".34576598876876867756765765")

me = D(".0001")

print bestFracForMinimumError(dec, me)


Traceback (most recent call last):
   File "fracSimple2-c.py", line 17, in 
 print bestFracForMinimumError(dec, me)
   File "fracSimple2-c.py", line 8, in bestFracForMinimumError
 error = abs(str((str(D(num) / D(str(denom))) - D(str(decimal))) 
/ str(D(str(
decimal)) * d("100"
   File "E:\Python25\lib\decimal.py", line 578, in __new__
 "First convert the float to a string")
TypeError: Cannot convert float to Decimal.  First convert the float 
to a string

I don't understand this TypeError. Seems to me that I've converted 
EVERYTHING in that line 8 to a string.

Dick

___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Dick Moores]

At 05:37 AM 1/1/2007, Kent Johnson wrote:
>Dick Moores wrote:
>>bestFracForMinimumError() is only a small part of a program I wrote 
>>long ago, called frac.py 
>>(). I'm trying to 
>>rewrite it so as to get more precision by using the Decimal module, 
>>but am getting nowhere. Errors all over the place. Can some kind 
>>and knowledgeable soul show me how?
>
>I don't see where you are actually using the decimal module here, am 
>I missing something?

You're missing the explanation that I didn't write. ;-)

I pasted what I had before even starting to use Decimal, but left in 
d(), thinking that my helper could use it..

Dick


>It doesn't help that you cal the argument to your function 'decimal'.
>
>>dec = .345765988765560057657654654
>
>Your Python interpreter probably doesn't have enough precision to 
>represent this number directly. Here is what I get:
>
>In [1]: dec = .345765988765560057657654654
>
>In [2]: dec
>Out[2]: 0.34576598876556008
>
>Kent
>
>


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Re: [Tutor] Need help with rewriting script to use Decimal module

[Kent Johnson]

Dick Moores wrote:
> bestFracForMinimumError() is only a small part of a program I wrote 
> long ago, called frac.py 
> (). I'm trying to rewrite 
> it so as to get more precision by using the Decimal module, but am 
> getting nowhere. Errors all over the place. Can some kind and 
> knowledgeable soul show me how?

I don't see where you are actually using the decimal module here, am I 
missing something?

It doesn't help that you cal the argument to your function 'decimal'.

> dec = .345765988765560057657654654

Your Python interpreter probably doesn't have enough precision to 
represent this number directly. Here is what I get:

In [1]: dec = .345765988765560057657654654

In [2]: dec
Out[2]: 0.34576598876556008

Kent


___
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor