QPI: convert decimal to p/q, ln(p/q), p/q*pi, e^(p/q), or sqrt(p/q)

[Han]

Jan. 31, 2018: This program now exists in the HP Prime firmware.

This program takes a decimal value and returns a one of the following expressions that is a "close" rational approximation of the specified decimal value:

\[ \frac{p}{q}, \quad \frac{a}{b}\cdot\sqrt{\frac{p}{q}},
\quad \frac{p}{q}\cdot \pi, \quad e^{\frac{p}{q}}, \quad
\text{ or } \quad \ln \left(\frac{p}{q}\right) \]

Also works for complex numbers and lists of real/complex numbers. The main algorithm basically finds the continued fraction representation of a decimal and the process either self-terminates (in the case of a rational value) or terminates due to reaching the accuracy limit. For example,
\[ \frac{47}{13} = 3 + \frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{2}}}}} \]
To get the continued fraction representation, note that
\[ \frac{47}{13} \approx 3.61538461538 \]
This decimal is converted to a rational expression
\[ \frac{361538461538}{10^{11}} \]
which is converted into the continued fraction as follows:
\[ \frac{361538461538}{10^{11}} = 3 +
\frac{61538461538}{10^{11}} =
3 + \frac{1}{\frac{10^{11}}{61538461538}}
= 3 + \frac{1}{1+ \frac{38461538462}{61538461538}} = \dotsm \]
and so on. While computing the continued fraction, the algorithm simultaneously reduces the partial continued fraction into a rational value of the form \(\frac{p}{q}\). The other forms are merely variations in which the decimal value is either squared, divided by \(\pi\), etc.

Source code below, and also included in attached zip file below.

Code:

// QPI by Han Duong
// ported from QPI 4.3 for the HP48G/GX by Mika Heiskanen & Andre Schoorl

export qpiEXPLN:=100; // max denom for exp(p/q) or ln(p/q)
export qpiMAXINT:=2^20; // largest n allowed for sqrt(n)=a*sqrt(b)
export qpiDIGITS:=10; // controls accuracy (best results at 9 or 10)

qpi_approx();
qpi_asqrtb();
qpi_out();
qpi_outsqrt();
qpi_real();
qpi_complex();
qpi_list();
qpi_root();
qpi_pi();
qpi_ln();
qpi_exp();


EXPORT QPI(r)
BEGIN

  case
    if TYPE(r)==0 then qpi_real(r); end;
    if TYPE(r)==1 then RETURN(r); end;
    if TYPE(r)==3 then qpi_complex(r); end;
    if TYPE(r)==6 then qpi_list(r); end;
    if TYPE(r)==8 then QPI(approx(r)); end;
    DEFAULT msgbox("Object type: " + TYPE(r) + " not supported!");
  end;

END;


qpi_real(r)
BEGIN
  local frac;

  if r then 
    frac:=qpi_approx(r);

    if frac(2)<100 then
      qpi_out(frac);
    else
      qpi_root(r,frac);
    end;
  else
    RETURN(0);
  end;
END;


qpi_complex(c)
BEGIN
  local rpart, ipart;

  rpart:=STRING(qpi_real(RE(c)));
  ipart:=STRING(qpi_real(abs(IM(c))));

  if IM(c)>0 then
    expr("'" + rpart + "+" + ipart + "*'"); // bold i symbol
  else
    expr("'" + rpart + "-" + ipart + "*'"); // bold i symbol
  end;

END;


qpi_list(l)
BEGIN
  local i,n;

  n:=SIZE(l);
  for i from 1 to n do
    l(i):=QPI(l(i));
  end;

  RETURN(l);
  
END;


qpi_root(r,frac)
BEGIN
  local frac1;

  if r^2<500001 then

    frac1:=qpi_approx(r^2);
    if r<0 then frac1(1):=-frac1(1); end;
    frac1(3):=1;
    if (frac1(2)<1000) AND (frac1(2)<=frac(2)) then
      if frac1(2)<10 then
        qpi_out(frac1);
      else
        qpi_pi(r,frac1);
      end;
    else // sqrt denom not smaller
      qpi_pi(r,frac);
    end;    

  else // r^2>500000

    qpi_pi(r,frac);

  end; // end_if r^2<500000
END;


