1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME COMPONENTS --
5 -- A D A . N U M E R I C S . A U X --
8 -- (Machine Version for x86) --
11 -- Copyright (C) 1998-2001 Free Software Foundation, Inc. --
13 -- GNAT is free software; you can redistribute it and/or modify it under --
14 -- terms of the GNU General Public License as published by the Free Soft- --
15 -- ware Foundation; either version 2, or (at your option) any later ver- --
16 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
17 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
18 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
19 -- for more details. You should have received a copy of the GNU General --
20 -- Public License distributed with GNAT; see file COPYING. If not, write --
21 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
22 -- MA 02111-1307, USA. --
24 -- As a special exception, if other files instantiate generics from this --
25 -- unit, or you link this unit with other files to produce an executable, --
26 -- this unit does not by itself cause the resulting executable to be --
27 -- covered by the GNU General Public License. This exception does not --
28 -- however invalidate any other reasons why the executable file might be --
29 -- covered by the GNU Public License. --
31 -- GNAT was originally developed by the GNAT team at New York University. --
32 -- Extensive contributions were provided by Ada Core Technologies Inc. --
34 ------------------------------------------------------------------------------
36 -- File a-numaux.adb <- 86numaux.adb
38 -- This version of Numerics.Aux is for the IEEE Double Extended floating
39 -- point format on x86.
41 with System
.Machine_Code
; use System
.Machine_Code
;
43 package body Ada
.Numerics
.Aux
is
45 NL
: constant String := ASCII
.LF
& ASCII
.HT
;
47 type FPU_Stack_Pointer
is range 0 .. 7;
48 for FPU_Stack_Pointer
'Size use 3;
50 type FPU_Status_Word
is record
51 B
: Boolean; -- FPU Busy (for 8087 compatibility only)
52 ES
: Boolean; -- Error Summary Status
53 SF
: Boolean; -- Stack Fault
55 Top
: FPU_Stack_Pointer
;
57 -- Condition Code Flags
59 -- C2 is set by FPREM and FPREM1 to indicate incomplete reduction.
60 -- In case of successfull recorction, C0, C3 and C1 are set to the
61 -- three least significant bits of the result (resp. Q2, Q1 and Q0).
63 -- C2 is used by FPTAN, FSIN, FCOS, and FSINCOS to indicate that
64 -- that source operand is beyond the allowable range of
65 -- -2.0**63 .. 2.0**63.
74 PE
: Boolean; -- Precision
75 UE
: Boolean; -- Underflow
76 OE
: Boolean; -- Overflow
77 ZE
: Boolean; -- Zero Divide
78 DE
: Boolean; -- Denormalized Operand
79 IE
: Boolean; -- Invalid Operation
82 for FPU_Status_Word
use record
83 B
at 0 range 15 .. 15;
84 C3
at 0 range 14 .. 14;
85 Top
at 0 range 11 .. 13;
86 C2
at 0 range 10 .. 10;
99 for FPU_Status_Word
'Size use 16;
101 -----------------------
102 -- Local subprograms --
103 -----------------------
105 function Is_Nan
(X
: Double
) return Boolean;
106 -- Return True iff X is a IEEE NaN value
108 function Logarithmic_Pow
(X
, Y
: Double
) return Double
;
109 -- Implementation of X**Y using Exp and Log functions (binary base)
110 -- to calculate the exponentiation. This is used by Pow for values
111 -- for values of Y in the open interval (-0.25, 0.25)
113 function Reduce
(X
: Double
) return Double
;
114 -- Implement partial reduction of X by Pi in the x86.
116 -- Note that for the Sin, Cos and Tan functions completely accurate
117 -- reduction of the argument is done for arguments in the range of
118 -- -2.0**63 .. 2.0**63, using a 66-bit approximation of Pi.
120 pragma Inline
(Is_Nan
);
121 pragma Inline
(Reduce
);
123 ---------------------------------
124 -- Basic Elementary Functions --
125 ---------------------------------
127 -- This section implements a few elementary functions that are
128 -- used to build the more complex ones. This ordering enables
135 function Atan
(X
: Double
) return Double
is
142 Outputs
=> Double
'Asm_Output ("=t", Result
),
143 Inputs
=> Double
'Asm_Input ("0", X
));
145 -- The result value is NaN iff input was invalid
147 if not (Result
= Result
) then
148 raise Argument_Error
;
158 function Exp
(X
: Double
) return Double
is
163 & "fmulp %%st, %%st(1)" & NL
-- X * log2 (E)
164 & "fld %%st(0) " & NL
165 & "frndint " & NL
-- Integer (X * Log2 (E))
166 & "fsubr %%st, %%st(1)" & NL
-- Fraction (X * Log2 (E))
168 & "f2xm1 " & NL
-- 2**(...) - 1
170 & "faddp %%st, %%st(1)" & NL
-- 2**(Fraction (X * Log2 (E)))
171 & "fscale " & NL
-- E ** X
173 Outputs
=> Double
'Asm_Output ("=t", Result
),
174 Inputs
=> Double
'Asm_Input ("0", X
));
182 function Is_Nan
(X
: Double
) return Boolean is
184 -- The IEEE NaN values are the only ones that do not equal themselves
193 function Log
(X
: Double
) return Double
is
201 Outputs
=> Double
'Asm_Output ("=t", Result
),
202 Inputs
=> Double
'Asm_Input ("0", X
));
210 function Reduce
(X
: Double
) return Double
is
215 -- Partial argument reduction
217 & "fadd %%st(0), %%st" & NL
218 & "fxch %%st(1) " & NL
221 Outputs
=> Double
'Asm_Output ("=t", Result
),
222 Inputs
=> Double
'Asm_Input ("0", X
));
230 function Sqrt
(X
: Double
) return Double
is
235 raise Argument_Error
;
238 Asm
(Template
=> "fsqrt",
239 Outputs
=> Double
'Asm_Output ("=t", Result
),
240 Inputs
=> Double
'Asm_Input ("0", X
));
245 ---------------------------------
246 -- Other Elementary Functions --
247 ---------------------------------
249 -- These are built using the previously implemented basic functions
255 function Acos
(X
: Double
) return Double
is
258 Result
:= 2.0 * Atan
(Sqrt
((1.0 - X
) / (1.0 + X
)));
260 -- The result value is NaN iff input was invalid
262 if Is_Nan
(Result
) then
263 raise Argument_Error
;
273 function Asin
(X
: Double
) return Double
is
277 Result
:= Atan
(X
/ Sqrt
((1.0 - X
) * (1.0 + X
)));
279 -- The result value is NaN iff input was invalid
281 if Is_Nan
(Result
) then
282 raise Argument_Error
;
292 function Cos
(X
: Double
) return Double
is
293 Reduced_X
: Double
:= X
;
295 Status
: FPU_Status_Word
;
303 & "xorl %%eax, %%eax " & NL
305 Outputs
=> (Double
'Asm_Output ("=t", Result
),
306 FPU_Status_Word
'Asm_Output ("=a", Status
)),
307 Inputs
=> Double
'Asm_Input ("0", Reduced_X
));
309 exit when not Status
.C2
;
311 -- Original argument was not in range and the result
312 -- is the unmodified argument.
314 Reduced_X
:= Reduce
(Result
);
320 ---------------------
321 -- Logarithmic_Pow --
322 ---------------------
324 function Logarithmic_Pow
(X
, Y
: Double
) return Double
is
328 Asm
(Template
=> "" -- X : Y
329 & "fyl2x " & NL
-- Y * Log2 (X)
330 & "fst %%st(1) " & NL
-- Y * Log2 (X) : Y * Log2 (X)
331 & "frndint " & NL
-- Int (...) : Y * Log2 (X)
332 & "fsubr %%st, %%st(1)" & NL
-- Int (...) : Fract (...)
333 & "fxch " & NL
-- Fract (...) : Int (...)
334 & "f2xm1 " & NL
-- 2**Fract (...) - 1 : Int (...)
335 & "fld1 " & NL
-- 1 : 2**Fract (...) - 1 : Int (...)
336 & "faddp %%st, %%st(1)" & NL
-- 2**Fract (...) : Int (...)
337 & "fscale " & NL
-- 2**(Fract (...) + Int (...))
339 Outputs
=> Double
'Asm_Output ("=t", Result
),
341 (Double
'Asm_Input ("0", X
),
342 Double
'Asm_Input ("u", Y
)));
351 function Pow
(X
, Y
: Double
) return Double
is
352 type Mantissa_Type
is mod 2**Double
'Machine_Mantissa;
353 -- Modular type that can hold all bits of the mantissa of Double
355 -- For negative exponents, a division is done
356 -- at the end of the processing.
358 Negative_Y
: constant Boolean := Y
< 0.0;
359 Abs_Y
: constant Double
:= abs Y
;
361 -- During this function the following invariant is kept:
362 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
366 Exp_High
: Double
:= Double
'Floor (Abs_Y
);
369 Exp_Int
: Mantissa_Type
;
371 Factor
: Double
:= 1.0;
374 -- Select algorithm for calculating Pow:
375 -- integer cases fall through
377 if Exp_High
>= 2.0**Double
'Machine_Mantissa then
379 -- In case of Y that is IEEE infinity, just raise constraint error
381 if Exp_High
> Double
'Safe_Last then
382 raise Constraint_Error
;
385 -- Large values of Y are even integers and will stay integer
386 -- after division by two.
389 -- Exp_Mid and Exp_Low are zero, so
390 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
392 Exp_High
:= Exp_High
/ 2.0;
394 exit when Exp_High
< 2.0**Double
'Machine_Mantissa;
397 elsif Exp_High
/= Abs_Y
then
398 Exp_Low
:= Abs_Y
- Exp_High
;
402 if Exp_Low
/= 0.0 then
404 -- Exp_Low now is in interval (0.0, 1.0)
405 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
408 Exp_Low
:= Exp_Low
- Exp_Mid
;
410 if Exp_Low
>= 0.5 then
412 Exp_Low
:= Exp_Low
- 0.5; -- exact
414 if Exp_Low
>= 0.25 then
415 Factor
:= Factor
* Sqrt
(Factor
);
416 Exp_Low
:= Exp_Low
- 0.25; -- exact
419 elsif Exp_Low
>= 0.25 then
420 Factor
:= Sqrt
(Sqrt
(X
));
421 Exp_Low
:= Exp_Low
- 0.25; -- exact
424 -- Exp_Low now is in interval (0.0, 0.25)
426 -- This means it is safe to call Logarithmic_Pow
427 -- for the remaining part.
429 Factor
:= Factor
* Logarithmic_Pow
(X
, Exp_Low
);
436 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
438 Exp_Int
:= Mantissa_Type
(Exp_High
);
440 -- Standard way for processing integer powers > 0
442 while Exp_Int
> 1 loop
443 if (Exp_Int
and 1) = 1 then
445 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
447 Factor
:= Factor
* Base
;
450 -- Exp_Int is even and Exp_Int > 0, so
451 -- Base**Y = (Base**2)**(Exp_Int / 2)
454 Exp_Int
:= Exp_Int
/ 2;
457 -- Exp_Int = 1 or Exp_Int = 0
460 Factor
:= Base
* Factor
;
464 Factor
:= 1.0 / Factor
;
474 function Sin
(X
: Double
) return Double
is
475 Reduced_X
: Double
:= X
;
477 Status
: FPU_Status_Word
;
485 & "xorl %%eax, %%eax " & NL
487 Outputs
=> (Double
'Asm_Output ("=t", Result
),
488 FPU_Status_Word
'Asm_Output ("=a", Status
)),
489 Inputs
=> Double
'Asm_Input ("0", Reduced_X
));
491 exit when not Status
.C2
;
493 -- Original argument was not in range and the result
494 -- is the unmodified argument.
496 Reduced_X
:= Reduce
(Result
);
506 function Tan
(X
: Double
) return Double
is
507 Reduced_X
: Double
:= X
;
509 Status
: FPU_Status_Word
;
517 & "xorl %%eax, %%eax " & NL
518 & "fnstsw %%ax " & NL
519 & "ffree %%st(0) " & NL
522 Outputs
=> (Double
'Asm_Output ("=t", Result
),
523 FPU_Status_Word
'Asm_Output ("=a", Status
)),
524 Inputs
=> Double
'Asm_Input ("0", Reduced_X
));
526 exit when not Status
.C2
;
528 -- Original argument was not in range and the result
529 -- is the unmodified argument.
531 Reduced_X
:= Reduce
(Result
);
541 function Sinh
(X
: Double
) return Double
is
543 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
546 return (Exp
(X
) - Exp
(-X
)) / 2.0;
549 return Exp
(X
) / 2.0;
558 function Cosh
(X
: Double
) return Double
is
560 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
563 return (Exp
(X
) + Exp
(-X
)) / 2.0;
566 return Exp
(X
) / 2.0;
575 function Tanh
(X
: Double
) return Double
is
577 -- Return the Hyperbolic Tangent of x
581 -- Tanh (X) is defined to be ----------- = --------
586 return Double
'Copy_Sign (1.0, X
);
589 return 1.0 / (1.0 + Exp
(-2.0 * X
)) - 1.0 / (1.0 + Exp
(2.0 * X
));
593 end Ada
.Numerics
.Aux
;