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: 228)
Search
Member List
Calendar
Help
/
QPI: convert decimal to p/q, ln(p/q), p/q*pi, e^(p/q), or sqrt(p/q)