qpi_pi(r,frac)
BEGIN
  local frac1;

  if abs(r/pi)<101 then

    frac1:=qpi_approx(r/pi);
    frac1(3):=2;
    if (frac1(2)<1000) AND (frac1(2)<=frac(2)) then
      if frac1(2)<10 then
        qpi_out(frac1);
      else
        qpi_ln(r,frac1);
      end;
    else // (r/pi) denom not smaller
      qpi_ln(r,frac);
    end;

  else // abs(r/pi)>100

    qpi_ln(r,frac);

  end; // end_if abs(r/pi)<101
END;


qpi_ln(r,frac)
BEGIN
  local frac1,tmp;

  tmp:=e^(r);

  if tmp<1001 then

    // check for LN(0)
    if tmp then
      frac1:=qpi_approx(tmp);
    else
      frac1:=qpi_approx(MINREAL);
    end;
    frac1(3):=3;

    if (frac1(1)*frac1(2)==1) OR (frac1(2)>qpiEXPLN) then

      qpi_exp(r,frac);

    else

      if (frac1(2)<=frac(2)) then
        if frac1(2)<10 then
          qpi_out(frac1);
        else
          qpi_exp(r,frac1);
        end;
      else
        qpi_exp(r,frac);
      end;

    end; // end_if p*q==1 or q>50

  else // e^(r)>1000

    qpi_exp(r,frac);

  end; // end_if e^(r)<1001
END;


qpi_exp(r,frac)
BEGIN
  local frac1;

  if r<0 then
    qpi_out(frac);
  else

    frac1:=qpi_approx(LN(r));
    frac1(3):=4;
    if frac1(2)>qpiEXPLN then
      qpi_out(frac);
    else
      if frac1(2)<=frac(2) then
        qpi_out(frac1);
      else
        qpi_out(frac);
      end;
    end;

  end;
END;


// returns frac(t+1) where
// frac:={p/q, sqrt(p/q), p/q*pi, ln(p/q), e^(p/q)}
// and list:={p,q,t}
qpi_out(list)
BEGIN
  local s0="(", s1=")'";

  if list(3)==1 then

    qpi_outsqrt(list);

  else

    if list(1)<0 then s0:="-" + s0; end;
    if list(3) then s1:=")" + s1; end;

    case
      if list(3)==2 then s1:=")*π'"; end;
      if list(3)==3 then s0:="LN(" + s0; end;
      if list(3)==4 then s0:="e^(" + s0; end;
    end;

    s0:="'" + s0;

    if list(2)==1 then
      expr(s0 + abs(list(1)) + s1);
    else
      expr(s0 + abs(list(1)) + "/" + list(2) + s1);
    end;
  end;
END;


qpi_outsqrt(list)
BEGIN
  local ab1, ab2;
  local s0="'";

  if list(1)<0 then s0:=s0+"-"; end;

  ab1:=qpi_asqrtb(abs(list(1)));
  ab2:=qpi_asqrtb(list(2));

  if ab1(1)<>ab2(1) then
    if ab2(1)==1 then
      s0:=s0 + ab1(1) + "*";
    else
      s0:=s0 + "(" + ab1(1) + "/" + ab2(1) + ")*";
    end;
  end;

  s0:=s0 + "√(";

  if ab2(2)==1 then
    s0:=s0 + ab1(2);
  else
    s0:=s0 + ab1(2) + "/" + ab2(2);
  end;

  expr(s0+")'");
END;


// returns {a,b} where n=a*sqrt(b)
qpi_asqrtb(n) 
BEGIN
  local div,quo,rem,num,den,nodd;

  if n>qpiMAXINT then RETURN({1,n}); end;

  div:=1;
  num:=n;
  den:=4;
  nodd:=3;

  repeat
    quo:=IP(num/den);
    rem:=num MOD den;

    if rem==0 then
      div:=div*(IP(nodd/2)+1);
      num:=quo; 
    else
      nodd:=nodd+2;
      den:=den+nodd;
    end;
  until (quo==0) OR (nodd>qpiMAXINT) end; 

  RETURN({div,num});
END;


