MODULE MathMitLef;
IMPORT NbrInt, NbrRe, NbrCplx, DataErrors, MathRe, MathGamma, MathCplx, MathCplxSeries, CalcGauss;
VAR
maxIterations: NbrInt.Integer; storeAlpha, storeBeta, tolerance: NbrRe.Real;
yk: ARRAY 2 OF NbrRe.Real;
yp: ARRAY 3 OF NbrRe.Real;
TYPE
MitLef = OBJECT (MathCplxSeries.Coefficient)
PROCEDURE Evaluate;
VAR x: NbrRe.Real;
BEGIN
x := storeBeta + n * storeAlpha;
IF GammaIsSingularAt( x ) THEN coef := 0 ELSE coef := 1 / MathGamma.Fn( x ) END;
IF n = maxIterations THEN eos := TRUE; DataErrors.ReWarning( x, "Did not converge - timed out." ) END
END Evaluate;
END MitLef;
AsympMitLef = OBJECT (MathCplxSeries.Coefficient)
PROCEDURE Evaluate;
VAR x: NbrRe.Real;
BEGIN
x := storeBeta - n * storeAlpha;
IF GammaIsSingularAt( x ) THEN coef := 0 ELSE coef := 1 / MathGamma.Fn( x ) END;
IF n = maxIterations THEN eos := TRUE; DataErrors.ReWarning( x, "Did not converge - timed out." ) END
END Evaluate;
END AsympMitLef;
DMitLef = OBJECT (MathCplxSeries.Coefficient)
PROCEDURE Evaluate;
VAR x: NbrRe.Real;
BEGIN
x := (1 + n) * storeAlpha + storeBeta;
IF GammaIsSingularAt( x ) THEN coef := 0 ELSE coef := (1 + n) / MathGamma.Fn( x ) END;
IF n = maxIterations THEN eos := TRUE; DataErrors.ReWarning( x, "Did not converge - timed out." ) END
END Evaluate;
END DMitLef;
PROCEDURE K( x: NbrRe.Real; z: NbrCplx.Complex ): NbrCplx.Complex;
VAR coef1, coef2, r1: NbrRe.Real; answer, denom, numer: NbrCplx.Complex;
BEGIN
coef1 := MathRe.Power( x, (1 - yk[1]) / yk[0] ) / (NbrRe.Pi * yk[0]);
coef2 := MathRe.Exp( -MathRe.Power( x, 1 / yk[0] ) ); r1 := NbrRe.Pi * (1 - yk[1]);
numer := x * MathRe.Sin( r1 ) - z * MathRe.Sin( r1 + NbrRe.Pi * yk[0] );
denom := x * x - 2 * x * z * MathRe.Cos( NbrRe.Pi * yk[0] ) + z * z; answer := coef1 * coef2 * numer / denom;
RETURN answer
END K;
PROCEDURE P( x: NbrRe.Real; z: NbrCplx.Complex ): NbrCplx.Complex;
VAR coef1, coef2, r1: NbrRe.Real; answer, denom, numer: NbrCplx.Complex;
BEGIN
coef1 := MathRe.Power( yp[2], 1 + (1 - yp[1]) / yp[0] ) / (2 * NbrRe.Pi * yp[0]);
coef2 := MathRe.Exp( MathRe.Power( yp[2], 1 / yp[0] ) * MathRe.Cos( x / yp[0] ) );
r1 := x * (1 + (1 - yp[1]) / yp[0]) + MathRe.Power( yp[2], 1 / yp[0] ) * MathRe.Sin( x / yp[0] );
numer := MathRe.Cos( r1 ) + NbrCplx.I * MathRe.Sin( r1 ); denom := yp[2] * MathCplx.Exp( NbrCplx.I * x ) - z;
answer := coef1 * coef2 * numer / denom; RETURN answer
END P;
PROCEDURE GammaIsSingularAt( x: NbrRe.Real ): BOOLEAN;
BEGIN
IF x > 0 THEN RETURN FALSE
ELSIF x < 0 THEN
IF x = NbrRe.Int( x ) + 1 THEN RETURN TRUE ELSE RETURN FALSE END
ELSE RETURN TRUE
END
END GammaIsSingularAt;
PROCEDURE Fn*( alpha, beta, x: NbrRe.Real ): NbrRe.Real;
VAR abs, answer, arg: NbrRe.Real; ml, z: NbrCplx.Complex;
BEGIN
z := x; ml := CplxFn( alpha, beta, z ); abs := NbrCplx.Abs( ml ); arg := NbrCplx.Arg( ml );
IF (arg < tolerance) & (-tolerance < arg) THEN answer := abs
ELSIF (arg > NbrRe.Pi - tolerance) OR (arg < -NbrRe.Pi + tolerance) THEN answer := -abs
ELSE DataErrors.Error( "The result is complex. Use CplxFn instead." )
END;
RETURN answer
END Fn;
PROCEDURE CplxFn*( alpha, beta: NbrRe.Real; z: NbrCplx.Complex ): NbrCplx.Complex;
VAR i1, i2, k, k0, res: NbrInt.Integer; a, absZ, argZ, b, r1, r2, x0: NbrRe.Real;
answer, c1, c2: NbrCplx.Complex; ml: MitLef; aml: AsympMitLef;
BEGIN
alpha := NbrRe.Abs( alpha ); absZ := NbrCplx.Abs( z ); argZ := NbrCplx.Arg( z );
IF absZ < 1 / NbrRe.MaxNbr THEN
IF GammaIsSingularAt( beta ) THEN answer := 0 ELSE answer := 1 / MathGamma.Fn( beta ) END
ELSIF (alpha = 1) & (beta = 1) THEN
answer := MathCplx.Exp( z )
ELSIF 1 < alpha THEN
k0 := NbrRe.Floor( alpha ) + 1; c1 := (2 * NbrRe.Pi / k0) * NbrCplx.I; r1 := 1 / k0; r2 := alpha / k0; answer := 0;
FOR k := 0 TO k0 - 1 DO
c2 := MathCplx.RealPower( z, r1 ) * MathCplx.Exp( c1 * k ); answer := answer + CplxFn( r2, beta, c2 )
END;
answer := answer / k0
ELSE
IF absZ < 0.9 THEN BEGIN {EXCLUSIVE}
storeAlpha := alpha; storeBeta := beta; i1 := NbrRe.Ceiling( (1 - beta) / alpha );
i2 := NbrRe.Ceiling( MathRe.Ln( NbrRe.Epsilon * (1 - absZ) ) / MathRe.Ln( absZ ) );
maxIterations := NbrInt.Max( i1, i2 ); NEW( ml ); answer := MathCplxSeries.PowerSeries( ml, z ) END
ELSIF absZ < NbrRe.Floor( 10 + 5 * alpha ) THEN
IF beta >= 0 THEN
r1 := 1; x0 := NbrRe.Max( r1, 2 * absZ );
x0 := NbrRe.Max( x0, MathRe.Power( -MathRe.Ln( NbrRe.Epsilon * NbrRe.Pi / 6 ), alpha ) )
ELSE
r1 := NbrRe.Abs( beta ); r2 := MathRe.Power( 1 + r1, alpha ); x0 := NbrRe.Max( r2, 2 * absZ );
r2 := 6 * (2 + r1) * MathRe.Power( 2 * r1, r1 );
x0 := NbrRe.Max( x0, MathRe.Power( -2 * MathRe.Ln( NbrRe.Epsilon * NbrRe.Pi / r2 ), alpha ) )
END;
IF (NbrRe.Abs( argZ ) > alpha * NbrRe.Pi) &
(NbrRe.Abs(NbrRe.Abs( argZ ) - alpha * NbrRe.Pi) > NbrRe.Epsilon) THEN
IF beta < 1 + alpha THEN
BEGIN {EXCLUSIVE}
a := 0; b := x0; yk[0] := alpha; yk[1] := beta;
answer := CalcGauss.SolveCplx( K, a, b, z, CalcGauss.Medium, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating K." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating K." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating K." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END
ELSE
BEGIN {EXCLUSIVE}
a := 1; b := x0; yk[0] := alpha; yk[1] := beta;
answer := CalcGauss.SolveCplx( K, a, b, z, CalcGauss.Medium, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating K." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating K." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating K." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END; BEGIN {EXCLUSIVE}
a := -alpha * NbrRe.Pi; b := -a; yp[0] := alpha; yp[1] := beta; yp[2] := 1;
answer := answer + CalcGauss.SolveCplx( P, a, b, z, CalcGauss.Fine, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating P." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating P." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating P." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END
END
ELSIF (NbrRe.Abs( argZ ) < alpha * NbrRe.Pi) &
(NbrRe.Abs(NbrRe.Abs( argZ ) - alpha * NbrRe.Pi) > NbrRe.Epsilon) THEN
IF beta < 1 + alpha THEN
BEGIN {EXCLUSIVE}
a := 0; b := x0; yk[0] := alpha; yk[1] := beta;
answer := CalcGauss.SolveCplx( K, a, b, z, CalcGauss.Medium, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating K." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating K." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating K." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END
ELSE
BEGIN {EXCLUSIVE}
a := absZ / 2; b := x0; yk[0] := alpha; yk[1] := beta;
answer := CalcGauss.SolveCplx( K, a, b, z, CalcGauss.Medium, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating K." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating K." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating K." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END; BEGIN {EXCLUSIVE}
a := -alpha * NbrRe.Pi; b := -a; yp[0] := alpha; yp[1] := beta; yp[2] := absZ / 2;
answer := answer + CalcGauss.SolveCplx( P, a, b, z, CalcGauss.Fine, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating P." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating P." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating P." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END
END;
answer :=
answer +
MathCplx.RealPower( z, (1 - beta) / alpha ) *
MathCplx.Exp( MathCplx.RealPower( z, 1 / alpha ) ) / alpha
ELSE
BEGIN {EXCLUSIVE}
a := (absZ + 1) / 2; b := x0; yk[0] := alpha; yk[1] := beta;
answer := CalcGauss.SolveCplx( K, a, b, z, CalcGauss.Medium, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating K." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating K." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating K." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END; BEGIN {EXCLUSIVE}
a := -alpha * NbrRe.Pi; b := -a; yp[0] := alpha; yp[1] := beta; yp[2] := (absZ + 1) / 2;
answer := answer + CalcGauss.SolveCplx( P, a, b, z, CalcGauss.Fine, tolerance, res ) END;
IF res # CalcGauss.OKay THEN
IF res = CalcGauss.MaxSubDivReached THEN
DataErrors.Warning( "Maximum subdivisions reached when integrating P." )
ELSIF res = CalcGauss.RoundoffError THEN
DataErrors.Warning( "Excessive roundoff error occurred when integrating P." )
ELSIF res = CalcGauss.RoughIntegrand THEN
DataErrors.Warning( "A rough integrand encountered when integrating P." )
ELSE DataErrors.Error( "Unknown error originating from CalcGauss.SolveCplx." )
END
END
END
ELSE
IF argZ < 3 * alpha * NbrRe.Pi / 4 THEN
answer :=
MathCplx.RealPower( z, (1 - beta) / alpha ) *
MathCplx.Exp( MathCplx.RealPower( z, 1 / alpha ) ) / alpha
ELSE answer := 0
END; BEGIN {EXCLUSIVE}
storeAlpha := alpha; storeBeta := beta;
maxIterations := 2 * NbrRe.Floor( -MathRe.Ln( NbrRe.Epsilon ) / MathRe.Ln( absZ ) ); NEW( aml );
answer := answer - MathCplxSeries.PowerSeries( aml, 1 / z );
IF ~GammaIsSingularAt( beta ) THEN answer := answer + 1 / MathGamma.Fn( beta ) END END
END
END;
RETURN answer
END CplxFn;
PROCEDURE DFn*( alpha, beta, x: NbrRe.Real ): NbrRe.Real;
VAR abs, answer, arg: NbrRe.Real; dml, z: NbrCplx.Complex;
BEGIN
z := x; dml := DCplxFn( alpha, beta, z ); abs := NbrCplx.Abs( dml ); arg := NbrCplx.Arg( dml );
IF (arg < tolerance) & (-tolerance < arg) THEN answer := abs
ELSIF (arg > NbrRe.Pi - tolerance) OR (arg < -NbrRe.Pi + tolerance) THEN answer := -abs
ELSE DataErrors.Error( "The result is complex. Use DCplxFn instead." )
END;
RETURN answer
END DFn;
PROCEDURE DCplxFn*( alpha, beta: NbrRe.Real; z: NbrCplx.Complex ): NbrCplx.Complex;
VAR answer: NbrCplx.Complex; absZ, d, k, k1, omega: NbrRe.Real; k0: NbrInt.Integer; dml: DMitLef;
BEGIN
alpha := NbrRe.Abs( alpha ); absZ := NbrCplx.Abs( z );
IF absZ < 1 / NbrRe.MaxNbr THEN
IF GammaIsSingularAt( alpha + beta ) THEN answer := 0 ELSE answer := 1 / MathGamma.Fn( alpha + beta ) END
ELSIF absZ < 0.9 THEN
IF alpha > 1 THEN k1 := 1 + (2 - alpha - beta) / (alpha - 1)
ELSE
d := 1 + alpha * (alpha - 4 * beta + 6); k := 1 + (3 - alpha - beta) / alpha;
IF d <= 1 THEN k1 := k
ELSE
omega := alpha + beta - 1.5;
k1 := NbrRe.Max( k, 1 + (1 - 2 * omega * alpha + MathRe.Sqrt( d )) / (2 * alpha * alpha) )
END
END; BEGIN {EXCLUSIVE}
k0 := NbrRe.Ceiling( NbrRe.Max( k1, MathRe.Ln( NbrRe.Epsilon * (1 - absZ) ) ) / MathRe.Ln( absZ ) );
maxIterations := k0; storeAlpha := alpha; storeBeta := beta; NEW( dml );
answer := MathCplxSeries.PowerSeries( dml, z ) END
ELSE
answer := CplxFn( alpha, beta - 1, z );
IF beta # 1 THEN answer := answer - (beta - 1) * CplxFn( alpha, beta, z ) END;
answer := answer / (alpha * z)
END;
RETURN answer
END DCplxFn;
BEGIN
tolerance := MathRe.Sqrt( NbrRe.Epsilon )
END MathMitLef.