decimal floating point

[srvaldez]

I am a beginner FreePascal tinkerer, here's a translation of my decimal math routines
note that the str2fp is especially horrible, it tries to parse a numeric string ignoring spaces or any non-number characters so it's incredibly clunky
criticism is welcome
the source code is about 2000 lines so you can download the zipped file from https://u.pcloud.link/publink/show?code=XZICyJ0ZuSSoJh415mSOWzyfB9O6UYzwJ4R7

[Bart]

Not exactly sure what his is about, but I saw multiple goto's, which made me shudder...

Bart

[srvaldez]

parsing a string is always difficult for me, did you run the code?
I tried to use the break statement but it failed to work in that convoluted code

[TRon]

I tried to use the break statement but it failed to work in that convoluted code

Instead of break (which will stop the while loop) you should be using continue.

[srvaldez]

thank you TRon for the tip 
code updated.

[srvaldez]

provided that you declare fpln2 after fppi at the top, then you can implement the log function as follows

Code: [Select]

function agm(const a_in:decfloat; const b_in:decfloat; const digits:longint=num_digits) : decfloat;
var
a, b, c, d, eps:decfloat;
s:string;
begin
a:=a_in;
b:=b_in;
str(digits-5, s);
eps:='1e-'+s;
repeat
c := (a + b)/2;
d := sqrt(a*b);
a := c;
b := d;
until abs(a-b)<eps;
result:=a;
end;

function calcln2(const digits:longint=num_digits) : decfloat;
begin
result:=decfloat(2);
result:=fpipow(result, -2*digits+1);
result:=(fppi/(2*agm(1, result, digits)))/(2 * digits + 1);
end;

function log_agm(const x:decfloat; const digits:longint=num_digits) : decfloat;
var
y:decfloat;
begin

if x <= 0 then
begin
writeln('can nott take fplog_agm of zero or a negative number');
exit(0);
end;

if x=1 then exit(0);

if x < 1 then
begin
y:=2;
y:=fpipow(y, 2 - 2 * digits);
y*=x;
y:=fppi/(2*agm(1, y, digits))-fpln2 * 2 * digits;
y.sign:=-1;
end
else
begin
y:=2;
y:=fpipow(y, 2 - 2 * digits);
y/=x;
y:=fppi/(2*agm(1, y, digits))-fpln2 * 2 * digits;
end;
result:=y
end;
you would first calculate fpln2 before calling the log_agm function
on my PC taking the log_agm of 2e1000 to 10000 digits (num_dwords=1250) of precision took less that a second

[TRon]

thank you TRon for the tip 

You're welcome 

code updated.

That was quick    (Yes, I am slow with this response but I noticed the refreshed upload directly after your reaction).

[srvaldez]

I need help, I use the geany IDE because of the multi-language support
something is tripping geany in that it won't show any functions/procedures after fpdiv in the source explorer pane
would you inspect the function fpdiv to see what could trip-up geany's source parser ?   I would be very grateful
the complete source is in the zip file linked in the OP, here's just the fpdiv function

Code: [Select]

function fpdiv(const x:decfloat; const y:decfloat; const dwords_in:longint = num_dwords) : decfloat;
const b = 10000;
var
dwords, i, is_power_of_ten, ubn, ubw:longint;
j, last, laststep, q, t, stp:longint;
xd, xn, rund:double;
resul: array of double;
n: array of double;
d: array of double;
w: array of double;
fac1, fac2:decfloat;
begin
dwords := dwords_in;
if (dwords > num_dwords) then dwords := num_dwords;
ubn := 2 * dwords + 4;
ubw := ubn + 5;
fac1 := x;
fac2 := y;

if (fac2.exponent = 0) then // if fac2 = 0, return
begin
// a divide-by-zero error and
// bail out.
for i := 0 to dwords do
fac1.mantissa[i] := 99999999;

fac1.exponent := 99999 + bias + 1;
//er=DIVZ_ERR;
exit(fac1);
end
else if (fac1.exponent = 0) then //fact1=0, just return
exit(fac1)
else
begin
//check to see if fac2 is a power of ten
is_power_of_ten := 0;
if (fac2.mantissa[0] = 10000000) then
begin
is_power_of_ten := 1;
for i := 1 to dwords do
begin
if (fac2.mantissa[i] <> 0) then
is_power_of_ten := 0;
break;
end;
end;
end;
//if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
if (is_power_of_ten = 1) then
begin
fac1.sign := fac1.sign xor fac2.sign;
fac1.exponent := fac1.exponent - fac2.exponent + bias + 1;
exit(fac1);
end;
setLength(w, ubw + 1);
setLength(resul, ubn + 5);
setLength(n, ubn + 5);
setLength(d, ubn + 5);

//double w[ubw];

for j := 0 to dwords do
begin
n[2 * j + 2] := fac1.mantissa[j] div 10000;
n[2 * j + 3] := fac1.mantissa[j] mod 10000;
d[2 * j + 2] := fac2.mantissa[j] div 10000;
d[2 * j + 3] := fac2.mantissa[j] mod 10000;
end;
n[1] := (fac1.exponent and $7fffffff) - bias - 1;
d[1] := (fac2.exponent and $7fffffff) - bias - 1;
for j := ubn to ubw do
w[j] := 0;
t := ubn - 1;
w[1] := n[1] - d[1] + 1;
w[2] := 0;
for j := 2 to ubn do
w[j + 1] := n[j];
xd := (d[2] * b + d[3]) * b + d[4] + d[5] / b;
laststep := t + 2;
for stp := 1 to laststep do
begin
xn := realw(w, (stp + 2), ubw);
q := round(xn / xd);
last := min(stp + t + 1, ubw);
subtract(w, q, d, (stp + 2), last);
normalize(w, (stp + 2), q);
end;
finalnorm(w, (laststep + 1));
If (w[2] <> 0) Then
dec(laststep);
rund := w[laststep + 1] / b;
w[laststep] := w[laststep];
If (rund >= 0.5) Then
w[laststep]+=1;
if (w[2] = 0) then
begin
for j := 1 to (t + 1) do
resul[j] := w[j + 1];
end
else
begin
for j := 1 to (t + 1) do
resul[j] := w[j];
end;
// left-shift out leading 0's if any to preserve accuracy
{if 1}
j := 0;
if (resul[2] < 10) then
j := 3
else if (resul[2] < 100) then
j := 2
else if (resul[2] < 1000) then
j := 1;

if (j > 0) then
LSHIFT_da(resul, 2 * dwords + 4, j);
{endif}
resul[1] := w[1];
If w[2] = 0 Then
resul[1]-=1;

for j := 0 to dwords do
fac1.mantissa[j] := round(resul[2 * j + 2] * 10000 + resul[2 * j + 3]);

j := norm_fac1(fac1, dwords);
fac1.exponent := round(resul[1] + bias);

fac1.sign := fac1.sign xor fac2.sign;
result := fac1;
end;