(FILECREATED "24-Feb-86 18:27:28" {DSK}<LISPFILES>PROBFNCTS.;2 30062
changes to: (FNS BETA)
previous date: "22-Feb-86 17:49:58"
{MITFS1-E40:SLOAN% SCHOOL:MASSINSTTECH}<EDAN% KABATCHNIK>PROBFNCTS.;7)
(* Copyright (c) 1986 by Massachusetts Institute of Technology. All rights reserved.)
(PRETTYCOMPRINT PROBFNCTSCOMS)
(RPAQQ PROBFNCTSCOMS ((FNS ALBETA ALNGAM ALNREL BETA CSEVL GAMLOG GAMMA HORNER INITS PBETA PCAUC
PCHIS PEXPT PFDIS PGAMM PLNORM PLOGI PNORM PNORMS PTDIS PUNIF R91GMC)
(CONSTANTS (\IEEE.SHORT.MACHEPS (EXPT 2.0 -23))
(\IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126))
(\IEEE.SHORT.BASE 2)
(\IEEE.SHORT.MANTISSA 23)
(\IEEE.SHORT.EXPONENT 8)
(\IEEE.SHORT.LEAST.EXPONENT -126)
(\IEEE.SHORT.GREATEST.EXPONENT 127)
(\IEEE.SHORT.EXPONENT.BIAS 127)
(\IEEE.SHORT.PI 3.141593))
(MACROS HORNER)))
(DEFINEQ
(ALBETA
[LAMBDA (A B) (* SCP "24-Feb-86 16:04")
(SETQ A (FLOAT A))
(SETQ B (FLOAT B))
(LET [(SmallerArgument (MIN A B))
(LargerArgument (MAX A B))
(NaturalLogSquareRootTwoPi (CONSTANT (LOG (SQRT (TIMES 2.0 \IEEE.SHORT.PI]
(if (LEQ SmallerArgument 0.0)
then (ERROR "ERROR in function ALBETA - Both arguments must be positive.")
elseif (GEQ SmallerArgument 10.0)
then [LET ((SumOfArguments (PLUS SmallerArgument LargerArgument)))
(PLUS (TIMES -.5 (LOG LargerArgument))
NaturalLogSquareRootTwoPi
(PLUS (R91GMC SmallerArgument)
(R91GMC LargerArgument)
(MINUS (R91GMC SumOfArguments)))
(TIMES (DIFFERENCE SmallerArgument .5)
(LOG (QUOTIENT SmallerArgument SumOfArguments)))
(TIMES LargerArgument (ALNREL (QUOTIENT (MINUS SmallerArgument)
SumOfArguments]
elseif (LESSP LargerArgument 10.0)
then [LOG (TIMES (GAMMA SmallerArgument)
(QUOTIENT (GAMMA LargerArgument)
(GAMMA (PLUS SmallerArgument LargerArgument]
else (LET ((SumOfArguments (PLUS SmallerArgument LargerArgument)))
(PLUS (ALNGAM SmallerArgument)
(DIFFERENCE (R91GMC LargerArgument)
(R91GMC SumOfArguments))
SmallerArgument
(TIMES (MINUS SmallerArgument)
(LOG SumOfArguments))
(TIMES (PLUS LargerArgument -.5)
(ALNREL (QUOTIENT (MINUS SmallerArgument)
SumOfArguments])
(ALNGAM
[LAMBDA (X) (* SCP "24-Feb-86 16:11")
(SETQ X (FLOAT X))
(LET [(AbsoluteX (ABS X))
[NaturalLogSquareRootPiOverTwo (CONSTANT (LOG (SQRT (QUOTIENT \IEEE.SHORT.PI 2.0]
[NaturalLogSquareRootTwoPi (CONSTANT (LOG (SQRT (TIMES 2.0 \IEEE.SHORT.PI]
[RelativeDX (CONSTANT (SQRT (TIMES 2.0 \IEEE.SHORT.MACHEPS]
(XMaximum (CONSTANT (QUOTIENT MAX.FLOAT (LOG MAX.FLOAT]
(if (LEQ AbsoluteX 10.0)
then (LOG (ABS (GAMMA X)))
elseif (GREATERP AbsoluteX XMaximum)
then (ERROR "ERROR in function ALNGAM - ABS(X) out of range.")
elseif (GREATERP X 0.0)
then (PLUS NaturalLogSquareRootTwoPi (TIMES (DIFFERENCE X .5)
(LOG X))
(MINUS X)
(R91GMC AbsoluteX))
else (LET [(SinePiTimesAbsoluteX (ABS (SIN (TIMES \IEEE.SHORT.PI AbsoluteX)
T]
(if (ZEROP SinePiTimesAbsoluteX)
then (ERROR "ERROR in function ALNGAM - X is a negative integer.")
elseif (LESSP (ABS (QUOTIENT (DIFFERENCE X (FIX (DIFFERENCE X .5)))
X))
RelativeDX)
then (ERROR "ERROR in function ALNGAM - Answer under half precision.")
else (PLUS NaturalLogSquareRootPiOverTwo (TIMES (DIFFERENCE X .5)
(LOG AbsoluteX))
(MINUS X)
(MINUS (LOG SinePiTimesAbsoluteX))
(MINUS (R91GMC AbsoluteX])
(ALNREL
[LAMBDA (X) (* SCP "24-Feb-86 16:13")
(SETQ X (FLOAT X))
(LET* [[Coefficients (CONSTANT (MAKEARRAY 23 (QUOTE INITIALCONTENTS)
(QUOTE (1.037869 -.133643 .1940825 -.003010755
.0004869461 -.00008105488 .00001377885
-2.380221E-6 4.164042E-7 -7.359583E-8
1.311761E-8 -2.354671E-9 4.252277E-10
-7.71909E-11 1.407575E-11
-2.576907E-12 4.734241E-13 -8.7249E-14
1.6124E-14 -2.987E-15 5.54E-16
-1.03E-16 1.9E-17]
[NumberOfChebyshevTerms (DEFERREDCONSTANT (INITS Coefficients 23 (TIMES .1 \IEEE.SHORT.MACHEPS]
(XMinimum (CONSTANT (PLUS -1.0 (SQRT (TIMES 2.0 \IEEE.SHORT.MACHEPS]
(if (LEQ X -1.0)
then (ERROR "ERROR in function ALNREL - X must be greater than -1.")
elseif (LESSP X XMinimum)
then (ERROR "ERROR in function ALNREL - Answer under half precision.")
elseif (LEQ (ABS X)
.375)
then [TIMES X (DIFFERENCE 1.0 (TIMES X (CSEVL (QUOTIENT X .375)
Coefficients NumberOfChebyshevTerms]
else (LOG (PLUS 1.0 X])
(BETA
[LAMBDA (X PInputValue QInputValue) (* SCP "24-Feb-86 18:22")
(SETQ X (FLOAT X))
(SETQ PInputValue (FLOAT PInputValue))
(SETQ QInputValue (FLOAT QInputValue))
(if (OR (MINUSP X)
(GREATERP X 1.0))
then (ERROR "ERROR in function BETA - Quantile must lie between 0 and 1.")
elseif (LEQ PInputValue 0.0)
then (ERROR "ERROR in function BETA - PInputValue must be positive.")
elseif (LEQ QInputValue 0.0)
then (ERROR "ERROR in function BETA - QInputValue must be positive.")
else
(LET* [(NaturalLogMachinePrecisionValue (CONSTANT (LOG \IEEE.SHORT.MACHEPS)))
(NaturalLogSmall (CONSTANT (LOG \IEEE.SHORT.LEAST.NORMAL)))
(PQXCondition1 (OR (AND (GREATERP QInputValue PInputValue)
(GEQ X .2))
(GEQ X .8)))
(Y (if PQXCondition1
then (DIFFERENCE 1.0 X)
else X))
(P (if PQXCondition1
then QInputValue
else PInputValue))
(Q (if PQXCondition1
then PInputValue
else QInputValue))
(PQXCondition2 (OR (NOT (EQP Y X))
(NOT (EQP P PInputValue]
(if (GREATERP \IEEE.SHORT.MACHEPS (QUOTIENT (TIMES (PLUS P Q)
Y)
(PLUS P 1.0)))
then (LET* ([XB (PLUS (TIMES P (LOG (FMAX Y Small)))
(MINUS (LOG P))
(MINUS (ALBETA P Q]
(A (if (AND (GREATERP XB NaturalLogSmall)
(NOT (EQP Y 0.0)))
then (ANTILOG XB)
else 0.0)))
(if PQXCondition2
then (DIFFERENCE 1.0 A)
else A))
else
(LET [(FinalResult
(LET* ((Betai (LET* [(F (LET [(B (DIFFERENCE Q (FIX Q]
(if (ZEROP B)
then 1.0
else B)))
(XB (PLUS (TIMES P (LOG Y))
(MINUS (ALBETA F P))
(MINUS (LOG P]
(if (GEQ XB NaturalLogSmall)
then (LET ((C (ANTILOG XB)))
(if (NOT (EQP F 1.0))
then (bind (TermValue ←(TIMES C P))
(LocalC ← C)
(LocalF ← F)
(LocalY ← Y)
(LocalP ← P)
(FloatI ← 1.0) for I from 1
to (FIX (FMAX (QUOTIENT
NaturalLogMachinePrecisionValue
(LOG Y))
4.0))
declare (TYPE FLOATP TermValue LocalC LocalF
LocalY LocalP FloatI)
do (SETQ TermValue (TIMES TermValue
(DIFFERENCE FloatI
LocalF)
(FQUOTIENT LocalY
FloatI)))
[SETQ LocalC
(PLUS LocalC (FQUOTIENT TermValue
(PLUS LocalP
FloatI]
(SETQ FloatI (PLUS FloatI 1.0))
finally (RETURN LocalC))
else C))
else 0.0)))
(TemporaryResult (if (GREATERP Q 1.0)
then [LET* [[XB (PLUS (TIMES P (LOG Y))
(TIMES Q (LOG (DIFFERENCE 1.0 Y)))
(MINUS (ALBETA P Q))
(MINUS (LOG Q]
(IB (FIX (FMAX (QUOTIENT XB NaturalLogSmall)
0.0)))
[TermValue (ANTILOG (DIFFERENCE XB (TIMES IB
NaturalLogSmall]
(D (QUOTIENT 1.0 (DIFFERENCE 1.0 Y)))
[ModifiedP (TIMES Q (QUOTIENT D (PLUS P Q
(MINUS 1.0]
(FiniteSum 0.0)
(E (if (EQP (FLOAT (FIX Q))
Q)
then (DIFFERENCE Q 1)
else (FIX Q]
(for I from 1 to E
do [SETQ TermValue
(TIMES (DIFFERENCE Q (DIFFERENCE I 1))
D
(QUOTIENT TermValue
(PLUS P Q (MINUS I]
(if (GREATERP TermValue 1.0)
then (SUB1 IB)
(SETQ TermValue (TIMES TermValue Small)
))
(if (ZEROP IB)
then (SETQ FiniteSum (PLUS FiniteSum
TermValue)))
repeatwhile (NOT (AND (LEQ ModifiedP 1.0)
(LEQ (QUOTIENT TermValue
MachinePrecisionValue)
FiniteSum)))
finally (RETURN (PLUS Betai FiniteSum]
else Betai)))
(if PQXCondition2
then (DIFFERENCE 1.0 TemporaryResult)
else TemporaryResult]
(FMAX (FMIN FinalResult 1.0)
0.0])
(CSEVL
[LAMBDA (X ChebyshevCoefs N) (* SCP "24-Feb-86 16:23")
(SETQ X (FLOAT X))
(SETQ N (FIX N))
(if (LESSP N 1)
then (ERROR "ERROR in function CSEVL - No terms to evaluate.")
elseif (GREATERP (ABS X)
1.0)
then (ERROR "ERROR in function CSEVL - X must lie between -1 and 1.")
else (bind (B0 ← 0.0)
(B1 ← 0.0)
B2
(TWOX ←(PLUS X X)) declare (TYPE FLOATP B0 B1 B2 TWOX) for I from (SUB1 N)
to 0 by -1
do (SETQ B2 B1)
(SETQ B1 B0)
(SETQ B0 (PLUS (TIMES TWOX B1)
(MINUS B2)
(AREF ChebyshevCoefs I)))
finally (RETURN (TIMES .5 (DIFFERENCE B0 B2])
(GAMLOG
[LAMBDA (Eta) (* SCP "24-Feb-86 16:28")
(SETQ Eta (FLOAT Eta))
(if (LEQ Eta 0.0)
then (ERROR "ERROR in function GAMLOG - Eta must be positive.")
elseif (LEQ Eta 33.0)
then (LOG (GAMMA Eta))
else (LET ((EtaMinusOne (DIFFERENCE Eta 1.0)))
(PLUS [CONSTANT (LOG (SQRT (TIMES 2.0 \IEEE.SHORT.PI]
(TIMES (PLUS EtaMinusOne .5)
(LOG EtaMinusOne))
(MINUS EtaMinusOne)
(QUOTIENT 1.0 (TIMES 12.0 EtaMinusOne))
(QUOTIENT 1.0 (TIMES 288.0 EtaMinusOne EtaMinusOne])
(GAMMA
[LAMBDA (Eta) (* SCP "24-Feb-86 15:57")
(SETQ Eta (FLOAT Eta))
(if (LEQ Eta 0.0)
then (ERROR "ERROR in function GAMMA - Eta must be positive.")
elseif (GREATERP Eta 33.0)
then (ERROR "ERROR in function GAMMA - Eta cannot be greater than 33.")
else (LET* ([Coefficients (CONSTANT (MAKEARRAY 9 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (.03586834 -.1935278 .4821994 -.7567041
.918207 -.897057 .988206
-.5771917 1.0]
(FEta (DIFFERENCE Eta (FIX Eta)))
(GLFEta (if (GREATERP FEta 0.0)
then (HORNER FEta Coefficients 0 8)
else 1.0)))
(if (LESSP Eta 1.0)
then (QUOTIENT GLFEta FEta)
elseif (LEQ Eta 2.0)
then GLFEta
else (LET ((Product 1.0)
(SeriesTerm (PLUS 1.0 FEta))
(TestValue (DIFFERENCE Eta .5)))
(do (SETQ Product (TIMES Product SeriesTerm))
(SETQ SeriesTerm (PLUS SeriesTerm 1.0)) repeatwhile (LESSP SeriesTerm
TestValue)
finally (RETURN (TIMES Product GLFEta])
(HORNER
[LAMBDA (X Polynomial Beginning Ending) (* SCP "24-Feb-86 15:59")
(if (EQP Beginning Ending)
then (AREF Polynomial Ending)
else (PLUS (AREF Polynomial Ending)
(TIMES X (HORNER X Polynomial Beginning (SUB1 Ending])
(INITS
[LAMBDA (OrthCoefs N Eta) (* SCP "24-Feb-86 16:40")
(SETQ N (FIX N))
(SETQ Eta (FLOAT Eta))
(if (LESSP N 1)
then (ERROR "ERROR in function INITS - No terms to evaluate.")
else (bind (Error ← 0.0) for I from N to 1 by -1 do [SETQ Error
(PLUS Error (ABS (AREF OrthCoefs
(SUB1 I]
repeatuntil (GREATERP Error Eta) finally (if (EQP I N)
then (ERROR
"ERROR in function INITS - Eta is too small.")
else (RETURN I])
(PBETA
[LAMBDA (X A B) (* SCP "24-Feb-86 15:57")
(SETQ X (FLOAT X))
(SETQ A (if (NULL A)
then 1.0
else (FLOAT A)))
(SETQ B (if (NULL B)
then 1.0
else (FLOAT B)))
(if (OR (LESSP X 0.0)
(GREATERP X 1.0))
then (ERROR "ERROR in function PBETA - Quantile must lie between 0 and 1.")
elseif (LEQ A 0.0)
then (ERROR "ERROR in function PBETA - A must be positive.")
elseif (LEQ B 0.0)
then (ERROR "ERROR in function PBETA - B must be positive.")
else (BETA X A B])
(PCAUC
[LAMBDA (X Mu Sigma) (* SCP "24-Feb-86 16:42")
(SETQ X (FLOAT X))
(SETQ Mu (if (NULL Mu)
then 0.0
else (FLOAT Mu)))
(SETQ Sigma (if (NULL Sigma)
then 1.0
else (FLOAT Sigma)))
(if (LEQ Sigma 0.0)
then (ERROR "ERROR in function PCAUC - Sigma must be positive.")
else (PLUS (QUOTIENT (ARCTAN2 (DIFFERENCE X Mu)
Sigma T)
\IEEE.SHORT.PI)
.5])
(PCHIS
[LAMBDA (X DegreesOfFreedom) (* SCP "24-Feb-86 15:57")
(SETQ X (FLOAT X))
(SETQ DegreesOfFreedom (if (NULL DegreesOfFreedom)
then 1.0
else (FLOAT DegreesOfFreedom)))
(if (MINUSP X)
then (ERROR "ERROR in function PCHIS - Quantile must be positive.")
elseif (LEQ DegreesOfFreedom 0.0)
then (ERROR "ERROR in function PCHIS - DegreesOfFreedom must be positive.")
elseif (EQP X 0.0)
then 0.0
else (PGAMM (QUOTIENT X 2.0)
(QUOTIENT DegreesOfFreedom 2.0])
(PEXPT
[LAMBDA (X Mu) (* SCP "24-Feb-86 15:56")
(SETQ X (FLOAT X))
(SETQ Mu (if (NULL Mu)
then 1.0
else (FLOAT Mu)))
(if (LEQ Mu 0.0)
then (ERROR "ERROR in function PEXPT - Mu must be positive.")
elseif (LEQ X 0.0)
then 0.0
else (DIFFERENCE 1.0 (ANTILOG (MINUS (QUOTIENT X Mu])
(PFDIS
[LAMBDA (X NumeratorDF DenominatorDF) (* SCP "24-Feb-86 15:57")
(SETQ X (FLOAT X))
(SETQ NumeratorDF (if (NULL NumeratorDF)
then 1.0
else (FLOAT NumeratorDF)))
(SETQ DenominatorDF (if (NULL DenominatorDF)
then 1.0
else (FLOAT DenominatorDF)))
(if (MINUSP X)
then (ERROR "ERROR in function PFDIS - Quantile must be positive.")
elseif (LEQ NumeratorDF 0.0)
then (ERROR "ERROR in function PFDIS - NumeratorDF must be positive.")
elseif (LEQ DenominatorDF 0.0)
then (ERROR "ERROR in function PFDIS - DenominatorDF must be positive.")
elseif (EQP X 0.0)
then 0.0
elseif (LET [(TemporaryResult (TIMES X (QUOTIENT NumeratorDF DenominatorDF]
(PBETA (QUOTIENT TemporaryResult (PLUS 1.0 TemporaryResult))
(QUOTIENT NumeratorDF 2.0)
(QUOTIENT DenominatorDF 2.0])
(PGAMM
[LAMBDA (X Eta) (* SCP "24-Feb-86 16:47")
(SETQ X (FLOAT X))
(SETQ Eta (if (NULL Eta)
then 1.0
else (FLOAT Eta)))
(if (LEQ X 0.0)
then (ERROR "ERROR in function PGAMM - Quantile must be positive.")
elseif (LEQ Eta 0.0)
then (ERROR "ERROR in function PGAMM - Eta must be positive.")
else (LET ((CopyOfEta Eta)
(SeriesTerm 1.0)
(TotalSum 1.0))
(if (AND (GEQ X Eta)
(GEQ X 7.0))
then [for I from 1 to (if (LESSP X 11.0)
then (FIX (PLUS X Eta -1))
else (FIX (PLUS Eta 10)))
do (SETQ CopyOfEta (DIFFERENCE CopyOfEta 1.0))
(SETQ SeriesTerm (QUOTIENT (TIMES SeriesTerm CopyOfEta)
X))
(SETQ TotalSum (PLUS TotalSum SeriesTerm))
finally (SETQ SeriesTerm (TIMES SeriesTerm (DIFFERENCE CopyOfEta 1.0)))
[SETQ TotalSum (PLUS TotalSum (QUOTIENT SeriesTerm
(PLUS X (MINUS CopyOfEta)
2.0]
(RETURN (DIFFERENCE 1.0 (ANTILOG (PLUS (MINUS X)
(TIMES (DIFFERENCE Eta 1.0)
(LOG X))
(LOG TotalSum)
(MINUS (GAMLOG Eta]
else (bind Scratch
do (SETQ CopyOfEta (PLUS CopyOfEta 1.0))
(SETQ Scratch (QUOTIENT X CopyOfEta))
(SETQ SeriesTerm (TIMES SeriesTerm Scratch))
(SETQ TotalSum (PLUS TotalSum SeriesTerm))
repeatwhile (OR (GREATERP Scratch .5)
(GREATERP (QUOTIENT SeriesTerm TotalSum)
(CONSTANT \IEEE.SHORT.MACHEPS)))
finally (RETURN (QUOTIENT [ANTILOG (PLUS (TIMES Eta (LOG X))
(MINUS X)
(LOG TotalSum)
(MINUS (GAMLOG Eta]
Eta])
(PLNORM
[LAMBDA (X Mu Sigma) (* SCP "24-Feb-86 15:57")
(SETQ X (FLOAT X))
(SETQ Mu (if (NULL Mu)
then 0.0
else (FLOAT Mu)))
(SETQ Sigma (if (NULL Sigma)
then 1.0
else (FLOAT Sigma)))
(if (LESSP X 0.0)
then (ERROR "ERROR in function PLNORM - Quantile must be positive.")
elseif (LEQ Sigma 0.0)
then (ERROR "ERROR in function PLNORM - Sigma must be positive.")
elseif (EQP X 0.0)
then 0.0
else (PNORMS (QUOTIENT (DIFFERENCE (LOG X)
Mu)
Sigma])
(PLOGI
[LAMBDA (X Mu Sigma) (* SCP "24-Feb-86 15:55")
(SETQ X (FLOAT X))
(SETQ Mu (if (NULL Mu)
then 0.0
else (FLOAT Mu)))
(SETQ Sigma (if (NULL Sigma)
then 1.0
else (FLOAT Sigma)))
(if (LEQ Sigma 0.0)
then (ERROR "ERROR in function PLOGI - Sigma must be positive.")
else (QUOTIENT 1.0 (PLUS 1.0 (ANTILOG (QUOTIENT (DIFFERENCE Mu X)
Sigma])
(PNORM
[LAMBDA (X Mu Sigma) (* SCP "24-Feb-86 15:55")
(SETQ X (FLOAT X))
(SETQ Mu (if (NULL Mu)
then 0.0
else (FLOAT Mu)))
(SETQ Sigma (if (NULL Sigma)
then 1.0
else (FLOAT Sigma)))
(if (LEQ Sigma 0.0)
then (ERROR "ERROR in function PNORM - Sigma must be positive.")
else (PNORMS (QUOTIENT (DIFFERENCE X Mu)
Sigma])
(PNORMS
[LAMBDA (X) (* SCP "24-Feb-86 16:49")
(SETQ X (FLOAT X))
(LET [[Bounds (CONSTANT (QUOTE (.625 1.25 2.0 2.45 3.5 4.62]
[Centers (CONSTANT (QUOTE (.3 .925 1.625 2.225 2.95 4.15]
(ListOfCoefficients (CONSTANT (LIST (MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (.000679824 .005102864 -.006016087
-.02272561 .03355577
.07270331 -.140939 -.1546891
.5156305 .1643134)))
(MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (-.001270975 -.001482265 .008771896
-.002020042 -.03468168
.0476424 .05685029
-.2218062 .2397904 .4045884)
))
(MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (.0006796452 -.0003458989 -.00308834
.007214109 -.001067704
-.02485975 .05742032
-.0653837 .04023613
.4892219)))
(MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (-.00008850032 .0005682608 -.0004501553
-.001387239 .005542219
-.01022113 .0118502
-.008886404 .003993889
.4991742)))
(MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (1.479155E-7 -.00006196501 .0002226671
-.0004439002 .0006297875
-.0006638263 .0005126568
-.0002765716 .00009375142
.4999849)))
(MAKEARRAY 10 (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)
(QUOTE INITIALCONTENTS)
(QUOTE (5.487908E-7 -8.665614E-7 6.687777E-7
-6.408465E-7 6.247178E-7
-4.210501E-7 2.074213E-7
-7.685604E-8 1.872202E-8 .5]
(if (EQP X 0.0)
then .5
else (LET [(XAdjustedBySqrtOf2 (ABS (QUOTIENT X (CONSTANT (SQRT 2]
(as Bound in Bounds as Center in Centers as Coefficients in ListOfCoefficients
thereis (LEQ XAdjustedBySqrtOf2 Bound)
finally (RETURN (if (NULL $$VAL)
then (if (LESSP X 0.0)
then 0.0
else 1.0)
else (LET ((FinalResult (HORNER (DIFFERENCE XAdjustedBySqrtOf2
Center)
Coefficients 0 9)))
(PLUS .5 (if (LESSP X 0.0)
then (MINUS FinalResult)
else FinalResult])
(PTDIS
[LAMBDA (X DegreesOfFreedom) (* SCP "24-Feb-86 16:55")
(SETQ X (FLOAT X))
(SETQ DegreesOfFreedom (if (NULL DegreesOfFreedom)
then 1.0
else (FLOAT DegreesOfFreedom)))
(if (LEQ DegreesOfFreedom 0.0)
then (ERROR "ERROR in function PTDIS - DegreesOfFreedom must be positive."))
(LET ((IntegerDegreesOfFreedom (FIX DegreesOfFreedom))
(XSquared (TIMES X X)))
(if (OR (GEQ (ABS (DIFFERENCE DegreesOfFreedom IntegerDegreesOfFreedom))
1.0E-6)
(AND (GEQ IntegerDegreesOfFreedom 20)
(LESSP XSquared DegreesOfFreedom))
(GREATERP IntegerDegreesOfFreedom 200))
then
(LET* [(A (QUOTIENT XSquared DegreesOfFreedom))
(B (DIFFERENCE DegreesOfFreedom .5))
(C (TIMES 48.0 B B))
(D (TIMES B (if (GREATERP A 1.0E-6)
then (LOG (PLUS 1.0 A))
else A)))
(FinalResult
(PNORMS
(MINUS
(TIMES
(PLUS (QUOTIENT (PLUS (QUOTIENT
(DIFFERENCE (TIMES (DIFFERENCE
(TIMES (DIFFERENCE (TIMES -.4 D)
3.3)
D)
24.0)
D)
85.5)
(PLUS (TIMES .8 D D)
100.0 C))
D 3.0)
C)
1.0)
(SQRT D]
(if (LESSP X 0.0)
then FinalResult
else (DIFFERENCE 1.0 FinalResult)))
else (LET* (A (B DegreesOfFreedom)
(C (QUOTIENT XSquared DegreesOfFreedom))
(D (PLUS 1.0 C))
(E 1.0))
[if (AND (LESSP IntegerDegreesOfFreedom 20)
(LESSP XSquared 4.0))
then (SETQ A (SQRT C))
(SETQ C A)
(if (EQP IntegerDegreesOfFreedom 1)
then (SETQ A 0.0))
else (SETQ A (SQRT D))
(SETQ C (TIMES A DegreesOfFreedom))
(bind ((F ← 0))
do (SETQ F (PLUS F 2))
(SETQ E A)
(SETQ C (QUOTIENT (TIMES C (DIFFERENCE F 1))
(TIMES D F)))
[SETQ A (PLUS A (QUOTIENT C (PLUS B F]
repeatwhile (GREATERP (ABS (DIFFERENCE A E))
Small)
finally (SETQ B (PLUS B 2.0))
(SETQ C 0.0)
(SETQ E 0.0)
(SETQ A (MINUS A]
(bind ((OldA ← A)
(OnePlusMachinePrecisionValue ←(PLUS 1.0 \IEEE.SHORT.MACHEPS)))
do (SETQ B (DIFFERENCE B 2.0))
(SETQ OldA A)
[SETQ A (PLUS C (TIMES A (QUOTIENT (DIFFERENCE B 1.0)
(TIMES D B]
repeatwhile (GREATERP B OnePlusMachinePrecisionValue) finally (SETQ A OldA))
(LET ((FinalResult (TIMES [DIFFERENCE E (if (GREATERP (ABS B)
\IEEE.SHORT.LEAST.NORMAL)
then (TIMES (PLUS (ARCTAN2 C 1 T)
(QUOTIENT A D))
(CONSTANT (QUOTIENT 2.0
\IEEE.SHORT.PI)))
else (QUOTIENT A (SQRT D]
.5)))
(if (LESSP X 0.0)
then FinalResult
else (DIFFERENCE 1.0 FinalResult])
(PUNIF
[LAMBDA (X BoundLow BoundHigh) (* SCP "24-Feb-86 15:56")
(SETQ X (FLOAT X))
(SETQ BoundLow (if (NULL BoundLow)
then 0.0
else (FLOAT BoundLow)))
(SETQ BoundHigh (if (NULL BoundHigh)
then 1.0
else (FLOAT BoundHigh)))
(if (LEQ BoundHigh BoundLow)
then (ERROR "ERROR in function PUNIF - No effective calculation range.")
elseif (OR (LEQ X BoundLow)
(LESSP BoundHigh X))
then 0.0
else (QUOTIENT (DIFFERENCE X BoundLow)
(DIFFERENCE BoundHigh BoundLow])
(R91GMC
[LAMBDA (X) (* SCP "24-Feb-86 16:56")
(SETQ X (FLOAT X))
(LET* [[Coefficients (CONSTANT (MAKEARRAY 6 (QUOTE INITIALCONTENTS)
(QUOTE (.1666389 -.00001384948 9.810826E-9 -1.80912E-11
6.22E-14 -3.0E-16]
[CutOffPoint (CONSTANT (QUOTIENT 1.0 (SQRT \IEEE.SHORT.MACHEPS]
(NumberOfChebyshevTerms (DEFERREDCONSTANT (INITS Coefficients 6 \IEEE.SHORT.MACHEPS)))
(XMaximum (CONSTANT (QUOTIENT MAX.FLOAT 12.0]
(if (LESSP X 10.0)
then (ERROR "ERROR in function R91GMC - X must be greater/equal to 10.")
elseif (GEQ X XMaximum)
then (ERROR "ERROR in function R91GMC - X out of range.")
elseif (LESSP X CutOffPoint)
then (LET ((TemporaryResult (QUOTIENT 10.0 X)))
(QUOTIENT (CSEVL (DIFFERENCE (TIMES TemporaryResult TemporaryResult 2.0)
1.0)
Coefficients NumberOfChebyshevTerms)
X))
else (QUOTIENT 1.0 (TIMES 12.0 X])
)
(DECLARE: EVAL@COMPILE
(RPAQ \IEEE.SHORT.MACHEPS (EXPT 2.0 -23))
(RPAQ \IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126))
(RPAQQ \IEEE.SHORT.BASE 2)
(RPAQQ \IEEE.SHORT.MANTISSA 23)
(RPAQQ \IEEE.SHORT.EXPONENT 8)
(RPAQQ \IEEE.SHORT.LEAST.EXPONENT -126)
(RPAQQ \IEEE.SHORT.GREATEST.EXPONENT 127)
(RPAQQ \IEEE.SHORT.EXPONENT.BIAS 127)
(RPAQQ \IEEE.SHORT.PI 3.141593)
(CONSTANTS (\IEEE.SHORT.MACHEPS (EXPT 2.0 -23))
(\IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126))
(\IEEE.SHORT.BASE 2)
(\IEEE.SHORT.MANTISSA 23)
(\IEEE.SHORT.EXPONENT 8)
(\IEEE.SHORT.LEAST.EXPONENT -126)
(\IEEE.SHORT.GREATEST.EXPONENT 127)
(\IEEE.SHORT.EXPONENT.BIAS 127)
(\IEEE.SHORT.PI 3.141593))
)
(DECLARE: EVAL@COMPILE
(PUTPROPS HORNER DMACRO [(X Polynomial Beginning Ending)
(\FLOATBOX ((OPCODES UBFLOAT3 0)
(\FLOATUNBOX X)
(\GETBASEPTR Polynomial Beginning)
(DIFFERENCE Ending Beginning])
)
(PUTPROPS PROBFNCTS COPYRIGHT ("Massachusetts Institute of Technology" 1986))
(DECLARE: DONTCOPY
(FILEMAP (NIL (874 29000 (ALBETA 884 . 2600) (ALNGAM 2602 . 4204) (ALNREL 4206 . 5473) (BETA 5475 .
10208) (CSEVL 10210 . 11006) (GAMLOG 11008 . 11679) (GAMMA 11681 . 12993) (HORNER 12995 . 13299) (
INITS 13301 . 13953) (PBETA 13955 . 14644) (PCAUC 14646 . 15191) (PCHIS 15193 . 15835) (PEXPT 15837 .
16285) (PFDIS 16287 . 17312) (PGAMM 17314 . 19324) (PLNORM 19326 . 20007) (PLOGI 20009 . 20549) (PNORM
20551 . 21054) (PNORMS 21056 . 23900) (PTDIS 23902 . 27235) (PUNIF 27237 . 27903) (R91GMC 27905 .
28998)))))
STOP