### Bookstore

 Computer Math and Games in Pascal (preview) Lazarus Handbook

### Author Topic: decimal floating point  (Read 937 times)

#### srvaldez

• New Member
• Posts: 36
##### decimal floating point
« on: February 17, 2024, 10:20:05 pm »
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

#### Bart

• Hero Member
• Posts: 5297
##### Re: decimal floating point
« Reply #1 on: February 17, 2024, 11:35:24 pm »
Not exactly sure what his is about, but I saw multiple goto's, which made me shudder...

Bart

#### srvaldez

• New Member
• Posts: 36
##### Re: decimal floating point
« Reply #2 on: February 17, 2024, 11:56:27 pm »
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
« Last Edit: February 17, 2024, 11:59:40 pm by srvaldez »

#### TRon

• Hero Member
• Posts: 2673
##### Re: decimal floating point
« Reply #3 on: February 18, 2024, 12:04:38 am »
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

• New Member
• Posts: 36
##### Re: decimal floating point
« Reply #4 on: February 18, 2024, 12:17:04 am »
thank you TRon for the tip
code updated.
« Last Edit: February 18, 2024, 12:35:04 am by srvaldez »

#### srvaldez

• New Member
• Posts: 36
##### Re: decimal floating point
« Reply #5 on: February 18, 2024, 03:15:27 am »
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
« Last Edit: February 18, 2024, 03:20:10 am by srvaldez »

#### TRon

• Hero Member
• Posts: 2673
##### Re: decimal floating point
« Reply #6 on: February 18, 2024, 08:07:23 am »
thank you TRon for the tip
You're welcome

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

#### srvaldez

• New Member
• Posts: 36
##### Re: decimal floating point
« Reply #7 on: February 19, 2024, 12:51:21 pm »
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;`