File : a-numaux-x86.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
5 -- A D A . N U M E R I C S . A U X --
6 -- --
7 -- B o d y --
8 -- (Machine Version for x86) --
9 -- --
10 -- Copyright (C) 1998-2014, Free Software Foundation, Inc. --
11 -- --
12 -- GNAT is free software; you can redistribute it and/or modify it under --
13 -- terms of the GNU General Public License as published by the Free Soft- --
14 -- ware Foundation; either version 3, or (at your option) any later ver- --
15 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
16 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
17 -- or FITNESS FOR A PARTICULAR PURPOSE. --
18 -- --
19 -- --
20 -- --
21 -- --
22 -- --
23 -- You should have received a copy of the GNU General Public License and --
24 -- a copy of the GCC Runtime Library Exception along with this program; --
25 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
26 -- <http://www.gnu.org/licenses/>. --
27 -- --
28 -- GNAT was originally developed by the GNAT team at New York University. --
29 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 -- --
31 ------------------------------------------------------------------------------
32
33 with System.Machine_Code; use System.Machine_Code;
34
35 package body Ada.Numerics.Aux is
36
37 NL : constant String := ASCII.LF & ASCII.HT;
38
39 -----------------------
40 -- Local subprograms --
41 -----------------------
42
43 function Is_Nan (X : Double) return Boolean;
44 -- Return True iff X is a IEEE NaN value
45
46 function Logarithmic_Pow (X, Y : Double) return Double;
47 -- Implementation of X**Y using Exp and Log functions (binary base)
48 -- to calculate the exponentiation. This is used by Pow for values
49 -- for values of Y in the open interval (-0.25, 0.25)
50
51 procedure Reduce (X : in out Double; Q : out Natural);
52 -- Implements reduction of X by Pi/2. Q is the quadrant of the final
53 -- result in the range 0 .. 3. The absolute value of X is at most Pi.
54
55 pragma Inline (Is_Nan);
56 pragma Inline (Reduce);
57
58 --------------------------------
59 -- Basic Elementary Functions --
60 --------------------------------
61
62 -- This section implements a few elementary functions that are used to
63 -- build the more complex ones. This ordering enables better inlining.
64
65 ----------
66 -- Atan --
67 ----------
68
69 function Atan (X : Double) return Double is
70 Result : Double;
71
72 begin
73 Asm (Template =>
74 "fld1" & NL
75 & "fpatan",
76 Outputs => Double'Asm_Output ("=t", Result),
77 Inputs => Double'Asm_Input ("0", X));
78
79 -- The result value is NaN iff input was invalid
80
81 if not (Result = Result) then
82 raise Argument_Error;
83 end if;
84
85 return Result;
86 end Atan;
87
88 ---------
89 -- Exp --
90 ---------
91
92 function Exp (X : Double) return Double is
93 Result : Double;
94 begin
95 Asm (Template =>
96 "fldl2e " & NL
97 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E)
98 & "fld %%st(0) " & NL
99 & "frndint " & NL -- Integer (X * Log2 (E))
100 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
101 & "fxch " & NL
102 & "f2xm1 " & NL -- 2**(...) - 1
103 & "fld1 " & NL
104 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
105 & "fscale " & NL -- E ** X
106 & "fstp %%st(1) ",
107 Outputs => Double'Asm_Output ("=t", Result),
108 Inputs => Double'Asm_Input ("0", X));
109 return Result;
110 end Exp;
111
112 ------------
113 -- Is_Nan --
114 ------------
115
116 function Is_Nan (X : Double) return Boolean is
117 begin
118 -- The IEEE NaN values are the only ones that do not equal themselves
119
120 return not (X = X);
121 end Is_Nan;
122
123 ---------
124 -- Log --
125 ---------
126
127 function Log (X : Double) return Double is
128 Result : Double;
129
130 begin
131 Asm (Template =>
132 "fldln2 " & NL
133 & "fxch " & NL
134 & "fyl2x " & NL,
135 Outputs => Double'Asm_Output ("=t", Result),
136 Inputs => Double'Asm_Input ("0", X));
137 return Result;
138 end Log;
139
140 ------------
141 -- Reduce --
142 ------------
143
144 procedure Reduce (X : in out Double; Q : out Natural) is
145 Half_Pi : constant := Pi / 2.0;
146 Two_Over_Pi : constant := 2.0 / Pi;
147
148 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
149 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
150 P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
151 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
152 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
153 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
154 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
155 - P4, HM);
156 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
157 K : Double := X * Two_Over_Pi;
158 begin
159 -- For X < 2.0**32, all products below are computed exactly.
160 -- Due to cancellation effects all subtractions are exact as well.
161 -- As no double extended floating-point number has more than 75
162 -- zeros after the binary point, the result will be the correctly
163 -- rounded result of X - K * (Pi / 2.0).
164
165 while abs K >= 2.0**HM loop
166 K := K * M - (K * M - K);
167 X := (((((X - K * P1) - K * P2) - K * P3)
168 - K * P4) - K * P5) - K * P6;
169 K := X * Two_Over_Pi;
170 end loop;
171
172 if K /= K then
173
174 -- K is not a number, because X was not finite
175
176 raise Constraint_Error;
177 end if;
178
179 K := Double'Rounding (K);
180 Q := Integer (K) mod 4;
181 X := (((((X - K * P1) - K * P2) - K * P3)
182 - K * P4) - K * P5) - K * P6;
183 end Reduce;
184
185 ----------
186 -- Sqrt --
187 ----------
188
189 function Sqrt (X : Double) return Double is
190 Result : Double;
191
192 begin
193 if X < 0.0 then
194 raise Argument_Error;
195 end if;
196
197 Asm (Template => "fsqrt",
198 Outputs => Double'Asm_Output ("=t", Result),
199 Inputs => Double'Asm_Input ("0", X));
200
201 return Result;
202 end Sqrt;
203
204 --------------------------------
205 -- Other Elementary Functions --
206 --------------------------------
207
208 -- These are built using the previously implemented basic functions
209
210 ----------
211 -- Acos --
212 ----------
213
214 function Acos (X : Double) return Double is
215 Result : Double;
216
217 begin
218 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
219
220 -- The result value is NaN iff input was invalid
221
222 if Is_Nan (Result) then
223 raise Argument_Error;
224 end if;
225
226 return Result;
227 end Acos;
228
229 ----------
230 -- Asin --
231 ----------
232
233 function Asin (X : Double) return Double is
234 Result : Double;
235
236 begin
237 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
238
239 -- The result value is NaN iff input was invalid
240
241 if Is_Nan (Result) then
242 raise Argument_Error;
243 end if;
244
245 return Result;
246 end Asin;
247
248 ---------
249 -- Cos --
250 ---------
251
252 function Cos (X : Double) return Double is
253 Reduced_X : Double := abs X;
254 Result : Double;
255 Quadrant : Natural range 0 .. 3;
256
257 begin
258 if Reduced_X > Pi / 4.0 then
259 Reduce (Reduced_X, Quadrant);
260
261 case Quadrant is
262 when 0 =>
263 Asm (Template => "fcos",
264 Outputs => Double'Asm_Output ("=t", Result),
265 Inputs => Double'Asm_Input ("0", Reduced_X));
266 when 1 =>
267 Asm (Template => "fsin",
268 Outputs => Double'Asm_Output ("=t", Result),
269 Inputs => Double'Asm_Input ("0", -Reduced_X));
270 when 2 =>
271 Asm (Template => "fcos ; fchs",
272 Outputs => Double'Asm_Output ("=t", Result),
273 Inputs => Double'Asm_Input ("0", Reduced_X));
274 when 3 =>
275 Asm (Template => "fsin",
276 Outputs => Double'Asm_Output ("=t", Result),
277 Inputs => Double'Asm_Input ("0", Reduced_X));
278 end case;
279
280 else
281 Asm (Template => "fcos",
282 Outputs => Double'Asm_Output ("=t", Result),
283 Inputs => Double'Asm_Input ("0", Reduced_X));
284 end if;
285
286 return Result;
287 end Cos;
288
289 ---------------------
290 -- Logarithmic_Pow --
291 ---------------------
292
293 function Logarithmic_Pow (X, Y : Double) return Double is
294 Result : Double;
295 begin
296 Asm (Template => "" -- X : Y
297 & "fyl2x " & NL -- Y * Log2 (X)
298 & "fld %%st(0) " & NL -- Y * Log2 (X) : Y * Log2 (X)
299 & "frndint " & NL -- Int (...) : Y * Log2 (X)
300 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...)
301 & "fxch " & NL -- Fract (...) : Int (...)
302 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...)
303 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...)
304 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...)
305 & "fscale ", -- 2**(Fract (...) + Int (...))
306 Outputs => Double'Asm_Output ("=t", Result),
307 Inputs =>
308 (Double'Asm_Input ("0", X),
309 Double'Asm_Input ("u", Y)));
310 return Result;
311 end Logarithmic_Pow;
312
313 ---------
314 -- Pow --
315 ---------
316
317 function Pow (X, Y : Double) return Double is
318 type Mantissa_Type is mod 2**Double'Machine_Mantissa;
319 -- Modular type that can hold all bits of the mantissa of Double
320
321 -- For negative exponents, do divide at the end of the processing
322
323 Negative_Y : constant Boolean := Y < 0.0;
324 Abs_Y : constant Double := abs Y;
325
326 -- During this function the following invariant is kept:
327 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
328
329 Base : Double := X;
330
331 Exp_High : Double := Double'Floor (Abs_Y);
332 Exp_Mid : Double;
333 Exp_Low : Double;
334 Exp_Int : Mantissa_Type;
335
336 Factor : Double := 1.0;
337
338 begin
339 -- Select algorithm for calculating Pow (integer cases fall through)
340
341 if Exp_High >= 2.0**Double'Machine_Mantissa then
342
343 -- In case of Y that is IEEE infinity, just raise constraint error
344
345 if Exp_High > Double'Safe_Last then
346 raise Constraint_Error;
347 end if;
348
349 -- Large values of Y are even integers and will stay integer
350 -- after division by two.
351
352 loop
353 -- Exp_Mid and Exp_Low are zero, so
354 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
355
356 Exp_High := Exp_High / 2.0;
357 Base := Base * Base;
358 exit when Exp_High < 2.0**Double'Machine_Mantissa;
359 end loop;
360
361 elsif Exp_High /= Abs_Y then
362 Exp_Low := Abs_Y - Exp_High;
363 Factor := 1.0;
364
365 if Exp_Low /= 0.0 then
366
367 -- Exp_Low now is in interval (0.0, 1.0)
368 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
369
370 Exp_Mid := 0.0;
371 Exp_Low := Exp_Low - Exp_Mid;
372
373 if Exp_Low >= 0.5 then
374 Factor := Sqrt (X);
375 Exp_Low := Exp_Low - 0.5; -- exact
376
377 if Exp_Low >= 0.25 then
378 Factor := Factor * Sqrt (Factor);
379 Exp_Low := Exp_Low - 0.25; -- exact
380 end if;
381
382 elsif Exp_Low >= 0.25 then
383 Factor := Sqrt (Sqrt (X));
384 Exp_Low := Exp_Low - 0.25; -- exact
385 end if;
386
387 -- Exp_Low now is in interval (0.0, 0.25)
388
389 -- This means it is safe to call Logarithmic_Pow
390 -- for the remaining part.
391
392 Factor := Factor * Logarithmic_Pow (X, Exp_Low);
393 end if;
394
395 elsif X = 0.0 then
396 return 0.0;
397 end if;
398
399 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
400
401 Exp_Int := Mantissa_Type (Exp_High);
402
403 -- Standard way for processing integer powers > 0
404
405 while Exp_Int > 1 loop
406 if (Exp_Int and 1) = 1 then
407
408 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
409
410 Factor := Factor * Base;
411 end if;
412
413 -- Exp_Int is even and Exp_Int > 0, so
414 -- Base**Y = (Base**2)**(Exp_Int / 2)
415
416 Base := Base * Base;
417 Exp_Int := Exp_Int / 2;
418 end loop;
419
420 -- Exp_Int = 1 or Exp_Int = 0
421
422 if Exp_Int = 1 then
423 Factor := Base * Factor;
424 end if;
425
426 if Negative_Y then
427 Factor := 1.0 / Factor;
428 end if;
429
430 return Factor;
431 end Pow;
432
433 ---------
434 -- Sin --
435 ---------
436
437 function Sin (X : Double) return Double is
438 Reduced_X : Double := X;
439 Result : Double;
440 Quadrant : Natural range 0 .. 3;
441
442 begin
443 if abs X > Pi / 4.0 then
444 Reduce (Reduced_X, Quadrant);
445
446 case Quadrant is
447 when 0 =>
448 Asm (Template => "fsin",
449 Outputs => Double'Asm_Output ("=t", Result),
450 Inputs => Double'Asm_Input ("0", Reduced_X));
451 when 1 =>
452 Asm (Template => "fcos",
453 Outputs => Double'Asm_Output ("=t", Result),
454 Inputs => Double'Asm_Input ("0", Reduced_X));
455 when 2 =>
456 Asm (Template => "fsin",
457 Outputs => Double'Asm_Output ("=t", Result),
458 Inputs => Double'Asm_Input ("0", -Reduced_X));
459 when 3 =>
460 Asm (Template => "fcos ; fchs",
461 Outputs => Double'Asm_Output ("=t", Result),
462 Inputs => Double'Asm_Input ("0", Reduced_X));
463 end case;
464
465 else
466 Asm (Template => "fsin",
467 Outputs => Double'Asm_Output ("=t", Result),
468 Inputs => Double'Asm_Input ("0", Reduced_X));
469 end if;
470
471 return Result;
472 end Sin;
473
474 ---------
475 -- Tan --
476 ---------
477
478 function Tan (X : Double) return Double is
479 Reduced_X : Double := X;
480 Result : Double;
481 Quadrant : Natural range 0 .. 3;
482
483 begin
484 if abs X > Pi / 4.0 then
485 Reduce (Reduced_X, Quadrant);
486
487 if Quadrant mod 2 = 0 then
488 Asm (Template => "fptan" & NL
489 & "ffree %%st(0)" & NL
490 & "fincstp",
491 Outputs => Double'Asm_Output ("=t", Result),
492 Inputs => Double'Asm_Input ("0", Reduced_X));
493 else
494 Asm (Template => "fsincos" & NL
495 & "fdivp %%st, %%st(1)" & NL
496 & "fchs",
497 Outputs => Double'Asm_Output ("=t", Result),
498 Inputs => Double'Asm_Input ("0", Reduced_X));
499 end if;
500
501 else
502 Asm (Template =>
503 "fptan " & NL
504 & "ffree %%st(0) " & NL
505 & "fincstp ",
506 Outputs => Double'Asm_Output ("=t", Result),
507 Inputs => Double'Asm_Input ("0", Reduced_X));
508 end if;
509
510 return Result;
511 end Tan;
512
513 ----------
514 -- Sinh --
515 ----------
516
517 function Sinh (X : Double) return Double is
518 begin
519 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
520
521 if abs X < 25.0 then
522 return (Exp (X) - Exp (-X)) / 2.0;
523 else
524 return Exp (X) / 2.0;
525 end if;
526 end Sinh;
527
528 ----------
529 -- Cosh --
530 ----------
531
532 function Cosh (X : Double) return Double is
533 begin
534 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
535
536 if abs X < 22.0 then
537 return (Exp (X) + Exp (-X)) / 2.0;
538 else
539 return Exp (X) / 2.0;
540 end if;
541 end Cosh;
542
543 ----------
544 -- Tanh --
545 ----------
546
547 function Tanh (X : Double) return Double is
548 begin
549 -- Return the Hyperbolic Tangent of x
550
551 -- x -x
552 -- e - e Sinh (X)
553 -- Tanh (X) is defined to be ----------- = --------
554 -- x -x Cosh (X)
555 -- e + e
556
557 if abs X > 23.0 then
558 return Double'Copy_Sign (1.0, X);
559 end if;
560
561 return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
562 end Tanh;
563
564 end Ada.Numerics.Aux;