← Back to branch summary

~ubuntu-branches/ubuntu/saucy/lazarus/saucy

« back to all changes in this revision

Viewing changes to components/tachart/numlib_fix/spe.pas

  • Committer: Package Import Robot
  • Author(s): Paul Gevers, Abou Al Montacir, Bart Martens, Paul Gevers
  • Date: 2013-06-08 14:12:17 UTC
  • mfrom: (1.1.9)
  • Revision ID: package-import@ubuntu.com-20130608141217-7k0cy9id8ifcnutc
Tags: 1.0.8+dfsg-1
http://bugs.debian.org/708695
[ Abou Al Montacir ]
* New upstream major release and multiple maintenace release offering many
  fixes and new features marking a new milestone for the Lazarus development
  and its stability level.
  - The detailed list of changes can be found here:
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_release_notes
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_fixes_branch
* LCL changes:
  - LCL is now a normal package.
      + Platform independent parts of the LCL are now in the package LCLBase
      + LCL is automatically recompiled when switching the target platform,
        unless pre-compiled binaries for this target are already installed.
      + No impact on existing projects.
      + Linker options needed by LCL are no more added to projects that do
        not use the LCL package.
  - Minor changes in LCL basic classes behaviour
      + TCustomForm.Create raises an exception if a form resource is not
        found.
      + TNotebook and TPage: a new implementation of these classes was added.
      + TDBNavigator: It is now possible to have focusable buttons by setting
        Options = [navFocusableButtons] and TabStop = True, useful for
        accessibility and for devices with neither mouse nor touch screen.
      + Names of TControlBorderSpacing.GetSideSpace and GetSpace were swapped
        and are now consistent. GetSideSpace = Around + GetSpace.
      + TForm.WindowState=wsFullscreen was added
      + TCanvas.TextFitInfo was added to calculate how many characters will
        fit into a specified Width. Useful for word-wrapping calculations.
      + TControl.GetColorResolvingParent and
        TControl.GetRGBColorResolvingParent were added, simplifying the work
        to obtain the final color of the control while resolving clDefault
        and the ParentColor.
      + LCLIntf.GetTextExtentExPoint now has a good default implementation
        which works in any platform not providing a specific implementation.
        However, Widgetset specific implementation is better, when available.
      + TTabControl was reorganized. Now it has the correct class hierarchy
        and inherits from TCustomTabControl as it should.
  - New unit in the LCL:
      + lazdialogs.pas: adds non-native versions of various native dialogs,
        for example TLazOpenDialog, TLazSaveDialog, TLazSelectDirectoryDialog.
        It is used by widgetsets which either do not have a native dialog, or
        do not wish to use it because it is limited. These dialogs can also be
        used by user applications directly.
      + lazdeviceapis.pas: offers an interface to more hardware devices such
        as the accelerometer, GPS, etc. See LazDeviceAPIs
      + lazcanvas.pas: provides a TFPImageCanvas descendent implementing
        drawing in a LCL-compatible way, but 100% in Pascal.
      + lazregions.pas. LazRegions is a wholly Pascal implementation of
        regions for canvas clipping, event clipping, finding in which control
        of a region tree one an event should reach, for drawing polygons, etc.
      + customdrawncontrols.pas, customdrawndrawers.pas,
        customdrawn_common.pas, customdrawn_android.pas and
        customdrawn_winxp.pas: are the Lazarus Custom Drawn Controls -controls
        which imitate the standard LCL ones, but with the difference that they
        are non-native and support skinning.
  - New APIs added to the LCL to improve support of accessibility software
    such as screen readers.
* IDE changes:
  - Many improvments.
  - The detailed list of changes can be found here:
    http://wiki.lazarus.freepascal.org/New_IDE_features_since#v1.0_.282012-08-29.29
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_release_notes#IDE_Changes
* Debugger / Editor changes:
  - Added pascal sources and breakpoints to the disassembler
  - Added threads dialog.
* Components changes:
  - TAChart: many fixes and new features
  - CodeTool: support Delphi style generics and new syntax extensions.
  - AggPas: removed to honor free licencing. (Closes: Bug#708695)
[Bart Martens]
* New debian/watch file fixing issues with upstream RC release.
[Abou Al Montacir]
* Avoid changing files in .pc hidden directory, these are used by quilt for
  internal purpose and could lead to surprises during build.
[Paul Gevers]
* Updated get-orig-source target and it compinion script orig-tar.sh so that they
  repack the source file, allowing bug 708695 to be fixed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
components/tachart/numlib_fix/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;
 
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
begin
 
1283
     ExitProc := ExitSave;
 
1284
     if (ExitCode>=400) AND (ExitCode<=408) then
 
1285
       begin
 
1286
         //write(ErrFil, 'critical error in ', ErrorS[ExitCode]);
 
1287
         ExitCode := 201
 
1288
       end;
 
1289
end;
 
1290
 
 
1291
begin
 
1292
   ExitSave := ExitProc;
 
1293
   ExitProc := @MyExit;
 
1294
end.