unit MathBase.SpecialFunctions;

{ Real Double special functions. The 2.1 development slices cover cylindrical
  Bessel J/Y and modified Bessel I/K of orders zero and one, plus bounded real
  Legendre elliptic integrals and Jacobi elliptic sn/cn/dn functions.

  Formula sources:
    https://dlmf.nist.gov/10.2  (J series)
    https://dlmf.nist.gov/10.8  (Y series)
    https://dlmf.nist.gov/10.17 (large-argument expansion)
    https://dlmf.nist.gov/10.25.2 (modified I series)
    https://dlmf.nist.gov/10.31 (modified K series)
    https://dlmf.nist.gov/10.40 (modified large-argument expansions)
    https://dlmf.nist.gov/19.2 (Legendre integral definitions)
    https://dlmf.nist.gov/19.25 (relations to Carlson forms)
    https://dlmf.nist.gov/19.36 (Carlson duplication algorithms)
    https://dlmf.nist.gov/6.2 (exponential-integral definitions)
    https://dlmf.nist.gov/6.6 (exponential-integral power series)
    https://dlmf.nist.gov/6.9 (E1 continued fraction)
    https://dlmf.nist.gov/6.12 (large-argument expansions)
    https://dlmf.nist.gov/15.2 (Gauss hypergeometric series)
    https://dlmf.nist.gov/15.19.i (Maclaurin computation guidance)
    https://dlmf.nist.gov/19.2 (Legendre integral definitions)
    https://dlmf.nist.gov/19.25 (relations to Carlson forms)
    https://dlmf.nist.gov/22.2 (Jacobi elliptic function definitions)
    https://dlmf.nist.gov/22.6 (Jacobi elliptic identities)
  The small-argument series loses accuracy through cancellation as X grows;
  the large-argument expansion becomes accurate only once X is sufficiently
  large. Middle-range Chebyshev coefficients bridge that gap and were generated
  from high-precision series by tools/generate_bessel_data.py. }

{$mode objfpc}{$H+}{$J-}

interface

{ J0 and J1 accept finite |X| <= 100. Outside this validated range they
  return NaN. J0(0)=1 and J1(0)=0; J0 is even and J1 is odd. }
function BesselJ0(const X: Double): Double;
function BesselJ1(const X: Double): Double;

{ Y0 and Y1 accept finite 0 < X <= 100. Both return -Infinity at zero; Y1
  also returns -Infinity when its pole exceeds Double range. Negative,
  nonfinite, or out-of-range arguments return NaN. }
function BesselY0(const X: Double): Double;
function BesselY1(const X: Double): Double;

{ I0 and I1 accept finite |X| <= 100. I0 is even and I1 is odd. }
function ModifiedBesselI0(const X: Double): Double;
function ModifiedBesselI1(const X: Double): Double;

{ K0 and K1 accept finite 0 < X <= 100. Both return +Infinity at zero;
  negative, nonfinite, or out-of-range inputs return NaN. }
function ModifiedBesselK0(const X: Double): Double;
function ModifiedBesselK1(const X: Double): Double;

{ Legendre elliptic integrals use parameter M = k^2 in [0,1]. The
  incomplete forms accept Phi in [-Pi/2, Pi/2]. Invalid or nonfinite inputs
  return NaN; K(1) and F(+-Pi/2, 1) are signed positive/negative infinities. }
function CompleteEllipticK(const M: Double): Double;
function CompleteEllipticE(const M: Double): Double;
function IncompleteEllipticF(const Phi, M: Double): Double;
function IncompleteEllipticE(const Phi, M: Double): Double;
{ Legendre third-kind integrals use N in [-16,1] and M=k^2 in [0,1].
  Incomplete Phi is restricted to [-Pi/2,Pi/2]. Incomplete Pi is odd in Phi;
  N=0 reduces to incomplete F. Singular incomplete endpoints return signed
  infinity and singular complete cases return positive infinity. Invalid or
  nonfinite inputs return NaN. }
function CompleteEllipticPi(const N, M: Double): Double;
function IncompleteEllipticPi(const Phi, N, M: Double): Double;

{ Real Jacobi elliptic functions using parameter M=k^2 in [0,1] and finite
  |U| <= 100. Nonfinite or out-of-domain inputs return NaN. }
function JacobiEllipticSN(const U, M: Double): Double;
function JacobiEllipticCN(const U, M: Double): Double;
function JacobiEllipticDN(const U, M: Double): Double;

{ Real exponential integrals. Ei accepts finite |X| <= 100 and returns
  -Infinity at zero. E1 accepts finite 0 <= X <= 100 and returns +Infinity
  at zero. Nonfinite and out-of-domain arguments return NaN. }
function ExponentialIntegralEi(const X: Double): Double;
function ExponentialIntegralE1(const X: Double): Double;

{ Real Gauss 2F1 for A,B in [-16,16], C in [0.5,32], and |X| <= 0.75.
  Nonfinite, out-of-range, and nonconvergent inputs return NaN. }
function GaussHypergeometric2F1(const A, B, C, X: Double): Double;

implementation

uses
  Math;

