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.

[Thomas_Sch]

(12-13-2013 05:09 AM)Joe Horn Wrote:  ...PDQ is a CAS program and must be created as such. ...

@Joe: Maybe it could be helpful to add some more words about entering a CAS program, since I didn't find hints about creating a CAS program in the manual(s).
My way (awkward?): First in CAS mode I created an empty program like

Code:

d2f:=(x)->x

then I could edit this by [Shift]+1 (Program), and enter your code. Works fine, besides a small typo ('END;;').

Are there other steps to create/enter a CAS program? Where did you find the operators like '/=' and '+=' ? The manuals are unhelpful ;-)

Many thanks for your programs!
Thomas

[Han]

(02-16-2014 02:01 PM)Thomas_Sch Wrote:  

(12-13-2013 05:09 AM)Joe Horn Wrote:  ...PDQ is a CAS program and must be created as such. ...

@Joe: Maybe it could be helpful to add some more words about entering a CAS program, since I didn't find hints about creating a CAS program in the manual(s).
My way (awkward?): First in CAS mode I created an empty program like

Code:

d2f:=(x)->x

then I could edit this by [Shift]+1 (Program), and enter your code. Works fine, besides a small typo ('END;;').

Are there other steps to create/enter a CAS program? Where did you find the operators like '/=' and '+=' ? The manuals are unhelpful ;-)

Many thanks for your programs!
Thomas

The CAS is xcas/giac. You can find relevant documentation here: http://www-fourier.ujf-grenoble.fr/~parisse/giac.html

As for creating CAS programs, your technique is essentially the only method on this current firmware.

[patrice]

(12-13-2013 05:09 AM)Joe Horn Wrote:  Example #1: What is the best fraction that's equal to pi plus or minus 1/800? Answer: 179/57, which is the same answer as given by Patrice's FareyDelta program. This is a difficult problem, because 179/57 is not a principal convergent of pi.

Not a big deal for a program that use Farey series because the Farey series try all the intermediate fractions until it find one within tolerance.
The counterpart is that Farey is slow to converge for values near an integer.

(12-13-2013 05:09 AM)Joe Horn Wrote:  Example #2: What is the best fraction that's equal to pi plus or minus 10^-14? FareyDelta can't handle that problem, because it is limited by the 12-digit accuracy of Home math. PDQ gets the right answer, though: 58466453/18610450.

That's Home math limitation.

[Joe Horn]

(03-26-2014 10:10 AM)patrice Wrote:  

(12-13-2013 05:09 AM)Joe Horn Wrote:  Example #1: What is the best fraction that's equal to pi plus or minus 1/800? Answer: 179/57, which is the same answer as given by Patrice's FareyDelta program. This is a difficult problem, because 179/57 is not a principal convergent of pi.

Not a big deal for a program that use Farey series because the Farey series try all the intermediate fractions until it find one within tolerance.

Yes, that was my point: the PDQ algorithm and Farey fractions find it, but the standard "continued fraction algorithm" (as used by everybody but HP) misses all intermediate convergents.

(03-26-2014 10:10 AM)patrice Wrote:  The counterpart is that Farey is slow to converge for values near an integer.

This is one of PDQ's unique strengths; it doesn't spend any time searching through intermediate convergents. As soon as it determines that the best fraction MIGHT lie between two principal convergents, it performs a single calculated jump to the best one. This is important when there are thousands of intermediate convergents in a row (which happens whenever there are very large partial quotients in the continued fraction expansion).

(03-26-2014 10:10 AM)patrice Wrote:  

(12-13-2013 05:09 AM)Joe Horn Wrote:  Example #2: What is the best fraction that's equal to pi plus or minus 10^-14? FareyDelta can't handle that problem, because it is limited by the 12-digit accuracy of Home math. PDQ gets the right answer, though: 58466453/18610450.

That's Home math limitation.

It would be very helpful (and fun!) to have a CAS version of your Farey fraction programs, so that they could use "infinite precision" integers to overcome Home's 12-digit limitation. Then we could really compare the performance of PDQ vs. Farey fractions on challenging inputs.

[Joe Horn]

Both programs in the original posting 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)

[dbbotkin]

Thank you Joe, for illustrating the format for a CAS function--I had searched in vain--and also for the excellent examples!

[Joe Horn]

Both listings in the original post have been updated to reflect the greatly improved handling of CAS programs in Prime rev 6940. CAS programs are now created the same way as regular programs; just wrap them in #cas / #end (see listings above). They will now appear in the Program Catalog, and not lose comments, and be handled by the Connectivity Kit correctly. Sweet!

[patrice]

Hi Joe,

You should also add the #pragma to ensure your programs compile no matter the user settings for decimal separator.

[Joe Horn]

(12-06-2014 10:47 AM)patrice Wrote:  Hi Joe,

You should also add the #pragma to ensure your programs compile no matter the user settings for decimal separator.

Does that matter when the program doesn't contain any decimal separators? Both D2F and PDQ use only integers.

[patrice]

(12-06-2014 02:38 PM)Joe Horn Wrote:  

