← Back to branch summary

~ubuntu-branches/ubuntu/lucid/fpc/lucid-proposed

« back to all changes in this revision

Viewing changes to fpcsrc/packages/extra/numlib/spe.pas

  • Committer: Bazaar Package Importer
  • Author(s): Mazen Neifer, Torsten Werner, Mazen Neifer
  • Date: 2008-10-09 23:29:00 UTC
  • mfrom: (4.1.1 sid)
  • Revision ID: james.westby@ubuntu.com-20081009232900-553f61m37jkp6upv
Tags: 2.2.2-4
http://bugs.debian.org/498882
[ Torsten Werner ]
* Update ABI version in fpc-depends automatically.
* Remove empty directories from binary package fpc-source.

[ Mazen Neifer ]
* Removed leading path when calling update-alternatives to remove a Linitian
  error.
* Fixed clean target.
* Improved description of packages. (Closes: #498882)

Show diffs side-by-side

added added

removed removed

Lines of Context:
fpcsrc/packages/extra/numlib/spe.pas
1
 
{
2
 
    This file is part of the Numlib package.
3
 
    Copyright (c) 1986-2000 by
4
 
     Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
5
 
     Computational centre of the Eindhoven University of Technology
6
 
 
7
 
    FPC port Code          by Marco van de Voort (marco@freepascal.org)
8
 
             documentation by Michael van Canneyt (Michael@freepascal.org)
9
 
 
10
 
    Special functions. (Currently, most of these functions only work for
11
 
            ArbFloat=REAL/DOUBLE, not for Extended(at least not at full
12
 
            precision, implement the tables in the diverse functions
13
 
            for extended)
14
 
 
15
 
    See the file COPYING.FPC, included in this distribution,
16
 
    for details about the copyright.
17
 
 
18
 
    This program is distributed in the hope that it will be useful,
19
 
    but WITHOUT ANY WARRANTY; without even the implied warranty of
20
 
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
21
 
 
22
 
 **********************************************************************}
23
 
unit spe;
24
 
{$I DIRECT.INC}
25
 
 
26
 
interface
27
 
 
28
 
uses typ;
29
 
 
30
 
{  Calculate modified Besselfunction "of the first kind" I0(x) }
31
 
function spebi0(x: ArbFloat): ArbFloat;
32
 
 
33
 
{  Calculate modified Besselfunction "of the first kind" I1(x) }
34
 
function spebi1(x: ArbFloat): ArbFloat;
35
 
 
36
 
{  Calculate Besselfunction "of the first kind" J0(x) }
37
 
function spebj0(x: ArbFloat): ArbFloat;
38
 
 
39
 
{  Calculate Besselfunction "of the first kind" J1(x) }
40
 
function spebj1(x: ArbFloat): ArbFloat;
41
 
 
42
 
{  Calculate modified Besselfunction "of the second kind" K0(x) }
43
 
function spebk0(x: ArbFloat): ArbFloat;
44
 
 
45
 
{  Calculate modified Besselfunction "of the second kind" K1(x) }
46
 
function spebk1(x: ArbFloat): ArbFloat;
47
 
 
48
 
{  Calculate Besselfunction "of the second kind" Y0(x) }
49
 
function speby0(x: ArbFloat): ArbFloat;
50
 
 
51
 
{  Calculate Besselfunction "of the second kind" Y1(x) }
52
 
function speby1(x: ArbFloat): ArbFloat;
53
 
 
54
 
{  Entier function, calculates first integer greater or equal than X}
55
 
function speent(x: ArbFloat): longint;
56
 
 
57
 
{  Errorfunction ( 2/sqrt(pi)* Int(t,0,pi,exp(sqr(t)) )}
58
 
function speerf(x: ArbFloat): ArbFloat;
59
 
 
60
 
{  Errorfunction's complement ( 2/sqrt(pi)* Int(t,pi,inf,exp(sqr(t)) )}
61
 
function speefc(x: ArbFloat): ArbFloat;
62
 
 
63
 
{  Function to calculate the Gamma function ( int(t,0,inf,t^(x-1)*exp(-t)) }
64
 
function spegam(x: ArbFloat): ArbFloat;
65
 
 
66
 
{  Function to calculate the natural logaritm of the Gamma function}
67
 
function spelga(x: ArbFloat): ArbFloat;
68
 
 
69
 
{  "Calculates" the maximum of two ArbFloat values     }
70
 
function spemax(a, b: ArbFloat): ArbFloat;
71
 
 
72
 
{  Calculates the functionvalue of a polynomalfunction with n coefficients in a
73
 
for variable X }
74
 
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
75
 
 
76
 
{ Calc a^b with a and b real numbers}
77
 
function spepow(a, b: ArbFloat): ArbFloat;
78
 
 
79
 
{ Returns sign of x (-1 for x<0, 0 for x=0 and 1 for x>0)  }
80
 
function spesgn(x: ArbFloat): ArbInt;
81
 
 
82
 
{  ArcSin(x) }
83
 
function spears(x: ArbFloat): ArbFloat;
84
 
 
85
 
{  ArcCos(x) }
86
 
function spearc(x: ArbFloat): ArbFloat;
87
 
 
88
 
{  Sinh(x) }
89
 
function spesih(x: ArbFloat): ArbFloat;
90
 
 
91
 
{  Cosh(x) }
92
 
function specoh(x: ArbFloat): ArbFloat;
93
 
 
94
 
{  Tanh(x) }
95
 
function spetah(x: ArbFloat): ArbFloat;
96
 
 
97
 
{  ArcSinH(x) }
98
 
function speash(x: ArbFloat): ArbFloat;
99
 
 
100
 
{  ArcCosh(x) }
101
 
function speach(x: ArbFloat): ArbFloat;
102
 
 
103
 
{  ArcTanH(x) }
104
 
function speath(x: ArbFloat): ArbFloat;
105
 
 
106
 
implementation
107
 
 
108
 
function spebi0(x: ArbFloat): ArbFloat;
109
 
 
110
 
const
111
 
 
112
 
     xvsmall = 3.2e-9;
113
 
          a1 : array[0..23] of ArbFloat =
114
 
               (  3.08508322553671039e-1, -1.86478066609466760e-1,
115
 
                  1.57686843969995904e-1, -1.28895621330524993e-1,
116
 
                  9.41616340200868389e-2, -6.04316795007737183e-2,
117
 
                  3.41505388391452157e-2, -1.71317947935716536e-2,
118
 
                  7.70061052263382555e-3, -3.12923286656374358e-3,
119
 
                  1.15888319775791686e-3, -3.93934532072526720e-4,
120
 
                  1.23682594989692688e-4, -3.60645571444886286e-5,
121
 
                  9.81395862769787105e-6, -2.50298975966588680e-6,
122
 
                  6.00566861079330132e-7, -1.36042013507151017e-7,
123
 
                  2.92096163521178835e-8, -5.94856273204259507e-9,
124
 
                  1.13415934215369209e-9, -2.10071360134551962e-10,
125
 
                  4.44484446637868974e-11,-7.48150165756234957e-12);
126
 
          a2 : array[0..26] of ArbFloat =
127
 
               (  1.43431781856850311e-1, -3.71571542566085323e-2,
128
 
                  1.44861237337359455e-2, -6.30121694459896307e-3,
129
 
                  2.89362046530968701e-3, -1.37638906941232170e-3,
130
 
                  6.72508592273773611e-4, -3.35833513200679384e-4,
131
 
                  1.70524543267970595e-4, -8.74354291104467762e-5,
132
 
                  4.48739019580173804e-5, -2.28278155280668483e-5,
133
 
                  1.14032404021741277e-5, -5.54917762110482949e-6,
134
 
                  2.61457634142262604e-6, -1.18752840689765504e-6,
135
 
                  5.18632519069546106e-7, -2.17653548816447667e-7,
136
 
                  8.75291839187305722e-8, -3.34900221934314738e-8,
137
 
                  1.24131668344616429e-8, -4.66215489983794905e-9,
138
 
                  1.58599776268172290e-9, -3.80370174256271589e-10,
139
 
                  1.23188158175419302e-10,-8.46900307934754898e-11,
140
 
                  2.45185252963941089e-11);
141
 
           a3: array[0..19] of ArbFloat =
142
 
               (  4.01071065066847416e-1,  2.18216817211694382e-3,
143
 
                  5.59848253337377763e-5,  2.79770701849785597e-6,
144
 
                  2.17160501061222148e-7,  2.36884434055843528e-8,
145
 
                  3.44345025431425567e-9,  6.47994117793472057e-10,
146
 
                  1.56147127476528831e-10, 4.82726630988879388e-11,
147
 
                  1.89599322920800794e-11, 1.05863621425699789e-11,
148
 
                  8.27719401266046976e-12, 2.82807056475555021e-12,
149
 
                 -4.34624739357691085e-12,-4.29417106720584499e-12,
150
 
                  4.30812577328136192e-13, 1.44572313799118029e-12,
151
 
                  4.73229306831831040e-14,-1.95679809047625728e-13);
152
 
 
153
 
 
154
 
var t : ArbFloat;
155
 
 
156
 
begin
157
 
  t:=abs(x);
158
 
  if t <=xvsmall
159
 
  then
160
 
    spebi0:=1
161
 
  else
162
 
  if t <= 4
163
 
  then
164
 
    spebi0 := exp(t)*spepol(t/2-1, a1[0], SizeOf(a1) div SizeOf(ArbFloat) -1)
165
 
  else
166
 
  if t <= 12
167
 
  then
168
 
    spebi0:=exp(t)*spepol(t/4-2, a2[0], SizeOf(a2) div SizeOf(ArbFloat) -1)
169
 
  else { t > 12}
170
 
    spebi0:=(exp(t)/sqrt(t))*
171
 
            spepol(24/t-1, a3[0], SizeOf(a3) div SizeOf(ArbFloat) -1)
172
 
end; {spebi0}
173
 
 
174
 
function spebi1(x: ArbFloat): ArbFloat;
175
 
 
176
 
 
177
 
const xvsmall = 3.2e-9;
178
 
      a1: array[0..11] of ArbFloat =
179
 
      ( 1.19741654963670236e+0, 9.28758890114609554e-1,
180
 
        2.68657659522092832e-1, 4.09286371827770484e-2,
181
 
        3.84763940423809498e-3, 2.45224314039278904e-4,
182
 
        1.12849795779951847e-5, 3.92368710996392755e-7,
183
 
        1.06662712314503955e-8, 2.32856921884663846e-10,
184
 
        4.17372709788222413e-12,6.24387910353848320e-14);
185
 
 
186
 
      a2: array[0..26] of ArbFloat =
187
 
      ( 1.34142493292698178e-1, -2.99140923897405570e-2,
188
 
        9.76021102528646704e-3, -3.40759647928956354e-3,
189
 
        1.17313412855965374e-3, -3.67626180992174570e-4,
190
 
        8.47999438119288094e-5,  5.21557319070236939e-6,
191
 
       -2.62051678511418163e-5,  2.47493270133518925e-5,
192
 
       -1.79026222757948636e-5,  1.13818992442463952e-5,
193
 
       -6.63144162982509821e-6,  3.60186151617732531e-6,
194
 
       -1.83910206626348772e-6,  8.86951515545183908e-7,
195
 
       -4.05456611578551130e-7,  1.76305222240064495e-7,
196
 
       -7.28978293484163628e-8,  2.84961041291017650e-8,
197
 
       -1.07563514207617768e-8,  4.11321223904934809e-9,
198
 
       -1.41575617446629553e-9,  3.38883570696523350e-10,
199
 
       -1.10970391104678003e-10, 7.79929176497056645e-11,
200
 
       -2.27061376122617856e-11);
201
 
 
202
 
       a3: array[0..19] of ArbFloat =
203
 
       ( 3.92624494204116555e-1, -6.40545360348237412e-3,
204
 
        -9.12475535508497109e-5, -3.82795135453556215e-6,
205
 
        -2.72684545741400871e-7, -2.82537120880041703e-8,
206
 
        -3.96757162863209348e-9, -7.28107961041827952e-10,
207
 
        -1.72060490748583241e-10,-5.23524129533553498e-11,
208
 
        -2.02947854602758139e-11,-1.11795516742222899e-11,
209
 
        -8.69631766630563635e-12,-3.05957293450420224e-12,
210
 
         4.42966462319664333e-12, 4.47735589657057690e-12,
211
 
        -3.95353303949377536e-13,-1.48765082315961139e-12,
212
 
        -5.77176811730370560e-14, 1.99448557598015488e-13);
213
 
 
214
 
var t : ArbFloat;
215
 
 
216
 
begin
217
 
  t:=abs(x);
218
 
  if t <= xvsmall
219
 
  then
220
 
    spebi1:=x/2
221
 
  else
222
 
  if t <= 4
223
 
  then
224
 
    spebi1:=x*spepol(sqr(t)/8-1, a1[0], sizeof(a1) div sizeof(ArbFloat)-1)
225
 
  else
226
 
  if t <= 12
227
 
  then
228
 
    spebi1:=
229
 
      exp(t)*spepol(t/4-2, a2[0], sizeof(a2) div sizeof(ArbFloat)-1)*spesgn(x)
230
 
  else { t > 12}
231
 
    spebi1:=
232
 
      (exp(t)/sqrt(t))*
233
 
      spepol(24/t-1, a3[0], sizeof(a3) div sizeof(ArbFloat)-1)*spesgn(x)
234
 
end; {spebi1}
235
 
 
236
 
function spebj0(x: ArbFloat): ArbFloat;
237
 
const
238
 
 
239
 
       xvsmall = 3.2e-9;
240
 
          tbpi = 6.36619772367581343e-1;
241
 
           a1 : array[0..5] of ArbFloat =
242
 
           ( 1.22200000000000000e-17, -1.94383469000000000e-12,
243
 
             7.60816359241900000e-8,  -4.60626166206275050e-4,
244
 
             1.58067102332097261e-1,  -8.72344235285222129e-3);
245
 
 
246
 
            b1 : array[0..5] of ArbFloat =
247
 
            ( - 7.58850000000000000e-16, 7.84869631400000000e-11,
248
 
              - 1.76194690776215000e-6,  4.81918006946760450e-3,
249
 
              - 3.70094993872649779e-1,  1.57727971474890120e-1);
250
 
 
251
 
            c1 : array[0..4] of ArbFloat =
252
 
            ( 4.12532100000000000e-14, - 2.67925353056000000e-9,
253
 
              3.24603288210050800e-5,  - 3.48937694114088852e-2,
254
 
              2.65178613203336810e-1);
255
 
 
256
 
            d1 : array[0..13] of ArbFloat =
257
 
            ( 9.99457275788251954e-1, -5.36367319213004570e-4,
258
 
              6.13741608010926000e-6, -2.05274481565160000e-7,
259
 
              1.28037614434400000e-8, -1.21211819632000000e-9,
260
 
              1.55005642880000000e-10,-2.48827276800000000e-11,
261
 
              4.78702080000000000e-12,-1.06365696000000000e-12,
262
 
              2.45294080000000000e-13,-6.41843200000000000e-14,
263
 
              3.34028800000000000e-14,-1.17964800000000000e-14);
264
 
 
265
 
             d2 : array[0..16] of ArbFloat =
266
 
             ( -1.55551138795135187e-2,  6.83314909934390000e-5,
267
 
               -1.47713883264594000e-6,  7.10621485930000000e-8,
268
 
               -5.66871613024000000e-9,  6.43278173280000000e-10,
269
 
               -9.47034774400000000e-11, 1.70330918400000000e-11,
270
 
               -3.59094272000000000e-12, 8.59855360000000000e-13,
271
 
               -2.28807680000000000e-13, 6.95193600000000000e-14,
272
 
               -2.27942400000000000e-14, 4.75136000000000000e-15,
273
 
               -1.14688000000000000e-15, 2.12992000000000000e-15,
274
 
               -9.83040000000000000e-16);
275
 
 
276
 
var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
277
 
    i, bov                       : ArbInt;
278
 
 
279
 
begin
280
 
  t:=abs(x);
281
 
  if t<=xvsmall
282
 
  then
283
 
    spebj0:=1
284
 
  else
285
 
  if t<=8
286
 
  then
287
 
    begin
288
 
      t:=0.03125*sqr(t)-1; t2:=2*t;
289
 
      b:=0; c:=0;
290
 
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
291
 
      for i:=0 to bov do
292
 
        begin
293
 
          a:=t2*c-b+a1[i];
294
 
          if i<5
295
 
          then
296
 
            b:=t2*a-c+b1[i]
297
 
          else
298
 
            spebj0:=t*a-c+b1[i];
299
 
          if i<bov
300
 
          then
301
 
            c:=t2*b-a+c1[i]
302
 
          else
303
 
            if i<5
304
 
            then
305
 
              spebj0:=t*b-a+c1[i]
306
 
        end {i}
307
 
    end {abs(x)<=8}
308
 
  else
309
 
    begin
310
 
      g:=t-1/(2*tbpi); y:=sqrt(tbpi/t);
311
 
      cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
312
 
      t:=128/sqr(t)-1;
313
 
      spebj0:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
314
 
              + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
315
 
    end {abs(x)>8}
316
 
 
317
 
end {spebj0};
318
 
 
319
 
function spebj1(x: ArbFloat): ArbFloat;
320
 
const
321
 
 
322
 
       xvsmall = 3.2e-9;
323
 
          tbpi = 6.36619772367581343e-1;
324
 
      a1 : array[0..5] of ArbFloat =
325
 
      ( 2.95000000000000000e-18, -5.77740420000000000e-13,
326
 
        2.94970700727800000e-8,  -2.60444389348580680e-4,
327
 
        1.77709117239728283e-1,  -1.19180116054121687e+0);
328
 
 
329
 
      b1 : array[0..5] of ArbFloat =
330
 
      ( -1.95540000000000000e-16, 2.52812366400000000e-11,
331
 
        -7.61758780540030000e-7,  3.24027018268385747e-3,
332
 
        -6.61443934134543253e-1,  6.48358770605264921e-1);
333
 
 
334
 
      c1 : array[0..4] of ArbFloat =
335
 
      ( 1.13857200000000000e-14, -9.42421298160000000e-10,
336
 
        1.58870192399321300e-5,  -2.91755248061542077e-2,
337
 
        1.28799409885767762e+0);
338
 
 
339
 
       d1 : array[0..13] of ArbFloat =
340
 
       ( 1.00090702627808217e+0,  8.98804941670557880e-4,
341
 
        -7.95969469843846000e-6,  2.45367662227560000e-7,
342
 
        -1.47085129889600000e-8,  1.36030580128000000e-9,
343
 
        -1.71310758400000000e-10, 2.72040729600000000e-11,
344
 
        -5.19113984000000000e-12, 1.14622464000000000e-12,
345
 
        -2.63372800000000000e-13, 6.86387200000000000e-14,
346
 
        -3.54508800000000000e-14, 1.24928000000000000e-14);
347
 
 
348
 
       d2 : array[0..15] of ArbFloat =
349
 
       ( 4.67768740274489776e-2,  -9.62145882205441600e-5,
350
 
         1.82120185123076000e-6,  -8.29196070929200000e-8,
351
 
         6.42013250344000000e-9,  -7.15110504800000000e-10,
352
 
         1.03950931840000000e-10, -1.85248000000000000e-11,
353
 
         3.87554432000000000e-12, -9.23228160000000000e-13,
354
 
         2.50224640000000000e-13, -7.43936000000000000e-14,
355
 
         1.75718400000000000e-14, -4.83328000000000000e-15,
356
 
         5.32480000000000000e-15, -2.29376000000000000e-15);
357
 
 
358
 
var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
359
 
    i, bov                       : ArbInt;
360
 
 
361
 
begin
362
 
  t:=abs(x);
363
 
  if t<xvsmall
364
 
  then
365
 
    spebj1:=x/2
366
 
  else
367
 
  if t<=8
368
 
  then
369
 
    begin
370
 
      t:=0.03125*sqr(t)-1; t2:=2*t;
371
 
      b:=0; c:=0;
372
 
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
373
 
      for i:=0 to bov do
374
 
        begin
375
 
          a:=t2*c-b+a1[i];
376
 
          if i<bov
377
 
          then
378
 
            begin
379
 
              b:=t2*a-c+b1[i];
380
 
              c:=t2*b-a+c1[i]
381
 
            end
382
 
          else
383
 
            spebj1:=(x/8)*(t*a-c+b1[i])
384
 
        end {i}
385
 
    end {abs(x)<=8}
386
 
  else
387
 
    begin
388
 
      g:=t-1.5/tbpi; y:=sqrt(tbpi/t)*spesgn(x);
389
 
      cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
390
 
      t:=128/sqr(t)-1;
391
 
      spebj1:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
392
 
              + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
393
 
    end {abs(x)>8}
394
 
end {spebj1};
395
 
 
396
 
function spebk0(x: ArbFloat): ArbFloat;
397
 
 
398
 
const
399
 
 
400
 
     egam = 0.57721566490153286;
401
 
     xvsmall = 3.2e-9;
402
 
     highexp = 745;
403
 
 
404
 
      a0: array[0..7] of ArbFloat =
405
 
      ( 1.12896092945412762e+0,  1.32976966478338191e-1,
406
 
        4.07157485171389048e-3,  5.59702338227915383e-5,
407
 
        4.34562671546158210e-7,  2.16382411824721532e-9,
408
 
        7.49110736894134794e-12, 1.90674197514561280e-14);
409
 
 
410
 
      a1: array[0..8] of ArbFloat =
411
 
      ( 2.61841879258687055e-1,  1.52436921799395196e-1,
412
 
        6.63513979313943827e-3,  1.09534292632401542e-4,
413
 
        9.57878493265929443e-7,  5.19906865800665633e-9,
414
 
        1.92405264219706684e-11, 5.16867886946332160e-14,
415
 
        1.05407718191360000e-16);
416
 
 
417
 
      a2: array[0..22] of ArbFloat =
418
 
      ( 9.58210053294896496e-1, -1.42477910128828254e-1,
419
 
        3.23582010649653009e-2, -8.27780350351692662e-3,
420
 
        2.24709729617770471e-3, -6.32678357460594866e-4,
421
 
        1.82652460089342789e-4, -5.37101208898441760e-5,
422
 
        1.60185974149720562e-5, -4.83134250336922161e-6,
423
 
        1.47055796078231691e-6, -4.51017292375200017e-7,
424
 
        1.39217270224614153e-7, -4.32185089841834127e-8,
425
 
        1.34790467361340101e-8, -4.20597329258249948e-9,
426
 
        1.32069362385968867e-9, -4.33326665618780914e-10,
427
 
        1.37999268074442719e-10, -3.19241059198852137e-11,
428
 
        9.74410152270679245e-12, -7.83738609108569293e-12,
429
 
        2.57466288575820595e-12);
430
 
 
431
 
      a3: array[0..22] of ArbFloat =
432
 
     ( 6.97761598043851776e-1, -1.08801882084935132e-1,
433
 
       2.56253646031960321e-2, -6.74459607940169198e-3,
434
 
       1.87292939725962385e-3, -5.37145622971910027e-4,
435
 
       1.57451516235860573e-4, -4.68936653814896712e-5,
436
 
       1.41376509343622727e-5, -4.30373871727268511e-6,
437
 
       1.32052261058932425e-6, -4.07851207862189007e-7,
438
 
       1.26672629417567360e-7, -3.95403255713518420e-8,
439
 
       1.23923137898346852e-8, -3.88349705250555658e-9,
440
 
       1.22424982779432970e-9, -4.03424607871960089e-10,
441
 
       1.28905587479980147e-10,-2.97787564633235128e-11,
442
 
       9.11109430833001267e-12,-7.39672783987933184e-12,
443
 
       2.43538242247537459e-12);
444
 
      a4: array[0..16] of ArbFloat =
445
 
      ( 1.23688664769425422e+0,  -1.72683652385321641e-2,
446
 
       -9.25551464765637133e-4,  -9.02553345187404564e-5,
447
 
       -6.31692398333746470e-6,  -7.69177622529272933e-7,
448
 
       -4.16044811174114579e-8,  -9.41555321137176073e-9,
449
 
        1.75359321273580603e-10, -2.22829582288833265e-10,
450
 
        3.49564293256545992e-11, -1.11391758572647639e-11,
451
 
        2.85481235167705907e-12, -7.31344482663931904e-13,
452
 
        2.06328892562554880e-13, -1.28108310826991616e-13,
453
 
        4.43741979886551040e-14);
454
 
 
455
 
 
456
 
var t: ArbFloat;
457
 
 
458
 
begin
459
 
  if x<=0
460
 
  then
461
 
    RunError(401);
462
 
  if x<=xvsmall
463
 
  then
464
 
    spebk0:=-(ln(x/2)+egam)
465
 
  else
466
 
  if x<=1
467
 
  then
468
 
    begin
469
 
      t:=2*sqr(x)-1;
470
 
      spebk0:=-ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
471
 
              + spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) - 1)
472
 
    end
473
 
  else
474
 
  if x<=2
475
 
  then
476
 
    spebk0:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
477
 
  else
478
 
  if x<=4
479
 
  then
480
 
    spebk0:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
481
 
  else
482
 
  if x <= highexp
483
 
  then
484
 
    spebk0:=exp(-x)*
485
 
            spepol(10/(1+x)-1, a4[0], sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
486
 
  else
487
 
    spebk0:=0
488
 
end; {spebk0}
489
 
 
490
 
function spebk1(x: ArbFloat): ArbFloat;
491
 
 
492
 
const
493
 
 
494
 
   xsmall = 7.9e-10;
495
 
  highexp = 745;
496
 
   a0: array[0..7] of ArbFloat =
497
 
   ( 5.31907865913352762e-1,  3.25725988137110495e-2,
498
 
     6.71642805873498653e-4,  6.95300274548206237e-6,
499
 
     4.32764823642997753e-8,  1.79784792380155752e-10,
500
 
     5.33888268665658944e-13, 1.18964962439910400e-15);
501
 
 
502
 
   a1: array[0..7] of ArbFloat =
503
 
   ( 3.51825828289325536e-1,  4.50490442966943726e-2,
504
 
     1.20333585658219028e-3,  1.44612432533006139e-5,
505
 
     9.96686689273781531e-8,  4.46828628435618679e-10,
506
 
     1.40917103024514301e-12, 3.29881058019865600e-15);
507
 
 
508
 
   a2: array[0..23] of ArbFloat =
509
 
   ( 1.24316587355255299e+0, -2.71910714388689413e-1,
510
 
     8.20250220860693888e-2, -2.62545818729427417e-2,
511
 
     8.57388087067410089e-3, -2.82450787841655951e-3,
512
 
     9.34594154387642940e-4, -3.10007681013626626e-4,
513
 
     1.02982746700060730e-4, -3.42424912211942134e-5,
514
 
     1.13930169202553526e-5, -3.79227698821142908e-6,
515
 
     1.26265578331941923e-6, -4.20507152338934956e-7,
516
 
     1.40138351985185509e-7, -4.66928912168020101e-8,
517
 
     1.54456653909012693e-8, -5.13783508140332214e-9,
518
 
     1.82808381381205361e-9, -6.15211416898895086e-10,
519
 
     1.28044023949946257e-10, -4.02591066627023831e-11,
520
 
     4.27404330568767242e-11, -1.46639291782948454e-11);
521
 
 
522
 
   a3: array[0..23] of ArbFloat =
523
 
   ( 8.06563480128786903e-1,  -1.60052611291327173e-1,
524
 
     4.58591528414023064e-2,  -1.42363136684423646e-2,
525
 
     4.55865751206724687e-3,  -1.48185472032688523e-3,
526
 
     4.85707174778663652e-4,  -1.59994873621599146e-4,
527
 
     5.28712919123131781e-5,  -1.75089594354079944e-5,
528
 
     5.80692311842296724e-6,  -1.92794586996432593e-6,
529
 
     6.40581814037398274e-7,  -2.12969229346310343e-7,
530
 
     7.08723366696569880e-8,  -2.35855618461025265e-8,
531
 
     7.79421651144832709e-9,  -2.59039399308009059e-9,
532
 
     9.20781685906110546e-10, -3.09667392343245062e-10,
533
 
     6.44913423545894175e-11, -2.02680401514735862e-11,
534
 
     2.14736751065133220e-11, -7.36478297050421658e-12);
535
 
 
536
 
    a4: array[0..16] of ArbFloat =
537
 
    ( 1.30387573604230402e+0,   5.44845254318931612e-2,
538
 
      4.31639434283445364e-3,   4.29973970898766831e-4,
539
 
      4.04720631528495020e-5,   4.32776409784235211e-6,
540
 
      4.07563856931843484e-7,   4.86651420008153956e-8,
541
 
      3.82717692121438315e-9,   6.77688943857588882e-10,
542
 
      6.97075379117731379e-12,  1.72026097285930936e-11,
543
 
     -2.60774502020271104e-12,  8.58211523713560576e-13,
544
 
     -2.19287104441802752e-13,  1.39321122940600320e-13,
545
 
     -4.77850238111580160e-14);
546
 
 
547
 
var t: ArbFloat;
548
 
 
549
 
begin
550
 
  if x<=0
551
 
  then
552
 
    RunError(402);
553
 
  if x<=xsmall
554
 
  then
555
 
    spebk1:=1/x
556
 
  else
557
 
  if x<=1
558
 
  then
559
 
    begin
560
 
      t:=2*sqr(x)-1;
561
 
      spebk1:=( ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
562
 
              -spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) -1) )*x + 1/x
563
 
    end
564
 
  else
565
 
  if x<=2
566
 
  then
567
 
    spebk1:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
568
 
  else
569
 
  if x<=4
570
 
  then
571
 
    spebk1:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
572
 
  else
573
 
  if x <= highexp
574
 
  then
575
 
    spebk1:=exp(-x)*spepol(10/(1+x)-1, a4[0],
576
 
            sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
577
 
  else
578
 
    spebk1:=0
579
 
end; {spebk1}
580
 
 
581
 
function speby0(x: ArbFloat): ArbFloat;
582
 
 
583
 
const
584
 
 
585
 
      tbpi = 6.36619772367581343e-1;
586
 
      egam = 5.77215664901532861e-1;
587
 
   xvsmall = 3.2e-9;
588
 
   a1 : array[0..5] of ArbFloat =
589
 
   ( 3.90000000000000000e-19, -8.74734100000000000e-14,
590
 
     5.24879478733000000e-9,  -5.63207914105698700e-5,
591
 
     4.71966895957633869e-2,   1.79034314077182663e-1);
592
 
 
593
 
   b1 : array[0..5] of ArbFloat =
594
 
   ( -2.69800000000000000e-17, 4.02633082000000000e-12,
595
 
     -1.44072332740190000e-7,  7.53113593257774230e-4,
596
 
     -1.77302012781143582e-1, -2.74474305529745265e-1);
597
 
 
598
 
   c1 : array[0..5] of ArbFloat =
599
 
   ( 1.64349000000000000e-15, -1.58375525420000000e-10,
600
 
     3.20653253765480000e-6,  -7.28796247955207918e-3,
601
 
     2.61567346255046637e-1,  -3.31461132032849417e-2);
602
 
 
603
 
    d1 : array[0..13] of ArbFloat =
604
 
    ( 9.99457275788251954e-1, -5.36367319213004570e-4,
605
 
      6.13741608010926000e-6, -2.05274481565160000e-7,
606
 
      1.28037614434400000e-8, -1.21211819632000000e-9,
607
 
      1.55005642880000000e-10,-2.48827276800000000e-11,
608
 
      4.78702080000000000e-12,-1.06365696000000000e-12,
609
 
      2.45294080000000000e-13,-6.41843200000000000e-14,
610
 
      3.34028800000000000e-14,-1.17964800000000000e-14);
611
 
 
612
 
    d2 : array[0..16] of ArbFloat =
613
 
    (-1.55551138795135187e-2,  6.83314909934390000e-5,
614
 
     -1.47713883264594000e-6,  7.10621485930000000e-8,
615
 
     -5.66871613024000000e-9,  6.43278173280000000e-10,
616
 
     -9.47034774400000000e-11, 1.70330918400000000e-11,
617
 
     -3.59094272000000000e-12, 8.59855360000000000e-13,
618
 
     -2.28807680000000000e-13, 6.95193600000000000e-14,
619
 
     -2.27942400000000000e-14, 4.75136000000000000e-15,
620
 
     -1.14688000000000000e-15, 2.12992000000000000e-15,
621
 
     -9.83040000000000000e-16);
622
 
 
623
 
var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
624
 
    i, bov                       : ArbInt;
625
 
 
626
 
begin
627
 
  if x<=0
628
 
  then
629
 
    RunError(403);
630
 
  if x<=xvsmall
631
 
  then
632
 
    speby0:=(ln(x/2)+egam)*tbpi
633
 
  else
634
 
  if x<=8
635
 
  then
636
 
    begin
637
 
      t:=0.03125*sqr(x)-1; t2:=2*t;
638
 
      b:=0; c:=0;
639
 
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
640
 
      for i:=0 to bov do
641
 
        begin
642
 
          a:=t2*c-b+a1[i];
643
 
          b:=t2*a-c+b1[i];
644
 
          if i<bov
645
 
          then
646
 
            c:=t2*b-a+c1[i]
647
 
          else
648
 
            speby0:=t*b-a+c1[i]+tbpi*spebj0(x)*ln(x)
649
 
        end {i}
650
 
    end {x<=8}
651
 
  else
652
 
    begin
653
 
      g:=x-1/(2*tbpi); y:=sqrt(tbpi/x);
654
 
      cx:=cos(g)*y*8/x; sx:=sin(g)*y;
655
 
      t:=128/sqr(x)-1;
656
 
      speby0:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
657
 
            + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
658
 
    end {x>8}
659
 
end {speby0};
660
 
 
661
 
function speby1(x: ArbFloat): ArbFloat;
662
 
 
663
 
const
664
 
    tbpi = 6.36619772367581343e-1;
665
 
    xsmall = 7.9e-10;
666
 
   a1 : array[0..5] of ArbFloat =
667
 
   (-6.58000000000000000e-18, 1.21143321000000000e-12,
668
 
    -5.68844003991900000e-8,  4.40478629867099510e-4,
669
 
    -2.26624991556754924e-1, -1.28697384381350001e-1);
670
 
 
671
 
   b1 : array[0..5] of ArbFloat =
672
 
   ( 4.27730000000000000e-16,-5.17212147300000000e-11,
673
 
     1.41662436449235000e-6, -5.13164116106108479e-3,
674
 
     6.75615780772187667e-1,  2.03041058859342538e-2);
675
 
 
676
 
   c1 : array[0..4] of ArbFloat =
677
 
   (-2.44094900000000000e-14, 1.87547032473000000e-9,
678
 
    -2.83046401495148000e-5,  4.23191803533369041e-2,
679
 
    -7.67296362886645940e-1);
680
 
 
681
 
   d1 : array[0..13] of ArbFloat =
682
 
   ( 1.00090702627808217e+0,  8.98804941670557880e-4,
683
 
    -7.95969469843846000e-6,  2.45367662227560000e-7,
684
 
    -1.47085129889600000e-8,  1.36030580128000000e-9,
685
 
    -1.71310758400000000e-10, 2.72040729600000000e-11,
686
 
    -5.19113984000000000e-12, 1.14622464000000000e-12,
687
 
    -2.63372800000000000e-13, 6.86387200000000000e-14,
688
 
    -3.54508800000000000e-14, 1.24928000000000000e-14);
689
 
 
690
 
    d2 : array[0..15] of ArbFloat =
691
 
    ( 4.67768740274489776e-2, -9.62145882205441600e-5,
692
 
      1.82120185123076000e-6, -8.29196070929200000e-8,
693
 
      6.42013250344000000e-9, -7.15110504800000000e-10,
694
 
      1.03950931840000000e-10,-1.85248000000000000e-11,
695
 
      3.87554432000000000e-12,-9.23228160000000000e-13,
696
 
      2.50224640000000000e-13,-7.43936000000000000e-14,
697
 
      1.75718400000000000e-14,-4.83328000000000000e-15,
698
 
      5.32480000000000000e-15,-2.29376000000000000e-15);
699
 
 
700
 
var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
701
 
    i, bov                       : ArbInt;
702
 
 
703
 
begin
704
 
  if x<=0
705
 
  then
706
 
    RunError(404);
707
 
  if x<=xsmall
708
 
  then
709
 
    speby1:=-tbpi/x
710
 
  else
711
 
  if x<=8
712
 
  then
713
 
    begin
714
 
      t:=0.03125*sqr(x)-1; t2:=2*t;
715
 
      b:=0; c:=0;
716
 
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
717
 
      for i:=0 to bov do
718
 
        begin
719
 
          a:=t2*c-b+a1[i];
720
 
          if i<bov
721
 
          then
722
 
            begin
723
 
              b:=t2*a-c+b1[i];
724
 
              c:=t2*b-a+c1[i]
725
 
            end
726
 
          else
727
 
          if bov=3   {single}
728
 
          then
729
 
            begin
730
 
              b:=t2*a-c+b1[i];
731
 
              speby1:=(t*b-a+c1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
732
 
            end
733
 
          else
734
 
            speby1:=(t*a-c+b1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
735
 
        end {i}
736
 
    end {x<=8}
737
 
  else
738
 
    begin
739
 
      g:=x-3/(2*tbpi); y:=sqrt(tbpi/x);
740
 
      cx:=cos(g)*y*8/x; sx:=sin(g)*y;
741
 
      t:=128/sqr(x)-1;
742
 
      speby1:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
743
 
            + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
744
 
    end {x>8}
745
 
end {speby1};
746
 
 
747
 
function speent(x : ArbFloat): longint;
748
 
 
749
 
var tx : longint;
750
 
 
751
 
begin
752
 
  tx:=trunc(x);
753
 
  if x>=0
754
 
  then
755
 
    speent:=tx
756
 
  else
757
 
    if x=tx
758
 
    then
759
 
      speent:=tx
760
 
    else
761
 
      speent:=tx-1
762
 
end; {speent}
763
 
 
764
 
function speerf(x : ArbFloat) : ArbFloat;
765
 
const
766
 
 
767
 
        xup = 6.25;
768
 
     sqrtpi = 1.7724538509055160;
769
 
     c : array[1..18] of ArbFloat =
770
 
     ( 1.9449071068178803e0,  4.20186582324414e-2, -1.86866103976769e-2,
771
 
       5.1281061839107e-3,   -1.0683107461726e-3,   1.744737872522e-4,
772
 
      -2.15642065714e-5,      1.7282657974e-6,     -2.00479241e-8,
773
 
      -1.64782105e-8,         2.0008475e-9,         2.57716e-11,
774
 
      -3.06343e-11,           1.9158e-12,           3.703e-13,
775
 
      -5.43e-14,             -4.0e-15,              1.2e-15);
776
 
 
777
 
     d : array[1..17] of ArbFloat =
778
 
     ( 1.4831105640848036e0, -3.010710733865950e-1, 6.89948306898316e-2,
779
 
      -1.39162712647222e-2,   2.4207995224335e-3,  -3.658639685849e-4,
780
 
       4.86209844323e-5,     -5.7492565580e-6,      6.113243578e-7,
781
 
      -5.89910153e-8,         5.2070091e-9,        -4.232976e-10,
782
 
       3.18811e-11,          -2.2361e-12,           1.467e-13,
783
 
      -9.0e-15,               5.0e-16);
784
 
 
785
 
  var t, s, s1, s2, x2: ArbFloat;
786
 
         bovc, bovd, j: ArbInt;
787
 
begin
788
 
  bovc:=sizeof(c) div sizeof(ArbFloat);
789
 
  bovd:=sizeof(d) div sizeof(ArbFloat);
790
 
  t:=abs(x);
791
 
  if t <= 2
792
 
  then
793
 
    begin
794
 
      x2:=sqr(x)-2;
795
 
      s1:=d[bovd]; s2:=0; j:=bovd-1;
796
 
      s:=x2*s1-s2+d[j];
797
 
      while j > 1 do
798
 
        begin
799
 
          s2:=s1; s1:=s; j:=j-1;
800
 
          s:=x2*s1-s2+d[j]
801
 
        end;
802
 
      speerf:=(s-s2)*x/2
803
 
    end
804
 
  else
805
 
    if t < xup
806
 
    then
807
 
      begin
808
 
        x2:=2-20/(t+3);
809
 
        s1:=c[bovc]; s2:=0; j:=bovc-1;
810
 
        s:=x2*s1-s2+c[j];
811
 
        while j > 1 do
812
 
          begin
813
 
            s2:=s1; s1:=s; j:=j-1;
814
 
            s:=x2*s1-s2+c[j]
815
 
          end;
816
 
        x2:=((s-s2)/(2*t))*exp(-sqr(x))/sqrtpi;
817
 
        speerf:=(1-x2)*spesgn(x)
818
 
      end
819
 
    else
820
 
      speerf:=spesgn(x)
821
 
end;  {speerf}
822
 
 
823
 
function spemax(a, b: ArbFloat): ArbFloat;
824
 
begin
825
 
  if a>b
826
 
  then
827
 
    spemax:=a
828
 
  else
829
 
    spemax:=b
830
 
end; {spemax}
831
 
 
832
 
function speefc(x : ArbFloat) : ArbFloat;
833
 
const
834
 
 
835
 
   xlow = -6.25;
836
 
  xhigh = 27.28;
837
 
      c : array[0..22] of ArbFloat =
838
 
      ( 1.455897212750385e-1, -2.734219314954260e-1,
839
 
        2.260080669166197e-1, -1.635718955239687e-1,
840
 
        1.026043120322792e-1, -5.480232669380236e-2,
841
 
        2.414322397093253e-2, -8.220621168415435e-3,
842
 
        1.802962431316418e-3, -2.553523453642242e-5,
843
 
       -1.524627476123466e-4,  4.789838226695987e-5,
844
 
        3.584014089915968e-6, -6.182369348098529e-6,
845
 
        7.478317101785790e-7,  6.575825478226343e-7,
846
 
       -1.822565715362025e-7, -7.043998994397452e-8,
847
 
        3.026547320064576e-8,  7.532536116142436e-9,
848
 
       -4.066088879757269e-9, -5.718639670776992e-10,
849
 
        3.328130055126039e-10);
850
 
 
851
 
  var t, s : ArbFloat;
852
 
begin
853
 
  if x <= xlow
854
 
  then
855
 
    speefc:=2
856
 
  else
857
 
  if x >= xhigh
858
 
  then
859
 
    speefc:=0
860
 
  else
861
 
    begin
862
 
      t:=1-7.5/(abs(x)+3.75);
863
 
      s:=exp(-x*x)*spepol(t, c[0], sizeof(c) div sizeof(ArbFloat) - 1);
864
 
      if x < 0
865
 
      then
866
 
        speefc:=2-s
867
 
      else
868
 
        speefc:=s
869
 
    end
870
 
end {speefc};
871
 
 
872
 
function spegam(x: ArbFloat): ArbFloat;
873
 
const
874
 
 
875
 
    tmax = 170;
876
 
    a: array[0..23] of ArbFloat =
877
 
    ( 8.86226925452758013e-1,  1.61691987244425092e-2,
878
 
      1.03703363422075456e-1, -1.34118505705967765e-2,
879
 
      9.04033494028101968e-3, -2.42259538436268176e-3,
880
 
      9.15785997288933120e-4, -2.96890121633200000e-4,
881
 
      1.00928148823365120e-4, -3.36375833240268800e-5,
882
 
      1.12524642975590400e-5, -3.75499034136576000e-6,
883
 
      1.25281466396672000e-6, -4.17808776355840000e-7,
884
 
      1.39383522590720000e-7, -4.64774927155200000e-8,
885
 
      1.53835215257600000e-8, -5.11961333760000000e-9,
886
 
      1.82243164160000000e-9, -6.13513953280000000e-10,
887
 
      1.27679856640000000e-10,-4.01499750400000000e-11,
888
 
      4.26560716800000000e-11,-1.46381209600000000e-11);
889
 
 
890
 
var tvsmall, t, g: ArbFloat;
891
 
             m, i: ArbInt;
892
 
begin
893
 
  if sizeof(ArbFloat) = 6
894
 
  then
895
 
    tvsmall:=2*midget
896
 
  else
897
 
    tvsmall:=midget;
898
 
  t:=abs(x);
899
 
  if t > tmax
900
 
  then
901
 
    RunError(407);
902
 
  if t < macheps
903
 
  then
904
 
    begin
905
 
      if t < tvsmall
906
 
      then
907
 
        RunError(407);
908
 
      spegam:=1/x
909
 
    end
910
 
  else  { abs(x) >= macheps }
911
 
    begin
912
 
      m:=trunc(x);
913
 
      if x > 0
914
 
      then
915
 
        begin
916
 
          t:=x-m; m:=m-1; g:=1;
917
 
          if m<0
918
 
          then
919
 
            g:=g/x
920
 
          else
921
 
            if m>0
922
 
            then
923
 
              for i:=1 to m do
924
 
                g:=(x-i)*g
925
 
        end
926
 
      else { x < 0 }
927
 
        begin
928
 
          t:=x-m+1;
929
 
          if t=1
930
 
          then
931
 
            RunError(407);
932
 
          m:=1-m;
933
 
          g:=x;
934
 
          for i:=1 to m do
935
 
            g:=(i+x)*g;
936
 
          g:=1/g
937
 
        end;
938
 
      spegam:=spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1)*g
939
 
    end { abs(x) >= macheps }
940
 
end; {spegam}
941
 
 
942
 
function spelga(x: ArbFloat): ArbFloat;
943
 
 
944
 
const
945
 
 
946
 
    xbig = 7.7e7;
947
 
    xmax = 2.559e305;
948
 
  lnr2pi = 9.18938533204672742e-1;
949
 
    a: array[0..23] of ArbFloat =
950
 
    ( 8.86226925452758013e-1,  1.61691987244425092e-2,
951
 
      1.03703363422075456e-1, -1.34118505705967765e-2,
952
 
      9.04033494028101968e-3, -2.42259538436268176e-3,
953
 
      9.15785997288933120e-4, -2.96890121633200000e-4,
954
 
      1.00928148823365120e-4, -3.36375833240268800e-5,
955
 
      1.12524642975590400e-5, -3.75499034136576000e-6,
956
 
      1.25281466396672000e-6, -4.17808776355840000e-7,
957
 
      1.39383522590720000e-7, -4.64774927155200000e-8,
958
 
      1.53835215257600000e-8, -5.11961333760000000e-9,
959
 
      1.82243164160000000e-9, -6.13513953280000000e-10,
960
 
      1.27679856640000000e-10,-4.01499750400000000e-11,
961
 
      4.26560716800000000e-11,-1.46381209600000000e-11);
962
 
    b: array[0..5] of ArbFloat =
963
 
    ( 8.33271644065786580e-2,  -6.16502049453716986e-6,
964
 
      3.89978899876484712e-9,  -6.45101975779653651e-12,
965
 
      2.00201927337982364e-14, -9.94561064728159347e-17);
966
 
 
967
 
 
968
 
var t, g : ArbFloat;
969
 
    m, i : ArbInt;
970
 
 
971
 
begin
972
 
  if x <= 0
973
 
  then
974
 
    RunError(408);
975
 
  if x <= macheps
976
 
  then
977
 
    spelga:=-ln(x)
978
 
  else
979
 
  if x <= 15
980
 
  then
981
 
    begin
982
 
      m:=trunc(x); t:=x-m; m:=m-1; g:=1;
983
 
      if m < 0
984
 
      then
985
 
        g:=g/x
986
 
      else
987
 
      if m > 0
988
 
      then
989
 
        for i:=1 to m do
990
 
          g:=(x-i)*g;
991
 
      spelga:=ln(g*spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1))
992
 
    end
993
 
  else { x > 15 }
994
 
  if x <= xbig
995
 
  then
996
 
    spelga:=(x-0.5)*ln(x) - x + lnr2pi
997
 
            + spepol(450/sqr(x)-1, b[0], sizeof(b) div sizeof(ArbFloat) - 1)/x
998
 
  else { x > xbig }
999
 
  if x <= xmax
1000
 
  then
1001
 
    spelga:=(x-0.5)*ln(x) - x + lnr2pi
1002
 
  else  { x > xmax => x*ln(x) > giant }
1003
 
    RunError(408)
1004
 
end; {spelga}
1005
 
 
1006
 
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
1007
 
var   pa : ^arfloat0;
1008
 
       i : ArbInt;
1009
 
    polx : ArbFloat;
1010
 
begin
1011
 
  pa:=@a;
1012
 
  polx:=0;
1013
 
  for i:=n downto 0 do
1014
 
    polx:=polx*x+pa^[i];
1015
 
  spepol:=polx
1016
 
end {spepol};
1017
 
 
1018
 
function spepow(a, b: ArbFloat): ArbFloat;
1019
 
 
1020
 
   function PowInt(a: double; n: longint): double;
1021
 
   var a1 : double;
1022
 
   begin
1023
 
     if n=0 then PowInt := 1 else
1024
 
     begin
1025
 
        a1 := 1;
1026
 
        if n<0 then begin a := 1/a; n := -n end;
1027
 
        while n>0
1028
 
        do if Odd(n)
1029
 
           then begin Dec(n); a1 := a1*a end
1030
 
           else begin n := n div 2; a := sqr(a) end;
1031
 
        PowInt := a1
1032
 
     end
1033
 
   end;
1034
 
 
1035
 
var tb : longint;
1036
 
    fb : double;
1037
 
begin
1038
 
 
1039
 
  { (a < 0, b niet geheel) of (a = 0, b <= 0), dan afbreken}
1040
 
  if (a=0) then if (b<=0) then RunError(400) else begin SpePow :=0; exit end;
1041
 
  tb := Trunc(b); fb := b-tb;
1042
 
  if (a<0) and (fb<>0) then RunError(400);
1043
 
 
1044
 
  if a>0
1045
 
  then if fb=0 then SpePow := PowInt(a, tb)
1046
 
               else SpePow := PowInt(a, tb)*exp(fb*ln(a))
1047
 
  else if odd(tb) then Spepow := -PowInt(-a, tb)
1048
 
                  else SpePow := PowInt(-a, tb)
1049
 
 
1050
 
end; {spepow}
1051
 
 
1052
 
function spesgn(x: ArbFloat): ArbInt;
1053
 
 
1054
 
begin
1055
 
  if x<0
1056
 
  then
1057
 
    spesgn:=-1
1058
 
  else
1059
 
    if x=0
1060
 
    then
1061
 
      spesgn:=0
1062
 
    else
1063
 
      spesgn:=1
1064
 
end; {spesgn}
1065
 
 
1066
 
function spears(x: ArbFloat): ArbFloat;
1067
 
const
1068
 
 
1069
 
    pi2 = 1.570796326794897;
1070
 
    a : array[0..17] of ArbFloat =
1071
 
    (  1.047197551196598e+0, 5.375149359132719e-2, 7.798902238957732e-3,
1072
 
       1.519668539582420e-3, 3.408637238430600e-4, 8.302317819598986e-5,
1073
 
       2.134554822576075e-5, 5.701781046148566e-6, 1.566985123962741e-6,
1074
 
       4.402076871418002e-7, 1.257811162594110e-7, 3.646577948300129e-8,
1075
 
       1.081021746966715e-8, 3.212744286269388e-9, 8.515014302985799e-10,
1076
 
       2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);
1077
 
 
1078
 
var    y, u, t, s : ArbFloat;
1079
 
    uprang        : boolean;
1080
 
begin
1081
 
  if abs(x) > 1
1082
 
  then
1083
 
    RunError(401);
1084
 
  u:=sqr(x); uprang:= u > 0.5;
1085
 
  if uprang
1086
 
  then
1087
 
    u:=1-u;
1088
 
  t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
1089
 
  if uprang
1090
 
  then
1091
 
    begin
1092
 
      s:=pi2-sqrt(u)*y;
1093
 
      if x < 0
1094
 
      then
1095
 
        s:=-s;
1096
 
      spears:=s
1097
 
    end
1098
 
  else
1099
 
    spears:=x*y
1100
 
end;  {spears}
1101
 
 
1102
 
function spearc(x: ArbFloat): ArbFloat;
1103
 
const
1104
 
 
1105
 
    pi2 = 1.570796326794897;
1106
 
    a : array[0..17] of ArbFloat =
1107
 
    ( 1.047197551196598e+0,  5.375149359132719e-2,  7.798902238957732e-3,
1108
 
      1.519668539582420e-3,  3.408637238430600e-4,  8.302317819598986e-5,
1109
 
      2.134554822576075e-5,  5.701781046148566e-6,  1.566985123962741e-6,
1110
 
      4.402076871418002e-7,  1.257811162594110e-7,  3.646577948300129e-8,
1111
 
      1.081021746966715e-8,  3.212744286269388e-9,  8.515014302985799e-10,
1112
 
      2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);
1113
 
 
1114
 
var u, t, y, s    : ArbFloat;
1115
 
    uprang        : boolean;
1116
 
begin
1117
 
  if abs(x)>1.0
1118
 
  then
1119
 
    RunError(402);
1120
 
  u:=sqr(x); uprang:=u>0.5;
1121
 
  if uprang
1122
 
  then
1123
 
    u:=1-u;
1124
 
  t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
1125
 
  if uprang
1126
 
  then
1127
 
    begin
1128
 
      s:=sqrt(u)*y;
1129
 
      if x<0
1130
 
      then
1131
 
        s:=2*pi2-s;
1132
 
      spearc:=s
1133
 
    end
1134
 
  else
1135
 
    spearc:=pi2-x*y
1136
 
end;  {spearc}
1137
 
 
1138
 
function spesih(x: ArbFloat): ArbFloat;
1139
 
const
1140
 
 
1141
 
    a : array[0..6] of ArbFloat =
1142
 
    ( 1.085441641272607e+0,  8.757509762437522e-2,  2.158779361257021e-3,
1143
 
      2.549839945498292e-5,  1.761854853281383e-7,  7.980704288665359e-10,
1144
 
      2.551377137317034e-12);
1145
 
 
1146
 
var u : ArbFloat;
1147
 
begin
1148
 
  if abs(x)<=1.0
1149
 
  then
1150
 
    begin
1151
 
      u:=2*sqr(x)-1;
1152
 
      spesih:=x*spepol(u, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
1153
 
    end
1154
 
  else
1155
 
  begin
1156
 
    u:=exp(x); spesih:=(u-1/u)/2
1157
 
  end
1158
 
end; {spesih}
1159
 
 
1160
 
function specoh(x: ArbFloat): ArbFloat;
1161
 
var u: ArbFloat;
1162
 
begin
1163
 
  u:=exp(x); specoh:=(u+1/u)/2
1164
 
end; {specoh}
1165
 
 
1166
 
function spetah(x: ArbFloat): ArbFloat;
1167
 
const
1168
 
    xhi = 18.50;
1169
 
    a : array[0..15] of ArbFloat =
1170
 
    ( 8.610571715805476e-1, -1.158834489728470e-1,  1.918072383973950e-2,
1171
 
     -3.225255180728459e-3,  5.433071386922689e-4, -9.154289983175165e-5,
1172
 
      1.542469328074432e-5, -2.599022539340038e-6,  4.379282308765732e-7,
1173
 
     -7.378980192173815e-8,  1.243517352745986e-8, -2.095373768837420e-9,
1174
 
      3.509758916273561e-10,-5.908745181531817e-11, 1.124199312776748e-11,
1175
 
     -1.907888434471600e-12);
1176
 
 
1177
 
var t, y: ArbFloat;
1178
 
 
1179
 
begin
1180
 
  t:=abs(x);
1181
 
  if t <= 1
1182
 
  then
1183
 
    begin
1184
 
      y:=2*sqr(x)-1;
1185
 
      spetah:=x*spepol(y, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
1186
 
    end
1187
 
  else
1188
 
  if t < xhi
1189
 
  then
1190
 
    begin
1191
 
      y:=exp(2*x); spetah:=(y-1)/(y+1)
1192
 
    end
1193
 
  else
1194
 
    spetah:=spesgn(x)
1195
 
end; {spetah}
1196
 
 
1197
 
function speash(x: ArbFloat): ArbFloat;
1198
 
const
1199
 
 
1200
 
    xhi = 1e9;
1201
 
    c : array[0..18] of ArbFloat =
1202
 
    (  9.312298594527122e-1,  -5.736663926249348e-2,
1203
 
       9.004288574881897e-3,  -1.833458667045431e-3,
1204
 
       4.230023450529706e-4,  -1.050715136470630e-4,
1205
 
       2.740790473603819e-5,  -7.402952157663977e-6,
1206
 
       2.052474396638805e-6,  -5.807433412373489e-7,
1207
 
       1.670117348345774e-7,  -4.863477336087045e-8,
1208
 
       1.432753532351304e-8,  -4.319978113584910e-9,
1209
 
       1.299779213740398e-9,  -3.394726871170490e-10,
1210
 
       1.008344962167889e-10, -5.731943029121004e-11,
1211
 
       1.810792296549804e-11);
1212
 
 
1213
 
 
1214
 
var t : ArbFloat;
1215
 
 
1216
 
begin
1217
 
  t:=abs(x);
1218
 
  if t <= 1 then
1219
 
    speash:=x*spepol(2*sqr(x)-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
1220
 
  else
1221
 
  if t < xhi then
1222
 
    speash:=ln(sqrt(sqr(x)+1)+t)*spesgn(x)
1223
 
  else
1224
 
    speash:=ln(2*t)*spesgn(x)
1225
 
end; {speash}
1226
 
 
1227
 
function speach(x: ArbFloat): ArbFloat;
1228
 
const
1229
 
 
1230
 
    xhi = 1e9;
1231
 
 
1232
 
begin
1233
 
  if x<1 then
1234
 
    RunError(405);
1235
 
  if x=1 then
1236
 
    speach:=0
1237
 
  else
1238
 
  if x<=xhi then
1239
 
    speach:=ln(x+sqrt(sqr(x)-1))
1240
 
  else
1241
 
    speach:=ln(2*x)
1242
 
end; {speach}
1243
 
 
1244
 
function speath(x: ArbFloat): ArbFloat;
1245
 
const
1246
 
 
1247
 
    c : array[0..19] of ArbFloat =
1248
 
    ( 1.098612288668110e+0,  1.173605223326117e-1,  2.309071936165689e-2,
1249
 
      5.449091889986991e-3,  1.404884102286929e-3,  3.816948426588058e-4,
1250
 
      1.073604335435426e-4,  3.095027782918129e-5,  9.088050814470148e-6,
1251
 
      2.706881064641104e-6,  8.155200644023077e-7,  2.479830612463254e-7,
1252
 
      7.588067811453948e-8,  2.339295963220429e-8,  7.408486568719348e-9,
1253
 
      2.319454882064018e-9,  5.960921368486746e-10, 1.820410351379402e-10,
1254
 
      1.184977617320312e-10, 3.856235316559190e-11);
1255
 
 
1256
 
var t, u: ArbFloat;
1257
 
begin
1258
 
  t:=abs(x);
1259
 
  if t >= 1 then
1260
 
    RunError(406);
1261
 
  u:=sqr(x);
1262
 
  if u < 0.5 then
1263
 
    speath:=x*spepol(4*u-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
1264
 
  else { 0.5 < x*x < 1 }
1265
 
    speath:=ln((1+t)/(1-t))/2*spesgn(x)
1266
 
end; {speath}
1267
 
 
1268
 
var exitsave : pointer;
1269
 
 
1270
 
procedure MyExit; Far;
1271
 
const ErrorS : array[400..408,1..6] of char =
1272
 
     ('spepow',
1273
 
      'spebk0',
1274
 
      'spebk1',
1275
 
      'speby0',
1276
 
      'speby1',
1277
 
      'speach',
1278
 
      'speath',
1279
 
      'spegam',
1280
 
      'spelga');
1281
 
 
1282
 
var ErrFil : text;
1283
 
 
1284
 
begin
1285
 
     ExitProc := ExitSave;
1286
 
     Assign(ErrFil, 'CON');
1287
 
     ReWrite(ErrFil);
1288
 
     if (ExitCode>=400) AND (ExitCode<=408) then
1289
 
       begin
1290
 
         write(ErrFil, 'critical error in ', ErrorS[ExitCode]);
1291
 
         ExitCode := 201
1292
 
       end;
1293
 
     Close(ErrFil)
1294
 
end;
1295
 
 
1296
 
begin
1297
 
   ExitSave := ExitProc;
1298
 
   ExitProc := @MyExit;
1299
 
end.