// returns {p,q,0} where r=p/q
qpi_approx(r) 
BEGIN
  local num,inum,den,iden;
  local p0,q0,p1,q1,p2,q2;
  local quo,rem;

  if NOT(r) then RETURN({0,1,0}); end;

  num:=abs(r);
  inum:=IP(num);
  den:=1;

  while num-inum do
    num:=num*10;
    den:=den*10;
    inum:=IP(num);
  end;

  iden:=den;

  rem:=den; den:=num; 
  p1:=0; p2:=1;
  q1:=1; q2:=0;

  repeat
    p0:=p1; p1:=p2;
    q0:=q1; q1:=q2;
    num:=den; den:=rem;
    quo:=IP(num/den);
    rem:=num MOD den;
    p2:=quo*p1+p0;
    q2:=quo*q1+q0;
  until 10^qpiDIGITS*abs(inum*q2-iden*p2)<iden*p2 end; 

  if (r>0) then
    RETURN({p2,q2,0});
  else
    RETURN({−p2,q2,0});
  end;
 
END;

  qpi.zip (Size: 18.24 KB / Downloads: 231)

[ramanvda]

Thank you very much!!!

[salvomic]

Han, I like very much your program!
I'd use it almost everywhere...

It, however, should treat also matrices...
As I need this, I use QPI with this ancillary program, after your advice:

Code:


#cas
qpimat(m):=
BEGIN
local s:=dim(m);
m:=mat2list(m);
m:=QPI(m)
m:=list2mat(m,s(2));
return m;
END;
#end

Do you think to include some routine to handle matrices in your original program or it couldn't work in Home mode with "extension"?

Thank you a lot,
Salvo

***
P.S. a part of π, there is a way to rationalize also √π, π^2...?

[Han]

(02-06-2015 02:57 PM)salvomic Wrote:  Han, I like very much your program!
I'd use it almost everywhere...

It, however, should treat also matrices...
Do you think to include some routine to handle matrices in your original program or it couldn't work in Home mode with "extension"?

Thank you a lot,
Salvo

***
P.S. a part of π, there is a way to rationalize also √π, π^2...?

If you want it to work in Home view, then a workaround is to convert the matrices into lists of lists. There is no getting around the fact that Home view forces all matrices to have numerical (non-symbolic) values. The drawback is that lists of lists are not displayed to look like matrices.

[salvomic]

(02-06-2015 07:54 PM)Han Wrote:  If you want it to work in Home view, then a workaround is to convert the matrices into lists of lists. There is no getting around the fact that Home view forces all matrices to have numerical (non-symbolic) values. The drawback is that lists of lists are not displayed to look like matrices.

yes, for that I prefer works almost always in CAS (I come from HP50g, with no difference between Home and CAS)...
I thought to list of lists, but I prefer real matrices.
For now my little program let me to use your QPI also with matrices (in CAS) and works well. Also I'd prefer that QPI could do the job also with matrices...

[salvomic]

please, Han, could you explain if with QPI would be possible to have rational approximation also for √π (1.77245385091), π^2 (9.86960440109), 2/π, and others "classic" irrational numbers?

Thank you!

Salvo

[Eddie W. Shore]

This is awesome, Han. Great program!

[compsystems]

Request#1: please HP-Prime team adhere this function to the firmware

[ggauny@live.fr]

Hi,
Yes a great good program !

[salvomic]

(11-05-2016 09:49 PM)compsystems Wrote:  Why this function is not built-in in the CAS?

I quote!
hoping soon...

Salvo

[compsystems]

Request#2: Please Han, QPI also for symbolic expressions

[Han]

For symbolic expressions you can use this:

Code:

#cas
qpif(f):=
begin
  local j,n,g,r;
  if (type(f) <> DOM_SYMBOLIC) then return(f); end;

  n:=dim(f)+1;

  for j from 2 to n do

    g:=f[j];
    if (type(g) == DOM_FLOAT) then
      r:=QPI(g);
      f[j]:=r;
    end;
    if (type(g) == DOM_SYMBOLIC) then
      r:=qpif(g);
      f[j]:=r;
    end;
  end;
  return(f);
end;
#end

When I have more time I may write a version that handles all known object types. But right now I am focusing on other projects.

   

[compsystems]

Thanks Han
Version QPI_4.4 that includes function with symbolic expressions, available below

Observing the source code, I see that the output could return in various formats,

Request#3:
Please Han adhere a new function called QPIRLNE( EXPR, FORMAT) where FORMAT = 0/1/.../6

case 0: (Default):
qpirlne( expr, 0) -> Expression as pi or root or ln or e


case 1: output only as a expression of QUOTIENT 1:
qpirlne( expr, 1) if it does not find the equivalent to quotient (1), but without
pi, root, ln, e, it returns the same value


