PDQ Algorithm: Infinite precision best fraction within tolerance

+- HP Forums (http://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: PDQ Algorithm: Infinite precision best fraction within tolerance (/thread-61.html)

PDQ Algorithm: Infinite precision best fraction within tolerance - Joe Horn - 12-13-2013 12:09 AM

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.

RE: PDQ Algorithm: Infinite precision best fraction within tolerance - pier4r - 03-15-2018 10:20 AM

If someone like me want to understand from where the algorithm come from, because the source code is nice but it is not always self explaining, here are some possible helpful links that one has to combine together to get the picture. Fortunately on this topic a lot was written and it helps.

http://www.hpmuseum.org/forum/thread-4738.html (50G) Fast Continued Fractions
http://www.hpmuseum.org/forum/thread-7984.html 42.86% responded “Yes” on an opinion poll. How many people responded?
http://www.hpmuseum.org/forum/thread-7900-post-69282.html The a b/c key
http://www.hpmuseum.org/forum/thread-3778-post-34316.html#pid34316 (50G) PDQ Algorithm in SRPL and URPL
https://groups.google.com/forum/#!topic/comp.sys.hp48/8NTK1yj3tJY/discussion PDQ Unleashed--update
https://groups.google.com/forum/#!msg/comp.sys.hp48/oH9AAn-43uU/discussion Re: 2 questions on 50G [decimal -> fraction]
Quote:The PDQ saga, published piecemeal in a plethora of past postings
True! But at least one has some chances to find the parts as long as good search engines exists and the content is public. Some other contributions are pretty silent or ready to disappear.
https://groups.google.com/d/topic/comp.sys.hp48/MrwIHmW_xwg/discussion 48/49: Best Fraction Challenge
https://groups.google.com/d/topic/comp.sys.hp48/7Wh5y_pzNCQ/discussion worst fraction
https://groups.google.com/d/topic/comp.sys.hp48/-10z8FXqmeA/discussion [49G]: Revival of the Wishlist

and by chance went to perusing holy joe (I did not know was a minister in the army!) that is pretty unrelated content. http://holyjoe.org/hp/flashback.htm

If I miss some other parts of the PDQ saga that may give insights about the idea behind the algorithm, don't hesitate to post them! (Maybe Joe knows better than anyone else)

RE: PDQ Algorithm: Infinite precision best fraction within tolerance - Joe Horn - 03-15-2018 09:41 PM

Thanks, Pier4r! PDQ was first "officially" presented to HP and friends at HHC 2003. Later, as PDQ developed, that initial version came to be called PDQ1. It was limited to the decimal accuracy of its host machine, because it was implemented on the HP-48 which didn't have infinite-precision integers built-in. It was also slightly slowed down by needing to perform a binary search for intermediate convergents. Here is the paper from the HHC 2003 Proceedings which accompanied that presentation. Needless to say, much of what it says is now obsolete.

http://hhuc.us/2003/files/PDQ1.pdf

PDQ2 was a great improvement, utilizing infinite-precision integers. It was the first version of PDQ which did all its work in the integer domain, thus avoiding all roundoff errors, and allowing PDQ to generate ALL the best fractions for any input, not just the ones up to the host machine's floating-point accuracy. However, it still used a binary search method to find intermediate convergents, because I couldn't figure out a way to avoid that. Here are the associated papers from the HHC 2004 Proceedings. There are many overlaps of information here.

http://hhuc.us/2004/files/PDQ2.pdf
http://hhuc.us/2004/files/PDQ2b.pdf

PDQ3 is the final version of PDQ which avoids the binary search by calculating a single jump to the proper intermediate convergent. That improvement was added by Rodger Rosenbaum. The HP 49/50 code for it was optimized by Tony Hutchins. Other details of PDQ's development history can be found in Pier4r's list above. Here are the final HHC 2012 write-ups about PDQ3, which has come to be known simply as PDQ.

http://hhuc.us/2012/PDQ.TXT <-- describes PDQ3 (AKA PDQ)
http://hhuc.us/2012/PDQ.USR.TXT <-- PDQ in User RPL for study purposes.
http://hhuc.us/2012/PDQ.hp <-- PDQ in binary form for HP49/50