lazarus Pascal and library about Mathematical Functions

[sedsil]

I wonder if some mathematical functionals exist in Lazarus Pascal such that Gamma and Hypergeometrical2F1 or , binomial coefficients. How can I learn what kind of mathematical routines it has.

[wp]

NumLib comes with Free Pascal and contains some special functions: https://wiki.lazarus.freepascal.org/NumLib_Documentation#Gamma_function

Other (external) libraries are:
* lmath (https://wiki.freepascal.org/LMath): e.g., special functions in folder lmGenMath
* similar: DMath (available in Online-Package-Manager)
* Wolfgang Ehrhard's lib: hypergeometric https://github.com/chadilukito/www.wolfgang-ehrhardt.de/blob/master/src/misc/amath/sfhyperg.pas
* jpmMath https://github.com/wp-xyz/jpmMath/blob/main/original/README.md (old code, requires adaption for general use)
* plus certainly some more...

[Thaddy]

If you are looking for the functions with complex numbers, there is no library available yet.

What a nonsense. Freepascal has very good support for complex numbers in the unit ucomplex (standard unit in the rtl).
That includes all basic operations, not only +-*/, for dealing with complex math.
You can simply write out code for any formulae that uses complex math using familiar operators.
And that support has been there for years, 26 years to be precise.... and well maintained.
I used it for my own fft implementation, but many users use it.

Shame on you.

May I ask why you did not know that?

[Thaddy]

Yes, all the basics.
So?
If you know your math it is very easy to write out any function relying on complex numbers.
Have you seen they are by now actual operators?
I mean I can write complex * complex, but also complex * real. etc. You missed that!

[LV]

For Gamma function, I wrote this function last year, however it's not a highly accuracy method.

Why reinvent the wheel?  Have you checked out the links provided by @wp?

[LV]

I didn't know! could you please share the link that works for the complex numbers?

For example, Lmath:

Declaration function CLnGamma(Z : Complex) : Complex;
Description Logarithm of Gamma function

If I'm not mistaken, you just need to take the complex exponent of the result Γ(z)=exp(CLnGamma(z)):

function CGamma(z: Complex): Complex;
begin
  CGamma := CExp(CLnGamma(z));
end;

[gues1]

If you want there is the GPMATH (from sourceforge), I attach the documentation (zipped 'cause limit dimension) so you can check.

[LV]

It uses Stirling’s approximation, enhanced with Bernoulli-derived coefficients, accurate to several digits. that's good when Re(Z) >9
for 1 < Re(z) < 9 or Re(z) < 0.5   has limited accuracy and when Im(z) > 10 Might cause instability or underflow in some cases.

So, that's not a really good and accurate and "universal" function.

This isn't my sport, but let's give it a shot.

Code: Pascal  [Select][+][-]

1. function CLnGamma(z: TComplex): TComplex;
2. const
3.   g = 7.0;
4.   lanczosCoef: array[0..8] of Double = (
5.     0.99999999999980993, 676.5203681218851, -1259.1392167224028,
6.     771.32342877765313, -176.61502916214059, 12.507343278686905,
7.     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
8.   );
9. var
10.   i: Integer;
11.   z1, sum, t, tmp: TComplex;
12. begin
13.   if z.Re < 0.5 then
14.   begin
15.     tmp := CMul(Complex(Pi, 0.0), z);
16.     Result := CSub(
17.                 CSub(
18.                   CLog(Complex(Pi, 0.0)),
19.                   CLog(CSin(tmp))
20.                 ),
21.                 CLnGamma(CSub(Complex(1.0, 0.0), z))
22.               );
23.     Exit;
24.   end;
25.  
26.   z1 := CSub(z, Complex(1.0, 0.0));
27.   sum := Complex(lanczosCoef[0], 0.0);
28.   for i := 1 to 8 do
29.     sum := CAdd(sum, CDiv(Complex(lanczosCoef[i], 0.0), CAdd(z1, Complex(i, 0.0))));
30.  
31.   t := CAdd(z1, Complex(g + 0.5, 0.0));
32.   Result := CAdd(
33.               Complex(Ln(Sqrt(2 * Pi)), 0.0),
34.               CSub(
35.                 CMul(CAdd(z1, Complex(0.5, 0.0)), CLog(t)),
36.                 t
37.               )
38.             );
39.   Result := CAdd(Result, CLog(sum));
40. end;
41.  

Some results:

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

1. Test 1: lnГ(1) = 0
2.   z = 1.000000 + 0.000000i
3.   Computed: 0.000000000000 + 0.000000000000i
4.   Expected: 0.000000000000 + 0.000000000000i
5.   Abs Error: 0.00E+000
6.   Passed
7.  
8. Test 2: lnГ(3) = ln(2)
9.   z = 3.000000 + 0.000000i
10.   Computed: 0.693147180560 + 0.000000000000i
11.   Expected: 0.693147180560 + 0.000000000000i
12.   Abs Error: 5.55E-016
13.   Passed
14.  
15. Test 3: lnГ(0.5) = ln(sqrt(pi))
16.   z = 0.500000 + 0.000000i
17.   Computed: 0.572364942925 + 0.000000000000i
18.   Expected: 0.572364942925 + 0.000000000000i
19.   Abs Error: 5.55E-016
20.   Passed
21.  
22. Test 4: lnГ(0.1 - 15i)
23.   z = 0.100000 + -15.000000i
24.   Computed: -23.726199767408 + -24.989878306376i
25.   Expected: -23.726199767408 + -24.989878306376i
26.   Abs Error: 8.29E-014
27.   Passed
28.  
29. Test 5: lnГ(0.1 + 30i)
30.   z = 0.100000 + 30.000000i
31.   Computed: -47.565423555699 + 71.406325063462i
32.   Expected: -47.565423555699 + 71.406325063462i
33.   Abs Error: 1.54E-013
34.   Passed
35.  

Expected values are obtained from Python

Code: Python  [Select][+][-]

1. from mpmath import loggamma
2. z = 0.1 + 30j
3. print(loggamma(z))
4.  

[srvaldez]

this implementation is based on https://www.researchgate.net/publication/302892889_Chebychev_interpolations_of_the_Gamma_and_Polygamma_Functions_and_their_analytical_properties
it took me some time to figure out how to use the tables because they are presented using the convolution operator which was new for me https://en.wikipedia.org/wiki/Convolution and https://youtu.be/KuXjwB4LzSA?t=346, I first converted the table to a power series and then used Maple to simplify the expression
to make this functions work with complex numbers simply change some variables and functions to complex

for some reason the accuracy is not as good as the FreeBasic version
in free Pascal the relative error for gamma(100) is 3.1859E-14, in FreeBasic the relative error is 5.0702E-16
see my FreeBasic code here https://www.freebasic.net/forum/viewtopic.php?p=303286&hilit=gamma#p303286

Code: [Select]

{$mode objfpc}{$H+}

program test_gamma_function;
uses
    Math;

function gamma( x : double) : double;
var
    f, sum, z : double;
    a : array [1..20] of double = (
        0.9999999999999999,
        0.08333333333338228,
        0.003472222216158396,
       -0.002681326870868177,
       -0.0002294790668608988,
        0.0007841331256329749,
        6.903992060449035e-05,
       -0.0005907032612269776,
       -2.877475047743023e-05,
        0.0005566293593820988,
        0.001799738210158344,
       -0.008767670094590723,
        0.01817828637250317,
       -0.02452583787937907,
        0.02361068245082701,
       -0.01654210549755366,
        0.008304315532029655,
       -0.00284326571576103,
        0.0005961678245858015,
       -5.783378931872318e-05);
begin
    z:=2.506628274631001*x**(x-0.5)*exp(-x);
    f:=x*x ; f:=f*f ; f:=f*f*x ; f:=f*f*x; //f := x^19
    sum:=z*(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)/f;

    result := sum;
end;

Function lngamma( x : double) : Double;
var
    f, sum : double;
    c : array [1..20] of double = (
        0.9189385332046727,
        0.08333333333338824,
       -6.875523375600131e-12,
       -0.002777777444245219,
       -8.229108205635117e-09,
        0.0007937666109114311,
       -9.417521476898452e-07,
       -0.0005917271548255839,
        1.10417689045145e-05,
        0.0006092585812576799,
        0.001574019576696141,
       -0.008450831265113272,
        0.01796959267198273,
       -0.02456297768317419,
        0.0238450482919148,
       -0.01680609543961501,
        0.008475041030242971,
       -0.002912109145921757,
        0.000612384508536461,
       -5.955145748126546e-05);
       
begin
    f:=x*x ; f:=f*f ; f:=f*f*x ; f:=f*f*x; //f := x^19
    sum:=(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]+(c[2]
    +(c[1]-0.5*ln(x)+(ln(x)-1)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f;

    Result := sum;
end;

var i : longint;

begin
    for i:=1 to 20 do writeln( i, gamma(i));
end.

the problem seems to be with freePascal's power function and possibly others, for example
x, y : double;
x := 100;
y:=x-0.5;
x=x**y;
result = 1.0000000000000318E+199
actual value should be 1E+199

[marcov]

Assuming that you are using 32-bit for both, it might be that some code stores a value to memory and reloads it, while in the other case it is kept in a register. If that rounding from 80-bit to 64-bit can cause small deviations.  You can try to mitigate this by making sure (floating point) optimisation is on. See also the discussion on the FB forum(with  Dodicat) about determining the macheps.

And as said it might be target (win32/win64) and FB backend dependent.

[LV]

@Ed78z. I don't know what you calculated. Here is the Pascal code and the results with a relative error of 1e-16...1e-15. The comparison was made with the reliable mpmath library, called from the Python code.

Code: Pascal  [Select][+][-]

1. program LnGammaComplexDemo;
2.  
3. {$codepage utf8}
4.  
5. uses
6.   Math, SysUtils;
7.  
8. type
9.   TComplex = record
10.     Re, Im: Double;
11.   end;
12.  
13.   TRefTest = record
14.     Z: TComplex;
15.     Expected: TComplex;
16.     LabelText: string;
17.   end;
18.  
19. function Complex(Re, Im: Double): TComplex;
20. begin
21.   Result.Re := Re;
22.   Result.Im := Im;
23. end;
24.  
25. function CAdd(a, b: TComplex): TComplex;
26. begin
27.   Result.Re := a.Re + b.Re;
28.   Result.Im := a.Im + b.Im;
29. end;
30.  
31. function CSub(a, b: TComplex): TComplex;
32. begin
33.   Result.Re := a.Re - b.Re;
34.   Result.Im := a.Im - b.Im;
35. end;
36.  
37. function CMul(a, b: TComplex): TComplex;
38. begin
39.   Result.Re := a.Re * b.Re - a.Im * b.Im;
40.   Result.Im := a.Re * b.Im + a.Im * b.Re;
41. end;
42.  
43. function CDiv(a, b: TComplex): TComplex;
44. var
45.   denom: Double;
46. begin
47.   denom := b.Re * b.Re + b.Im * b.Im;
48.   Result.Re := (a.Re * b.Re + a.Im * b.Im) / denom;
49.   Result.Im := (a.Im * b.Re - a.Re * b.Im) / denom;
50. end;
51.  
52. function CLog(z: TComplex): TComplex;
53. begin
54.   Result.Re := Ln(Hypot(z.Re, z.Im));
55.   Result.Im := ArcTan2(z.Im, z.Re);
56. end;
57.  
58. function CExp(z: TComplex): TComplex;
59. var
60.   expRe: Double;
61. begin
62.   expRe := Exp(z.Re);
63.   Result.Re := expRe * Cos(z.Im);
64.   Result.Im := expRe * Sin(z.Im);
65. end;
66.  
67. function CSin(z: TComplex): TComplex;
68. begin
69.   Result.Re := Sin(z.Re) * Cosh(z.Im);
70.   Result.Im := Cos(z.Re) * Sinh(z.Im);
71. end;
72.  
73. function CLnGamma(z: TComplex): TComplex;
74. const
75.   g = 7.0;
76.   lanczosCoef: array[0..8] of Double = (
77.     0.99999999999980993, 676.5203681218851, -1259.1392167224028,
78.     771.32342877765313, -176.61502916214059, 12.507343278686905,
79.     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
80.   );
81. var
82.   i: Integer;
83.   z1, sum, t, tmp: TComplex;
84. begin
85.   if z.Re < 0.5 then
86.   begin
87.     tmp := CMul(Complex(Pi, 0.0), z);
88.     Result := CSub(
89.                 CSub(
90.                   CLog(Complex(Pi, 0.0)),
91.                   CLog(CSin(tmp))
92.                 ),
93.                 CLnGamma(CSub(Complex(1.0, 0.0), z))
94.               );
95.     Exit;
96.   end;
97.  
98.   z1 := CSub(z, Complex(1.0, 0.0));
99.   sum := Complex(lanczosCoef[0], 0.0);
100.   for i := 1 to 8 do
101.     sum := CAdd(sum, CDiv(Complex(lanczosCoef[i], 0.0), CAdd(z1, Complex(i, 0.0))));
102.  
103.   t := CAdd(z1, Complex(g + 0.5, 0.0));
104.   Result := CAdd(
105.               Complex(Ln(Sqrt(2 * Pi)), 0.0),
106.               CSub(
107.                 CMul(CAdd(z1, Complex(0.5, 0.0)), CLog(t)),
108.                 t
109.               )
110.             );
111.   Result := CAdd(Result, CLog(sum));
112. end;
113.  
114. procedure RunTests;
115. const
116.   Eps = 1e-14;
117. var
118.   res, delta: TComplex;
119.   err, normExp: Double;
120.   i: Integer;
121.   tests: array[1..7] of TRefTest = (
122.     (Z: (Re: 1.0; Im: 0.0); Expected: (Re: 0.0; Im: 0.0); LabelText: 'lnГ(1) = 0'),
123.     (Z: (Re: 0.5; Im: 0.0); Expected: (Re: 0.572364942924700087; Im: 0.0); LabelText: 'lnГ(0.5) = ln(sqrt(pi))'),
124.     (Z: (Re: 30.0; Im: 0.0); Expected: (Re: 71.257038967168009; Im: 0.0); LabelText: 'lnГ(30)'),
125.     (Z: (Re: 0.1; Im: -30.0); Expected: (Re: -47.5654235556991727; Im: -71.4063250634621394); LabelText: 'lnГ(0.1 - 30i)'),
126.     (Z: (Re: 0.1; Im: 30); Expected: (Re: -47.5654235556991727; Im: 71.4063250634621394); LabelText: 'lnГ(0.1 + 30i)'),
127.     (Z: (Re: 0.1; Im: -100.0); Expected: (Re: -158.002761620672605; Im: -359.888316732655004); LabelText: 'lnГ(0.1 - 100i)'),
128.     (Z: (Re: 0.1; Im: 100.0); Expected: (Re: -158.002761620672605; Im: 359.888316732655004); LabelText: 'lnГ(0.1 + 100i)')
129.   );
130.  
131. begin
132.   for i := 1 to Length(tests) do
133.   begin
134.     WriteLn('Test ', i, ': ', tests[i].LabelText);
135.     WriteLn(Format('  z = %.6f + %.6fi', [tests[i].Z.Re, tests[i].Z.Im]));
136.     res := CLnGamma(tests[i].Z);
137.     WriteLn(Format('  Computed: %.15f + %.15fi', [res.Re, res.Im]));
138.     WriteLn(Format('  Expected: %.15f + %.15fi', [tests[i].Expected.Re, tests[i].Expected.Im]));
139.     delta := CSub(res, tests[i].Expected);
140.     normExp := Hypot(tests[i].Expected.Re, tests[i].Expected.Im);
141.     if normExp = 0 then
142.       err := Hypot(res.Re, res.Im)
143.     else
144.       err := Hypot(delta.Re, delta.Im) / normExp;
145.     WriteLn(Format('  Rel Error: %.3e', [err]));
146.     if err < Eps then
147.       WriteLn('  Passed')
148.     else
149.       WriteLn('  Failed');
150.     WriteLn;
151.   end;
152. end;
153.  
154. begin
155.   RunTests;
156.   Readln;
157. end.
158.  

Results:

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

1. Test 1: lnГ(1) = 0
2.   z = 1,000000 + 0,000000i
3.   Computed: 0,000000000000000 + 0,000000000000000i
4.   Expected: 0,000000000000000 + 0,000000000000000i
5.   Rel Error: 0,00E+000
6.   Passed
7.  
8. Test 2: lnГ(0.5) = ln(sqrt(pi))
9.   z = 0,500000 + 0,000000i
10.   Computed: 0,572364942924700 + 0,000000000000000i
11.   Expected: 0,572364942924700 + 0,000000000000000i
12.   Rel Error: 9,70E-016
13.   Passed
14.  
15. Test 3: lnГ(30)
16.   z = 30,000000 + 0,000000i
17.   Computed: 71,257038967168015 + 0,000000000000000i
18.   Expected: 71,257038967168015 + 0,000000000000000i
19.   Rel Error: 0,00E+000
20.   Passed
21.  
22. Test 4: lnГ(0.1 - 30i)
23.   z = 0,100000 + -30,000000i
24.   Computed: -47,565423555699326 + -71,406325063462191i
25.   Expected: -47,565423555699169 + -71,406325063462134i
26.   Rel Error: 1,94E-015
27.   Passed
28.  
29. Test 5: lnГ(0.1 + 30i)
30.   z = 0,100000 + 30,000000i
31.   Computed: -47,565423555699326 + 71,406325063462191i
32.   Expected: -47,565423555699169 + 71,406325063462134i
33.   Rel Error: 1,94E-015
34.   Passed
35.  
36. Test 6: lnГ(0.1 - 100i)
37.   z = 0,100000 + -100,000000i
38.   Computed: -158,002761620672520 + -359,888316732654800i
39.   Expected: -158,002761620672600 + -359,888316732655030i
40.   Rel Error: 6,18E-016
41.   Passed
42.  
43. Test 7: lnГ(0.1 + 100i)
44.   z = 0,100000 + 100,000000i
45.   Computed: -158,002761620672520 + 359,888316732654800i
46.   Expected: -158,002761620672600 + 359,888316732655030i
47.   Rel Error: 6,18E-016
48.   Passed
49.  

[LV]

Code: Python  [Select][+][-]

1. from mpmath import mp, loggamma
2. mp.dps = 17
3. z = 0.1 + 30j
4. print(loggamma(z))
5.  

output:

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

1. (-47.565423555699173 + 71.406325063462139j)
2.  

Good night!

[wp]

Instead of lamenting that there is no advanced math library contribute to this community project: write a bug report/feature request and attach a patch with the complex gamma function for inclusion in the ucomplex unit.

[srvaldez]

I tried my implementation, Mathematica, Maple and gp/pari
they all give the same answer
log(Gamma(0.1, -30) = -47.565423555699173 - 2.291286684486688 I
but
lnGamma(0.1, -30) = -47.565423555699173 - 71.406325063462139 I
similar results for 0.1+30*I
why are the imaginary values different ?

[LV]

The imaginary parts differ because of a difference in the choice of branch of the complex logarithm when taking log(Γ(0.1,−30)). The principal branch takes arg(z)∈(−π,π]. The analytic continuation is outside this range (i.e., multiple of 2π): 71.4063−2.2913≈69.115≈11⋅2π.

Code: Python  [Select][+][-]

1. from mpmath import arg, log, gamma, mp, loggamma
2. mp.dps = 17
3.  
4. z = 0.1 - 30j
5.  
6. print('The principal branch')
7. print(log(gamma(z)))  # arg in (-π, π]
8.  
9. print('The analytic continuation')
10. print(loggamma(z))
11.  

output:

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

1. The principal branch
2. (-47.565423555699173 - 2.2912866844866882j)
3. The analytic continuation
4. (-47.565423555699173 - 71.406325063462139j)
5.  

Code: Pascal  [Select][+][-]

1. // The analytic continuation
2. function CLnGamma(z: TComplex): TComplex;
3. const
4.   g = 7.0;
5.   lanczosCoef: array[0..8] of Double = (
6.     0.99999999999980993, 676.5203681218851, -1259.1392167224028,
7.     771.32342877765313, -176.61502916214059, 12.507343278686905,
8.     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
9.   );
10. var
11.   i: Integer;
12.   z1, sum, t, tmp: TComplex;
13. begin
14.   if z.Re < 0.5 then
15.   begin
16.     tmp := CMul(Complex(Pi, 0.0), z);
17.     Result := CSub(
18.                 CSub(
19.                   CLog(Complex(Pi, 0.0)),
20.                   CLog(CSin(tmp))
21.                 ),
22.                 CLnGamma(CSub(Complex(1.0, 0.0), z))
23.               );
24.     Exit;
25.   end;
26.  
27.   z1 := CSub(z, Complex(1.0, 0.0));
28.   sum := Complex(lanczosCoef[0], 0.0);
29.   for i := 1 to 8 do
30.     sum := CAdd(sum, CDiv(Complex(lanczosCoef[i], 0.0), CAdd(z1, Complex(i, 0.0))));
31.  
32.   t := CAdd(z1, Complex(g + 0.5, 0.0));
33.   Result := CAdd(
34.               Complex(Ln(Sqrt(2 * Pi)), 0.0),
35.               CSub(
36.                 CMul(CAdd(z1, Complex(0.5, 0.0)), CLog(t)),
37.                 t
38.               )
39.             );
40.   Result := CAdd(Result, CLog(sum));
41. end;
42.  
43. // The principal branch
44. function CLnGammaPrincipal(z: TComplex): TComplex;
45. begin
46.   Result := CLnGamma(z);
47.   Result.Im := ArcTan2(Sin(Result.Im), Cos(Result.Im));
48. end;    
49.  

output:

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

1. Test 1: The principal branch
2.   z = 0,10 + -30,00i
3.   Computed: -47,5654235556993 + -2,29128668448674i
4.   Expected: -47,5654235556992 + -2,29128668448669i
5.   Rel Error: 3,31E-015
6.   Passed
7.  
8. Test 2: The analytic continuation
9.   z = 0,10 + -30,00i
10.   Computed: -47,5654235556993 + -71,4063250634622i
11.   Expected: -47,5654235556992 + -71,4063250634621i
12.   Rel Error: 1,86E-015
13.   Passed
14.