{ BEGIN GENERATED BESSEL COEFFICIENTS }
{ Generated by tools/generate_bessel_data.py v1, Decimal precision 120. }
{ Chebyshev interpolation on [8, 16], 48 nodes, 40 coefficients. }
{ Formula source: https://dlmf.nist.gov/10.2 and https://dlmf.nist.gov/10.8. }
const
  J0Chebyshev: array[0..39] of Double = (
    -3.32757071513399108E-2,
    -2.72973786768856034E-2,
    -1.57547574209600824E-2,
    -1.97105061267964570E-1,
    3.88396594214053526E-2,
    5.68424062917470252E-2,
    -8.85904074702571111E-3,
    -6.04655909671513851E-3,
    8.29361866667398610E-4,
    3.43673179110411996E-4,
    -4.28910276251385126E-5,
    -1.23218453758903828E-5,
    1.42040890242856417E-6,
    3.06395666882318228E-7,
    -3.29068080930476561E-8,
    -5.61862016456436979E-9,
    5.65513025239583560E-10,
    7.93277603549392638E-11,
    -7.51587219232992200E-12,
    -8.90339593932641688E-13,
    7.96940119783780622E-14,
    8.14261483651086774E-15,
    -6.90703524584695219E-16,
    -6.18908327277302780E-17,
    4.98870721749267065E-18,
    3.97338806008743696E-19,
    -3.05074661704394298E-20,
    -2.18386171967241720E-21,
    1.60066991190402584E-22,
    1.03941009841955886E-23,
    -7.28716714305077179E-25,
    -4.32630739244591004E-26,
    2.90653883823670948E-27,
    1.58829092178254153E-28,
    -1.02424834432349241E-29,
    -5.18185680418853587E-31,
    3.21256461231234324E-32,
    1.51241218185090387E-33,
    -9.02722251696131202E-35,
    -3.97245544049885705E-36
  );
  J1Chebyshev: array[0..39] of Double = (
    1.86882513094999593E-1,
    -3.84587280723260655E-2,
    1.73233823756556791E-1,
    -5.42134854932861479E-2,
    -1.22423768145390063E-1,
    2.34658333495245574E-2,
    1.96822475839774996E-2,
    -3.11128889155257597E-3,
    -1.48070925452548524E-3,
    2.06158575117018475E-4,
    6.58200514713687463E-5,
    -8.29656300867408827E-6,
    -1.95009809602835894E-6,
    2.25890405897296761E-7,
    4.14737387067095398E-8,
    -4.45725075403683213E-9,
    -6.65912527523233630E-10,
    6.68534478798363459E-11,
    8.37343549375011243E-12,
    -7.89401851132952127E-13,
    -8.47906486099836051E-14,
    7.53826865082849405E-15,
    7.06807173380506213E-16,
    -5.94701196031533621E-17,
    -4.93740298839198325E-18,
    3.94367006758685783E-19,
    2.93320867173129515E-20,
    -2.23005345702680429E-21,
    -1.50046498264680808E-22,
    1.08844196388318867E-23,
    6.67966006155226631E-25,
    -4.63310757442709645E-26,
    -2.61163967388942535E-27,
    1.73545667516387154E-28,
    9.04034705176817874E-30,
    -5.76551018606555500E-31,
    -2.79023555617590325E-32,
    1.71061160966628412E-33,
    7.72698024826891389E-35,
    -4.56066855636516391E-36
  );
  Y0Chebyshev: array[0..39] of Double = (
    1.87644242402624852E-1,
    -3.71589118317390903E-2,
    1.72026720050605922E-1,
    -4.53264276025220173E-2,
    -1.25066336101071474E-1,
    2.13946409536613331E-2,
    2.03565636051884068E-2,
    -2.95745213330277308E-3,
    -1.53929399678976296E-3,
    2.00882916572971138E-4,
    6.84512997383048775E-5,
    -8.19919749168724856E-6,
    -2.02511515982470492E-6,
    2.25139725549853246E-7,
    4.29700838860601874E-8,
    -4.46472895888127956E-9,
    -6.88315157150652544E-10,
    6.71848429907763621E-11,
    8.63312767074907685E-12,
    -7.94622689982501646E-13,
    -8.72628678823227373E-14,
    7.60311604309732353E-15,
    7.24987900154830262E-16,
    -5.98682296626861837E-17,
    -5.08034030226358219E-18,
    4.01308507406739519E-19,
    2.94572472473407211E-20,
    -2.15855508303162524E-21,
    -1.68658188575764795E-22,
    1.35446098637761640E-23,
    2.55185694145745240E-25,
    2.43414393533746038E-26,
    -1.45279659152391041E-26,
    2.15676123454664144E-27,
    -3.21914773337179426E-28,
    5.47790598815838128E-29,
    -9.28965381846884574E-30,
    1.55189134760701413E-30,
    -2.59528098549073973E-31,
    4.34940469725697590E-32
  );
  Y1Chebyshev: array[0..39] of Double = (
    4.25732693773345723E-2,
    2.28630361544989012E-2,
    2.39938134614650272E-2,
    1.94889756205104824E-1,
    -4.39958279423179987E-2,
    -5.52429159970381245E-2,
    9.49077444183533399E-3,
    5.82677481852709584E-3,
    -8.60308024724371779E-4,
    -3.30401168631956015E-4,
    4.36650998539983426E-5,
    1.18553300595683721E-5,
    -1.43048635028152449E-6,
    -2.95360899379857415E-7,
    3.29218657925216073E-8,
    5.42968782256389722E-9,
    -5.63601399087989373E-10,
    -7.68334346413231295E-11,
    7.46976633360970496E-12,
    8.64714395418562135E-13,
    -7.91492212240606794E-14,
    -7.91428340466523837E-15,
    6.83497228461217683E-16,
    6.05834970378945064E-17,
    -4.98741265967342965E-18,
    -3.80586589268479974E-19,
    2.89436829108143321E-20,
    2.35762494694939991E-21,
    -1.96810710112608699E-22,
    -3.58969311130721905E-24,
    -4.13867087854321410E-25,
    2.38092300878959554E-25,
    -3.65747778770150511E-26,
    5.64484623513388914E-27,
    -9.88217506995467316E-28,
    1.72295088401838898E-28,
    -2.95839590677505914E-29,
    5.08131966939967476E-30,
    -8.73969137020829406E-31,
    1.50285796967265180E-31
  );
{ END GENERATED BESSEL COEFFICIENTS }

{ BEGIN GENERATED MODIFIED BESSEL COEFFICIENTS }
{ Generated by tools/generate_modified_bessel_data.py v1, Decimal precision 120. }
{ Chebyshev interpolation on [2, 16], 80 nodes, 64 coefficients. }
{ Formula source: https://dlmf.nist.gov/10.31.E1. }
const
  K0MiddleChebyshev: array[0..63] of Double = (
    3.20981414927557261E-2,
    -3.00339430583817502E-2,
    2.46636356956477503E-2,
    -1.78940826272889021E-2,
    1.15826110830163099E-2,
    -6.76802954722511712E-3,
    3.61715725188412871E-3,
    -1.79321995850520785E-3,
    8.36950278459322237E-4,
    -3.73419623678653047E-4,
    1.61689844416410546E-4,
    -6.88977967767444157E-5,
    2.92272275676857661E-5,
    -1.24456175845301355E-5,
    5.34485906457103108E-6,
    -2.31893637354987328E-6,
    1.01614919884093960E-6,
    -4.49169639455581102E-7,
    1.99992172294991216E-7,
    -8.95754582184474243E-8,
    4.03156478934965082E-8,
    -1.82184242517088234E-8,
    8.26098009848187187E-9,
    -3.75690994625044413E-9,
    1.71296554665549222E-9,
    -7.82812134404905603E-10,
    3.58470475865327346E-10,
    -1.64455895478253799E-10,
    7.55739939108896734E-11,
    -3.47823839610056263E-11,
    1.60308281884453023E-11,
    -7.39800632102476165E-12,
    3.41817058621490469E-12,
    -1.58108785145019324E-12,
    7.32096741309588593E-13,
    -3.39313825646950301E-13,
    1.57408655845501386E-13,
    -7.30845468980125441E-14,
    3.39602243247750122E-14,
    -1.57922381859688063E-14,
    7.34897522294677642E-15,
    -3.42218261430878265E-15,
    1.59462131509264180E-15,
    -7.43491189357236998E-16,
    3.46852582077965818E-16,
    -1.61902224463034867E-16,
    7.56115748656006981E-17,
    -3.53297796687933638E-17,
    1.65158548795016551E-17,
    -7.72430878293329869E-18,
    3.61416669993674106E-18,
    -1.69176020137006462E-18,
    7.92216852982754741E-19,
    -3.71122282522920211E-19,
    1.73920634842798162E-19,
    -8.15342632991259367E-20,
    3.82365176447048196E-20,
    -1.79374277385540103E-20,
    8.41744833428575848E-21,
    -3.95124817904679968E-21,
    1.85531247134981692E-21,
    -8.71413568841904041E-22,
    4.09403818757013392E-22,
    -1.92395896591356982E-22
  );
  K1MiddleChebyshev: array[0..63] of Double = (
    3.84096971335351832E-2,
    -3.60341754671531532E-2,
    2.98285705454261117E-2,
    -2.19406693553544387E-2,
    1.44907854363213385E-2,
    -8.70339954619294163E-3,
    4.82217179742831404E-3,
    -2.50255854296300670E-3,
    1.23573188041789835E-3,
    -5.89529335055984448E-4,
    2.75509990958504798E-4,
    -1.27558351009097172E-4,
    5.89740582315937769E-5,
    -2.73507136341745454E-5,
    1.27474786319104166E-5,
    -5.97127737589042109E-6,
    2.80917988812524542E-6,
    -1.32602389547469721E-6,
    6.27498782198137207E-7,
    -2.97492723671885241E-7,
    1.41232009012279761E-7,
    -6.71175928519051946E-8,
    3.19214635020268207E-8,
    -1.51914322328762857E-8,
    7.23319814095247351E-9,
    -3.44538277009576767E-9,
    1.64168289520314778E-9,
    -7.82459235096193099E-10,
    3.73023130085189904E-10,
    -1.77867283809075712E-10,
    8.48262344082861435E-11,
    -4.04601850509731207E-11,
    1.93010355649239689E-11,
    -9.20833969129399203E-12,
    4.39363582267928971E-12,
    -2.09654277571513141E-12,
    1.00049756620978670E-12,
    -4.77482315589974298E-13,
    2.27889499002225525E-13,
    -1.08771308635274166E-13,
    5.19188449300016829E-14,
    -2.47830203730252923E-14,
    1.18304200195273719E-14,
    -5.64756459191359072E-15,
    2.69609969313846018E-15,
    -1.28713213150487757E-15,
    6.14499664328011892E-16,
    -2.93380004699839828E-16,
    1.40071194489929578E-16,
    -6.68768520666742730E-17,
    3.19308715288633961E-17,
    -1.52458992104351149E-17,
    7.27950859461388304E-18,
    -3.47582025183418734E-18,
    1.65965688783823408E-18,
    -7.92473314259587040E-19,
    3.78404178851969190E-19,
    -1.80689031944309557E-19,
    8.62803556812316714E-20,
    -4.11998881190010542E-20,
    1.96736006630107056E-20,
    -9.39453146728253573E-21,
    4.48610703460764465E-21,
    -2.14223524926817724E-21
  );
{ END GENERATED MODIFIED BESSEL COEFFICIENTS }

const
  EulerGamma = 0.57721566490153286061;
  TwoOverPi = 0.63661977236758134308;
  MaxBesselArgument = 100.0;
  MaxExponentialIntegralArgument = 100.0;
  MaxJacobiEllipticArgument = 100.0;
  MaxHypergeometricParameter = 16.0;
  MaxHypergeometricDenominator = 32.0;
  MinHypergeometricDenominator = 0.5;
  MaxGaussHypergeometricArgument = 0.75;
  HalfPiValue: Double = 1.57079632679489661923;

function ChebyshevValue(const X, Center, HalfWidth: Double;
  const Coefficients: array of Double): Double;
var
  I: Integer;
  T, B0, B1, B2: Double;
begin
  T := (X - Center) / HalfWidth;
  B1 := 0.0;
  B2 := 0.0;
  for I := High(Coefficients) downto 1 do
  begin
    B0 := 2.0 * T * B1 - B2 + Coefficients[I];
    B2 := B1;
    B1 := B0;
  end;
  Result := T * B1 - B2 + 0.5 * Coefficients[0];
end;

function ModifiedBesselISeries(const X: Double; const Order: Integer): Double;
var
  K: Integer;
  Q, Term: Double;
begin
  Q := Sqr(0.5 * X);
  if Order = 0 then
    Term := 1.0
  else
    Term := 0.5 * X;
  Result := Term;
  for K := 1 to 256 do
  begin
    Term := Term * Q / (K * (K + Order));
    Result := Result + Term;
    if Term <= Abs(Result) * 1E-17 then
      Break;
  end;
end;

function ModifiedBesselK0Series(const X: Double): Double;
var
  K: Integer;
  Q, Term, Harmonic, SeriesSum: Double;
begin
  Q := Sqr(0.5 * X);
  Term := 1.0;
  Harmonic := 0.0;
  SeriesSum := 0.0;
  for K := 1 to 80 do
  begin
    Harmonic := Harmonic + 1.0 / Double(K);
    Term := Term * Q / (K * K);
    SeriesSum := SeriesSum + Harmonic * Term;
    if Term < 1E-18 then
      Break;
  end;
  Result := -(Ln(X) - Ln(2.0) + EulerGamma) *
    ModifiedBesselISeries(X, 0) + SeriesSum;
end;

function ModifiedBesselK1Series(const X: Double): Double;
var
  K: Integer;
  Q, Term, Harmonic, SeriesSum, PsiSum: Double;
begin
  Q := Sqr(0.5 * X);
  Term := 0.5 * X;
  Harmonic := 0.0;
  SeriesSum := -Term * (1.0 - 2.0 * EulerGamma) * 0.5;
  for K := 1 to 80 do
  begin
    Harmonic := Harmonic + 1.0 / Double(K);
    Term := Term * Q / (K * (K + 1));
    PsiSum := 2.0 * Harmonic + 1.0 / Double(K + 1) -
      2.0 * EulerGamma;
    SeriesSum := SeriesSum - 0.5 * PsiSum * Term;
    if Term < 1E-18 then
      Break;
  end;
  Result := 1.0 / X + (Ln(X) - Ln(2.0)) * ModifiedBesselISeries(X, 1) +
    SeriesSum;
end;

function ModifiedBesselKAsymptotic(const X: Double;
  const Order: Integer): Double;
var
  K: Integer;
  Term, NextTerm, SeriesSum: Double;
begin
  Term := 1.0;
  SeriesSum := Term;
  for K := 1 to 96 do
  begin
    NextTerm := Term * (4.0 * Order * Order - Sqr(2 * K - 1)) /
      (8.0 * K * X);
    if Abs(NextTerm) > Abs(Term) then
      Break;
    Term := NextTerm;
    SeriesSum := SeriesSum + Term;
    if Abs(Term) < Abs(SeriesSum) * 1E-17 then
      Break;
  end;
  Result := Sqrt(Pi / (2.0 * X)) * Exp(-X) * SeriesSum;
end;

function JSeries(const X: Double; const Order: Integer): Double;
var
  K: Integer;
  Q, Term: Double;
begin
  Q := Sqr(0.5 * X);
  if Order = 0 then
    Term := 1.0
  else
    Term := 0.5 * X;
  Result := Term;
  for K := 1 to 80 do
  begin
    Term := -Term * Q / (K * (K + Order));
    Result := Result + Term;
    if Abs(Term) < 1E-17 then
      Break;
  end;
end;

function Y0Series(const X: Double): Double;
var
  K: Integer;
  Q, Term, Harmonic, SeriesSum: Double;
begin
  Q := Sqr(0.5 * X);
  Term := 1.0;
  Harmonic := 0.0;
  SeriesSum := 0.0;
  for K := 1 to 80 do
  begin
    Harmonic := Harmonic + 1.0 / Double(K);
    Term := -Term * Q / (K * K);
    SeriesSum := SeriesSum - Harmonic * Term;
    if Abs(Term) < 1E-17 then
      Break;
  end;
  Result := TwoOverPi * ((Ln(X) - Ln(2.0) + EulerGamma) *
    JSeries(X, 0) + SeriesSum);
end;

function Y1Series(const X: Double): Double;
var
  K: Integer;
  Q, Term, Harmonic, SeriesSum: Double;
begin
  Q := Sqr(0.5 * X);
  Term := 0.5 * X;
  Harmonic := 0.0;
  SeriesSum := Term;
  for K := 1 to 80 do
  begin
    Harmonic := Harmonic + 1.0 / Double(K);
    Term := -Term * Q / (K * (K + 1));
    SeriesSum := SeriesSum + (2.0 * Harmonic + 1.0 / Double(K + 1)) * Term;
    if Abs(Term) < 1E-17 then
      Break;
  end;
  Result := -TwoOverPi / X + TwoOverPi * (Ln(X) - Ln(2.0) +
    EulerGamma) * JSeries(X, 1) - SeriesSum / Pi;
end;

procedure AsymptoticJY(const X: Double; const Order: Integer;
  out JValue, YValue: Double);
var
  K: Integer;
  Term, NextTerm, P, Q, Phase, Amplitude: Double;
begin
  P := 1.0;
  Q := 0.0;
  Term := 1.0;
  for K := 1 to 96 do
  begin
    NextTerm := Term * (4.0 * Order * Order - Sqr(2 * K - 1)) /
      (8.0 * K * X);
    if Abs(NextTerm) > Abs(Term) then
      Break;
    Term := NextTerm;
    if Odd(K) then
    begin
      if Odd((K - 1) div 2) then
        Q := Q - Term
      else
        Q := Q + Term;
    end
    else if Odd(K div 2) then
      P := P - Term
    else
      P := P + Term;
    if Abs(Term) < 1E-18 then
      Break;
  end;
  Phase := X - (Order * 0.5 + 0.25) * Pi;
  Amplitude := Sqrt(2.0 / (Pi * X));
  JValue := Amplitude * (P * Cos(Phase) - Q * Sin(Phase));
  YValue := Amplitude * (P * Sin(Phase) + Q * Cos(Phase));
end;

function BesselJ0(const X: Double): Double;
var
  AX, UnusedY: Double;
begin
  if IsNan(X) or IsInfinite(X) or (Abs(X) > MaxBesselArgument) then
    Exit(NaN);
  AX := Abs(X);
  if AX <= 8.0 then
    Result := JSeries(AX, 0)
  else if AX < 16.0 then
    Result := ChebyshevValue(AX, 12.0, 4.0, J0Chebyshev)
  else
    AsymptoticJY(AX, 0, Result, UnusedY);
end;

function BesselJ1(const X: Double): Double;
var
  AX, UnusedY: Double;
begin
  if IsNan(X) or IsInfinite(X) or (Abs(X) > MaxBesselArgument) then
    Exit(NaN);
  AX := Abs(X);
  if AX <= 8.0 then
    Result := JSeries(AX, 1)
  else if AX < 16.0 then
    Result := ChebyshevValue(AX, 12.0, 4.0, J1Chebyshev)
  else
    AsymptoticJY(AX, 1, Result, UnusedY);
  if X < 0.0 then
    Result := -Result;
end;

function BesselY0(const X: Double): Double;
var
  UnusedJ: Double;
begin
  if IsNan(X) or IsInfinite(X) or (X < 0.0) or
    (X > MaxBesselArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(-Infinity);
  if X <= 8.0 then
    Result := Y0Series(X)
  else if X < 16.0 then
    Result := ChebyshevValue(X, 12.0, 4.0, Y0Chebyshev)
  else
    AsymptoticJY(X, 0, UnusedJ, Result);
end;

function BesselY1(const X: Double): Double;
var
  UnusedJ: Double;
begin
  if IsNan(X) or IsInfinite(X) or (X < 0.0) or
    (X > MaxBesselArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(-Infinity);
  { Guard the leading 2/(Pi*X) pole before division: FPC raises EOverflow
    for a finite quotient outside Double instead of returning infinity. }
  if (X < 1.0) and (X * MaxDouble < TwoOverPi) then
    Exit(-Infinity);
  if X <= 8.0 then
    Result := Y1Series(X)
  else if X < 16.0 then
    Result := ChebyshevValue(X, 12.0, 4.0, Y1Chebyshev)
  else
    AsymptoticJY(X, 1, UnusedJ, Result);
end;

function ModifiedBesselI0(const X: Double): Double;
begin
  if IsNan(X) or IsInfinite(X) or (Abs(X) > MaxBesselArgument) then
    Exit(NaN);
  Result := ModifiedBesselISeries(Abs(X), 0);
end;

function ModifiedBesselI1(const X: Double): Double;
begin
  if IsNan(X) or IsInfinite(X) or (Abs(X) > MaxBesselArgument) then
    Exit(NaN);
  Result := ModifiedBesselISeries(Abs(X), 1);
  if X < 0.0 then
    Result := -Result;
end;

function ModifiedBesselK0(const X: Double): Double;
begin
  if IsNan(X) or IsInfinite(X) or (X < 0.0) or
    (X > MaxBesselArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(Infinity);
  if X <= 2.0 then
    Result := ModifiedBesselK0Series(X)
  else if X < 16.0 then
    Result := ChebyshevValue(X, 9.0, 7.0, K0MiddleChebyshev)
  else
    Result := ModifiedBesselKAsymptotic(X, 0);
end;

function ModifiedBesselK1(const X: Double): Double;
begin
  if IsNan(X) or IsInfinite(X) or (X < 0.0) or
    (X > MaxBesselArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(Infinity);
  { Avoid a hardware overflow in the 1/X leading pole. }
  if (X < 1.0) and (X * MaxDouble < 1.0) then
    Exit(Infinity);
  if X <= 2.0 then
    Result := ModifiedBesselK1Series(X)
  else if X < 16.0 then
    Result := ChebyshevValue(X, 9.0, 7.0, K1MiddleChebyshev)
  else
    Result := ModifiedBesselKAsymptotic(X, 1);
end;

function CarlsonRF(XIn, YIn, ZIn: Double): Double;
const
  ErrorTolerance = 0.0015;
var
  X, Y, Z, Average, DX, DY, DZ: Double;
  SqrtX, SqrtY, SqrtZ, Lambda, E2, E3: Double;
  Iteration: Integer;
begin
  X := XIn;
  Y := YIn;
  Z := ZIn;
  for Iteration := 1 to 64 do
  begin
    Average := (X + Y + Z) / 3.0;
    DX := (Average - X) / Average;
    DY := (Average - Y) / Average;
    DZ := (Average - Z) / Average;
    if Max(Abs(DX), Max(Abs(DY), Abs(DZ))) <= ErrorTolerance then
    begin
      E2 := DX * DY - DZ * DZ;
      E3 := DX * DY * DZ;
      Result := (1.0 + (E2 / 24.0 - 0.1 - 3.0 * E3 / 44.0) * E2 +
        E3 / 14.0) / Sqrt(Average);
      Exit;
    end;
    SqrtX := Sqrt(X);
    SqrtY := Sqrt(Y);
    SqrtZ := Sqrt(Z);
    Lambda := SqrtX * (SqrtY + SqrtZ) + SqrtY * SqrtZ;
    X := 0.25 * (X + Lambda);
    Y := 0.25 * (Y + Lambda);
    Z := 0.25 * (Z + Lambda);
  end;
  Result := NaN;
end;

function CarlsonRD(XIn, YIn, ZIn: Double): Double;
const
  ErrorTolerance = 0.0012;
var
  X, Y, Z, Average, DX, DY, DZ: Double;
  SqrtX, SqrtY, SqrtZ, Lambda, Sum, Factor: Double;
  EA, EB, EC, ED, EE, Correction: Double;
  Iteration: Integer;
begin
  X := XIn;
  Y := YIn;
  Z := ZIn;
  Sum := 0.0;
  Factor := 1.0;
  for Iteration := 1 to 64 do
  begin
    SqrtX := Sqrt(X);
    SqrtY := Sqrt(Y);
    SqrtZ := Sqrt(Z);
    Lambda := SqrtX * (SqrtY + SqrtZ) + SqrtY * SqrtZ;
    Sum := Sum + Factor / (SqrtZ * (Z + Lambda));
    Factor := Factor * 0.25;
    X := 0.25 * (X + Lambda);
    Y := 0.25 * (Y + Lambda);
    Z := 0.25 * (Z + Lambda);

    Average := (X + Y + 3.0 * Z) / 5.0;
    DX := (Average - X) / Average;
    DY := (Average - Y) / Average;
    DZ := (Average - Z) / Average;
    if Max(Abs(DX), Max(Abs(DY), Abs(DZ))) <= ErrorTolerance then
    begin
      EA := DX * DY;
      EB := DZ * DZ;
      EC := EA - EB;
      ED := EA - 6.0 * EB;
      EE := ED + 2.0 * EC;
      Correction := 1.0 + ED * (-3.0 / 14.0 + 9.0 * ED / 88.0 -
        9.0 * DZ * EE / 52.0) + DZ * (EE / 6.0 + DZ *
        (-9.0 * EC / 22.0 + 3.0 * DZ * EA / 26.0));
      Result := 3.0 * Sum + Factor * Correction /
        (Average * Sqrt(Average));
      Exit;
    end;
  end;
  Result := NaN;
end;

function CarlsonRC(XIn, YIn: Double): Double;
var
  RootRatio, Difference: Double;
begin
  if (XIn < 0.0) or (YIn <= 0.0) then
    Exit(NaN);
  if XIn = YIn then
    Exit(1.0 / Sqrt(XIn));
  if XIn < YIn then
  begin
    Difference := YIn - XIn;
    if XIn = 0.0 then
      Exit(HalfPiValue / Sqrt(Difference));
    RootRatio := Sqrt(Difference / XIn);
    Result := ArcTan(RootRatio) / Sqrt(Difference);
  end
  else
  begin
    Difference := XIn - YIn;
    RootRatio := Sqrt(Difference / XIn);
    if RootRatio < 0.5 then
      Result := ArcTanH(RootRatio) / Sqrt(Difference)
    else
      Result := Ln((Sqrt(XIn) + Sqrt(Difference)) / Sqrt(YIn)) /
        Sqrt(Difference);
  end;
end;

function CarlsonRJ(XIn, YIn, ZIn, PIn: Double): Double;
const
  ErrorTolerance = 1E-5;
var
  X, Y, Z, P, Average, DX, DY, DZ, DP: Double;
  SqrtX, SqrtY, SqrtZ, Lambda, Alpha, Beta, Sum, Factor: Double;
  EA, EB, EC, ED, EE, Correction: Double;
  Iteration: Integer;
begin
  if (XIn < 0.0) or (YIn < 0.0) or (ZIn <= 0.0) or (PIn <= 0.0) then
    Exit(NaN);
  X := XIn;
  Y := YIn;
  Z := ZIn;
  P := PIn;
  Sum := 0.0;
  Factor := 1.0;
  for Iteration := 1 to 64 do
  begin
    Average := (X + Y + Z + 2.0 * P) / 5.0;
    DX := (Average - X) / Average;
    DY := (Average - Y) / Average;
    DZ := (Average - Z) / Average;
    DP := (Average - P) / Average;
    if Max(Max(Abs(DX), Abs(DY)), Max(Abs(DZ), Abs(DP))) <=
      ErrorTolerance then
    begin
      EA := DX * (DY + DZ) + DY * DZ;
      EB := DX * DY * DZ;
      EC := DP * DP;
      ED := EA - 3.0 * EC;
      EE := EB + 2.0 * DP * (EA - EC);
      Correction := 1.0 + ED * (-3.0 / 14.0 + 9.0 * ED / 88.0 -
        9.0 * EE / 52.0) + EB * (1.0 / 6.0 + DP * (-3.0 / 11.0 +
        3.0 * DP / 26.0)) + DP * EA * (1.0 / 3.0 - 3.0 * DP / 22.0) -
        DP * EC / 3.0;
      Result := 3.0 * Sum + Factor * Correction /
        (Average * Sqrt(Average));
      Exit;
    end;
    SqrtX := Sqrt(X);
    SqrtY := Sqrt(Y);
    SqrtZ := Sqrt(Z);
    Lambda := SqrtX * (SqrtY + SqrtZ) + SqrtY * SqrtZ;
    Alpha := Sqr(P * (SqrtX + SqrtY + SqrtZ) + SqrtX * SqrtY * SqrtZ);
    Beta := P * Sqr(P + Lambda);
    Sum := Sum + Factor * CarlsonRC(Alpha, Beta);
    Factor := Factor * 0.25;
    X := 0.25 * (X + Lambda);
    Y := 0.25 * (Y + Lambda);
    Z := 0.25 * (Z + Lambda);
    P := 0.25 * (P + Lambda);
  end;
  Result := NaN;
end;

function CompleteEllipticK(const M: Double): Double;
begin
  if IsNan(M) or IsInfinite(M) or (M < 0.0) or (M > 1.0) then
    Exit(NaN);
  if M = 1.0 then
    Exit(Infinity);
  if M = 0.0 then
    Exit(Pi / 2.0);
  Result := CarlsonRF(0.0, 1.0 - M, 1.0);
end;

function CompleteEllipticE(const M: Double): Double;
var
  RFValue, RDValue: Double;
begin
  if IsNan(M) or IsInfinite(M) or (M < 0.0) or (M > 1.0) then
    Exit(NaN);
  if M = 1.0 then
    Exit(1.0);
  if M = 0.0 then
    Exit(Pi / 2.0);
  RFValue := CarlsonRF(0.0, 1.0 - M, 1.0);
  RDValue := CarlsonRD(0.0, 1.0 - M, 1.0);
  Result := RFValue - M * RDValue / 3.0;
end;

function IncompleteEllipticF(const Phi, M: Double): Double;
var
  S, C, C2, Y: Double;
begin
  if IsNan(Phi) or IsInfinite(Phi) or IsNan(M) or IsInfinite(M) or
    (M < 0.0) or (M > 1.0) or (Abs(Phi) > HalfPiValue) then
    Exit(NaN);
  if Phi = 0.0 then
    Exit(Phi);
  if M = 0.0 then
    Exit(Phi);
  if M = 1.0 then
  begin
    if Abs(Phi) = HalfPiValue then
    begin
      if Phi < 0.0 then
        Exit(-Infinity)
      else
        Exit(Infinity);
    end;
    S := Sin(Phi);
    C := Cos(Phi);
    Result := Sign(S) * Ln((1.0 + Abs(S)) / Abs(C));
    Exit;
  end;
  S := Sin(Phi);
  C := Cos(Phi);
  C2 := C * C;
  Y := C2 + (1.0 - M) * S * S;
  Result := S * CarlsonRF(C2, Y, 1.0);
end;

function IncompleteEllipticE(const Phi, M: Double): Double;
var
  S, C, C2, Y, RFValue, RDValue: Double;
begin
  if IsNan(Phi) or IsInfinite(Phi) or IsNan(M) or IsInfinite(M) or
    (M < 0.0) or (M > 1.0) or (Abs(Phi) > HalfPiValue) then
    Exit(NaN);
  S := Sin(Phi);
  if M = 0.0 then
    Exit(Phi);
  if (M = 1.0) or (Phi = 0.0) then
    Exit(S);
  C := Cos(Phi);
  C2 := C * C;
  Y := C2 + (1.0 - M) * S * S;
  RFValue := CarlsonRF(C2, Y, 1.0);
  RDValue := CarlsonRD(C2, Y, 1.0);
  Result := S * RFValue - M * S * S * S * RDValue / 3.0;
end;

function IncompleteEllipticPi(const Phi, N, M: Double): Double;
var
  S, C, C2, Y, P, RFValue, RJValue: Double;
begin
  if IsNan(Phi) or IsInfinite(Phi) or IsNan(N) or IsInfinite(N) or
    IsNan(M) or IsInfinite(M) or (N < -16.0) or (N > 1.0) or
    (M < 0.0) or (M > 1.0) or (Abs(Phi) > HalfPiValue) then
    Exit(NaN);
  if Phi = 0.0 then
    Exit(Phi);
  if N = 0.0 then
    Exit(IncompleteEllipticF(Phi, M));
  if Abs(Phi) = HalfPiValue then
  begin
    if (N = 1.0) or (M = 1.0) then
    begin
      if Phi < 0.0 then
        Exit(-Infinity)
      else
        Exit(Infinity);
    end;
  end;
  S := Sin(Phi);
  C := Cos(Phi);
  C2 := C * C;
  Y := C2 + (1.0 - M) * S * S;
  P := C2 + (1.0 - N) * S * S;
  RFValue := CarlsonRF(C2, Y, 1.0);
  RJValue := CarlsonRJ(C2, Y, 1.0, P);
  Result := S * RFValue + N * S * S * S * RJValue / 3.0;
end;

function CompleteEllipticPi(const N, M: Double): Double;
begin
  if IsNan(N) or IsInfinite(N) or IsNan(M) or IsInfinite(M) or
    (N < -16.0) or (N > 1.0) or (M < 0.0) or (M > 1.0) then
    Exit(NaN);
  if (N = 1.0) or (M = 1.0) then
    Exit(Infinity);
  Result := IncompleteEllipticPi(HalfPiValue, N, M);
end;

procedure JacobiEllipticValues(const U, M: Double;
  out SNValue, CNValue, DNValue: Double);
const
  MaxIterations = 32;
  MaxBisectionIterations = 96;
  FunctionTolerance = 1E-14;
  AmplitudeTolerance = 2E-14;
var
  PhiValue, FunctionValue, Difference, LowPhi, HighPhi, CandidatePhi,
    InverseDerivative, ExponentialArgument, QuarterPeriod, ReducedU,
    SignValue: Double;
  I: Integer;
  PeriodCount: Int64;
  Converged: Boolean;
begin
  if IsNan(U) or IsInfinite(U) or
    (Abs(U) > MaxJacobiEllipticArgument) or IsNan(M) or IsInfinite(M) or
    (M < 0.0) or (M > 1.0) then
  begin
    SNValue := NaN;
    CNValue := NaN;
    DNValue := NaN;
    Exit;
  end;

  if M = 0.0 then
  begin
    SNValue := Sin(U);
    CNValue := Cos(U);
    DNValue := 1.0;
    Exit;
  end;
  if M = 1.0 then
  begin
    SNValue := Tanh(U);
    ExponentialArgument := Exp(-Abs(U));
    CNValue := 2.0 * ExponentialArgument /
      (1.0 + ExponentialArgument * ExponentialArgument);
    DNValue := CNValue;
    Exit;
  end;

  { Invert F(phi|M)=U on the principal interval. The inverse derivative is
    sqrt(1-M*sin(phi)^2); safeguarding Newton with bisection keeps convergence
    reliable near the quarter period. Reduce by 2K first, where K is evaluated
    by the qualified Carlson RF implementation. }
  QuarterPeriod := CompleteEllipticK(M);
  PeriodCount := Round(U / (2.0 * QuarterPeriod));
  ReducedU := U - PeriodCount * (2.0 * QuarterPeriod);
  if ReducedU > QuarterPeriod then
    ReducedU := QuarterPeriod
  else if ReducedU < -QuarterPeriod then
    ReducedU := -QuarterPeriod;
  if Abs(PeriodCount mod 2) = 1 then
    SignValue := -1.0
  else
    SignValue := 1.0;
  LowPhi := -Pi / 2.0;
  HighPhi := Pi / 2.0;
  PhiValue := ReducedU * (Pi / (2.0 * QuarterPeriod));
  if PhiValue > HighPhi then
    PhiValue := HighPhi
  else if PhiValue < LowPhi then
    PhiValue := LowPhi;
  Converged := False;
  for I := 1 to MaxIterations do
  begin
    FunctionValue := IncompleteEllipticF(PhiValue, M);
    Difference := FunctionValue - ReducedU;
    if Abs(Difference) <= FunctionTolerance then
    begin
      Converged := True;
      Break;
    end;
    if Difference < 0.0 then
      LowPhi := PhiValue
    else
      HighPhi := PhiValue;
    InverseDerivative := Sqrt(Max(0.0,
      1.0 - M * Sqr(Sin(PhiValue))));
    CandidatePhi := PhiValue - Difference * InverseDerivative;
    if (CandidatePhi <= LowPhi) or (CandidatePhi >= HighPhi) or
      (CandidatePhi < LowPhi + 0.1 * (HighPhi - LowPhi)) or
      (CandidatePhi > HighPhi - 0.1 * (HighPhi - LowPhi)) or
      IsNan(CandidatePhi) or IsInfinite(CandidatePhi) then
      CandidatePhi := (LowPhi + HighPhi) * 0.5;
    PhiValue := CandidatePhi;
  end;
  if not Converged then
    for I := 1 to MaxBisectionIterations do
    begin
      PhiValue := (LowPhi + HighPhi) * 0.5;
      FunctionValue := IncompleteEllipticF(PhiValue, M);
      Difference := FunctionValue - ReducedU;
      if IsNan(Difference) or IsInfinite(Difference) then
        Break;
      if Abs(Difference) <= FunctionTolerance then
      begin
        Converged := True;
        Break;
      end;
      if Difference < 0.0 then
        LowPhi := PhiValue
      else
        HighPhi := PhiValue;
      if HighPhi - LowPhi <= AmplitudeTolerance then
      begin
        PhiValue := (LowPhi + HighPhi) * 0.5;
        Converged := True;
        Break;
      end;
    end;
  if not Converged then
  begin
    SNValue := NaN;
    CNValue := NaN;
    DNValue := NaN;
    Exit;
  end;
  SNValue := SignValue * Sin(PhiValue);
  CNValue := SignValue * Cos(PhiValue);
  { For real U and M, dn is positive. This identity remains stable near
    quarter-period values. }
  DNValue := Sqrt(Max(0.0, 1.0 - M * SNValue * SNValue));
end;

function JacobiEllipticSN(const U, M: Double): Double;
var
  CNValue, DNValue: Double;
begin
  JacobiEllipticValues(U, M, Result, CNValue, DNValue);
end;

function JacobiEllipticCN(const U, M: Double): Double;
var
  SNValue, DNValue: Double;
begin
  JacobiEllipticValues(U, M, SNValue, Result, DNValue);
end;

function JacobiEllipticDN(const U, M: Double): Double;
var
  SNValue, CNValue: Double;
begin
  JacobiEllipticValues(U, M, SNValue, CNValue, Result);
end;

function ExponentialIntegralE1Series(const X: Double): Double;
var
  K: Integer;
  Term, Addend, SeriesSum: Double;
begin
  Term := 1.0;
  SeriesSum := 0.0;
  for K := 1 to 256 do
  begin
    Term := -Term * X / K;
    Addend := Term / K;
    SeriesSum := SeriesSum - Addend;
    if Abs(Addend) <= Abs(SeriesSum) * 1E-17 then
      Break;
  end;
  Result := -EulerGamma - Ln(X) + SeriesSum;
end;

function ExponentialIntegralE1ContinuedFraction(const X: Double): Double;
const
  Tiny = 1E-300;
  RelativeTolerance = 2E-16;
var
  F, C, D, A, B, Delta: Double;
  N, Numerator: Integer;
begin
  { DLMF 6.9.1 in modified Lentz form. The even/odd denominators alternate
    between X and 1; convergence is rapid once X is bounded away from zero. }
  F := X;
  C := F;
  D := 0.0;
  for N := 1 to 512 do
  begin
    Numerator := (N + 1) div 2;
    A := Numerator;
    if Odd(N) then
      B := 1.0
    else
      B := X;

    D := B + A * D;
    if Abs(D) < Tiny then
      D := Tiny;
    D := 1.0 / D;
    C := B + A / C;
    if Abs(C) < Tiny then
      C := Tiny;
    Delta := C * D;
    F := F * Delta;
    if Abs(Delta - 1.0) <= RelativeTolerance then
    begin
      Result := Exp(-X) / F;
      Exit;
    end;
  end;
  Result := NaN;
end;

function ExponentialIntegralE1(const X: Double): Double;
begin
  if IsNan(X) or IsInfinite(X) or (X < 0.0) or
    (X > MaxExponentialIntegralArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(Infinity);
  if X <= 2.0 then
    Result := ExponentialIntegralE1Series(X)
  else
    Result := ExponentialIntegralE1ContinuedFraction(X);
end;

function ExponentialIntegralEi(const X: Double): Double;
var
  K: Integer;
  Term, Addend, SeriesSum, AX: Double;
begin
  if IsNan(X) or IsInfinite(X) or
    (Abs(X) > MaxExponentialIntegralArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(-Infinity);
  if X < 0.0 then
    Exit(-ExponentialIntegralE1(-X));

  { The positive-axis DLMF series has no cancellation in this bounded range. }
  AX := X;
  Term := AX;
  SeriesSum := 0.0;
  for K := 1 to 512 do
  begin
    if K > 1 then
      Term := Term * AX / K;
    Addend := Term / K;
    SeriesSum := SeriesSum + Addend;
    if Abs(Addend) <= Abs(SeriesSum) * 1E-17 then
      Break;
  end;
  Result := EulerGamma + Ln(AX) + SeriesSum;
end;

function GaussHypergeometric2F1(const A, B, C, X: Double): Double;
const
  MaxIterations = 10000;
  MinimumStoppingIteration = 128;
  RelativeTermTolerance = 2E-16;
var
  N: Integer;
  Term, SumValue, Compensation, AdjustedTerm, UpdatedSum: Double;
  NumeratorA, NumeratorB: Double;
begin
  if IsNan(A) or IsInfinite(A) or (Abs(A) > MaxHypergeometricParameter) or
    IsNan(B) or IsInfinite(B) or (Abs(B) > MaxHypergeometricParameter) or
    IsNan(C) or IsInfinite(C) or
    (C < MinHypergeometricDenominator) or
    (C > MaxHypergeometricDenominator) or
    IsNan(X) or IsInfinite(X) or
    (Abs(X) > MaxGaussHypergeometricArgument) then
    Exit(NaN);
  if X = 0.0 then
    Exit(1.0);

  { DLMF 15.2.1 Gauss series. Kahan compensation limits loss when the
    alternating terms nearly cancel; the compact argument disk gives a
    convergent series without analytic continuation or branch handling. }
  Term := 1.0;
  SumValue := 1.0;
  Compensation := 0.0;
  for N := 1 to MaxIterations do
  begin
    NumeratorA := A + N - 1;
    NumeratorB := B + N - 1;
    if (NumeratorA = 0.0) or (NumeratorB = 0.0) then
      Exit(SumValue);
    Term := Term * NumeratorA * NumeratorB * X /
      ((C + N - 1) * N);
    if IsNan(Term) or IsInfinite(Term) then
      Exit(NaN);
    if Term = 0.0 then
      Exit(SumValue);

    AdjustedTerm := Term - Compensation;
    UpdatedSum := SumValue + AdjustedTerm;
    Compensation := (UpdatedSum - SumValue) - AdjustedTerm;
    SumValue := UpdatedSum;
    if IsNan(SumValue) or IsInfinite(SumValue) then
      Exit(NaN);

    { For N >= 128, the largest possible next-term ratio is below 0.94
      under the public parameter bounds, so the omitted tail is below 17 terms
      of the current magnitude. }
    if (N >= MinimumStoppingIteration) and
      (Abs(Term) <= RelativeTermTolerance * Max(1.0, Abs(SumValue))) then
      Exit(SumValue);
  end;
  Result := NaN;
end;

end.