summaryrefslogtreecommitdiff
path: root/packages/numlib/src/spe.pas
diff options
context:
space:
mode:
Diffstat (limited to 'packages/numlib/src/spe.pas')
-rw-r--r--packages/numlib/src/spe.pas1299
1 files changed, 1299 insertions, 0 deletions
diff --git a/packages/numlib/src/spe.pas b/packages/numlib/src/spe.pas
new file mode 100644
index 0000000000..3058525eb4
--- /dev/null
+++ b/packages/numlib/src/spe.pas
@@ -0,0 +1,1299 @@
+{
+ This file is part of the Numlib package.
+ Copyright (c) 1986-2000 by
+ Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
+ Computational centre of the Eindhoven University of Technology
+
+ FPC port Code by Marco van de Voort (marco@freepascal.org)
+ documentation by Michael van Canneyt (Michael@freepascal.org)
+
+ Special functions. (Currently, most of these functions only work for
+ ArbFloat=REAL/DOUBLE, not for Extended(at least not at full
+ precision, implement the tables in the diverse functions
+ for extended)
+
+ See the file COPYING.FPC, included in this distribution,
+ for details about the copyright.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+ **********************************************************************}
+unit spe;
+{$I DIRECT.INC}
+
+interface
+
+uses typ;
+
+{ Calculate modified Besselfunction "of the first kind" I0(x) }
+function spebi0(x: ArbFloat): ArbFloat;
+
+{ Calculate modified Besselfunction "of the first kind" I1(x) }
+function spebi1(x: ArbFloat): ArbFloat;
+
+{ Calculate Besselfunction "of the first kind" J0(x) }
+function spebj0(x: ArbFloat): ArbFloat;
+
+{ Calculate Besselfunction "of the first kind" J1(x) }
+function spebj1(x: ArbFloat): ArbFloat;
+
+{ Calculate modified Besselfunction "of the second kind" K0(x) }
+function spebk0(x: ArbFloat): ArbFloat;
+
+{ Calculate modified Besselfunction "of the second kind" K1(x) }
+function spebk1(x: ArbFloat): ArbFloat;
+
+{ Calculate Besselfunction "of the second kind" Y0(x) }
+function speby0(x: ArbFloat): ArbFloat;
+
+{ Calculate Besselfunction "of the second kind" Y1(x) }
+function speby1(x: ArbFloat): ArbFloat;
+
+{ Entier function, calculates first integer greater or equal than X}
+function speent(x: ArbFloat): longint;
+
+{ Errorfunction ( 2/sqrt(pi)* Int(t,0,pi,exp(sqr(t)) )}
+function speerf(x: ArbFloat): ArbFloat;
+
+{ Errorfunction's complement ( 2/sqrt(pi)* Int(t,pi,inf,exp(sqr(t)) )}
+function speefc(x: ArbFloat): ArbFloat;
+
+{ Function to calculate the Gamma function ( int(t,0,inf,t^(x-1)*exp(-t)) }
+function spegam(x: ArbFloat): ArbFloat;
+
+{ Function to calculate the natural logaritm of the Gamma function}
+function spelga(x: ArbFloat): ArbFloat;
+
+{ "Calculates" the maximum of two ArbFloat values }
+function spemax(a, b: ArbFloat): ArbFloat;
+
+{ Calculates the functionvalue of a polynomalfunction with n coefficients in a
+for variable X }
+function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
+
+{ Calc a^b with a and b real numbers}
+function spepow(a, b: ArbFloat): ArbFloat;
+
+{ Returns sign of x (-1 for x<0, 0 for x=0 and 1 for x>0) }
+function spesgn(x: ArbFloat): ArbInt;
+
+{ ArcSin(x) }
+function spears(x: ArbFloat): ArbFloat;
+
+{ ArcCos(x) }
+function spearc(x: ArbFloat): ArbFloat;
+
+{ Sinh(x) }
+function spesih(x: ArbFloat): ArbFloat;
+
+{ Cosh(x) }
+function specoh(x: ArbFloat): ArbFloat;
+
+{ Tanh(x) }
+function spetah(x: ArbFloat): ArbFloat;
+
+{ ArcSinH(x) }
+function speash(x: ArbFloat): ArbFloat;
+
+{ ArcCosh(x) }
+function speach(x: ArbFloat): ArbFloat;
+
+{ ArcTanH(x) }
+function speath(x: ArbFloat): ArbFloat;
+
+implementation
+
+function spebi0(x: ArbFloat): ArbFloat;
+
+const
+
+ xvsmall = 3.2e-9;
+ a1 : array[0..23] of ArbFloat =
+ ( 3.08508322553671039e-1, -1.86478066609466760e-1,
+ 1.57686843969995904e-1, -1.28895621330524993e-1,
+ 9.41616340200868389e-2, -6.04316795007737183e-2,
+ 3.41505388391452157e-2, -1.71317947935716536e-2,
+ 7.70061052263382555e-3, -3.12923286656374358e-3,
+ 1.15888319775791686e-3, -3.93934532072526720e-4,
+ 1.23682594989692688e-4, -3.60645571444886286e-5,
+ 9.81395862769787105e-6, -2.50298975966588680e-6,
+ 6.00566861079330132e-7, -1.36042013507151017e-7,
+ 2.92096163521178835e-8, -5.94856273204259507e-9,
+ 1.13415934215369209e-9, -2.10071360134551962e-10,
+ 4.44484446637868974e-11,-7.48150165756234957e-12);
+ a2 : array[0..26] of ArbFloat =
+ ( 1.43431781856850311e-1, -3.71571542566085323e-2,
+ 1.44861237337359455e-2, -6.30121694459896307e-3,
+ 2.89362046530968701e-3, -1.37638906941232170e-3,
+ 6.72508592273773611e-4, -3.35833513200679384e-4,
+ 1.70524543267970595e-4, -8.74354291104467762e-5,
+ 4.48739019580173804e-5, -2.28278155280668483e-5,
+ 1.14032404021741277e-5, -5.54917762110482949e-6,
+ 2.61457634142262604e-6, -1.18752840689765504e-6,
+ 5.18632519069546106e-7, -2.17653548816447667e-7,
+ 8.75291839187305722e-8, -3.34900221934314738e-8,
+ 1.24131668344616429e-8, -4.66215489983794905e-9,
+ 1.58599776268172290e-9, -3.80370174256271589e-10,
+ 1.23188158175419302e-10,-8.46900307934754898e-11,
+ 2.45185252963941089e-11);
+ a3: array[0..19] of ArbFloat =
+ ( 4.01071065066847416e-1, 2.18216817211694382e-3,
+ 5.59848253337377763e-5, 2.79770701849785597e-6,
+ 2.17160501061222148e-7, 2.36884434055843528e-8,
+ 3.44345025431425567e-9, 6.47994117793472057e-10,
+ 1.56147127476528831e-10, 4.82726630988879388e-11,
+ 1.89599322920800794e-11, 1.05863621425699789e-11,
+ 8.27719401266046976e-12, 2.82807056475555021e-12,
+ -4.34624739357691085e-12,-4.29417106720584499e-12,
+ 4.30812577328136192e-13, 1.44572313799118029e-12,
+ 4.73229306831831040e-14,-1.95679809047625728e-13);
+
+
+var t : ArbFloat;
+
+begin
+ t:=abs(x);
+ if t <=xvsmall
+ then
+ spebi0:=1
+ else
+ if t <= 4
+ then
+ spebi0 := exp(t)*spepol(t/2-1, a1[0], SizeOf(a1) div SizeOf(ArbFloat) -1)
+ else
+ if t <= 12
+ then
+ spebi0:=exp(t)*spepol(t/4-2, a2[0], SizeOf(a2) div SizeOf(ArbFloat) -1)
+ else { t > 12}
+ spebi0:=(exp(t)/sqrt(t))*
+ spepol(24/t-1, a3[0], SizeOf(a3) div SizeOf(ArbFloat) -1)
+end; {spebi0}
+
+function spebi1(x: ArbFloat): ArbFloat;
+
+
+const xvsmall = 3.2e-9;
+ a1: array[0..11] of ArbFloat =
+ ( 1.19741654963670236e+0, 9.28758890114609554e-1,
+ 2.68657659522092832e-1, 4.09286371827770484e-2,
+ 3.84763940423809498e-3, 2.45224314039278904e-4,
+ 1.12849795779951847e-5, 3.92368710996392755e-7,
+ 1.06662712314503955e-8, 2.32856921884663846e-10,
+ 4.17372709788222413e-12,6.24387910353848320e-14);
+
+ a2: array[0..26] of ArbFloat =
+ ( 1.34142493292698178e-1, -2.99140923897405570e-2,
+ 9.76021102528646704e-3, -3.40759647928956354e-3,
+ 1.17313412855965374e-3, -3.67626180992174570e-4,
+ 8.47999438119288094e-5, 5.21557319070236939e-6,
+ -2.62051678511418163e-5, 2.47493270133518925e-5,
+ -1.79026222757948636e-5, 1.13818992442463952e-5,
+ -6.63144162982509821e-6, 3.60186151617732531e-6,
+ -1.83910206626348772e-6, 8.86951515545183908e-7,
+ -4.05456611578551130e-7, 1.76305222240064495e-7,
+ -7.28978293484163628e-8, 2.84961041291017650e-8,
+ -1.07563514207617768e-8, 4.11321223904934809e-9,
+ -1.41575617446629553e-9, 3.38883570696523350e-10,
+ -1.10970391104678003e-10, 7.79929176497056645e-11,
+ -2.27061376122617856e-11);
+
+ a3: array[0..19] of ArbFloat =
+ ( 3.92624494204116555e-1, -6.40545360348237412e-3,
+ -9.12475535508497109e-5, -3.82795135453556215e-6,
+ -2.72684545741400871e-7, -2.82537120880041703e-8,
+ -3.96757162863209348e-9, -7.28107961041827952e-10,
+ -1.72060490748583241e-10,-5.23524129533553498e-11,
+ -2.02947854602758139e-11,-1.11795516742222899e-11,
+ -8.69631766630563635e-12,-3.05957293450420224e-12,
+ 4.42966462319664333e-12, 4.47735589657057690e-12,
+ -3.95353303949377536e-13,-1.48765082315961139e-12,
+ -5.77176811730370560e-14, 1.99448557598015488e-13);
+
+var t : ArbFloat;
+
+begin
+ t:=abs(x);
+ if t <= xvsmall
+ then
+ spebi1:=x/2
+ else
+ if t <= 4
+ then
+ spebi1:=x*spepol(sqr(t)/8-1, a1[0], sizeof(a1) div sizeof(ArbFloat)-1)
+ else
+ if t <= 12
+ then
+ spebi1:=
+ exp(t)*spepol(t/4-2, a2[0], sizeof(a2) div sizeof(ArbFloat)-1)*spesgn(x)
+ else { t > 12}
+ spebi1:=
+ (exp(t)/sqrt(t))*
+ spepol(24/t-1, a3[0], sizeof(a3) div sizeof(ArbFloat)-1)*spesgn(x)
+end; {spebi1}
+
+function spebj0(x: ArbFloat): ArbFloat;
+const
+
+ xvsmall = 3.2e-9;
+ tbpi = 6.36619772367581343e-1;
+ a1 : array[0..5] of ArbFloat =
+ ( 1.22200000000000000e-17, -1.94383469000000000e-12,
+ 7.60816359241900000e-8, -4.60626166206275050e-4,
+ 1.58067102332097261e-1, -8.72344235285222129e-3);
+
+ b1 : array[0..5] of ArbFloat =
+ ( - 7.58850000000000000e-16, 7.84869631400000000e-11,
+ - 1.76194690776215000e-6, 4.81918006946760450e-3,
+ - 3.70094993872649779e-1, 1.57727971474890120e-1);
+
+ c1 : array[0..4] of ArbFloat =
+ ( 4.12532100000000000e-14, - 2.67925353056000000e-9,
+ 3.24603288210050800e-5, - 3.48937694114088852e-2,
+ 2.65178613203336810e-1);
+
+ d1 : array[0..13] of ArbFloat =
+ ( 9.99457275788251954e-1, -5.36367319213004570e-4,
+ 6.13741608010926000e-6, -2.05274481565160000e-7,
+ 1.28037614434400000e-8, -1.21211819632000000e-9,
+ 1.55005642880000000e-10,-2.48827276800000000e-11,
+ 4.78702080000000000e-12,-1.06365696000000000e-12,
+ 2.45294080000000000e-13,-6.41843200000000000e-14,
+ 3.34028800000000000e-14,-1.17964800000000000e-14);
+
+ d2 : array[0..16] of ArbFloat =
+ ( -1.55551138795135187e-2, 6.83314909934390000e-5,
+ -1.47713883264594000e-6, 7.10621485930000000e-8,
+ -5.66871613024000000e-9, 6.43278173280000000e-10,
+ -9.47034774400000000e-11, 1.70330918400000000e-11,
+ -3.59094272000000000e-12, 8.59855360000000000e-13,
+ -2.28807680000000000e-13, 6.95193600000000000e-14,
+ -2.27942400000000000e-14, 4.75136000000000000e-15,
+ -1.14688000000000000e-15, 2.12992000000000000e-15,
+ -9.83040000000000000e-16);
+
+var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
+ i, bov : ArbInt;
+
+begin
+ t:=abs(x);
+ if t<=xvsmall
+ then
+ spebj0:=1
+ else
+ if t<=8
+ then
+ begin
+ t:=0.03125*sqr(t)-1; t2:=2*t;
+ b:=0; c:=0;
+ bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
+ for i:=0 to bov do
+ begin
+ a:=t2*c-b+a1[i];
+ if i<5
+ then
+ b:=t2*a-c+b1[i]
+ else
+ spebj0:=t*a-c+b1[i];
+ if i<bov
+ then
+ c:=t2*b-a+c1[i]
+ else
+ if i<5
+ then
+ spebj0:=t*b-a+c1[i]
+ end {i}
+ end {abs(x)<=8}
+ else
+ begin
+ g:=t-1/(2*tbpi); y:=sqrt(tbpi/t);
+ cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
+ t:=128/sqr(t)-1;
+ spebj0:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
+ + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
+ end {abs(x)>8}
+
+end {spebj0};
+
+function spebj1(x: ArbFloat): ArbFloat;
+const
+
+ xvsmall = 3.2e-9;
+ tbpi = 6.36619772367581343e-1;
+ a1 : array[0..5] of ArbFloat =
+ ( 2.95000000000000000e-18, -5.77740420000000000e-13,
+ 2.94970700727800000e-8, -2.60444389348580680e-4,
+ 1.77709117239728283e-1, -1.19180116054121687e+0);
+
+ b1 : array[0..5] of ArbFloat =
+ ( -1.95540000000000000e-16, 2.52812366400000000e-11,
+ -7.61758780540030000e-7, 3.24027018268385747e-3,
+ -6.61443934134543253e-1, 6.48358770605264921e-1);
+
+ c1 : array[0..4] of ArbFloat =
+ ( 1.13857200000000000e-14, -9.42421298160000000e-10,
+ 1.58870192399321300e-5, -2.91755248061542077e-2,
+ 1.28799409885767762e+0);
+
+ d1 : array[0..13] of ArbFloat =
+ ( 1.00090702627808217e+0, 8.98804941670557880e-4,
+ -7.95969469843846000e-6, 2.45367662227560000e-7,
+ -1.47085129889600000e-8, 1.36030580128000000e-9,
+ -1.71310758400000000e-10, 2.72040729600000000e-11,
+ -5.19113984000000000e-12, 1.14622464000000000e-12,
+ -2.63372800000000000e-13, 6.86387200000000000e-14,
+ -3.54508800000000000e-14, 1.24928000000000000e-14);
+
+ d2 : array[0..15] of ArbFloat =
+ ( 4.67768740274489776e-2, -9.62145882205441600e-5,
+ 1.82120185123076000e-6, -8.29196070929200000e-8,
+ 6.42013250344000000e-9, -7.15110504800000000e-10,
+ 1.03950931840000000e-10, -1.85248000000000000e-11,
+ 3.87554432000000000e-12, -9.23228160000000000e-13,
+ 2.50224640000000000e-13, -7.43936000000000000e-14,
+ 1.75718400000000000e-14, -4.83328000000000000e-15,
+ 5.32480000000000000e-15, -2.29376000000000000e-15);
+
+var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
+ i, bov : ArbInt;
+
+begin
+ t:=abs(x);
+ if t<xvsmall
+ then
+ spebj1:=x/2
+ else
+ if t<=8
+ then
+ begin
+ t:=0.03125*sqr(t)-1; t2:=2*t;
+ b:=0; c:=0;
+ bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
+ for i:=0 to bov do
+ begin
+ a:=t2*c-b+a1[i];
+ if i<bov
+ then
+ begin
+ b:=t2*a-c+b1[i];
+ c:=t2*b-a+c1[i]
+ end
+ else
+ spebj1:=(x/8)*(t*a-c+b1[i])
+ end {i}
+ end {abs(x)<=8}
+ else
+ begin
+ g:=t-1.5/tbpi; y:=sqrt(tbpi/t)*spesgn(x);
+ cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
+ t:=128/sqr(t)-1;
+ spebj1:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
+ + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
+ end {abs(x)>8}
+end {spebj1};
+
+function spebk0(x: ArbFloat): ArbFloat;
+
+const
+
+ egam = 0.57721566490153286;
+ xvsmall = 3.2e-9;
+ highexp = 745;
+
+ a0: array[0..7] of ArbFloat =
+ ( 1.12896092945412762e+0, 1.32976966478338191e-1,
+ 4.07157485171389048e-3, 5.59702338227915383e-5,
+ 4.34562671546158210e-7, 2.16382411824721532e-9,
+ 7.49110736894134794e-12, 1.90674197514561280e-14);
+
+ a1: array[0..8] of ArbFloat =
+ ( 2.61841879258687055e-1, 1.52436921799395196e-1,
+ 6.63513979313943827e-3, 1.09534292632401542e-4,
+ 9.57878493265929443e-7, 5.19906865800665633e-9,
+ 1.92405264219706684e-11, 5.16867886946332160e-14,
+ 1.05407718191360000e-16);
+
+ a2: array[0..22] of ArbFloat =
+ ( 9.58210053294896496e-1, -1.42477910128828254e-1,
+ 3.23582010649653009e-2, -8.27780350351692662e-3,
+ 2.24709729617770471e-3, -6.32678357460594866e-4,
+ 1.82652460089342789e-4, -5.37101208898441760e-5,
+ 1.60185974149720562e-5, -4.83134250336922161e-6,
+ 1.47055796078231691e-6, -4.51017292375200017e-7,
+ 1.39217270224614153e-7, -4.32185089841834127e-8,
+ 1.34790467361340101e-8, -4.20597329258249948e-9,
+ 1.32069362385968867e-9, -4.33326665618780914e-10,
+ 1.37999268074442719e-10, -3.19241059198852137e-11,
+ 9.74410152270679245e-12, -7.83738609108569293e-12,
+ 2.57466288575820595e-12);
+
+ a3: array[0..22] of ArbFloat =
+ ( 6.97761598043851776e-1, -1.08801882084935132e-1,
+ 2.56253646031960321e-2, -6.74459607940169198e-3,
+ 1.87292939725962385e-3, -5.37145622971910027e-4,
+ 1.57451516235860573e-4, -4.68936653814896712e-5,
+ 1.41376509343622727e-5, -4.30373871727268511e-6,
+ 1.32052261058932425e-6, -4.07851207862189007e-7,
+ 1.26672629417567360e-7, -3.95403255713518420e-8,
+ 1.23923137898346852e-8, -3.88349705250555658e-9,
+ 1.22424982779432970e-9, -4.03424607871960089e-10,
+ 1.28905587479980147e-10,-2.97787564633235128e-11,
+ 9.11109430833001267e-12,-7.39672783987933184e-12,
+ 2.43538242247537459e-12);
+ a4: array[0..16] of ArbFloat =
+ ( 1.23688664769425422e+0, -1.72683652385321641e-2,
+ -9.25551464765637133e-4, -9.02553345187404564e-5,
+ -6.31692398333746470e-6, -7.69177622529272933e-7,
+ -4.16044811174114579e-8, -9.41555321137176073e-9,
+ 1.75359321273580603e-10, -2.22829582288833265e-10,
+ 3.49564293256545992e-11, -1.11391758572647639e-11,
+ 2.85481235167705907e-12, -7.31344482663931904e-13,
+ 2.06328892562554880e-13, -1.28108310826991616e-13,
+ 4.43741979886551040e-14);
+
+
+var t: ArbFloat;
+
+begin
+ if x<=0
+ then
+ RunError(401);
+ if x<=xvsmall
+ then
+ spebk0:=-(ln(x/2)+egam)
+ else
+ if x<=1
+ then
+ begin
+ t:=2*sqr(x)-1;
+ spebk0:=-ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
+ + spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) - 1)
+ end
+ else
+ if x<=2
+ then
+ spebk0:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
+ else
+ if x<=4
+ then
+ spebk0:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
+ else
+ if x <= highexp
+ then
+ spebk0:=exp(-x)*
+ spepol(10/(1+x)-1, a4[0], sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
+ else
+ spebk0:=0
+end; {spebk0}
+
+function spebk1(x: ArbFloat): ArbFloat;
+
+const
+
+ xsmall = 7.9e-10;
+ highexp = 745;
+ a0: array[0..7] of ArbFloat =
+ ( 5.31907865913352762e-1, 3.25725988137110495e-2,
+ 6.71642805873498653e-4, 6.95300274548206237e-6,
+ 4.32764823642997753e-8, 1.79784792380155752e-10,
+ 5.33888268665658944e-13, 1.18964962439910400e-15);
+
+ a1: array[0..7] of ArbFloat =
+ ( 3.51825828289325536e-1, 4.50490442966943726e-2,
+ 1.20333585658219028e-3, 1.44612432533006139e-5,
+ 9.96686689273781531e-8, 4.46828628435618679e-10,
+ 1.40917103024514301e-12, 3.29881058019865600e-15);
+
+ a2: array[0..23] of ArbFloat =
+ ( 1.24316587355255299e+0, -2.71910714388689413e-1,
+ 8.20250220860693888e-2, -2.62545818729427417e-2,
+ 8.57388087067410089e-3, -2.82450787841655951e-3,
+ 9.34594154387642940e-4, -3.10007681013626626e-4,
+ 1.02982746700060730e-4, -3.42424912211942134e-5,
+ 1.13930169202553526e-5, -3.79227698821142908e-6,
+ 1.26265578331941923e-6, -4.20507152338934956e-7,
+ 1.40138351985185509e-7, -4.66928912168020101e-8,
+ 1.54456653909012693e-8, -5.13783508140332214e-9,
+ 1.82808381381205361e-9, -6.15211416898895086e-10,
+ 1.28044023949946257e-10, -4.02591066627023831e-11,
+ 4.27404330568767242e-11, -1.46639291782948454e-11);
+
+ a3: array[0..23] of ArbFloat =
+ ( 8.06563480128786903e-1, -1.60052611291327173e-1,
+ 4.58591528414023064e-2, -1.42363136684423646e-2,
+ 4.55865751206724687e-3, -1.48185472032688523e-3,
+ 4.85707174778663652e-4, -1.59994873621599146e-4,
+ 5.28712919123131781e-5, -1.75089594354079944e-5,
+ 5.80692311842296724e-6, -1.92794586996432593e-6,
+ 6.40581814037398274e-7, -2.12969229346310343e-7,
+ 7.08723366696569880e-8, -2.35855618461025265e-8,
+ 7.79421651144832709e-9, -2.59039399308009059e-9,
+ 9.20781685906110546e-10, -3.09667392343245062e-10,
+ 6.44913423545894175e-11, -2.02680401514735862e-11,
+ 2.14736751065133220e-11, -7.36478297050421658e-12);
+
+ a4: array[0..16] of ArbFloat =
+ ( 1.30387573604230402e+0, 5.44845254318931612e-2,
+ 4.31639434283445364e-3, 4.29973970898766831e-4,
+ 4.04720631528495020e-5, 4.32776409784235211e-6,
+ 4.07563856931843484e-7, 4.86651420008153956e-8,
+ 3.82717692121438315e-9, 6.77688943857588882e-10,
+ 6.97075379117731379e-12, 1.72026097285930936e-11,
+ -2.60774502020271104e-12, 8.58211523713560576e-13,
+ -2.19287104441802752e-13, 1.39321122940600320e-13,
+ -4.77850238111580160e-14);
+
+var t: ArbFloat;
+
+begin
+ if x<=0
+ then
+ RunError(402);
+ if x<=xsmall
+ then
+ spebk1:=1/x
+ else
+ if x<=1
+ then
+ begin
+ t:=2*sqr(x)-1;
+ spebk1:=( ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
+ -spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) -1) )*x + 1/x
+ end
+ else
+ if x<=2
+ then
+ spebk1:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
+ else
+ if x<=4
+ then
+ spebk1:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
+ else
+ if x <= highexp
+ then
+ spebk1:=exp(-x)*spepol(10/(1+x)-1, a4[0],
+ sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
+ else
+ spebk1:=0
+end; {spebk1}
+
+function speby0(x: ArbFloat): ArbFloat;
+
+const
+
+ tbpi = 6.36619772367581343e-1;
+ egam = 5.77215664901532861e-1;
+ xvsmall = 3.2e-9;
+ a1 : array[0..5] of ArbFloat =
+ ( 3.90000000000000000e-19, -8.74734100000000000e-14,
+ 5.24879478733000000e-9, -5.63207914105698700e-5,
+ 4.71966895957633869e-2, 1.79034314077182663e-1);
+
+ b1 : array[0..5] of ArbFloat =
+ ( -2.69800000000000000e-17, 4.02633082000000000e-12,
+ -1.44072332740190000e-7, 7.53113593257774230e-4,
+ -1.77302012781143582e-1, -2.74474305529745265e-1);
+
+ c1 : array[0..5] of ArbFloat =
+ ( 1.64349000000000000e-15, -1.58375525420000000e-10,
+ 3.20653253765480000e-6, -7.28796247955207918e-3,
+ 2.61567346255046637e-1, -3.31461132032849417e-2);
+
+ d1 : array[0..13] of ArbFloat =
+ ( 9.99457275788251954e-1, -5.36367319213004570e-4,
+ 6.13741608010926000e-6, -2.05274481565160000e-7,
+ 1.28037614434400000e-8, -1.21211819632000000e-9,
+ 1.55005642880000000e-10,-2.48827276800000000e-11,
+ 4.78702080000000000e-12,-1.06365696000000000e-12,
+ 2.45294080000000000e-13,-6.41843200000000000e-14,
+ 3.34028800000000000e-14,-1.17964800000000000e-14);
+
+ d2 : array[0..16] of ArbFloat =
+ (-1.55551138795135187e-2, 6.83314909934390000e-5,
+ -1.47713883264594000e-6, 7.10621485930000000e-8,
+ -5.66871613024000000e-9, 6.43278173280000000e-10,
+ -9.47034774400000000e-11, 1.70330918400000000e-11,
+ -3.59094272000000000e-12, 8.59855360000000000e-13,
+ -2.28807680000000000e-13, 6.95193600000000000e-14,
+ -2.27942400000000000e-14, 4.75136000000000000e-15,
+ -1.14688000000000000e-15, 2.12992000000000000e-15,
+ -9.83040000000000000e-16);
+
+var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
+ i, bov : ArbInt;
+
+begin
+ if x<=0
+ then
+ RunError(403);
+ if x<=xvsmall
+ then
+ speby0:=(ln(x/2)+egam)*tbpi
+ else
+ if x<=8
+ then
+ begin
+ t:=0.03125*sqr(x)-1; t2:=2*t;
+ b:=0; c:=0;
+ bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
+ for i:=0 to bov do
+ begin
+ a:=t2*c-b+a1[i];
+ b:=t2*a-c+b1[i];
+ if i<bov
+ then
+ c:=t2*b-a+c1[i]
+ else
+ speby0:=t*b-a+c1[i]+tbpi*spebj0(x)*ln(x)
+ end {i}
+ end {x<=8}
+ else
+ begin
+ g:=x-1/(2*tbpi); y:=sqrt(tbpi/x);
+ cx:=cos(g)*y*8/x; sx:=sin(g)*y;
+ t:=128/sqr(x)-1;
+ speby0:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
+ + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
+ end {x>8}
+end {speby0};
+
+function speby1(x: ArbFloat): ArbFloat;
+
+const
+ tbpi = 6.36619772367581343e-1;
+ xsmall = 7.9e-10;
+ a1 : array[0..5] of ArbFloat =
+ (-6.58000000000000000e-18, 1.21143321000000000e-12,
+ -5.68844003991900000e-8, 4.40478629867099510e-4,
+ -2.26624991556754924e-1, -1.28697384381350001e-1);
+
+ b1 : array[0..5] of ArbFloat =
+ ( 4.27730000000000000e-16,-5.17212147300000000e-11,
+ 1.41662436449235000e-6, -5.13164116106108479e-3,
+ 6.75615780772187667e-1, 2.03041058859342538e-2);
+
+ c1 : array[0..4] of ArbFloat =
+ (-2.44094900000000000e-14, 1.87547032473000000e-9,
+ -2.83046401495148000e-5, 4.23191803533369041e-2,
+ -7.67296362886645940e-1);
+
+ d1 : array[0..13] of ArbFloat =
+ ( 1.00090702627808217e+0, 8.98804941670557880e-4,
+ -7.95969469843846000e-6, 2.45367662227560000e-7,
+ -1.47085129889600000e-8, 1.36030580128000000e-9,
+ -1.71310758400000000e-10, 2.72040729600000000e-11,
+ -5.19113984000000000e-12, 1.14622464000000000e-12,
+ -2.63372800000000000e-13, 6.86387200000000000e-14,
+ -3.54508800000000000e-14, 1.24928000000000000e-14);
+
+ d2 : array[0..15] of ArbFloat =
+ ( 4.67768740274489776e-2, -9.62145882205441600e-5,
+ 1.82120185123076000e-6, -8.29196070929200000e-8,
+ 6.42013250344000000e-9, -7.15110504800000000e-10,
+ 1.03950931840000000e-10,-1.85248000000000000e-11,
+ 3.87554432000000000e-12,-9.23228160000000000e-13,
+ 2.50224640000000000e-13,-7.43936000000000000e-14,
+ 1.75718400000000000e-14,-4.83328000000000000e-15,
+ 5.32480000000000000e-15,-2.29376000000000000e-15);
+
+var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
+ i, bov : ArbInt;
+
+begin
+ if x<=0
+ then
+ RunError(404);
+ if x<=xsmall
+ then
+ speby1:=-tbpi/x
+ else
+ if x<=8
+ then
+ begin
+ t:=0.03125*sqr(x)-1; t2:=2*t;
+ b:=0; c:=0;
+ bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
+ for i:=0 to bov do
+ begin
+ a:=t2*c-b+a1[i];
+ if i<bov
+ then
+ begin
+ b:=t2*a-c+b1[i];
+ c:=t2*b-a+c1[i]
+ end
+ else
+ if bov=3 {single}
+ then
+ begin
+ b:=t2*a-c+b1[i];
+ speby1:=(t*b-a+c1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
+ end
+ else
+ speby1:=(t*a-c+b1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
+ end {i}
+ end {x<=8}
+ else
+ begin
+ g:=x-3/(2*tbpi); y:=sqrt(tbpi/x);
+ cx:=cos(g)*y*8/x; sx:=sin(g)*y;
+ t:=128/sqr(x)-1;
+ speby1:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
+ + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
+ end {x>8}
+end {speby1};
+
+function speent(x : ArbFloat): longint;
+
+var tx : longint;
+
+begin
+ tx:=trunc(x);
+ if x>=0
+ then
+ speent:=tx
+ else
+ if x=tx
+ then
+ speent:=tx
+ else
+ speent:=tx-1
+end; {speent}
+
+function speerf(x : ArbFloat) : ArbFloat;
+const
+
+ xup = 6.25;
+ sqrtpi = 1.7724538509055160;
+ c : array[1..18] of ArbFloat =
+ ( 1.9449071068178803e0, 4.20186582324414e-2, -1.86866103976769e-2,
+ 5.1281061839107e-3, -1.0683107461726e-3, 1.744737872522e-4,
+ -2.15642065714e-5, 1.7282657974e-6, -2.00479241e-8,
+ -1.64782105e-8, 2.0008475e-9, 2.57716e-11,
+ -3.06343e-11, 1.9158e-12, 3.703e-13,
+ -5.43e-14, -4.0e-15, 1.2e-15);
+
+ d : array[1..17] of ArbFloat =
+ ( 1.4831105640848036e0, -3.010710733865950e-1, 6.89948306898316e-2,
+ -1.39162712647222e-2, 2.4207995224335e-3, -3.658639685849e-4,
+ 4.86209844323e-5, -5.7492565580e-6, 6.113243578e-7,
+ -5.89910153e-8, 5.2070091e-9, -4.232976e-10,
+ 3.18811e-11, -2.2361e-12, 1.467e-13,
+ -9.0e-15, 5.0e-16);
+
+ var t, s, s1, s2, x2: ArbFloat;
+ bovc, bovd, j: ArbInt;
+begin
+ bovc:=sizeof(c) div sizeof(ArbFloat);
+ bovd:=sizeof(d) div sizeof(ArbFloat);
+ t:=abs(x);
+ if t <= 2
+ then
+ begin
+ x2:=sqr(x)-2;
+ s1:=d[bovd]; s2:=0; j:=bovd-1;
+ s:=x2*s1-s2+d[j];
+ while j > 1 do
+ begin
+ s2:=s1; s1:=s; j:=j-1;
+ s:=x2*s1-s2+d[j]
+ end;
+ speerf:=(s-s2)*x/2
+ end
+ else
+ if t < xup
+ then
+ begin
+ x2:=2-20/(t+3);
+ s1:=c[bovc]; s2:=0; j:=bovc-1;
+ s:=x2*s1-s2+c[j];
+ while j > 1 do
+ begin
+ s2:=s1; s1:=s; j:=j-1;
+ s:=x2*s1-s2+c[j]
+ end;
+ x2:=((s-s2)/(2*t))*exp(-sqr(x))/sqrtpi;
+ speerf:=(1-x2)*spesgn(x)
+ end
+ else
+ speerf:=spesgn(x)
+end; {speerf}
+
+function spemax(a, b: ArbFloat): ArbFloat;
+begin
+ if a>b
+ then
+ spemax:=a
+ else
+ spemax:=b
+end; {spemax}
+
+function speefc(x : ArbFloat) : ArbFloat;
+const
+
+ xlow = -6.25;
+ xhigh = 27.28;
+ c : array[0..22] of ArbFloat =
+ ( 1.455897212750385e-1, -2.734219314954260e-1,
+ 2.260080669166197e-1, -1.635718955239687e-1,
+ 1.026043120322792e-1, -5.480232669380236e-2,
+ 2.414322397093253e-2, -8.220621168415435e-3,
+ 1.802962431316418e-3, -2.553523453642242e-5,
+ -1.524627476123466e-4, 4.789838226695987e-5,
+ 3.584014089915968e-6, -6.182369348098529e-6,
+ 7.478317101785790e-7, 6.575825478226343e-7,
+ -1.822565715362025e-7, -7.043998994397452e-8,
+ 3.026547320064576e-8, 7.532536116142436e-9,
+ -4.066088879757269e-9, -5.718639670776992e-10,
+ 3.328130055126039e-10);
+
+ var t, s : ArbFloat;
+begin
+ if x <= xlow
+ then
+ speefc:=2
+ else
+ if x >= xhigh
+ then
+ speefc:=0
+ else
+ begin
+ t:=1-7.5/(abs(x)+3.75);
+ s:=exp(-x*x)*spepol(t, c[0], sizeof(c) div sizeof(ArbFloat) - 1);
+ if x < 0
+ then
+ speefc:=2-s
+ else
+ speefc:=s
+ end
+end {speefc};
+
+function spegam(x: ArbFloat): ArbFloat;
+const
+
+ tmax = 170;
+ a: array[0..23] of ArbFloat =
+ ( 8.86226925452758013e-1, 1.61691987244425092e-2,
+ 1.03703363422075456e-1, -1.34118505705967765e-2,
+ 9.04033494028101968e-3, -2.42259538436268176e-3,
+ 9.15785997288933120e-4, -2.96890121633200000e-4,
+ 1.00928148823365120e-4, -3.36375833240268800e-5,
+ 1.12524642975590400e-5, -3.75499034136576000e-6,
+ 1.25281466396672000e-6, -4.17808776355840000e-7,
+ 1.39383522590720000e-7, -4.64774927155200000e-8,
+ 1.53835215257600000e-8, -5.11961333760000000e-9,
+ 1.82243164160000000e-9, -6.13513953280000000e-10,
+ 1.27679856640000000e-10,-4.01499750400000000e-11,
+ 4.26560716800000000e-11,-1.46381209600000000e-11);
+
+var tvsmall, t, g: ArbFloat;
+ m, i: ArbInt;
+begin
+ if sizeof(ArbFloat) = 6
+ then
+ tvsmall:=2*midget
+ else
+ tvsmall:=midget;
+ t:=abs(x);
+ if t > tmax
+ then
+ RunError(407);
+ if t < macheps
+ then
+ begin
+ if t < tvsmall
+ then
+ RunError(407);
+ spegam:=1/x
+ end
+ else { abs(x) >= macheps }
+ begin
+ m:=trunc(x);
+ if x > 0
+ then
+ begin
+ t:=x-m; m:=m-1; g:=1;
+ if m<0
+ then
+ g:=g/x
+ else
+ if m>0
+ then
+ for i:=1 to m do
+ g:=(x-i)*g
+ end
+ else { x < 0 }
+ begin
+ t:=x-m+1;
+ if t=1
+ then
+ RunError(407);
+ m:=1-m;
+ g:=x;
+ for i:=1 to m do
+ g:=(i+x)*g;
+ g:=1/g
+ end;
+ spegam:=spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1)*g
+ end { abs(x) >= macheps }
+end; {spegam}
+
+function spelga(x: ArbFloat): ArbFloat;
+
+const
+
+ xbig = 7.7e7;
+ xmax = 2.559e305;
+ lnr2pi = 9.18938533204672742e-1;
+ a: array[0..23] of ArbFloat =
+ ( 8.86226925452758013e-1, 1.61691987244425092e-2,
+ 1.03703363422075456e-1, -1.34118505705967765e-2,
+ 9.04033494028101968e-3, -2.42259538436268176e-3,
+ 9.15785997288933120e-4, -2.96890121633200000e-4,
+ 1.00928148823365120e-4, -3.36375833240268800e-5,
+ 1.12524642975590400e-5, -3.75499034136576000e-6,
+ 1.25281466396672000e-6, -4.17808776355840000e-7,
+ 1.39383522590720000e-7, -4.64774927155200000e-8,
+ 1.53835215257600000e-8, -5.11961333760000000e-9,
+ 1.82243164160000000e-9, -6.13513953280000000e-10,
+ 1.27679856640000000e-10,-4.01499750400000000e-11,
+ 4.26560716800000000e-11,-1.46381209600000000e-11);
+ b: array[0..5] of ArbFloat =
+ ( 8.33271644065786580e-2, -6.16502049453716986e-6,
+ 3.89978899876484712e-9, -6.45101975779653651e-12,
+ 2.00201927337982364e-14, -9.94561064728159347e-17);
+
+
+var t, g : ArbFloat;
+ m, i : ArbInt;
+
+begin
+ if x <= 0
+ then
+ RunError(408);
+ if x <= macheps
+ then
+ spelga:=-ln(x)
+ else
+ if x <= 15
+ then
+ begin
+ m:=trunc(x); t:=x-m; m:=m-1; g:=1;
+ if m < 0
+ then
+ g:=g/x
+ else
+ if m > 0
+ then
+ for i:=1 to m do
+ g:=(x-i)*g;
+ spelga:=ln(g*spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1))
+ end
+ else { x > 15 }
+ if x <= xbig
+ then
+ spelga:=(x-0.5)*ln(x) - x + lnr2pi
+ + spepol(450/sqr(x)-1, b[0], sizeof(b) div sizeof(ArbFloat) - 1)/x
+ else { x > xbig }
+ if x <= xmax
+ then
+ spelga:=(x-0.5)*ln(x) - x + lnr2pi
+ else { x > xmax => x*ln(x) > giant }
+ RunError(408)
+end; {spelga}
+
+function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
+var pa : ^arfloat0;
+ i : ArbInt;
+ polx : ArbFloat;
+begin
+ pa:=@a;
+ polx:=0;
+ for i:=n downto 0 do
+ polx:=polx*x+pa^[i];
+ spepol:=polx
+end {spepol};
+
+function spepow(a, b: ArbFloat): ArbFloat;
+
+ function PowInt(a: double; n: longint): double;
+ var a1 : double;
+ begin
+ if n=0 then PowInt := 1 else
+ begin
+ a1 := 1;
+ if n<0 then begin a := 1/a; n := -n end;
+ while n>0
+ do if Odd(n)
+ then begin Dec(n); a1 := a1*a end
+ else begin n := n div 2; a := sqr(a) end;
+ PowInt := a1
+ end
+ end;
+
+var tb : longint;
+ fb : double;
+begin
+
+ { (a < 0, b niet geheel) of (a = 0, b <= 0), dan afbreken}
+ if (a=0) then if (b<=0) then RunError(400) else begin SpePow :=0; exit end;
+ tb := Trunc(b); fb := b-tb;
+ if (a<0) and (fb<>0) then RunError(400);
+
+ if a>0
+ then if fb=0 then SpePow := PowInt(a, tb)
+ else SpePow := PowInt(a, tb)*exp(fb*ln(a))
+ else if odd(tb) then Spepow := -PowInt(-a, tb)
+ else SpePow := PowInt(-a, tb)
+
+end; {spepow}
+
+function spesgn(x: ArbFloat): ArbInt;
+
+begin
+ if x<0
+ then
+ spesgn:=-1
+ else
+ if x=0
+ then
+ spesgn:=0
+ else
+ spesgn:=1
+end; {spesgn}
+
+function spears(x: ArbFloat): ArbFloat;
+const
+
+ pi2 = 1.570796326794897;
+ a : array[0..17] of ArbFloat =
+ ( 1.047197551196598e+0, 5.375149359132719e-2, 7.798902238957732e-3,
+ 1.519668539582420e-3, 3.408637238430600e-4, 8.302317819598986e-5,
+ 2.134554822576075e-5, 5.701781046148566e-6, 1.566985123962741e-6,
+ 4.402076871418002e-7, 1.257811162594110e-7, 3.646577948300129e-8,
+ 1.081021746966715e-8, 3.212744286269388e-9, 8.515014302985799e-10,
+ 2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);
+
+var y, u, t, s : ArbFloat;
+ uprang : boolean;
+begin
+ if abs(x) > 1
+ then
+ RunError(401);
+ u:=sqr(x); uprang:= u > 0.5;
+ if uprang
+ then
+ u:=1-u;
+ t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
+ if uprang
+ then
+ begin
+ s:=pi2-sqrt(u)*y;
+ if x < 0
+ then
+ s:=-s;
+ spears:=s
+ end
+ else
+ spears:=x*y
+end; {spears}
+
+function spearc(x: ArbFloat): ArbFloat;
+const
+
+ pi2 = 1.570796326794897;
+ a : array[0..17] of ArbFloat =
+ ( 1.047197551196598e+0, 5.375149359132719e-2, 7.798902238957732e-3,
+ 1.519668539582420e-3, 3.408637238430600e-4, 8.302317819598986e-5,
+ 2.134554822576075e-5, 5.701781046148566e-6, 1.566985123962741e-6,
+ 4.402076871418002e-7, 1.257811162594110e-7, 3.646577948300129e-8,
+ 1.081021746966715e-8, 3.212744286269388e-9, 8.515014302985799e-10,
+ 2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);
+
+var u, t, y, s : ArbFloat;
+ uprang : boolean;
+begin
+ if abs(x)>1.0
+ then
+ RunError(402);
+ u:=sqr(x); uprang:=u>0.5;
+ if uprang
+ then
+ u:=1-u;
+ t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
+ if uprang
+ then
+ begin
+ s:=sqrt(u)*y;
+ if x<0
+ then
+ s:=2*pi2-s;
+ spearc:=s
+ end
+ else
+ spearc:=pi2-x*y
+end; {spearc}
+
+function spesih(x: ArbFloat): ArbFloat;
+const
+
+ a : array[0..6] of ArbFloat =
+ ( 1.085441641272607e+0, 8.757509762437522e-2, 2.158779361257021e-3,
+ 2.549839945498292e-5, 1.761854853281383e-7, 7.980704288665359e-10,
+ 2.551377137317034e-12);
+
+var u : ArbFloat;
+begin
+ if abs(x)<=1.0
+ then
+ begin
+ u:=2*sqr(x)-1;
+ spesih:=x*spepol(u, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
+ end
+ else
+ begin
+ u:=exp(x); spesih:=(u-1/u)/2
+ end
+end; {spesih}
+
+function specoh(x: ArbFloat): ArbFloat;
+var u: ArbFloat;
+begin
+ u:=exp(x); specoh:=(u+1/u)/2
+end; {specoh}
+
+function spetah(x: ArbFloat): ArbFloat;
+const
+ xhi = 18.50;
+ a : array[0..15] of ArbFloat =
+ ( 8.610571715805476e-1, -1.158834489728470e-1, 1.918072383973950e-2,
+ -3.225255180728459e-3, 5.433071386922689e-4, -9.154289983175165e-5,
+ 1.542469328074432e-5, -2.599022539340038e-6, 4.379282308765732e-7,
+ -7.378980192173815e-8, 1.243517352745986e-8, -2.095373768837420e-9,
+ 3.509758916273561e-10,-5.908745181531817e-11, 1.124199312776748e-11,
+ -1.907888434471600e-12);
+
+var t, y: ArbFloat;
+
+begin
+ t:=abs(x);
+ if t <= 1
+ then
+ begin
+ y:=2*sqr(x)-1;
+ spetah:=x*spepol(y, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
+ end
+ else
+ if t < xhi
+ then
+ begin
+ y:=exp(2*x); spetah:=(y-1)/(y+1)
+ end
+ else
+ spetah:=spesgn(x)
+end; {spetah}
+
+function speash(x: ArbFloat): ArbFloat;
+const
+
+ xhi = 1e9;
+ c : array[0..18] of ArbFloat =
+ ( 9.312298594527122e-1, -5.736663926249348e-2,
+ 9.004288574881897e-3, -1.833458667045431e-3,
+ 4.230023450529706e-4, -1.050715136470630e-4,
+ 2.740790473603819e-5, -7.402952157663977e-6,
+ 2.052474396638805e-6, -5.807433412373489e-7,
+ 1.670117348345774e-7, -4.863477336087045e-8,
+ 1.432753532351304e-8, -4.319978113584910e-9,
+ 1.299779213740398e-9, -3.394726871170490e-10,
+ 1.008344962167889e-10, -5.731943029121004e-11,
+ 1.810792296549804e-11);
+
+
+var t : ArbFloat;
+
+begin
+ t:=abs(x);
+ if t <= 1 then
+ speash:=x*spepol(2*sqr(x)-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
+ else
+ if t < xhi then
+ speash:=ln(sqrt(sqr(x)+1)+t)*spesgn(x)
+ else
+ speash:=ln(2*t)*spesgn(x)
+end; {speash}
+
+function speach(x: ArbFloat): ArbFloat;
+const
+
+ xhi = 1e9;
+
+begin
+ if x<1 then
+ RunError(405);
+ if x=1 then
+ speach:=0
+ else
+ if x<=xhi then
+ speach:=ln(x+sqrt(sqr(x)-1))
+ else
+ speach:=ln(2*x)
+end; {speach}
+
+function speath(x: ArbFloat): ArbFloat;
+const
+
+ c : array[0..19] of ArbFloat =
+ ( 1.098612288668110e+0, 1.173605223326117e-1, 2.309071936165689e-2,
+ 5.449091889986991e-3, 1.404884102286929e-3, 3.816948426588058e-4,
+ 1.073604335435426e-4, 3.095027782918129e-5, 9.088050814470148e-6,
+ 2.706881064641104e-6, 8.155200644023077e-7, 2.479830612463254e-7,
+ 7.588067811453948e-8, 2.339295963220429e-8, 7.408486568719348e-9,
+ 2.319454882064018e-9, 5.960921368486746e-10, 1.820410351379402e-10,
+ 1.184977617320312e-10, 3.856235316559190e-11);
+
+var t, u: ArbFloat;
+begin
+ t:=abs(x);
+ if t >= 1 then
+ RunError(406);
+ u:=sqr(x);
+ if u < 0.5 then
+ speath:=x*spepol(4*u-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
+ else { 0.5 < x*x < 1 }
+ speath:=ln((1+t)/(1-t))/2*spesgn(x)
+end; {speath}
+
+var exitsave : pointer;
+
+procedure MyExit;
+const ErrorS : array[400..408,1..6] of char =
+ ('spepow',
+ 'spebk0',
+ 'spebk1',
+ 'speby0',
+ 'speby1',
+ 'speach',
+ 'speath',
+ 'spegam',
+ 'spelga');
+
+var ErrFil : text;
+
+begin
+ ExitProc := ExitSave;
+ Assign(ErrFil, 'CON');
+ ReWrite(ErrFil);
+ if (ExitCode>=400) AND (ExitCode<=408) then
+ begin
+ write(ErrFil, 'critical error in ', ErrorS[ExitCode]);
+ ExitCode := 201
+ end;
+ Close(ErrFil)
+end;
+
+begin
+ ExitSave := ExitProc;
+ ExitProc := @MyExit;
+end.
View Lorry Controller Status