case 2: output only as a expression of QUOTIENT 2:
qpirlne( expr, 2) if it does not find the equivalent to quotient (2), but without
pi, root, ln, e, it returns the same value


case 3: output only as a expression of PI:
qpirlne( expr, 3) if it does not find the equivalent to PI, it returns the same value


case 4: output only as a expression of ROOT
qpirlne( expr, 4) if it does not find the equivalent to ROOT, it returns the same value


case 5: output only as a expression of LN
qpirlne( expr, 5) if it does not find the equivalent to LN, it returns the same value


case 6: output only as a expression of EXPR
qpirlne( expr, 6) if it does not find the equivalent to EXP, it returns the same value

PHP Code:

ex#0:
 qpirlne( (2*π/3)+(3*π/4) , 0)  -> 17/12*π 

 qpirlne( (2*π/3)+(3*π/4) , 1)  -> 1137949/255685 // Q1

 qpirlne( (2*π/3)+(3*π/4) , 2)  -> 4+(115209/255685) // Q2

 qpirlne( (2*π/3)+(3*π/4) , 3)  ->  17/12*π // PI

 qpirlne( (2*π/3)+(3*π/4) , 4)  ->  (2*π/3)+(3*π/4) // ROOT

 qpirlne( (2*π/3)+(3*π/4) , 5)  ->  (2*π/3)+(3*π/4) // LN

 qpirlne( (2*π/3)+(3*π/4) , 6)  ->  (2*π/3)+(3*π/4) // e

ex#1:
 qpirlne( LN(3*π)-LN(√(5)), 0)  -> LN( (3*π*√(5)/5) ) 
 qpirlne( LN(3*π)-LN(√(5)), 1)  -> 55715/38728 // Q1
 qpirlne( LN(3*π)-LN(√(5)), 2)  -> 1+(16987/38728) // Q2
 qpirlne( LN(3*π)-LN(√(5)), 3)  -> LN( (3*π*√(5)/5) )  // PI
 qpirlne( LN(3*π)-LN(√(5)), 4)  -> LN( (3*π*√(5)/5) )  // ROOT
 qpirlne( LN(3*π)-LN(√(5)), 5)  -> LN( (3*π*√(5)/5) ) // LN
 qpirlne( LN(3*π)-LN(√(5)), 6)  -> LN(3*π)-LN(√(5)) // e

ex#2:
 qpirlne( LN((2/5))-LN(√(2)), 0)  -> -LN((25/2))/2
 qpirlne( LN((2/5))-LN(√(2)), 1)  -> -116599/92329 // Q1
 qpirlne( LN((2/5))-LN(√(2)), 2)  -> -1+(-24270/92329) // Q2
 qpirlne( LN((2/5))-LN(√(2)), 3)  -> LN((2/5))-LN(√(2)) // PI
 qpirlne( LN((2/5))-LN(√(2)), 4)  -> LN((2/5))-LN(√(2)) // ROOT
 qpirlne( LN((2/5))-LN(√(2)), 5) -> -LN((25/2))/2 // LN
 qpirlne( LN((2/5))-LN(√(2)), 6)  -> LN((2/5))-LN(√(2)) // e

ex#3:

 qpirlne( e^(2*π/(3*√(7))), 0)  -> e^((2*π*√(7)/21))
 qpirlne( e^(2*π/(3*√(7))), 1)  -> 224192/101585 // Q1
 qpirlne( e^(2*π/(3*√(7))), 2)  -> 1+(21022/101585) // Q2
 qpirlne( e^(2*π/(3*√(7))), 3)  -> e^(2*π/(3*√(7))) // PI
 qpirlne( e^(2*π/(3*√(7))), 4)  -> e^(2*π/(3*√(7))) // ROOT
 qpirlne( e^(2*π/(3*√(7))), 5)  -> e^(2*π/(3*√(7))) // LN
 qpirlne( e^(2*π/(3*√(7))), 6)  -> e^(2*π/(3*√(7))) // e


ex#4:
 qpirlne( 7*π/√(90), 0)  -> 7*π*√(10)/30 

 qpirlne( 7*π/√(90), 1)  -> 171470/73971 

 qpirlne( 7*π/√(90), 2)  -> 260521/112387

 qpirlne( 7*π/√(90), 3)  -> 7*π/√(90)


