|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
131
|
9 -- Copyright (C) 1992-2018, Free Software Foundation, Inc. --
|
|
111
|
10 -- --
|
|
|
11 -- GNAT is free software; you can redistribute it and/or modify it under --
|
|
|
12 -- terms of the GNU General Public License as published by the Free Soft- --
|
|
|
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
|
|
|
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
|
|
|
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
|
|
|
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
|
|
|
17 -- --
|
|
|
18 -- As a special exception under Section 7 of GPL version 3, you are granted --
|
|
|
19 -- additional permissions described in the GCC Runtime Library Exception, --
|
|
|
20 -- version 3.1, as published by the Free Software Foundation. --
|
|
|
21 -- --
|
|
|
22 -- You should have received a copy of the GNU General Public License and --
|
|
|
23 -- a copy of the GCC Runtime Library Exception along with this program; --
|
|
|
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
|
|
|
25 -- <http://www.gnu.org/licenses/>. --
|
|
|
26 -- --
|
|
|
27 -- GNAT was originally developed by the GNAT team at New York University. --
|
|
|
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
|
|
|
29 -- --
|
|
|
30 ------------------------------------------------------------------------------
|
|
|
31
|
|
|
32 -- This body is specifically for using an Ada interface to C math.h to get
|
|
|
33 -- the computation engine. Many special cases are handled locally to avoid
|
|
|
34 -- unnecessary calls or to meet Annex G strict mode requirements.
|
|
|
35
|
|
|
36 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, sinh,
|
|
|
37 -- cosh, tanh from C library via math.h
|
|
|
38
|
|
|
39 with Ada.Numerics.Aux;
|
|
|
40
|
|
|
41 package body Ada.Numerics.Generic_Elementary_Functions with
|
|
|
42 SPARK_Mode => Off
|
|
|
43 is
|
|
|
44
|
|
|
45 use type Ada.Numerics.Aux.Double;
|
|
|
46
|
|
|
47 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
|
|
|
48 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
|
|
|
49
|
|
|
50 Half_Log_Two : constant := Log_Two / 2;
|
|
|
51
|
|
|
52 subtype T is Float_Type'Base;
|
|
|
53 subtype Double is Aux.Double;
|
|
|
54
|
|
|
55 Two_Pi : constant T := 2.0 * Pi;
|
|
|
56 Half_Pi : constant T := Pi / 2.0;
|
|
|
57
|
|
|
58 Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two;
|
|
|
59 Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two;
|
|
|
60 Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
|
|
|
61
|
|
|
62 -----------------------
|
|
|
63 -- Local Subprograms --
|
|
|
64 -----------------------
|
|
|
65
|
|
|
66 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base;
|
|
|
67 -- Cody/Waite routine, supposedly more precise than the library version.
|
|
|
68 -- Currently only needed for Sinh/Cosh on X86 with the largest FP type.
|
|
|
69
|
|
|
70 function Local_Atan
|
|
|
71 (Y : Float_Type'Base;
|
|
|
72 X : Float_Type'Base := 1.0) return Float_Type'Base;
|
|
|
73 -- Common code for arc tangent after cycle reduction
|
|
|
74
|
|
|
75 ----------
|
|
|
76 -- "**" --
|
|
|
77 ----------
|
|
|
78
|
|
|
79 function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
|
|
|
80 A_Right : Float_Type'Base;
|
|
|
81 Int_Part : Integer;
|
|
|
82 Result : Float_Type'Base;
|
|
|
83 R1 : Float_Type'Base;
|
|
|
84 Rest : Float_Type'Base;
|
|
|
85
|
|
|
86 begin
|
|
|
87 if Left = 0.0
|
|
|
88 and then Right = 0.0
|
|
|
89 then
|
|
|
90 raise Argument_Error;
|
|
|
91
|
|
|
92 elsif Left < 0.0 then
|
|
|
93 raise Argument_Error;
|
|
|
94
|
|
|
95 elsif Right = 0.0 then
|
|
|
96 return 1.0;
|
|
|
97
|
|
|
98 elsif Left = 0.0 then
|
|
|
99 if Right < 0.0 then
|
|
|
100 raise Constraint_Error;
|
|
|
101 else
|
|
|
102 return 0.0;
|
|
|
103 end if;
|
|
|
104
|
|
|
105 elsif Left = 1.0 then
|
|
|
106 return 1.0;
|
|
|
107
|
|
|
108 elsif Right = 1.0 then
|
|
|
109 return Left;
|
|
|
110
|
|
|
111 else
|
|
|
112 begin
|
|
|
113 if Right = 2.0 then
|
|
|
114 return Left * Left;
|
|
|
115
|
|
|
116 elsif Right = 0.5 then
|
|
|
117 return Sqrt (Left);
|
|
|
118
|
|
|
119 else
|
|
|
120 A_Right := abs (Right);
|
|
|
121
|
|
|
122 -- If exponent is larger than one, compute integer exponen-
|
|
|
123 -- tiation if possible, and evaluate fractional part with more
|
|
|
124 -- precision. The relative error is now proportional to the
|
|
|
125 -- fractional part of the exponent only.
|
|
|
126
|
|
|
127 if A_Right > 1.0
|
|
|
128 and then A_Right < Float_Type'Base (Integer'Last)
|
|
|
129 then
|
|
|
130 Int_Part := Integer (Float_Type'Base'Truncation (A_Right));
|
|
|
131 Result := Left ** Int_Part;
|
|
|
132 Rest := A_Right - Float_Type'Base (Int_Part);
|
|
|
133
|
|
|
134 -- Compute with two leading bits of the mantissa using
|
|
|
135 -- square roots. Bound to be better than logarithms, and
|
|
|
136 -- easily extended to greater precision.
|
|
|
137
|
|
|
138 if Rest >= 0.5 then
|
|
|
139 R1 := Sqrt (Left);
|
|
|
140 Result := Result * R1;
|
|
|
141 Rest := Rest - 0.5;
|
|
|
142
|
|
|
143 if Rest >= 0.25 then
|
|
|
144 Result := Result * Sqrt (R1);
|
|
|
145 Rest := Rest - 0.25;
|
|
|
146 end if;
|
|
|
147
|
|
|
148 elsif Rest >= 0.25 then
|
|
|
149 Result := Result * Sqrt (Sqrt (Left));
|
|
|
150 Rest := Rest - 0.25;
|
|
|
151 end if;
|
|
|
152
|
|
|
153 Result := Result *
|
|
|
154 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));
|
|
|
155
|
|
|
156 if Right >= 0.0 then
|
|
|
157 return Result;
|
|
|
158 else
|
|
|
159 return (1.0 / Result);
|
|
|
160 end if;
|
|
|
161 else
|
|
|
162 return
|
|
|
163 Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
|
|
|
164 end if;
|
|
|
165 end if;
|
|
|
166
|
|
|
167 exception
|
|
|
168 when others =>
|
|
|
169 raise Constraint_Error;
|
|
|
170 end;
|
|
|
171 end if;
|
|
|
172 end "**";
|
|
|
173
|
|
|
174 ------------
|
|
|
175 -- Arccos --
|
|
|
176 ------------
|
|
|
177
|
|
|
178 -- Natural cycle
|
|
|
179
|
|
|
180 function Arccos (X : Float_Type'Base) return Float_Type'Base is
|
|
|
181 Temp : Float_Type'Base;
|
|
|
182
|
|
|
183 begin
|
|
|
184 if abs X > 1.0 then
|
|
|
185 raise Argument_Error;
|
|
|
186
|
|
|
187 elsif abs X < Sqrt_Epsilon then
|
|
|
188 return Pi / 2.0 - X;
|
|
|
189
|
|
|
190 elsif X = 1.0 then
|
|
|
191 return 0.0;
|
|
|
192
|
|
|
193 elsif X = -1.0 then
|
|
|
194 return Pi;
|
|
|
195 end if;
|
|
|
196
|
|
|
197 Temp := Float_Type'Base (Aux.Acos (Double (X)));
|
|
|
198
|
|
|
199 if Temp < 0.0 then
|
|
|
200 Temp := Pi + Temp;
|
|
|
201 end if;
|
|
|
202
|
|
|
203 return Temp;
|
|
|
204 end Arccos;
|
|
|
205
|
|
|
206 -- Arbitrary cycle
|
|
|
207
|
|
|
208 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
209 Temp : Float_Type'Base;
|
|
|
210
|
|
|
211 begin
|
|
|
212 if Cycle <= 0.0 then
|
|
|
213 raise Argument_Error;
|
|
|
214
|
|
|
215 elsif abs X > 1.0 then
|
|
|
216 raise Argument_Error;
|
|
|
217
|
|
|
218 elsif abs X < Sqrt_Epsilon then
|
|
|
219 return Cycle / 4.0;
|
|
|
220
|
|
|
221 elsif X = 1.0 then
|
|
|
222 return 0.0;
|
|
|
223
|
|
|
224 elsif X = -1.0 then
|
|
|
225 return Cycle / 2.0;
|
|
|
226 end if;
|
|
|
227
|
|
|
228 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);
|
|
|
229
|
|
|
230 if Temp < 0.0 then
|
|
|
231 Temp := Cycle / 2.0 + Temp;
|
|
|
232 end if;
|
|
|
233
|
|
|
234 return Temp;
|
|
|
235 end Arccos;
|
|
|
236
|
|
|
237 -------------
|
|
|
238 -- Arccosh --
|
|
|
239 -------------
|
|
|
240
|
|
|
241 function Arccosh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
242 begin
|
|
|
243 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or the proper
|
|
|
244 -- approximation for X close to 1 or >> 1.
|
|
|
245
|
|
|
246 if X < 1.0 then
|
|
|
247 raise Argument_Error;
|
|
|
248
|
|
|
249 elsif X < 1.0 + Sqrt_Epsilon then
|
|
|
250 return Sqrt (2.0 * (X - 1.0));
|
|
|
251
|
|
|
252 elsif X > 1.0 / Sqrt_Epsilon then
|
|
|
253 return Log (X) + Log_Two;
|
|
|
254
|
|
|
255 else
|
|
|
256 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
|
|
|
257 end if;
|
|
|
258 end Arccosh;
|
|
|
259
|
|
|
260 ------------
|
|
|
261 -- Arccot --
|
|
|
262 ------------
|
|
|
263
|
|
|
264 -- Natural cycle
|
|
|
265
|
|
|
266 function Arccot
|
|
|
267 (X : Float_Type'Base;
|
|
|
268 Y : Float_Type'Base := 1.0)
|
|
|
269 return Float_Type'Base
|
|
|
270 is
|
|
|
271 begin
|
|
|
272 -- Just reverse arguments
|
|
|
273
|
|
|
274 return Arctan (Y, X);
|
|
|
275 end Arccot;
|
|
|
276
|
|
|
277 -- Arbitrary cycle
|
|
|
278
|
|
|
279 function Arccot
|
|
|
280 (X : Float_Type'Base;
|
|
|
281 Y : Float_Type'Base := 1.0;
|
|
|
282 Cycle : Float_Type'Base)
|
|
|
283 return Float_Type'Base
|
|
|
284 is
|
|
|
285 begin
|
|
|
286 -- Just reverse arguments
|
|
|
287
|
|
|
288 return Arctan (Y, X, Cycle);
|
|
|
289 end Arccot;
|
|
|
290
|
|
|
291 -------------
|
|
|
292 -- Arccoth --
|
|
|
293 -------------
|
|
|
294
|
|
|
295 function Arccoth (X : Float_Type'Base) return Float_Type'Base is
|
|
|
296 begin
|
|
|
297 if abs X > 2.0 then
|
|
|
298 return Arctanh (1.0 / X);
|
|
|
299
|
|
|
300 elsif abs X = 1.0 then
|
|
|
301 raise Constraint_Error;
|
|
|
302
|
|
|
303 elsif abs X < 1.0 then
|
|
|
304 raise Argument_Error;
|
|
|
305
|
|
|
306 else
|
|
|
307 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the other
|
|
|
308 -- has error 0 or Epsilon.
|
|
|
309
|
|
|
310 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
|
|
|
311 end if;
|
|
|
312 end Arccoth;
|
|
|
313
|
|
|
314 ------------
|
|
|
315 -- Arcsin --
|
|
|
316 ------------
|
|
|
317
|
|
|
318 -- Natural cycle
|
|
|
319
|
|
|
320 function Arcsin (X : Float_Type'Base) return Float_Type'Base is
|
|
|
321 begin
|
|
|
322 if abs X > 1.0 then
|
|
|
323 raise Argument_Error;
|
|
|
324
|
|
|
325 elsif abs X < Sqrt_Epsilon then
|
|
|
326 return X;
|
|
|
327
|
|
|
328 elsif X = 1.0 then
|
|
|
329 return Pi / 2.0;
|
|
|
330
|
|
|
331 elsif X = -1.0 then
|
|
|
332 return -(Pi / 2.0);
|
|
|
333 end if;
|
|
|
334
|
|
|
335 return Float_Type'Base (Aux.Asin (Double (X)));
|
|
|
336 end Arcsin;
|
|
|
337
|
|
|
338 -- Arbitrary cycle
|
|
|
339
|
|
|
340 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
341 begin
|
|
|
342 if Cycle <= 0.0 then
|
|
|
343 raise Argument_Error;
|
|
|
344
|
|
|
345 elsif abs X > 1.0 then
|
|
|
346 raise Argument_Error;
|
|
|
347
|
|
|
348 elsif X = 0.0 then
|
|
|
349 return X;
|
|
|
350
|
|
|
351 elsif X = 1.0 then
|
|
|
352 return Cycle / 4.0;
|
|
|
353
|
|
|
354 elsif X = -1.0 then
|
|
|
355 return -(Cycle / 4.0);
|
|
|
356 end if;
|
|
|
357
|
|
|
358 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
|
|
|
359 end Arcsin;
|
|
|
360
|
|
|
361 -------------
|
|
|
362 -- Arcsinh --
|
|
|
363 -------------
|
|
|
364
|
|
|
365 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
366 begin
|
|
|
367 if abs X < Sqrt_Epsilon then
|
|
|
368 return X;
|
|
|
369
|
|
|
370 elsif X > 1.0 / Sqrt_Epsilon then
|
|
|
371 return Log (X) + Log_Two;
|
|
|
372
|
|
|
373 elsif X < -(1.0 / Sqrt_Epsilon) then
|
|
|
374 return -(Log (-X) + Log_Two);
|
|
|
375
|
|
|
376 elsif X < 0.0 then
|
|
|
377 return -Log (abs X + Sqrt (X * X + 1.0));
|
|
|
378
|
|
|
379 else
|
|
|
380 return Log (X + Sqrt (X * X + 1.0));
|
|
|
381 end if;
|
|
|
382 end Arcsinh;
|
|
|
383
|
|
|
384 ------------
|
|
|
385 -- Arctan --
|
|
|
386 ------------
|
|
|
387
|
|
|
388 -- Natural cycle
|
|
|
389
|
|
|
390 function Arctan
|
|
|
391 (Y : Float_Type'Base;
|
|
|
392 X : Float_Type'Base := 1.0)
|
|
|
393 return Float_Type'Base
|
|
|
394 is
|
|
|
395 begin
|
|
|
396 if X = 0.0 and then Y = 0.0 then
|
|
|
397 raise Argument_Error;
|
|
|
398
|
|
|
399 elsif Y = 0.0 then
|
|
|
400 if X > 0.0 then
|
|
|
401 return 0.0;
|
|
|
402 else -- X < 0.0
|
|
|
403 return Pi * Float_Type'Copy_Sign (1.0, Y);
|
|
|
404 end if;
|
|
|
405
|
|
|
406 elsif X = 0.0 then
|
|
|
407 return Float_Type'Copy_Sign (Half_Pi, Y);
|
|
|
408
|
|
|
409 else
|
|
|
410 return Local_Atan (Y, X);
|
|
|
411 end if;
|
|
|
412 end Arctan;
|
|
|
413
|
|
|
414 -- Arbitrary cycle
|
|
|
415
|
|
|
416 function Arctan
|
|
|
417 (Y : Float_Type'Base;
|
|
|
418 X : Float_Type'Base := 1.0;
|
|
|
419 Cycle : Float_Type'Base)
|
|
|
420 return Float_Type'Base
|
|
|
421 is
|
|
|
422 begin
|
|
|
423 if Cycle <= 0.0 then
|
|
|
424 raise Argument_Error;
|
|
|
425
|
|
|
426 elsif X = 0.0 and then Y = 0.0 then
|
|
|
427 raise Argument_Error;
|
|
|
428
|
|
|
429 elsif Y = 0.0 then
|
|
|
430 if X > 0.0 then
|
|
|
431 return 0.0;
|
|
|
432 else -- X < 0.0
|
|
|
433 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
|
|
|
434 end if;
|
|
|
435
|
|
|
436 elsif X = 0.0 then
|
|
|
437 return Float_Type'Copy_Sign (Cycle / 4.0, Y);
|
|
|
438
|
|
|
439 else
|
|
|
440 return Local_Atan (Y, X) * Cycle / Two_Pi;
|
|
|
441 end if;
|
|
|
442 end Arctan;
|
|
|
443
|
|
|
444 -------------
|
|
|
445 -- Arctanh --
|
|
|
446 -------------
|
|
|
447
|
|
|
448 function Arctanh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
449 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;
|
|
|
450
|
|
|
451 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa;
|
|
|
452
|
|
|
453 begin
|
|
|
454 -- The naive formula:
|
|
|
455
|
|
|
456 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X)
|
|
|
457
|
|
|
458 -- is not well-behaved numerically when X < 0.5 and when X is close
|
|
|
459 -- to one. The following is accurate but probably not optimal.
|
|
|
460
|
|
|
461 if abs X = 1.0 then
|
|
|
462 raise Constraint_Error;
|
|
|
463
|
|
|
464 elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then
|
|
|
465
|
|
|
466 if abs X >= 1.0 then
|
|
|
467 raise Argument_Error;
|
|
|
468 else
|
|
|
469
|
|
|
470 -- The one case that overflows if put through the method below:
|
|
|
471 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is
|
|
|
472 -- accurate. This simplifies to:
|
|
|
473
|
|
|
474 return Float_Type'Copy_Sign (
|
|
|
475 Half_Log_Two * Float_Type'Base (Mantissa + 1), X);
|
|
|
476 end if;
|
|
|
477
|
|
|
478 -- elsif abs X <= 0.5 then
|
|
|
479 -- why is above line commented out ???
|
|
|
480
|
|
|
481 else
|
|
|
482 -- Use several piecewise linear approximations. A is close to X,
|
|
|
483 -- chosen so 1.0 + A, 1.0 - A, and X - A are exact. The two scalings
|
|
|
484 -- remove the low-order bits of X.
|
|
|
485
|
|
|
486 A := Float_Type'Base'Scaling (
|
|
|
487 Float_Type'Base (Long_Long_Integer
|
|
|
488 (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa);
|
|
|
489
|
|
|
490 B := X - A; -- This is exact; abs B <= 2**(-Mantissa).
|
|
|
491 A_Plus_1 := 1.0 + A; -- This is exact.
|
|
|
492 A_From_1 := 1.0 - A; -- Ditto.
|
|
|
493 D := A_Plus_1 * A_From_1; -- 1 - A*A.
|
|
|
494
|
|
|
495 -- use one term of the series expansion:
|
|
|
496
|
|
|
497 -- f (x + e) = f(x) + e * f'(x) + ..
|
|
|
498
|
|
|
499 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is
|
|
|
500 -- A*(B/D)**2 (if a quadratic approximation is ever needed).
|
|
|
501
|
|
|
502 return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D;
|
|
|
503 end if;
|
|
|
504 end Arctanh;
|
|
|
505
|
|
|
506 ---------
|
|
|
507 -- Cos --
|
|
|
508 ---------
|
|
|
509
|
|
|
510 -- Natural cycle
|
|
|
511
|
|
|
512 function Cos (X : Float_Type'Base) return Float_Type'Base is
|
|
|
513 begin
|
|
|
514 if abs X < Sqrt_Epsilon then
|
|
|
515 return 1.0;
|
|
|
516 end if;
|
|
|
517
|
|
|
518 return Float_Type'Base (Aux.Cos (Double (X)));
|
|
|
519 end Cos;
|
|
|
520
|
|
|
521 -- Arbitrary cycle
|
|
|
522
|
|
|
523 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
524 begin
|
|
|
525 -- Just reuse the code for Sin. The potential small loss of speed is
|
|
|
526 -- negligible with proper (front-end) inlining.
|
|
|
527
|
|
|
528 return -Sin (abs X - Cycle * 0.25, Cycle);
|
|
|
529 end Cos;
|
|
|
530
|
|
|
531 ----------
|
|
|
532 -- Cosh --
|
|
|
533 ----------
|
|
|
534
|
|
|
535 function Cosh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
536 Lnv : constant Float_Type'Base := 8#0.542714#;
|
|
|
537 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
|
|
|
538 Y : constant Float_Type'Base := abs X;
|
|
|
539 Z : Float_Type'Base;
|
|
|
540
|
|
|
541 begin
|
|
|
542 if Y < Sqrt_Epsilon then
|
|
|
543 return 1.0;
|
|
|
544
|
|
|
545 elsif Y > Log_Inverse_Epsilon then
|
|
|
546 Z := Exp_Strict (Y - Lnv);
|
|
|
547 return (Z + V2minus1 * Z);
|
|
|
548
|
|
|
549 else
|
|
|
550 Z := Exp_Strict (Y);
|
|
|
551 return 0.5 * (Z + 1.0 / Z);
|
|
|
552 end if;
|
|
|
553
|
|
|
554 end Cosh;
|
|
|
555
|
|
|
556 ---------
|
|
|
557 -- Cot --
|
|
|
558 ---------
|
|
|
559
|
|
|
560 -- Natural cycle
|
|
|
561
|
|
|
562 function Cot (X : Float_Type'Base) return Float_Type'Base is
|
|
|
563 begin
|
|
|
564 if X = 0.0 then
|
|
|
565 raise Constraint_Error;
|
|
|
566
|
|
|
567 elsif abs X < Sqrt_Epsilon then
|
|
|
568 return 1.0 / X;
|
|
|
569 end if;
|
|
|
570
|
|
|
571 return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
|
|
|
572 end Cot;
|
|
|
573
|
|
|
574 -- Arbitrary cycle
|
|
|
575
|
|
|
576 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
577 T : Float_Type'Base;
|
|
|
578
|
|
|
579 begin
|
|
|
580 if Cycle <= 0.0 then
|
|
|
581 raise Argument_Error;
|
|
|
582 end if;
|
|
|
583
|
|
|
584 T := Float_Type'Base'Remainder (X, Cycle);
|
|
|
585
|
|
|
586 if T = 0.0 or else abs T = 0.5 * Cycle then
|
|
|
587 raise Constraint_Error;
|
|
|
588
|
|
|
589 elsif abs T < Sqrt_Epsilon then
|
|
|
590 return 1.0 / T;
|
|
|
591
|
|
|
592 elsif abs T = 0.25 * Cycle then
|
|
|
593 return 0.0;
|
|
|
594
|
|
|
595 else
|
|
|
596 T := T / Cycle * Two_Pi;
|
|
|
597 return Cos (T) / Sin (T);
|
|
|
598 end if;
|
|
|
599 end Cot;
|
|
|
600
|
|
|
601 ----------
|
|
|
602 -- Coth --
|
|
|
603 ----------
|
|
|
604
|
|
|
605 function Coth (X : Float_Type'Base) return Float_Type'Base is
|
|
|
606 begin
|
|
|
607 if X = 0.0 then
|
|
|
608 raise Constraint_Error;
|
|
|
609
|
|
|
610 elsif X < Half_Log_Epsilon then
|
|
|
611 return -1.0;
|
|
|
612
|
|
|
613 elsif X > -Half_Log_Epsilon then
|
|
|
614 return 1.0;
|
|
|
615
|
|
|
616 elsif abs X < Sqrt_Epsilon then
|
|
|
617 return 1.0 / X;
|
|
|
618 end if;
|
|
|
619
|
|
|
620 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
|
|
|
621 end Coth;
|
|
|
622
|
|
|
623 ---------
|
|
|
624 -- Exp --
|
|
|
625 ---------
|
|
|
626
|
|
|
627 function Exp (X : Float_Type'Base) return Float_Type'Base is
|
|
|
628 Result : Float_Type'Base;
|
|
|
629
|
|
|
630 begin
|
|
|
631 if X = 0.0 then
|
|
|
632 return 1.0;
|
|
|
633 end if;
|
|
|
634
|
|
|
635 Result := Float_Type'Base (Aux.Exp (Double (X)));
|
|
|
636
|
|
|
637 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
|
|
|
638 -- is False, then we can just leave it as an infinity (and indeed we
|
|
|
639 -- prefer to do so). But if Machine_Overflows is True, then we have
|
|
|
640 -- to raise a Constraint_Error exception as required by the RM.
|
|
|
641
|
|
|
642 if Float_Type'Machine_Overflows and then not Result'Valid then
|
|
|
643 raise Constraint_Error;
|
|
|
644 end if;
|
|
|
645
|
|
|
646 return Result;
|
|
|
647 end Exp;
|
|
|
648
|
|
|
649 ----------------
|
|
|
650 -- Exp_Strict --
|
|
|
651 ----------------
|
|
|
652
|
|
|
653 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
|
|
|
654 G : Float_Type'Base;
|
|
|
655 Z : Float_Type'Base;
|
|
|
656
|
|
|
657 P0 : constant := 0.25000_00000_00000_00000;
|
|
|
658 P1 : constant := 0.75753_18015_94227_76666E-2;
|
|
|
659 P2 : constant := 0.31555_19276_56846_46356E-4;
|
|
|
660
|
|
|
661 Q0 : constant := 0.5;
|
|
|
662 Q1 : constant := 0.56817_30269_85512_21787E-1;
|
|
|
663 Q2 : constant := 0.63121_89437_43985_02557E-3;
|
|
|
664 Q3 : constant := 0.75104_02839_98700_46114E-6;
|
|
|
665
|
|
|
666 C1 : constant := 8#0.543#;
|
|
|
667 C2 : constant := -2.1219_44400_54690_58277E-4;
|
|
|
668 Le : constant := 1.4426_95040_88896_34074;
|
|
|
669
|
|
|
670 XN : Float_Type'Base;
|
|
|
671 P, Q, R : Float_Type'Base;
|
|
|
672
|
|
|
673 begin
|
|
|
674 if X = 0.0 then
|
|
|
675 return 1.0;
|
|
|
676 end if;
|
|
|
677
|
|
|
678 XN := Float_Type'Base'Rounding (X * Le);
|
|
|
679 G := (X - XN * C1) - XN * C2;
|
|
|
680 Z := G * G;
|
|
|
681 P := G * ((P2 * Z + P1) * Z + P0);
|
|
|
682 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
|
|
|
683 R := 0.5 + P / (Q - P);
|
|
|
684
|
|
|
685 R := Float_Type'Base'Scaling (R, Integer (XN) + 1);
|
|
|
686
|
|
|
687 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
|
|
|
688 -- is False, then we can just leave it as an infinity (and indeed we
|
|
|
689 -- prefer to do so). But if Machine_Overflows is True, then we have to
|
|
|
690 -- raise a Constraint_Error exception as required by the RM.
|
|
|
691
|
|
|
692 if Float_Type'Machine_Overflows and then not R'Valid then
|
|
|
693 raise Constraint_Error;
|
|
|
694 else
|
|
|
695 return R;
|
|
|
696 end if;
|
|
|
697
|
|
|
698 end Exp_Strict;
|
|
|
699
|
|
|
700 ----------------
|
|
|
701 -- Local_Atan --
|
|
|
702 ----------------
|
|
|
703
|
|
|
704 function Local_Atan
|
|
|
705 (Y : Float_Type'Base;
|
|
|
706 X : Float_Type'Base := 1.0) return Float_Type'Base
|
|
|
707 is
|
|
|
708 Z : Float_Type'Base;
|
|
|
709 Raw_Atan : Float_Type'Base;
|
|
|
710
|
|
|
711 begin
|
|
|
712 Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));
|
|
|
713
|
|
|
714 Raw_Atan :=
|
|
|
715 (if Z < Sqrt_Epsilon then Z
|
|
|
716 elsif Z = 1.0 then Pi / 4.0
|
|
|
717 else Float_Type'Base (Aux.Atan (Double (Z))));
|
|
|
718
|
|
|
719 if abs Y > abs X then
|
|
|
720 Raw_Atan := Half_Pi - Raw_Atan;
|
|
|
721 end if;
|
|
|
722
|
|
|
723 if X > 0.0 then
|
|
|
724 return Float_Type'Copy_Sign (Raw_Atan, Y);
|
|
|
725 else
|
|
|
726 return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
|
|
|
727 end if;
|
|
|
728 end Local_Atan;
|
|
|
729
|
|
|
730 ---------
|
|
|
731 -- Log --
|
|
|
732 ---------
|
|
|
733
|
|
|
734 -- Natural base
|
|
|
735
|
|
|
736 function Log (X : Float_Type'Base) return Float_Type'Base is
|
|
|
737 begin
|
|
|
738 if X < 0.0 then
|
|
|
739 raise Argument_Error;
|
|
|
740
|
|
|
741 elsif X = 0.0 then
|
|
|
742 raise Constraint_Error;
|
|
|
743
|
|
|
744 elsif X = 1.0 then
|
|
|
745 return 0.0;
|
|
|
746 end if;
|
|
|
747
|
|
|
748 return Float_Type'Base (Aux.Log (Double (X)));
|
|
|
749 end Log;
|
|
|
750
|
|
|
751 -- Arbitrary base
|
|
|
752
|
|
|
753 function Log (X, Base : Float_Type'Base) return Float_Type'Base is
|
|
|
754 begin
|
|
|
755 if X < 0.0 then
|
|
|
756 raise Argument_Error;
|
|
|
757
|
|
|
758 elsif Base <= 0.0 or else Base = 1.0 then
|
|
|
759 raise Argument_Error;
|
|
|
760
|
|
|
761 elsif X = 0.0 then
|
|
|
762 raise Constraint_Error;
|
|
|
763
|
|
|
764 elsif X = 1.0 then
|
|
|
765 return 0.0;
|
|
|
766 end if;
|
|
|
767
|
|
|
768 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
|
|
|
769 end Log;
|
|
|
770
|
|
|
771 ---------
|
|
|
772 -- Sin --
|
|
|
773 ---------
|
|
|
774
|
|
|
775 -- Natural cycle
|
|
|
776
|
|
|
777 function Sin (X : Float_Type'Base) return Float_Type'Base is
|
|
|
778 begin
|
|
|
779 if abs X < Sqrt_Epsilon then
|
|
|
780 return X;
|
|
|
781 end if;
|
|
|
782
|
|
|
783 return Float_Type'Base (Aux.Sin (Double (X)));
|
|
|
784 end Sin;
|
|
|
785
|
|
|
786 -- Arbitrary cycle
|
|
|
787
|
|
|
788 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
789 T : Float_Type'Base;
|
|
|
790
|
|
|
791 begin
|
|
|
792 if Cycle <= 0.0 then
|
|
|
793 raise Argument_Error;
|
|
|
794
|
|
|
795 -- If X is zero, return it as the result, preserving the argument sign.
|
|
|
796 -- Is this test really needed on any machine ???
|
|
|
797
|
|
|
798 elsif X = 0.0 then
|
|
|
799 return X;
|
|
|
800 end if;
|
|
|
801
|
|
|
802 T := Float_Type'Base'Remainder (X, Cycle);
|
|
|
803
|
|
|
804 -- The following two reductions reduce the argument to the interval
|
|
|
805 -- [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed
|
|
|
806 -- to prevent inaccuracy that may result if the sine function uses a
|
|
|
807 -- different (more accurate) value of Pi in its reduction than is used
|
|
|
808 -- in the multiplication with Two_Pi.
|
|
|
809
|
|
|
810 if abs T > 0.25 * Cycle then
|
|
|
811 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
|
|
|
812 end if;
|
|
|
813
|
|
|
814 -- Could test for 12.0 * abs T = Cycle, and return an exact value in
|
|
|
815 -- those cases. It is not clear this is worth the extra test though.
|
|
|
816
|
|
|
817 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
|
|
|
818 end Sin;
|
|
|
819
|
|
|
820 ----------
|
|
|
821 -- Sinh --
|
|
|
822 ----------
|
|
|
823
|
|
|
824 function Sinh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
825 Lnv : constant Float_Type'Base := 8#0.542714#;
|
|
|
826 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
|
|
|
827 Y : constant Float_Type'Base := abs X;
|
|
|
828 F : constant Float_Type'Base := Y * Y;
|
|
|
829 Z : Float_Type'Base;
|
|
|
830
|
|
|
831 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;
|
|
|
832
|
|
|
833 begin
|
|
|
834 if Y < Sqrt_Epsilon then
|
|
|
835 return X;
|
|
|
836
|
|
|
837 elsif Y > Log_Inverse_Epsilon then
|
|
|
838 Z := Exp_Strict (Y - Lnv);
|
|
|
839 Z := Z + V2minus1 * Z;
|
|
|
840
|
|
|
841 elsif Y < 1.0 then
|
|
|
842
|
|
|
843 if Float_Digits_1_6 then
|
|
|
844
|
|
|
845 -- Use expansion provided by Cody and Waite, p. 226. Note that
|
|
|
846 -- leading term of the polynomial in Q is exactly 1.0.
|
|
|
847
|
|
|
848 declare
|
|
|
849 P0 : constant := -0.71379_3159E+1;
|
|
|
850 P1 : constant := -0.19033_3399E+0;
|
|
|
851 Q0 : constant := -0.42827_7109E+2;
|
|
|
852
|
|
|
853 begin
|
|
|
854 Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
|
|
|
855 end;
|
|
|
856
|
|
|
857 else
|
|
|
858 declare
|
|
|
859 P0 : constant := -0.35181_28343_01771_17881E+6;
|
|
|
860 P1 : constant := -0.11563_52119_68517_68270E+5;
|
|
|
861 P2 : constant := -0.16375_79820_26307_51372E+3;
|
|
|
862 P3 : constant := -0.78966_12741_73570_99479E+0;
|
|
|
863 Q0 : constant := -0.21108_77005_81062_71242E+7;
|
|
|
864 Q1 : constant := 0.36162_72310_94218_36460E+5;
|
|
|
865 Q2 : constant := -0.27773_52311_96507_01667E+3;
|
|
|
866
|
|
|
867 begin
|
|
|
868 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
|
|
|
869 / (((F + Q2) * F + Q1) * F + Q0);
|
|
|
870 end;
|
|
|
871 end if;
|
|
|
872
|
|
|
873 else
|
|
|
874 Z := Exp_Strict (Y);
|
|
|
875 Z := 0.5 * (Z - 1.0 / Z);
|
|
|
876 end if;
|
|
|
877
|
|
|
878 if X > 0.0 then
|
|
|
879 return Z;
|
|
|
880 else
|
|
|
881 return -Z;
|
|
|
882 end if;
|
|
|
883 end Sinh;
|
|
|
884
|
|
|
885 ----------
|
|
|
886 -- Sqrt --
|
|
|
887 ----------
|
|
|
888
|
|
|
889 function Sqrt (X : Float_Type'Base) return Float_Type'Base is
|
|
|
890 begin
|
|
|
891 if X < 0.0 then
|
|
|
892 raise Argument_Error;
|
|
|
893
|
|
|
894 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
|
|
|
895
|
|
|
896 elsif X = 0.0 then
|
|
|
897 return X;
|
|
|
898 end if;
|
|
|
899
|
|
|
900 return Float_Type'Base (Aux.Sqrt (Double (X)));
|
|
|
901 end Sqrt;
|
|
|
902
|
|
|
903 ---------
|
|
|
904 -- Tan --
|
|
|
905 ---------
|
|
|
906
|
|
|
907 -- Natural cycle
|
|
|
908
|
|
|
909 function Tan (X : Float_Type'Base) return Float_Type'Base is
|
|
|
910 begin
|
|
|
911 if abs X < Sqrt_Epsilon then
|
|
|
912 return X;
|
|
|
913 end if;
|
|
|
914
|
|
|
915 -- Note: if X is exactly pi/2, then we should raise an exception, since
|
|
|
916 -- the result would overflow. But for all floating-point formats we deal
|
|
|
917 -- with, it is impossible for X to be exactly pi/2, and the result is
|
|
|
918 -- always in range.
|
|
|
919
|
|
|
920 return Float_Type'Base (Aux.Tan (Double (X)));
|
|
|
921 end Tan;
|
|
|
922
|
|
|
923 -- Arbitrary cycle
|
|
|
924
|
|
|
925 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
|
|
|
926 T : Float_Type'Base;
|
|
|
927
|
|
|
928 begin
|
|
|
929 if Cycle <= 0.0 then
|
|
|
930 raise Argument_Error;
|
|
|
931
|
|
|
932 elsif X = 0.0 then
|
|
|
933 return X;
|
|
|
934 end if;
|
|
|
935
|
|
|
936 T := Float_Type'Base'Remainder (X, Cycle);
|
|
|
937
|
|
|
938 if abs T = 0.25 * Cycle then
|
|
|
939 raise Constraint_Error;
|
|
|
940
|
|
|
941 elsif abs T = 0.5 * Cycle then
|
|
|
942 return 0.0;
|
|
|
943
|
|
|
944 else
|
|
|
945 T := T / Cycle * Two_Pi;
|
|
|
946 return Sin (T) / Cos (T);
|
|
|
947 end if;
|
|
|
948
|
|
|
949 end Tan;
|
|
|
950
|
|
|
951 ----------
|
|
|
952 -- Tanh --
|
|
|
953 ----------
|
|
|
954
|
|
|
955 function Tanh (X : Float_Type'Base) return Float_Type'Base is
|
|
|
956 P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
|
|
|
957 P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
|
|
|
958 P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;
|
|
|
959
|
|
|
960 Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4;
|
|
|
961 Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4;
|
|
|
962 Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3;
|
|
|
963 Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1;
|
|
|
964
|
|
|
965 Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
|
|
|
966
|
|
|
967 P, Q, R : Float_Type'Base;
|
|
|
968 Y : constant Float_Type'Base := abs X;
|
|
|
969 G : constant Float_Type'Base := Y * Y;
|
|
|
970
|
|
|
971 Float_Type_Digits_15_Or_More : constant Boolean :=
|
|
|
972 Float_Type'Digits > 14;
|
|
|
973
|
|
|
974 begin
|
|
|
975 if X < Half_Log_Epsilon then
|
|
|
976 return -1.0;
|
|
|
977
|
|
|
978 elsif X > -Half_Log_Epsilon then
|
|
|
979 return 1.0;
|
|
|
980
|
|
|
981 elsif Y < Sqrt_Epsilon then
|
|
|
982 return X;
|
|
|
983
|
|
|
984 elsif Y < Half_Ln3
|
|
|
985 and then Float_Type_Digits_15_Or_More
|
|
|
986 then
|
|
|
987 P := (P2 * G + P1) * G + P0;
|
|
|
988 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
|
|
|
989 R := G * (P / Q);
|
|
|
990 return X + X * R;
|
|
|
991
|
|
|
992 else
|
|
|
993 return Float_Type'Base (Aux.Tanh (Double (X)));
|
|
|
994 end if;
|
|
|
995 end Tanh;
|
|
|
996
|
|
|
997 end Ada.Numerics.Generic_Elementary_Functions;
|
