PDQ Algorithm: Infinite precision best fraction within tolerance

[Joe Horn]

PDQ Algorithm for HP Prime, by Joe Horn

PDQ finds best rational approximations, with infinite precision. This means it finds the two smallest integers whose ratio is equal to some target real number plus or minus some desired tolerance. In other words, it finds the simplest fraction in any given interval. Unlike other methods, it always finds the unique best answer, and uses the infinite precision of CAS long integers.

Example #1: What is the best fraction that's equal to \(\pi\pm\dfrac{1}{800}\)? Answer: \(\dfrac{179}{57}\), which is the same answer as given by Patrice's FareyDelta program. This is a difficult problem, because \(\dfrac{179}{57}\) is not a principal convergent of pi.

Example #2: What is the best fraction that's equal to \(\pi\pm\dfrac{1}{10^{14}}\)? Patrice's FareyDelta program can't handle that problem, because it is limited by the 12-digit accuracy of Home math. PDQ gets the right answer, though: \(\dfrac{58466453}{18610450}\).

Example #3: What is the best fraction for \(\pi\pm\dfrac{13131}{10^{440}}\)? PDQ returns the correct ratio of two huge integers (221 digits over 220 digits) in less than one second. (It takes the HP 50g over two minutes using System RPL). This is an extreme case, since the numbers are so huge, and once again the answer is not one of pi's principal convergents. Piece of cake for Prime+PDQ, which yields the unique correct answer:
\(\frac{197589170636854062408454380413327813798855733721902369198118555167226743​04730662906703593620215835931889230827416036013979330716090096564056017111952129​356153172850632330284830147063755110178945173800035059898203820427519}{628945864​16566610363809944329675177193853240546238323897010336439848926113959966464032961​05227342931817856832425836690143454522392566443183092738965162183777760340854050​065970057153850182984658932709234874407013579584688}\)

PDQ calls another program called D2F (Decimal To Fraction) which returns the EXACT internal CAS representation of any real number.

Example: d2f(.5) returns \(\dfrac{1}{2}\), but d2f(.6) returns \(\dfrac{168884986026393}{2^{48}}\), because that's how Prime's CAS stores 0.6 internally.

Here's the source code for D2F (be sure to name it "d2f" in lowercase):

Code:

#cas
d2f(x):=BEGIN  
 LOCAL a,p,s,t,j; 
  IF x==0 THEN return(0);  END ;  
  IF type(x)==DOM_RAT THEN RETURN(abs(x));  END ;  
  IF type(x)==DOM_INT THEN RETURN(abs(x));  END ;  
  x:=abs(approx(x));  
  p:=floor(logb(x,2));  
  x/=2^p;  
  t:=-p;  
  a:=IP(x);  
  FOR j FROM 1 TO 3 DO  
  x:=65536*FP(x); 
  a:=65536*a+IP(x); 
  t+=16; 
  END;;  
  RETURN(a/2^t);  
END;
#end

And here's the code for PDQ (be sure to name it "pdq" in lowercase):

Code:

#cas
pdq(j,t):=BEGIN  
 LOCAL n,n0,d,d0,c,p,s,t1,t2,t3,a,b; 
  IF t==0 THEN  
    err:=0; 
    ic:=0;
    RETURN(d2f(j))  END ;  
  IF FP(t) AND (ABS(t)>1) THEN RETURN("Illegal Tolerance");  END ;  
  IF FP(t)==0 THEN t:=1/10^(exact(ABS(t)));  END ;  
  j:=d2f(j);  
  n0:=numer(j);  
  d0:=denom(j);  
  t:=d2f(t);  
  a:=numer(t);  
  b:=denom(t);  
  n:=n0;  
  d:=d0;  
  c:=0;  
  p:=1;  
  s:=1;  
  REPEAT  
  t1:=c; 
  t2:=d; 
  t3:=irem(n,d); 
  c:=c*iquo(n,d)+p; 
  p:=t1; 
  n:=t2; 
  d:=t3; 
  s:=-s 
  UNTIL (b*d)<=(c*a*d0);  
  t1:=iquo(c*a*d0-b*d,p*a*d0+b*n);  
  t2:=(n*t1+d)*s;  
  t3:=c-p*t1;  
  err:=t2/(d0*t3);  
  ic:=t1;
  RETURN((t2+t3*n0)/(d0*t3));  
END;
#end

Instructions for PDQ:

Syntax: pdq(target, tolerance) where you want to find the best fraction for target +/- tolerance.

The target can be input as either a floating-point real number, or as a fraction (ratio of two integers). In the latter case, you may use as many digits as you wish, for extended precision.

The tolerance can also be a floating-point real number or a fraction, in which case it is used as-is. The tolerance can also be an integer n > 1, in which case it is converted to the exact ratio 1/10^n. This shortcut is handy for testing purposes.

PDQ not only returns its result, but also stores the exact error in 'err'. The output of PDQ, minus err, exactly equals PDQ's input target.

If the global variable 'ic' is non-zero when PDQ exits, then it is the number of intermediate convergents between the result and the next convergent. Example: pdq(pi,3) returns \(\dfrac{201}{64}\), and stores 6 into 'ic' because there are 6 intermediate convergents between \(\dfrac{201}{64}\) (inclusive) and the next convergent \(\dfrac{333}{106}\). "Intermediate convergents" are also known as "semiconvergents" and "secondary convergents".

Here is a handy variable for PDQ experimentation.

PI500 (This fraction's decimal expansion has the same first 500 decimal places as pi. Remove all carriage returns):

Code:

27530008686166622188536681168621832641085194972343166639705257535483379211746872​24521381642611856603178539596529812288248903337810098177795117288227409717155741​87957420619251445521692137166819636595557228499775776315464391353285480273592327​83581546654/87630739315324660697093180818659895483560383602191807997668834668010919518358106​20316848615678705592355956178010775004819317000458201135712222333217497308440528​56473605433480720429471645084774186874516644432633187107318767999674836818181677​0785368793

Examples (performed in CAS, not Home, for perfect accuracy):

• pdq(PI500,14) --> \(\dfrac{58466453}{18610450}\) (example #2 above).

• pdq(pi,14) --> \(\dfrac{47627751}{15160384}\) (different, because pi is not equal to PI500).

• d2f(pi) --> \(\dfrac{27633741218861}{2^{43}}\) (that's what pi is equal to in CAS).

• pdq(pi,14)-err --> same as d2f(pi).

• time(pdq(PI500,13131/alog10(440))) --> approximately 0.6 seconds.

Note: The exact( ) function in CAS has three severe shortcomings: (1) it only finds answers which are principal convergents; (2) it only allows the tolerance (epsilon) to lie between 10^-6 and 10^-14; and (3) it sometimes returns incorrect answers (outside the specified tolerance). PDQ has none of these shortcomings. It always finds the unique best answer for any target and tolerance.

Edit: added the warnings to name "d2f" and "pdq" in lowercase.

Edit: Changed example #3 at the beginning to a non-convergent. Thanks to Rodger Rosenbaum for pointing out this great example.

Edit: Both programs above have been revised slightly to reflect two minor changes in the Prime update released today, rev 6030:

exact(iPart(x)) is now IP(x)
frac(x) is now FP(x)

Edit (5 December 2014): Updated for compatibility with the greatly improved CAS program handling of Prime rev 6940. Also added the 'ic' variable for tracking intermediate convergents.