View Issue Details

IDProjectCategoryView StatusLast Update
0031534FPCRTLpublic2017-04-22 16:48
ReporterwpAssigned ToMichael Van Canneyt 
PrioritynormalSeverityminorReproducibilityhave not tried
Status closedResolutionfixed 
Product Version3.1.1Product Build 
Target Version3.2.0Fixed in Version3.1.1 
Summary0031534: [Feature Request] Add calculation of incomplete gamma function to numlib
DescriptionI am planning to write a wrapper for a more general fitting routine than the polynomial fit found in numlib (when finished it will be posted as an additional procedure for numlib).

For calculation of the goodness-of-fit parameter I need the incomplete gamma function. The attached patch adds gammap and gammaq to the unit spe of numlib.
Additional InformationCode is based on formulas in wikipedia. Results were checked against the procedures of Numerical Recipes.
TagsNo tags attached.
Fixed in Revision35594, 35687
FPCOldBugId
FPCTarget
Attached Files
  • spe.pas.patch (4,055 bytes)
    Index: numlib/src/spe.pas
    ===================================================================
    --- numlib/src/spe.pas	(revision 35574)
    +++ numlib/src/spe.pas	(working copy)
    @@ -20,6 +20,9 @@
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
     
      **********************************************************************}
    +{$mode objfpc}{$H+}
    +{$modeswitch nestedprocvars}
    +
     unit spe;
     {$I DIRECT.INC}
     
    @@ -66,6 +69,15 @@
     {  Function to calculate the natural logaritm of the Gamma function}
     function spelga(x: ArbFloat): ArbFloat;
     
    +{  Function to calculate the lower incomplete Gamma function
    +     int(t,0,x,exp(-t)t^(s-1)) / spegam(s)  (s > 0) }
    +function gammap(s, x: ArbFloat): ArbFloat;
    +
    +{  Function to calculate the upper incomplete Gamma function
    +     int(t,x,inf,exp(-t)t^(s-1)) / spegam(s)  (s > 0)
    +     gammaq(s,x) = 1 - gammap(s,x) }
    +function gammaq(s, x: ArbFloat): ArbFloat;
    +
     {  "Calculates" the maximum of two ArbFloat values     }
     function spemax(a, b: ArbFloat): ArbFloat;
     
    @@ -1003,6 +1015,115 @@
         RunError(408)
     end; {spelga}
     
    +{ Implements power series expansion for lower incomplete gamma function
    +  according to
    +  https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae
    +    gamma(s, x) = sum {k = 0 to inf, x^s exp(-x) x^k / (s (s+1) ... (s+k) ) }
    +  Converges rapidly for x < s + 1 }
    +function gammaser(s, x: ArbFloat): ArbFloat;
    +const
    +  MAX_IT = 100;
    +  EPS = 1E-7;
    +var
    +  delta: Arbfloat;
    +  sum: ArbFloat;
    +  k: Integer;
    +  lngamma: ArbFloat;
    +begin
    +  delta := 1 / s;
    +  sum := delta;
    +  for k := 1 to MAX_IT do begin
    +    delta := delta * x / (s + k);
    +    sum := sum + delta;
    +    if delta < EPS then break;
    +  end;
    +  lngamma := spelga(s);  // log of complete gamma(s)
    +  Result := exp(s * ln(x) - x + ln(sum) - lngamma);
    +end;
    +
    +type
    +  TCFFunc = function(n: Integer): ArbFloat is nested;
    +
    +{ Calculates the continued fraction a0 + (b1 / (a1 + b2 / (a2 + b3 / (a3 + b4 /...))))
    +  Ref.: https://rosettacode.org/wiki/Continued_fraction#C}
    +function CalcCF(funca, funcb: TCfFunc; NumIt: Integer): ArbFloat;
    +var
    +  r: ArbFloat;
    +  i: Integer;
    +begin
    +  r := 0;
    +  for i := NumIt downTo 1 do
    +    r := funcb(i) / (funca(i) + r);
    +  Result := funca(0) + r;
    +end;
    +
    +{ calculates the upper incomplete gamma function using its continued
    +  fraction expansion
    +  https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function }
    +function gammacf(s, x: ArbFloat): ArbFloat;
    +
    +  function funca(i: Integer): ArbFloat;
    +  begin
    +    if i = 0 then
    +      Result := 0
    +    else
    +    if odd(i) then
    +      Result := x
    +    else
    +      Result := 1;
    +  end;
    +
    +  function funcb(i: Integer): ArbFloat;
    +  begin
    +    if i = 1 then
    +      Result := 1
    +    else
    +    if odd(i) then
    +      Result := (i-1) div 2
    +    else
    +      Result := i div 2 - s;
    +  end;
    +
    +const
    +  MAX_IT = 100;
    +  EPS = 1E-7;
    +var
    +  gamma, prevgamma, lngamma: ArbFloat;
    +  i: Integer;
    +begin
    +  prevgamma := giant;
    +  i := 0;
    +  while i < MAX_IT do begin
    +    gamma := CalcCF(@funca, @funcb, i);
    +    if (abs(gamma - prevgamma) < EPS) then
    +      break;
    +    prevgamma := gamma;
    +    inc(i);
    +  end;
    +  lngamma := spelga(s);  // logarithm of complete gamma(s)
    +  Result := exp(-x + s*ln(x) - lngamma) * gamma;
    +end;
    +
    +function gammap(s, x: ArbFloat): ArbFloat;
    +begin
    +  if (x < 0.0) or (s <= 0.0) then
    +    RunError(401);                  // Invalid argument of gammap
    +  if x < s + 1 then
    +    Result := gammaser(s, x)        // Use series expansion
    +  else
    +    Result := 1.0 - gammacf(s, x);  // Use continued fraction
    +end;
    +
    +function gammaq(s, x: ArbFloat): ArbFloat;
    +begin
    +  if (x < 0.0) or (s <= 0.0) then
    +    RunError(401);                  // Invalid argument of gammaq
    +  if x < s + 1 then
    +    Result := 1.0 - gammaser(s, x)  // Use series expansion
    +  else
    +    Result := gammacf(s, x);        // Use continued fraction
    +end;
    +
     function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
     var   pa : ^arfloat0;
            i : ArbInt;
    
    spe.pas.patch (4,055 bytes)
  • spe-v2.patch (1,850 bytes)
    Index: numlib/src/spe.pas
    ===================================================================
    --- numlib/src/spe.pas	(revision 35650)
    +++ numlib/src/spe.pas	(working copy)
    @@ -115,6 +115,21 @@
     {  ArcTanH(x) }
     function speath(x: ArbFloat): ArbFloat;
     
    +{ Error numbers used within this unit:
    +
    +  401 - function spebk0(x) is not defined for x <= 0.
    +  402 - function spebk1(x) is not defined for x <= 0.
    +  403 - function speby0(x) is not defined for x <= 0.
    +  404 - function speby1(x) is not defined for x <= 0.
    +  405 - function speach(x) is not defined for x < 1
    +  406 - function speath(x) is not defined for x <= -1 or x >= 1
    +  407 - function spgam(x): x is too small or too large.
    +  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
    +  409 - function spears(s, x) is not defined for x < -1 or x > 1
    +  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
    +  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
    +}
    +
     implementation
     
     function spebi0(x: ArbFloat): ArbFloat;
    @@ -1107,7 +1122,7 @@
     function gammap(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammap
    +    RunError(410);                  // Invalid argument of gammap
       if x < s + 1 then
         Result := gammaser(s, x)        // Use series expansion
       else
    @@ -1117,7 +1132,7 @@
     function gammaq(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammaq
    +    RunError(411);                  // Invalid argument of gammaq
       if x < s + 1 then
         Result := 1.0 - gammaser(s, x)  // Use series expansion
       else
    @@ -1201,7 +1216,7 @@
     begin
       if abs(x) > 1
       then
    -    RunError(401);
    +    RunError(411);
       u:=sqr(x); uprang:= u > 0.5;
       if uprang
       then
    
    spe-v2.patch (1,850 bytes)
  • spe-v3.patch (1,847 bytes)
    Index: numlib/src/spe.pas
    ===================================================================
    --- numlib/src/spe.pas	(revision 35654)
    +++ numlib/src/spe.pas	(working copy)
    @@ -115,6 +115,21 @@
     {  ArcTanH(x) }
     function speath(x: ArbFloat): ArbFloat;
     
    +{ Error numbers used within this unit:
    +
    +  401 - function spebk0(x) is not defined for x <= 0.
    +  402 - function spebk1(x) is not defined for x <= 0.
    +  403 - function speby0(x) is not defined for x <= 0.
    +  404 - function speby1(x) is not defined for x <= 0.
    +  405 - function speach(x) is not defined for x < 1
    +  406 - function speath(x) is not defined for x <= -1 or x >= 1
    +  407 - function spgam(x): x is too small or too large.
    +  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
    +  409 - function spears(x) is not defined for x < -1 or x > 1
    +  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
    +  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
    +}
    +
     implementation
     
     function spebi0(x: ArbFloat): ArbFloat;
    @@ -1107,7 +1122,7 @@
     function gammap(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammap
    +    RunError(410);                  // Invalid argument of gammap
       if x < s + 1 then
         Result := gammaser(s, x)        // Use series expansion
       else
    @@ -1117,7 +1132,7 @@
     function gammaq(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammaq
    +    RunError(411);                  // Invalid argument of gammaq
       if x < s + 1 then
         Result := 1.0 - gammaser(s, x)  // Use series expansion
       else
    @@ -1201,7 +1216,7 @@
     begin
       if abs(x) > 1
       then
    -    RunError(401);
    +    RunError(409);
       u:=sqr(x); uprang:= u > 0.5;
       if uprang
       then
    
    spe-v3.patch (1,847 bytes)
  • spe-v5.patch (4,146 bytes)
    Index: numlib/src/spe.pas
    ===================================================================
    --- numlib/src/spe.pas	(revision 35654)
    +++ numlib/src/spe.pas	(working copy)
    @@ -115,6 +115,21 @@
     {  ArcTanH(x) }
     function speath(x: ArbFloat): ArbFloat;
     
    +{ Error numbers used within this unit:
    +
    +  401 - function spebk0(x) is not defined for x <= 0.
    +  402 - function spebk1(x) is not defined for x <= 0.
    +  403 - function speby0(x) is not defined for x <= 0.
    +  404 - function speby1(x) is not defined for x <= 0.
    +  405 - function speach(x) is not defined for x < 1
    +  406 - function speath(x) is not defined for x <= -1 or x >= 1
    +  407 - function spgam(x): x is too small or too large.
    +  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
    +  409 - function spears(s, x) is not defined for x < -1 or x > 1
    +  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
    +  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
    +}
    +
     implementation
     
     function spebi0(x: ArbFloat): ArbFloat;
    @@ -1045,18 +1060,42 @@
       TCFFunc = function(n: Integer): ArbFloat is nested;
     
     { Calculates the continued fraction a0 + (b1 / (a1 + b2 / (a2 + b3 / (a3 + b4 /...))))
    -  Ref.: https://rosettacode.org/wiki/Continued_fraction#C}
    -function CalcCF(funca, funcb: TCfFunc; NumIt: Integer): ArbFloat;
    +  using convergents.
    +  Ref.: http://lib.dr.iastate.edu/cgi/viewcontent.cgi?article=8639&context=rtd
    +    nth convergent: wn = P(n)/Q(n).
    +    P(n) = a(n) P(n-1) + b(n) P(n-2)
    +    Q(n) = a(n) Q(n-1) + b(n) Q(n-2)
    +    P(-1) = 1, P(0) = a(0), Q(-1) = 0, Q(0) = 1 }
    +function CalcCF(funca, funcb: TCfFunc; MaxIt: Integer; Eps: ArbFloat): ArbFloat;
     var
    -  r: ArbFloat;
    -  i: Integer;
    +  Pn, Pn1, Pn2: ArbFloat;
    +  Qn, Qn1, Qn2: ArbFloat;
    +  it: Integer;
    +  prev: ArbFloat;
    +  a, b: ArbFloat;
     begin
    -  r := 0;
    -  for i := NumIt downTo 1 do
    -    r := funcb(i) / (funca(i) + r);
    -  Result := funca(0) + r;
    +  Pn2 := 1.0;
    +  Pn1 := funca(0);
    +  Qn2 := 0.0;
    +  Qn1 := 1.0;
    +  prev := Giant;
    +  for it := 1 to MaxIt do begin
    +    a := funca(it);
    +    b := funcb(it);
    +    Pn := a * Pn1 + b * Pn2;
    +    Qn := a * Qn1 + b * Qn2;
    +    Result := Pn/Qn;
    +    if abs(Result - prev) < EPS * abs(Result) then
    +      exit;
    +    prev := Result;
    +    Pn2 := Pn1;
    +    Pn1 := Pn;
    +    Qn2 := Qn1;
    +    Qn1 := Qn;
    +  end;
     end;
     
    +
     { calculates the upper incomplete gamma function using its continued
       fraction expansion
       https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function }
    @@ -1087,28 +1126,17 @@
     const
       MAX_IT = 100;
       EPS = 1E-7;
    -var
    -  gamma, prevgamma, lngamma: ArbFloat;
    -  i: Integer;
     begin
    -  prevgamma := giant;
    -  i := 0;
    -  while i < MAX_IT do begin
    -    gamma := CalcCF(@funca, @funcb, i);
    -    if (abs(gamma - prevgamma) < EPS) then
    -      break;
    -    prevgamma := gamma;
    -    inc(i);
    -  end;
    -  lngamma := spelga(s);  // logarithm of complete gamma(s)
    -  Result := exp(-x + s*ln(x) - lngamma) * gamma;
    +  Result := exp(-x + s*ln(x) - spelga(s)) * CalcCF(@funca, @funcb, MAX_IT, EPS);
     end;
     
     function gammap(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammap
    -  if x < s + 1 then
    +    RunError(410);                  // Invalid argument of gammap
    +  if x = 0.0 then
    +    Result := 0.0
    +  else if x < s + 1 then
         Result := gammaser(s, x)        // Use series expansion
       else
         Result := 1.0 - gammacf(s, x);  // Use continued fraction
    @@ -1117,8 +1145,10 @@
     function gammaq(s, x: ArbFloat): ArbFloat;
     begin
       if (x < 0.0) or (s <= 0.0) then
    -    RunError(401);                  // Invalid argument of gammaq
    -  if x < s + 1 then
    +    RunError(411);                  // Invalid argument of gammaq
    +  if x = 0.0 then
    +    Result := 1.0
    +  else if x < s + 1 then
         Result := 1.0 - gammaser(s, x)  // Use series expansion
       else
         Result := gammacf(s, x);        // Use continued fraction
    @@ -1201,7 +1231,7 @@
     begin
       if abs(x) > 1
       then
    -    RunError(401);
    +    RunError(411);
       u:=sqr(x); uprang:= u > 0.5;
       if uprang
       then
    
    spe-v5.patch (4,146 bytes)

Activities

wp

2017-03-13 13:43

reporter  

spe.pas.patch (4,055 bytes)
Index: numlib/src/spe.pas
===================================================================
--- numlib/src/spe.pas	(revision 35574)
+++ numlib/src/spe.pas	(working copy)
@@ -20,6 +20,9 @@
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
  **********************************************************************}
+{$mode objfpc}{$H+}
+{$modeswitch nestedprocvars}
+
 unit spe;
 {$I DIRECT.INC}
 
@@ -66,6 +69,15 @@
 {  Function to calculate the natural logaritm of the Gamma function}
 function spelga(x: ArbFloat): ArbFloat;
 
+{  Function to calculate the lower incomplete Gamma function
+     int(t,0,x,exp(-t)t^(s-1)) / spegam(s)  (s > 0) }
+function gammap(s, x: ArbFloat): ArbFloat;
+
+{  Function to calculate the upper incomplete Gamma function
+     int(t,x,inf,exp(-t)t^(s-1)) / spegam(s)  (s > 0)
+     gammaq(s,x) = 1 - gammap(s,x) }
+function gammaq(s, x: ArbFloat): ArbFloat;
+
 {  "Calculates" the maximum of two ArbFloat values     }
 function spemax(a, b: ArbFloat): ArbFloat;
 
@@ -1003,6 +1015,115 @@
     RunError(408)
 end; {spelga}
 
+{ Implements power series expansion for lower incomplete gamma function
+  according to
+  https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae
+    gamma(s, x) = sum {k = 0 to inf, x^s exp(-x) x^k / (s (s+1) ... (s+k) ) }
+  Converges rapidly for x < s + 1 }
+function gammaser(s, x: ArbFloat): ArbFloat;
+const
+  MAX_IT = 100;
+  EPS = 1E-7;
+var
+  delta: Arbfloat;
+  sum: ArbFloat;
+  k: Integer;
+  lngamma: ArbFloat;
+begin
+  delta := 1 / s;
+  sum := delta;
+  for k := 1 to MAX_IT do begin
+    delta := delta * x / (s + k);
+    sum := sum + delta;
+    if delta < EPS then break;
+  end;
+  lngamma := spelga(s);  // log of complete gamma(s)
+  Result := exp(s * ln(x) - x + ln(sum) - lngamma);
+end;
+
+type
+  TCFFunc = function(n: Integer): ArbFloat is nested;
+
+{ Calculates the continued fraction a0 + (b1 / (a1 + b2 / (a2 + b3 / (a3 + b4 /...))))
+  Ref.: https://rosettacode.org/wiki/Continued_fraction#C}
+function CalcCF(funca, funcb: TCfFunc; NumIt: Integer): ArbFloat;
+var
+  r: ArbFloat;
+  i: Integer;
+begin
+  r := 0;
+  for i := NumIt downTo 1 do
+    r := funcb(i) / (funca(i) + r);
+  Result := funca(0) + r;
+end;
+
+{ calculates the upper incomplete gamma function using its continued
+  fraction expansion
+  https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function }
+function gammacf(s, x: ArbFloat): ArbFloat;
+
+  function funca(i: Integer): ArbFloat;
+  begin
+    if i = 0 then
+      Result := 0
+    else
+    if odd(i) then
+      Result := x
+    else
+      Result := 1;
+  end;
+
+  function funcb(i: Integer): ArbFloat;
+  begin
+    if i = 1 then
+      Result := 1
+    else
+    if odd(i) then
+      Result := (i-1) div 2
+    else
+      Result := i div 2 - s;
+  end;
+
+const
+  MAX_IT = 100;
+  EPS = 1E-7;
+var
+  gamma, prevgamma, lngamma: ArbFloat;
+  i: Integer;
+begin
+  prevgamma := giant;
+  i := 0;
+  while i < MAX_IT do begin
+    gamma := CalcCF(@funca, @funcb, i);
+    if (abs(gamma - prevgamma) < EPS) then
+      break;
+    prevgamma := gamma;
+    inc(i);
+  end;
+  lngamma := spelga(s);  // logarithm of complete gamma(s)
+  Result := exp(-x + s*ln(x) - lngamma) * gamma;
+end;
+
+function gammap(s, x: ArbFloat): ArbFloat;
+begin
+  if (x < 0.0) or (s <= 0.0) then
+    RunError(401);                  // Invalid argument of gammap
+  if x < s + 1 then
+    Result := gammaser(s, x)        // Use series expansion
+  else
+    Result := 1.0 - gammacf(s, x);  // Use continued fraction
+end;
+
+function gammaq(s, x: ArbFloat): ArbFloat;
+begin
+  if (x < 0.0) or (s <= 0.0) then
+    RunError(401);                  // Invalid argument of gammaq
+  if x < s + 1 then
+    Result := 1.0 - gammaser(s, x)  // Use series expansion
+  else
+    Result := gammacf(s, x);        // Use continued fraction
+end;
+
 function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
 var   pa : ^arfloat0;
        i : ArbInt;
spe.pas.patch (4,055 bytes)

Michael Van Canneyt

2017-03-15 16:05

administrator   ~0098922

Applied. Thanks for the contribution!

wp

2017-03-24 21:34

reporter   ~0099205

Last edited: 2017-03-25 11:05

View 2 revisions

I did not notice that I reused one of the error numbers. The new patch assigns a new error number to gammap and gammaq (and spears which has used the same number). I also added a comment with a list of all error numbers to the header of the unit.

Sorry for the inconvenience of having to touch this again.

wp

2017-03-24 21:34

reporter  

spe-v2.patch (1,850 bytes)
Index: numlib/src/spe.pas
===================================================================
--- numlib/src/spe.pas	(revision 35650)
+++ numlib/src/spe.pas	(working copy)
@@ -115,6 +115,21 @@
 {  ArcTanH(x) }
 function speath(x: ArbFloat): ArbFloat;
 
+{ Error numbers used within this unit:
+
+  401 - function spebk0(x) is not defined for x <= 0.
+  402 - function spebk1(x) is not defined for x <= 0.
+  403 - function speby0(x) is not defined for x <= 0.
+  404 - function speby1(x) is not defined for x <= 0.
+  405 - function speach(x) is not defined for x < 1
+  406 - function speath(x) is not defined for x <= -1 or x >= 1
+  407 - function spgam(x): x is too small or too large.
+  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
+  409 - function spears(s, x) is not defined for x < -1 or x > 1
+  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
+  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
+}
+
 implementation
 
 function spebi0(x: ArbFloat): ArbFloat;
@@ -1107,7 +1122,7 @@
 function gammap(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammap
+    RunError(410);                  // Invalid argument of gammap
   if x < s + 1 then
     Result := gammaser(s, x)        // Use series expansion
   else
@@ -1117,7 +1132,7 @@
 function gammaq(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammaq
+    RunError(411);                  // Invalid argument of gammaq
   if x < s + 1 then
     Result := 1.0 - gammaser(s, x)  // Use series expansion
   else
@@ -1201,7 +1216,7 @@
 begin
   if abs(x) > 1
   then
-    RunError(401);
+    RunError(411);
   u:=sqr(x); uprang:= u > 0.5;
   if uprang
   then
spe-v2.patch (1,850 bytes)

wp

2017-03-24 22:50

reporter  

spe-v3.patch (1,847 bytes)
Index: numlib/src/spe.pas
===================================================================
--- numlib/src/spe.pas	(revision 35654)
+++ numlib/src/spe.pas	(working copy)
@@ -115,6 +115,21 @@
 {  ArcTanH(x) }
 function speath(x: ArbFloat): ArbFloat;
 
+{ Error numbers used within this unit:
+
+  401 - function spebk0(x) is not defined for x <= 0.
+  402 - function spebk1(x) is not defined for x <= 0.
+  403 - function speby0(x) is not defined for x <= 0.
+  404 - function speby1(x) is not defined for x <= 0.
+  405 - function speach(x) is not defined for x < 1
+  406 - function speath(x) is not defined for x <= -1 or x >= 1
+  407 - function spgam(x): x is too small or too large.
+  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
+  409 - function spears(x) is not defined for x < -1 or x > 1
+  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
+  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
+}
+
 implementation
 
 function spebi0(x: ArbFloat): ArbFloat;
@@ -1107,7 +1122,7 @@
 function gammap(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammap
+    RunError(410);                  // Invalid argument of gammap
   if x < s + 1 then
     Result := gammaser(s, x)        // Use series expansion
   else
@@ -1117,7 +1132,7 @@
 function gammaq(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammaq
+    RunError(411);                  // Invalid argument of gammaq
   if x < s + 1 then
     Result := 1.0 - gammaser(s, x)  // Use series expansion
   else
@@ -1201,7 +1216,7 @@
 begin
   if abs(x) > 1
   then
-    RunError(401);
+    RunError(409);
   u:=sqr(x); uprang:= u > 0.5;
   if uprang
   then
spe-v3.patch (1,847 bytes)

wp

2017-03-24 22:51

reporter   ~0099208

Last edited: 2017-03-26 10:21

View 4 revisions

Oh no, patch v2 is non-sense! Please use v4.(v3 crashes for x=0 although gammap and gammaq are defined here).

wp

2017-03-25 23:55

reporter  

spe-v5.patch (4,146 bytes)
Index: numlib/src/spe.pas
===================================================================
--- numlib/src/spe.pas	(revision 35654)
+++ numlib/src/spe.pas	(working copy)
@@ -115,6 +115,21 @@
 {  ArcTanH(x) }
 function speath(x: ArbFloat): ArbFloat;
 
+{ Error numbers used within this unit:
+
+  401 - function spebk0(x) is not defined for x <= 0.
+  402 - function spebk1(x) is not defined for x <= 0.
+  403 - function speby0(x) is not defined for x <= 0.
+  404 - function speby1(x) is not defined for x <= 0.
+  405 - function speach(x) is not defined for x < 1
+  406 - function speath(x) is not defined for x <= -1 or x >= 1
+  407 - function spgam(x): x is too small or too large.
+  408 - function splga(x) cannot be calculated for x <= 0, or x is too large.
+  409 - function spears(s, x) is not defined for x < -1 or x > 1
+  410 - function gammap(s, x) is not defined for s <= 0 or x < 0
+  411 - function gammaq(s, x) is not defined for s <= 0 or x < 0
+}
+
 implementation
 
 function spebi0(x: ArbFloat): ArbFloat;
@@ -1045,18 +1060,42 @@
   TCFFunc = function(n: Integer): ArbFloat is nested;
 
 { Calculates the continued fraction a0 + (b1 / (a1 + b2 / (a2 + b3 / (a3 + b4 /...))))
-  Ref.: https://rosettacode.org/wiki/Continued_fraction#C}
-function CalcCF(funca, funcb: TCfFunc; NumIt: Integer): ArbFloat;
+  using convergents.
+  Ref.: http://lib.dr.iastate.edu/cgi/viewcontent.cgi?article=8639&context=rtd
+    nth convergent: wn = P(n)/Q(n).
+    P(n) = a(n) P(n-1) + b(n) P(n-2)
+    Q(n) = a(n) Q(n-1) + b(n) Q(n-2)
+    P(-1) = 1, P(0) = a(0), Q(-1) = 0, Q(0) = 1 }
+function CalcCF(funca, funcb: TCfFunc; MaxIt: Integer; Eps: ArbFloat): ArbFloat;
 var
-  r: ArbFloat;
-  i: Integer;
+  Pn, Pn1, Pn2: ArbFloat;
+  Qn, Qn1, Qn2: ArbFloat;
+  it: Integer;
+  prev: ArbFloat;
+  a, b: ArbFloat;
 begin
-  r := 0;
-  for i := NumIt downTo 1 do
-    r := funcb(i) / (funca(i) + r);
-  Result := funca(0) + r;
+  Pn2 := 1.0;
+  Pn1 := funca(0);
+  Qn2 := 0.0;
+  Qn1 := 1.0;
+  prev := Giant;
+  for it := 1 to MaxIt do begin
+    a := funca(it);
+    b := funcb(it);
+    Pn := a * Pn1 + b * Pn2;
+    Qn := a * Qn1 + b * Qn2;
+    Result := Pn/Qn;
+    if abs(Result - prev) < EPS * abs(Result) then
+      exit;
+    prev := Result;
+    Pn2 := Pn1;
+    Pn1 := Pn;
+    Qn2 := Qn1;
+    Qn1 := Qn;
+  end;
 end;
 
+
 { calculates the upper incomplete gamma function using its continued
   fraction expansion
   https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function }
@@ -1087,28 +1126,17 @@
 const
   MAX_IT = 100;
   EPS = 1E-7;
-var
-  gamma, prevgamma, lngamma: ArbFloat;
-  i: Integer;
 begin
-  prevgamma := giant;
-  i := 0;
-  while i < MAX_IT do begin
-    gamma := CalcCF(@funca, @funcb, i);
-    if (abs(gamma - prevgamma) < EPS) then
-      break;
-    prevgamma := gamma;
-    inc(i);
-  end;
-  lngamma := spelga(s);  // logarithm of complete gamma(s)
-  Result := exp(-x + s*ln(x) - lngamma) * gamma;
+  Result := exp(-x + s*ln(x) - spelga(s)) * CalcCF(@funca, @funcb, MAX_IT, EPS);
 end;
 
 function gammap(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammap
-  if x < s + 1 then
+    RunError(410);                  // Invalid argument of gammap
+  if x = 0.0 then
+    Result := 0.0
+  else if x < s + 1 then
     Result := gammaser(s, x)        // Use series expansion
   else
     Result := 1.0 - gammacf(s, x);  // Use continued fraction
@@ -1117,8 +1145,10 @@
 function gammaq(s, x: ArbFloat): ArbFloat;
 begin
   if (x < 0.0) or (s <= 0.0) then
-    RunError(401);                  // Invalid argument of gammaq
-  if x < s + 1 then
+    RunError(411);                  // Invalid argument of gammaq
+  if x = 0.0 then
+    Result := 1.0
+  else if x < s + 1 then
     Result := 1.0 - gammaser(s, x)  // Use series expansion
   else
     Result := gammacf(s, x);        // Use continued fraction
@@ -1201,7 +1231,7 @@
 begin
   if abs(x) > 1
   then
-    RunError(401);
+    RunError(411);
   u:=sqr(x); uprang:= u > 0.5;
   if uprang
   then
spe-v5.patch (4,146 bytes)

wp

2017-03-25 23:57

reporter   ~0099225

Last edited: 2017-03-26 10:20

View 2 revisions

Sorry for the noise, hopefully this is the last modification. Patch v5 replaces the inefficient backward calculation of the continued fraction by a forward approach, should be much faster.

Marco van de Voort

2017-03-30 16:31

manager   ~0099279

committed v5, thanks

wp

2017-04-22 16:48

reporter   ~0099783

Thank you

Issue History

Date Modified Username Field Change
2017-03-13 13:41 wp New Issue
2017-03-13 13:43 wp File Added: spe.pas.patch
2017-03-15 16:05 Michael Van Canneyt Fixed in Revision => 35594
2017-03-15 16:05 Michael Van Canneyt Note Added: 0098922
2017-03-15 16:05 Michael Van Canneyt Status new => resolved
2017-03-15 16:05 Michael Van Canneyt Fixed in Version => 3.1.1
2017-03-15 16:05 Michael Van Canneyt Resolution open => fixed
2017-03-15 16:05 Michael Van Canneyt Assigned To => Michael Van Canneyt
2017-03-15 16:05 Michael Van Canneyt Target Version => 3.2.0
2017-03-24 21:34 wp Note Added: 0099205
2017-03-24 21:34 wp Status resolved => feedback
2017-03-24 21:34 wp Resolution fixed => reopened
2017-03-24 21:34 wp File Added: spe-v2.patch
2017-03-24 22:50 wp File Added: spe-v3.patch
2017-03-24 22:51 wp Note Added: 0099208
2017-03-24 22:51 wp Status feedback => assigned
2017-03-25 11:05 wp Note Edited: 0099205 View Revisions
2017-03-25 12:02 wp Note Edited: 0099208 View Revisions
2017-03-25 23:55 wp File Added: spe-v5.patch
2017-03-25 23:57 wp Note Added: 0099225
2017-03-26 10:20 wp Note Edited: 0099225 View Revisions
2017-03-26 10:20 wp Note Edited: 0099208 View Revisions
2017-03-26 10:21 wp Note Edited: 0099208 View Revisions
2017-03-30 16:31 Marco van de Voort Fixed in Revision 35594 => 35594, 35687
2017-03-30 16:31 Marco van de Voort Note Added: 0099279
2017-03-30 16:31 Marco van de Voort Status assigned => resolved
2017-03-30 16:31 Marco van de Voort Resolution reopened => fixed
2017-04-22 16:48 wp Note Added: 0099783
2017-04-22 16:48 wp Status resolved => closed