Re: What Algorithm(s) Do(es) 'Rationalize[]' Use?
- To: mathgroup@smc.vnet.net
- Subject: [mg12427] Re: What Algorithm(s) Do(es) 'Rationalize[]' Use?
- From: Daniel Lichtblau <danl@wolfram.com>
- Date: Thu, 14 May 1998 11:15:36 -0400
Benjamin Tubb wrote: > > What algorithm(s) do(es) the function 'Rationalize[]' use? > > ------------- > Benjamin Tubb > AIM: brtubb > brtubb@cybertron.com > http://home.cybertron.com/~brtubb [This reply was written in part by Mark Sofroniou. I hope the merging of two responses does not cause confusion. dl] Rationalize uses the following algorithm. Given an inexact number x Rationalize computes the partial quotients in the continued fraction expansion of the number. At each step of the approximation Rationalize then finds a corresponding continued fraction convergent, which is a rational approximation p/q to x. Continued fraction convergents have several important properties: 1) p and q are relatively prime, so that p/q is expressed in lowest terms. 2) The approximant p/q is the closest to x amongst all approximants with denominator <= q. When you use the form Rationalize[x, tol] then Rationalize stops when the convergent p/q satisfies |x-p/q| < tol, or the absolute error is within the specified tolerance. When tol is 0, the most accurate approximant p/q, as governed by the accuracy of x, is chosen. When you use the form Rationalize[x] then Rationalize finds the approximant q with the property that |p/q - x| < c/q^2, where c is chosen to be 10^-4, provided it exists. Here we show explicitly what happens when we rationalize Pi. In[11]:= Rationalize[Pi, 10^-6] 355 Out[11]= --- 113 We check that the result is within 10^(-6) of Pi: In[12]:= N[Pi]-% < 10^(-6) Out[12]= True Let's look at the continued fraction representation of Pi. This uses functionality from our development version. In[13]:= ContinuedFraction[Pi,5] Out[13]= {3, 7, 15, 1, 292} So the continued fraction approximation, to first four denominators, is 3 + 1/(7 + 1/(15 + 1/(1 + 1/292)))) By trial and error (or recognizing that 292 is relatively large and so constitutes a small error term) we find that a sufficiently close approximation is In[14]:= 3 + 1/(7+1/(15+1)) 355 Out[14]= --- 113 If we give a slightly smaller tolerance, we'll need more terms. In[15]:= Rationalize[Pi, 10^(-7)] 104348 Out[15]= ------ 33215 In[16]:= 3 + 1/(7+1/(15+1/(1+1/292))) 103993 Out[16]= ------ 33102 Not quite the same. Using N we learn that either is a sufficiently close approximation. ContinuedFraction tells us we actually computed one more term than needed. This might be a small bug in the implementation. In[22]:= ContinuedFraction[N[Pi], 6] Out[22]= {3, 7, 15, 1, 292, 1} In[23]:= ff = 3 + 1/(7+1/(15+1/(1+1/293))) 104348 Out[23]= ------ 33215 If you read the fine print of the documentation notes for Rationalize you will note a failure mode that arises when the denominator of the result gets too large. While this can happen, it apparently occurred for spurious reasons in the past. Alot of this has been improved for our next release. Here is some related functionality from our development version. It uses continued fraction approximations to get successive rational approximations to Pi. In[1]:= <<NumberTheory`ContinuedFractions`; In[2]:= Convergents[Pi,8] 22 333 355 103993 104348 208341 312689 Out[2]= {3, --, ---, ---, ------, ------, ------, ------} 7 106 113 33102 33215 66317 99532 In[3]:= N[Pi-%,20] Out[3]= {0.14159265358979323846, -0.0012644892673496186802, -7 > 0.000083219627529087519247, -2.6676418906242231237 10 , -10 -10 > 5.7789063439038188851 10 , -3.3162780624607255831 10 , -10 -11 > 1.2235653294218859793 10 , -2.9143384934856918131 10 } Daniel Lichtblau Wolfram Research