Gamma and lnGamma using the FPU

[srvaldez]

inspired by this thread https://forum.lazarus.freepascal.org/index.php/topic,71919.0/topicseen.html and https://forum.lazarus.freepascal.org/index.php/topic,62198.msg470193.html#msg470193
I translated my gamma implementation to FPC but since fpc x64 seems to be using sse instead of the FPU there was severe precision loss due to no guard digits, so here are versions of the gamma and lnGamma functions using the FPU, they are more than likely not the fastest by a long shot but I think that the accuracy is very good

Code: [Select]

    {$MODE OBJFPC}{$H+}
    {$MINFPCONSTPREC 64}
    {$MODESWITCH ADVANCEDRECORDS}

    program gamma_using_FPU_test1;
    uses ctypes, sysutils, windows, math;

    type extfloat = record
        ext_number : array[0..4] of uint16;
    end;

    type extfloatp = ^extfloat;

    function parseExtfloat(str: PChar): extfloat; forward;

    function extsqrt(x: extfloat):extfloat;
    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldt     x
            fsqrt
            fstpt y
        end;
        extsqrt:=y;
    end;

    function extpow(x, y: extfloat):extfloat;
    var
        z:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld tbyte ptr [y]
            fld tbyte ptr [x]
            fyl2x
            fld st(0)
            frndint
            fsub st(1), st(0)
            fld1
            fscale
            fxch
            fxch st(2)
            f2xm1
            fld1
            faddp st(1), st(0)
            fmulp st(1), st(0)
            fstp st(1)
            fstp tbyte ptr [z] // z
        end;
        extpow:=z;
    end;

    function extlog(x: extfloat):extfloat;
    var
        y:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fldln2 //load loge(2)
            fld tbyte ptr [x]
            fyl2x //st(1)*log2(x)
            fstp tbyte ptr [y]
        end;
        extlog:=y;
    end;

    function extexp(x: extfloat):extfloat;
    var
        y:extfloat;
        x_e:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        x_e.ext_number[0]:=$4A9B; // the number e
        x_e.ext_number[1]:=$A2BB;
        x_e.ext_number[2]:=$5458;
        x_e.ext_number[3]:=$ADF8;
        x_e.ext_number[4]:=$4000;
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld tbyte ptr [x]
            fld tbyte ptr [x_e]
            fyl2x
            fld st(0)
            frndint
            fsub st(1), st(0)
            fld1
            fscale
            fxch
            fxch st(2)
            f2xm1
            fld1
            faddp st(1), st(0)
            fmulp st(1), st(0)
            fstp st(1)
            fstp tbyte ptr [y]
        end;
        extexp:=y;
    end;

    function extexp10(x: extfloat):extfloat;
    var
        y:extfloat;
        x_e:int32=10;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld tbyte ptr [x]
            fild dword ptr [x_e]
            fyl2x
            fld st(0)
            frndint
            fsub st(1), st(0)
            fld1
            fscale
            fxch
            fxch st(2)
            f2xm1
            fld1
            faddp st(1), st(0)
            fmulp st(1), st(0)
            fstp st(1)
            fstp tbyte ptr [y]
        end;
        extexp10:=y;
    end;

    function dbl2ext(x: double):extfloat;
    var
        y:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld qword ptr [x]
            fstp tbyte ptr [y]
        end;
        dbl2ext:=y;
    end;

    function int64ext(x: int64):extfloat;
    var
        y:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fild qword ptr [x]
            fstp tbyte ptr [y]
        end;
        int64ext:=y;
    end;

    function ext2dbl(x: extfloat):double;
    var
        y:double;

    begin
        {$ASMMODE intel}
        asm
            fld tbyte ptr [x]
            fstp qword ptr [y]
        end;
        ext2dbl:=y;
    end;

    function extadd(a: extfloat; b: extfloat):extfloat;

    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldt     a
            fldt     b
            faddp   %st, %st(1)
            fstpt y
        end;
        extadd:=y;
    end;

    function extsub(a: extfloat; b: extfloat):extfloat;

    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldt     a
            fldt     b
            fsubrp   %st, %st(1)
            fstpt y
        end;
        extsub:=y;
    end;

    function extmul(a: extfloat; b: extfloat):extfloat;

    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldt     a
            fldt     b
            fmulp   %st, %st(1)
            fstpt y
        end;
        extmul:=y;
    end;

    function extdiv(a: extfloat; b: extfloat):extfloat;

    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldt     a
            fldt     b
            fdivrp   %st, %st(1)
            fstpt y
        end;
        extdiv:=y;
    end;

    operator := (r : ansistring) z : extfloat;

    begin
        z:=parseExtfloat(PChar(r));
    end;

    operator := (r : int32) z : extfloat;
    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fild  r
            fstpt y
        end;
        z:=y;
    end;

    operator := (r : int64) z : extfloat;
    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fildl    r
            fstpt y
        end;
        z:=y;
    end;

    operator := (r : double) z : extfloat;
    var
        y:extfloat;
    begin
        {$ASMMODE att}
        asm
            fldl    r
            fstpt y
        end;
        z:=y;
    end;

    operator + (x : extfloat; y : extfloat) z : extfloat;

    begin
        z:=extadd(x, y);
    end;

    operator + (x : extfloat; y : int32) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extadd(x, w);
    end;

    operator + (x : extfloat; y : int64) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extadd(x, w);
    end;

    operator + (x : extfloat; y : double) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extadd(x, w);
    end;

    operator - (x : extfloat) z : extfloat;
    var c:extfloat;
    begin
        {$ASMMODE intel}
        asm
            fld tbyte ptr [x]
            fchs
            fstp tbyte ptr [c]
        end;
        z:=c;
    end;

    operator - (x : extfloat; y : extfloat) z : extfloat;

    begin
        z:=extsub(x, y);
    end;

    operator - (x : extfloat; y : int32) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extsub(x, w);
    end;

    operator - (x : extfloat; y : int64) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extsub(x, w);
    end;

    operator - (x : extfloat; y : double) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extsub(x, w);
    end;

    operator * (x : extfloat; y : extfloat) z : extfloat;

    begin
        z:=extmul(x, y);
    end;

    operator * (x : extfloat; y : int32) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extmul(x, w);
    end;

    operator * (x : extfloat; y : int64) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extmul(x, w);
    end;

    operator * (x : extfloat; y : double) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extmul(x, w);
    end;

    operator / (x : extfloat; y : extfloat) z : extfloat;

    begin
        z:=extdiv(x, y);
    end;

    operator / (x : extfloat; y : int32) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extdiv(x, w);
    end;

    operator / (x : extfloat; y : int64) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extdiv(x, w);
    end;

    operator / (x : extfloat; y : double) z : extfloat;
    var
        w:extfloat;
    begin
        w:=y;
        z:=extdiv(x, w);
    end;

    function exp(x : extfloat): extfloat;Overload;
    begin
        result:=extexp(x);
    end;

    function ln(x : extfloat): extfloat;Overload;
    begin
        result:=extlog(x);
    end;

    function sin(x : extfloat): extfloat;Overload;
    var
        y:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd]
            fld tbyte ptr [x]
            fsin
            fstp tbyte ptr [y]
        End;
        result:=y;
    End;

    function Abs(x : extfloat): extfloat;Overload;
    var
        y:extfloat;
        ctrlwrd:Uint16 = 4927;
    begin
        {$ASMMODE intel}
        asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld tbyte ptr [x]
            fabs
            fstp tbyte ptr [y]
        End;
        result:=y;
    End;

    function extFrac(x : extfloat) : extfloat;
        var oldcw, newcw : uint16;
            y:extfloat;
        begin
        {$ASMMODE intel}
        Asm
            fstcw word ptr [oldcw]
            mov ax,word ptr [oldcw]
            or ax,110000000000b
            mov word ptr [newcw],ax
            fldcw word ptr [newcw]
            fld tbyte ptr [x]
            fld st(0)
            frndint
            fsubp st(1),st(0)
            mov rax,[y]
            fstp tbyte ptr [y]
            fldcw word ptr [oldcw]
        End;
        extFrac:=y;
    End;

    function extTrunc(x : extfloat) : extfloat;
        var oldcw, newcw : uint16;
            y:extfloat;
        begin
        {$ASMMODE intel}
        Asm
            fstcw word ptr [oldcw]
            mov ax,word ptr [oldcw]
            or ax,110000000000b
            mov word ptr [newcw],ax
            fldcw word ptr [newcw]
            fld tbyte ptr [x]
            fld st(0)
            frndint
            fsubp st(1),st(0)
            mov rax,[y]
            fstp tbyte ptr [y]
            fldcw word ptr [oldcw]
        End;
        extTrunc:=x-y;
    End;

    function fcmp(x : extfloat; y : extfloat) : int32;
        var flag:int32;
        ctrlwrd:Uint16 = 4927;
        begin
        {$ASMMODE intel}
        Asm
            fldcw word ptr [ctrlwrd] // this guarantees extended precision
            fld tbyte ptr [x]
            fld tbyte ptr [y]
            fsubp st(1),st(0)
            ftst
            fstp st(0)
            fstsw ax
            mov al,ah
            shr al,6
            Xor ah,1
            Xor ah,al
            shl ah,1
            Or al,ah
            And rax,3
            dec rax
            mov [flag],rax
        End;
        fcmp := flag; // returns -1 if x<y, 0 if x=y and 1 if x>y
    end;

    operator ** (x : extfloat; y : extfloat) z : extfloat;

    begin
        z:=extpow(x, y);
    end;

    operator ** (x : extfloat; y : int64) z : extfloat;
    var c:extfloat;
    begin
    {$ASMMODE intel}
        asm
            fild qword ptr [y]
            fstp tbyte ptr [c]
        end;
        z:=extpow(x, c);
    end;

    operator ** (x : extfloat; y : int32) z : extfloat;
    var c:extfloat;
    begin
    {$ASMMODE intel}
        asm
            fild dword ptr [y]
            fstp tbyte ptr [c]
        end;
        z:=extpow(x, c);
    end;

    Operator <> (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop<>0;
    End;

    Operator < (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop=-1;
    End;

    Operator <= (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop<=0;
    End;

    Operator = (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop=0;
    End;

    Operator >= (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop>=0;
    End;

    Operator > (x : extfloat; y : extfloat) z : Boolean;
    var relop:int32;
    begin
        relop:=fcmp ( x, y );
        z:= relop=1;
    End;

    function parseExtfloat(str: PChar): extfloat;
    var
        parse_result: extfloat;
        fraction: extfloat;
        sign: int32;
        frac_divisor: extfloat;
        seen_decimal: int32;
        digit: int32;
        exp_sign: int32;
        exponent: int32;
    begin
        parse_result := 0;
        fraction := 0;
        sign := 1;
        frac_divisor := 1;
        seen_decimal := 0;

        // Skip leading whitespace
        while (str^ <> #0) and (str^ in [' ', #9, #10, #11, #12, #13]) do
            Inc(str);

        // Handle optional sign
        if str^ = '-' then
        begin
            sign := -1;
            Inc(str);
        end
        else if str^ = '+' then
            Inc(str);

        // Parse main number
        while str^ <> #0 do
        begin
            if str^ in ['0'..'9'] then
            begin
                digit := Ord(str^) - Ord('0');
                if seen_decimal = 0 then
                    parse_result := parse_result * 10 + digit
                else
                begin
                    fraction := fraction * 10 + digit;
                    frac_divisor := frac_divisor * 10;
                end;
            end
            else if str^ = '.' then
            begin
                if seen_decimal <> 0 then
                    Break;  // Second decimal? stop
                seen_decimal := 1;
            end
            else
                Break;
            Inc(str);
        end;

        parse_result := parse_result + fraction / frac_divisor;

        // Handle exponent part
        if (str^ = 'e') or (str^ = 'E') then
        begin
            Inc(str);
            exp_sign := 1;
            exponent := 0;

            // Optional + or - sign
            if str^ = '-' then
            begin
                exp_sign := -1;
                Inc(str);
            end
            else if str^ = '+' then
                Inc(str);

            // Parse exponent digits
            while (str^ <> #0) and (str^ in ['0'..'9']) do
            begin
                exponent := exponent * 10 + (Ord(str^) - Ord('0'));
                Inc(str);
            end;

            //parse_result := parse_result * (dbl2ext(10) ** (exp_sign * exponent));
            parse_result := parse_result * extexp10(dbl2ext((exp_sign * exponent)));
        end;

        parseExtfloat := parse_result * sign;
    end;

function gamma( x : extfloat) : double;
var
    f, sum, z, pi : extfloat;
    a : array [1..21] of extfloat;
begin
    Pi:='3.141592653589793238';
    if x<0 then begin
        if extFrac(x)<>0 then
            exit(ext2dbl(pi / (sin(pi * x) * gamma(1 - x))))
        else exit(0.0/0.0);
    end;
    if x<1 then
        exit(ext2dbl(gamma(1+x)/x));
    a[ 1]:='1.000000000000000013';
    a[ 2]:='0.0833333333333206093';
    a[ 3]:='0.003472222224360703039';
    a[ 4]:='-0.002681327303950003052';
    a[ 5]:='-0.0002294669715041962899';
    a[ 6]:='0.0007839266982119188621';
    a[ 7]:='7.13856867528538703e-05';
    a[ 8]:='-0.000609469390413885729';
    a[ 9]:='8.10071052659824655e-05';
    a[10]:='7.445022827377632428e-05';
    a[11]:='0.003418844976722079676';
    a[12]:='-0.01297314221553548962';
    a[13]:='0.02668065392136976901';
    a[14]:='-0.03792033691020890777';
    a[15]:='0.03998173682184125536';
    a[16]:='-0.03189702545933253455';
    a[17]:='0.01913944050505910382';
    a[18]:='-0.008405682814214650224';
    a[19]:='0.002556257659278982075';
    a[20]:='-0.0004816369968199514169';
    a[21]:='4.238032075012282382e-05';

    z:='2.506628274631000502'*x**(x-0.5)*exp(-x);
    f:=x*x ; sum:=f*f; sum:=sum*sum*f; f:=sum*sum; //f := x^20
    sum:=z*(a[21]+(a[20]+(a[19]+(a[18]+(a[17]+(a[16]+(a[15]+(a[14]+(a[13]
    +(a[12]+(a[11]+(a[10]+(a[9]+(a[8]+(a[7]+(a[6]+(a[5]+(a[4]+(a[3]+(x*a[1]
    +a[2])*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f;

    result := ext2dbl(sum);
end;

Function lngamma( x : extfloat) : double;
var
    f, sum : extfloat;
    c : array [1..21] of extfloat;

begin
    if x<0 then exit(0.0/0.0);
    if (x=1) or (x=2) then exit(0.0);
    if x<1 then exit(ext2dbl(lngamma(1+x)-ln(x)));
    c[ 1]:='0.9189385332046727528';
    c[ 2]:='0.08333333333332207063';
    c[ 3]:='1.925089361519014712e-12';
    c[ 4]:='-0.002777777908917571328';
    c[ 5]:='4.748526782598978435e-09';
    c[ 6]:='0.0007935451259409652588';
    c[ 7]:='1.575122516694949084e-06';
    c[ 8]:='-0.0006118621521406622235';
    c[ 9]:='0.0001288315031977228723';
    c[10]:='9.190759142084308022e-05';
    c[11]:='0.003311229742569308571';
    c[12]:='-0.01296306546218643507';
    c[13]:='0.02709215311389151191';
    c[14]:='-0.03893451905627356963';
    c[15]:='0.04141026552570293094';
    c[16]:='-0.033281057810616288';
    c[17]:='0.02010051850978016536';
    c[18]:='-0.008880268707609300334';
    c[19]:='0.002715450258845404893';
    c[20]:='-0.000514268376466983066';
    c[21]:='4.547169189857176013e-05';

    f:=x*x ; sum:=f*f; sum:=sum*sum*f; f:=sum*sum; //f := x^20
    sum:=(x-0.5)*ln(x)-x;
    sum:=sum+(c[21]+(c[20]+(c[19]+(c[18]+(c[17]+(c[16]+(c[15]+(c[14]+(c[13]
    +(c[12]+(c[11]+(c[10]+(c[9]+(c[8]+(c[7]+(c[6]+(c[5]+(c[4]+(c[3]+(x*c[1]
    +c[2])*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f;

    Result := ext2dbl(sum);
end;
(**********************************************************************)
    //const frmt=' %.20Lg';
    var x, y, z:extfloat; //try with double
        i : longint;
        d : double;
    begin
        writeln('x       gamma(x)');
        for i:=1 to 10 do
        begin;
            writeln( i, '   ',gamma(i));
        end;
        writeln;
        writeln('   x               gamma(x)');
        z:='0.1';
        x:=z;
        for i:=1 to 10 do
        begin;
            writeln( ext2dbl(x):10:10, '   ',gamma(x));
            x:=x+z;
        end;
        writeln;
        writeln('x       lngamma(x)');
        for i:=1 to 10 do
        begin;
            writeln( i, '   ',lngamma(i));
        end;
        writeln;
        writeln('   x               lngamma(x)');
        z:='0.1';
        x:=z;
        for i:=1 to 10 do
        begin;
            writeln( ext2dbl(x):10:10, '   ',lngamma(x));
            x:=x+z;
        end;
        readln;
    end.


[LV]

Great job!
Just out of curiosity, I decided to check the accuracy of the Pascal library https://github.com/chadilukito/www.wolfgang-ehrhardt.de/blob/master/src/misc/damath/sdgamma.pas

Code: Text  [Select][+][-]

1. Gamma(x) comparison for x = 1..10:
2. x     Pascal                  Scipy.special           RelError
3. 1,0   1,000000000000000E+00   1,000000000000000E+00   0,000000000000000E+00
4. 2,0   1,000000000000000E+00   1,000000000000000E+00   0,000000000000000E+00
5. 3,0   2,000000000000000E+00   2,000000000000000E+00   0,000000000000000E+00
6. 4,0   6,000000000000000E+00   6,000000000000000E+00   0,000000000000000E+00
7. 5,0   2,400000000000000E+01   2,400000000000000E+01   0,000000000000000E+00
8. 6,0   1,200000000000000E+02   1,200000000000000E+02   0,000000000000000E+00
9. 7,0   7,200000000000000E+02   7,200000000000000E+02   0,000000000000000E+00
10. 8,0   5,040000000000000E+03   5,040000000000000E+03   0,000000000000000E+00
11. 9,0   4,032000000000000E+04   4,032000000000000E+04   0,000000000000000E+00
12. 10,0   3,628800000000000E+05   3,628800000000000E+05   0,000000000000000E+00
13.  
14. Gamma(x) comparison for x = 0.1..1.0:
15. x     Pascal                  Scipy.special           RelError
16. 0,1   9,513507698668731E+00   9,513507698668732E+00   1,867194409953359E-16
17. 0,2   4,590843711998804E+00   4,590843711998803E+00   1,934673614304814E-16
18. 0,3   2,991568987687590E+00   2,991568987687590E+00   1,484469225606369E-16
19. 0,4   2,218159543757688E+00   2,218159543757688E+00   0,000000000000000E+00
20. 0,5   1,772453850905516E+00   1,772453850905516E+00   1,252752531816795E-16
21. 0,6   1,489192248812817E+00   1,489192248812817E+00   1,491040563110942E-16
22. 0,7   1,298055332647558E+00   1,298055332647558E+00   1,710594297025394E-16
23. 0,8   1,164229713725303E+00   1,164229713725303E+00   1,907223310892254E-16
24. 0,9   1,068628702119320E+00   1,068628702119319E+00   4,155692327647001E-16
25. 1,0   1,000000000000000E+00   1,000000000000000E+00   0,000000000000000E+00
26.  
27. lnGamma(x) comparison for x = 1..10:
28. x     Pascal                  Scipy.special           RelError
29. 1,0   0,000000000000000E+00   0,000000000000000E+00   0,000000000000000E+00
30. 2,0   0,000000000000000E+00   0,000000000000000E+00   0,000000000000000E+00
31. 3,0   6,931471805599453E-01   6,931471805599453E-01   0,000000000000000E+00
32. 4,0   1,791759469228055E+00   1,791759469228055E+00   0,000000000000000E+00
33. 5,0   3,178053830347946E+00   3,178053830347946E+00   1,397362139084478E-16
34. 6,0   4,787491742782046E+00   4,787491742782046E+00   0,000000000000000E+00
35. 7,0   6,579251212010101E+00   6,579251212010101E+00   0,000000000000000E+00
36. 8,0   8,525161361065415E+00   8,525161361065415E+00   0,000000000000000E+00
37. 9,0   1,060460290274525E+01   1,060460290274525E+01   0,000000000000000E+00
38. 10,0   1,280182748008147E+01   1,280182748008147E+01   1,387580673278176E-16
39.  
40. lnGamma(x) comparison for x = 0.1..1.0:
41. x     Pascal                  Scipy.special           RelError
42. 0,1   2,252712651734206E+00   2,252712651734206E+00   1,971353112915619E-16
43. 0,2   1,524063822430784E+00   1,524063822430785E+00   4,370773749570755E-16
44. 0,3   1,095797994818076E+00   1,095797994818075E+00   6,078983698867654E-16
45. 0,4   7,966778177017837E-01   7,966778177017838E-01   1,393565880656590E-16
46. 0,5   5,723649429247001E-01   5,723649429247000E-01   1,939711784149605E-16
47. 0,6   3,982338580692350E-01   3,982338580692350E-01   1,393933491752651E-16
48. 0,7   2,608672465316666E-01   2,608672465316665E-01   4,255892755361250E-16
49. 0,8   1,520596783998376E-01   1,520596783998377E-01   9,126540285928469E-16
50. 0,9   6,637623973474296E-02   6,637623973474301E-02   6,272329313896251E-16
51. 1,0   0,000000000000000E+00   2,220446049250313E-16   2,220446049250313E-16
52.  

P.S. By the way, when I followed the first link provided, https://forum.lazarus.freepascal.org/index.php/topic,71919.0/topicseen.html, I noticed that @Ed78z's posts had mysteriously disappeared. This may make it difficult to understand the context of the discussion in this thread

[srvaldez]

LV
thanks for testing