(12-06-2014 10:47 AM)patrice Wrote:  Hi Joe,

You should also add the #pragma to ensure your programs compile no matter the user settings for decimal separator.

Does that matter when the program doesn't contain any decimal separators? Both D2F and PDQ use only integers.

YES, when someone use comma as decimal separator, the semicolon is the parameter separator.

[Joe Horn]

(12-06-2014 05:06 PM)patrice Wrote:  YES, when someone use comma as decimal separator, the semicolon is the parameter separator.

Hmmm... adding this line to the beginning of the D2F program:
#pragma mode( separator(.,;) )
causes a syntax error. What exactly must be added, and where? Thanx in advance.

[patrice]

(12-07-2014 04:30 AM)Joe Horn Wrote:  

(12-06-2014 05:06 PM)patrice Wrote:  YES, when someone use comma as decimal separator, the semicolon is the parameter separator.

Hmmm... adding this line to the beginning of the D2F program:
#pragma mode( separator(.,;) )
causes a syntax error. What exactly must be added, and where? Thanx in advance.

Code:

#pragma mode( separator(.,;) integer(h64) )

As first line of code
As done with Menu > insert pragma while in program editor
... unless an incompatibility have been introduced lately.

[Joe Horn]

(12-07-2014 08:51 AM)patrice Wrote:  

Code:

#pragma mode( separator(.,;) integer(h64) )

As first line of code
As done with Menu > insert pragma while in program editor
... unless an incompatibility have been introduced lately.

Bummer... Apparently an incompatibility has in fact been introduced, because what you are recommending always causes a syntax error. Please try it yourself and see what I mean, using either the D2F or PDQ program above. It seems that #pragma and #cas do not play well together in rev 6940.

[Joe Horn]

(05-01-2015 01:54 AM)compsystems Wrote:  Great program! which means the acronym PDQ?

Thanks. Two meanings: (1) Since the math literature often uses P/Q when talking about continued fractions, I wanted to call the original HP48 function 'P/Q' but the parser didn't like that name, so I called it 'PDQ' with "D" standing for "Divide". (2) It's a common English abbreviation for "Pretty Darn Quick" (or equivalent words for "D"), and is used as a synonym for "immediately" or "at once", as in "Get that project started PDQ!". It's ASAP on steroids. Since the PDQ Algorithm is so fast, it seemed an appropriate name. No relation to PDQ Bach, unfortunately.

(05-01-2015 01:54 AM)compsystems Wrote:  A question, should be called arbitrary presicion, there tolerance. The tolerance parameter makes it finite.

Counterintuitively, in English "infinite precision" and "arbitrary precision" are more or less synonyms. That's why "Infinite Precision" in Wikipedia defaults to the "Arbitrary Precision" page: http://en.wikipedia.org/w/index.php?titl...edirect=no

I prefer "infinite precision" in this case because the user-specified precision is not limited by the algorithm or by the machine's floating-point bandwidth, but only by available memory. "Infinite precision" is never really infinite of course; it's just what it's called.

(05-01-2015 01:54 AM)compsystems Wrote:  Please can port this code into the HP50G

Done, long ago, in System RPL for speed. It was the subject of my HHC 2003 presentation (on the 49G at the time, I think). I'll post the 50g program to the "General Software Library" ASAP, if not PDQ.

EDIT: Done. The 50g version is HERE.

[pier4r]

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-790...69282.html The a b/c key
http://www.hpmuseum.org/forum/thread-377...l#pid34316 (50G) PDQ Algorithm in SRPL and URPL
https://groups.google.com/forum/#!topic/...discussion PDQ Unleashed--update
https://groups.google.com/forum/#!msg/co...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.s...discussion 48/49: Best Fraction Challenge
https://groups.google.com/d/topic/comp.s...discussion worst fraction
https://groups.google.com/d/topic/comp.s...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)

[Luigi Vampa]

@Pier4r, thanks for sharing that comprehensive summary!

[Joe Horn]

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

[pier4r]

Discovered by chance another bit

http://www.hpmuseum.org/forum/thread-202...l#pid17844

http://www.hpmuseum.org/forum/thread-202...l#pid17846

[smartin]

(12-13-2013 05:09 AM)Joe Horn Wrote:  PDQ Algorithm for HP Prime, by Joe Horn

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}\)

Finally inspired to try out PDQ on the Prime, but I could not get all the examples to work out. I'm using PDQ from hpcalc.org (https://www.hpcalc.org/details/7477) on a Prime with CAS ver 1.4.9 and ROM 2.1.14181.

Example #1: works as advertised
but,
Example #2: pdq(\(\pi\),14) = \(\dfrac{111513555}{35495867}\)

Example #3: pdq(\(\pi,\dfrac{13131}{10^{440}}\)) = \(\dfrac{27633741218861}{8796093022208}\)

As with trying any new program my first thought is operator error (or some setting on the Prime).

Another related question:
For example #3, even though it appears that the return fraction accurately represents pi to about 436 decimal places, is it true that one is not actually getting pi to that accuracy, just a lot of zeros after the precision of pi within the Prime is reached (12 digits)?

Thanks,
Steve