ex#5:
 qpirlne( 1/(3+i*√(3)), 0)  -> (1/4)-i*((√(3)/12))
 qpirlne( 1/(3+i*√(3)), 1)  -> (1/4)-(1/4)*i*√(1/3) // Q1
 qpirlne( 1/(3+i*√(3)), 2)  -> (1/4)-(i*37829/262087) // Q1
 qpirlne( 1/(3+i*√(3)), 3)  -> 1/(3+i*√(3)) // PI
 qpirlne( 1/(3+i*√(3)), 4)  -> (1/4)-i*((√(3)/12)) // ROOT
 qpirlne( 1/(3+i*√(3)), 5)  -> 1/(3+i*√(3)) // LN
 qpirlne( 1/(3+i*√(3)), 6)  -> 1/(3+i*√(3)) // e

ex#6:
 qpirlne( ACOS((-1/2)), 0)  -> 2/3*PI
 qpirlne( ACOS((-1/2)), 1)  -> 138894/66317 // Q1
 qpirlne( ACOS((-1/2)), 2)  -> 2*(6260/66317) // Q2
 qpirlne( ACOS((-1/2)), 3)  -> 2/3*PI // PI
 qpirlne( ACOS((-1/2)), 4)  -> ACOS((-1/2) // ROOT
 qpirlne( ACOS((-1/2)), 5)  -> ACOS((-1/2) // LN
 qpirlne( ACOS((-1/2)), 6)  -> ACOS((-1/2) // e

ex#7:
 qpirlne( COS((3*π/4)), 0)  -> -√(-2)/2
 qpirlne( COS((3*π/4)), 1)  -> -195025/275807 // Q1
 qpirlne( COS((3*π/4)), 2)  -> -195025/275807  // Q2
 qpirlne( COS((3*π/4)), 3)  -> COS((3*π/4)) // PI
 qpirlne( COS((3*π/4)), 4)  -> -√(1/2) // ROOT
 qpirlne( COS((3*π/4)), 5)  -> COS((3*π/4)) // LN
 qpirlne( COS((3*π/4)), 6)  -> COS((3*π/4)) // e

ex#8:
 qpirlne( COS(π/12), 0)  -> (√(3)+1)*(√(2)/4)
 qpirlne( COS(π/12), 1)  -> 129209/133767 // Q1
 qpirlne( COS(π/12), 2)  -> 272847/282472  // Q2
 qpirlne( COS(π/12), 3)  -> COS((3*π/4)) // PI
 qpirlne( COS(π/12), 4)  -> (√(3)+1)*(√(2)/4) // ROOT
 qpirlne( COS(π/12), 5)  -> COS(π/12) // LN
 qpirlne( COS(π/12), 6)  -> COS(π/12) // e

ex#9:
 qpirlne( SIN(π/10), 0)  -> (-1+√((5)))/4
 qpirlne( SIN(π/10), 1)  -> 98209/317811 // Q1
 qpirlne( SIN(π/10), 2)  -> 98209/317811  // Q2
 qpirlne( SIN(π/10), 3)  -> SIN(π/10) // PI
 qpirlne( SIN(π/10), 4)  -> (-1+√((5)))/4 // ROOT
 qpirlne( SIN(π/10), 5)  -> SIN(π/10) // LN
 qpirlne( SIN(π/10), 6)  -> SIN(π/10) // e

ex#10:
 qpirlne( SIN(π/8), 0)  -> √(2-√(2))/2
 qpirlne( SIN(π/8), 1)  -> 69237/180925 // Q1
 qpirlne( SIN(π/8), 2)  -> 69237/180925  // Q2
 qpirlne( SIN(π/8), 3)  -> SIN(π/8) // PI
 qpirlne( SIN(π/8), 4)  -> √(2-√(2))/2 // ROOT
 qpirlne( SIN(π/8), 5)  -> SIN(π/8) // LN
 qpirlne( SIN(π/8), 6)  -> SIN(π/8) // e

ex#11:
 qpirlne( COS(π/5), 0)  -> (1+(√(5)))/4
 qpirlne( COS(π/5), 1)  -> 98209/121393 // Q1
 qpirlne( COS(π/5), 2)  -> 317811/392836  // Q2
 qpirlne( COS(π/5), 3)  -> COS(π/5) // PI
 qpirlne( COS(π/5), 4)  ->  (1+(√(5)))/4 // ROOT
 qpirlne( COS(π/5), 5)  -> COS(π/5) // LN
 qpirlne( COS(π/5),, 6)  -> COS(π/5) // e