Mercurial > hg > CbC > CbC_gcc

comparison gcc/ada/libgnat/a-numaux__x86.adb @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents
children 84e7813d76e9
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
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-2017, 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 -- As a special exception under Section 7 of GPL version 3, you are granted --
20 -- additional permissions described in the GCC Runtime Library Exception, --
21 -- version 3.1, as published by the Free Software Foundation. --
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 -- Implement 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/4.
54 -- It is needed to avoid a loss of accuracy for sin near Pi and cos
55 -- near Pi/2 due to the use of an insufficiently precise value of Pi
56 -- in the range reduction.
57
58 pragma Inline (Is_Nan);
59 pragma Inline (Reduce);
60
61 --------------------------------
62 -- Basic Elementary Functions --
63 --------------------------------
64
65 -- This section implements a few elementary functions that are used to
66 -- build the more complex ones. This ordering enables better inlining.
67
68 ----------
69 -- Atan --
70 ----------
71
72 function Atan (X : Double) return Double is
73 Result : Double;
74
75 begin
76 Asm (Template =>
77 "fld1" & NL
78 & "fpatan",
79 Outputs => Double'Asm_Output ("=t", Result),
80 Inputs => Double'Asm_Input ("0", X));
81
82 -- The result value is NaN iff input was invalid
83
84 if not (Result = Result) then
85 raise Argument_Error;
86 end if;
87
88 return Result;
89 end Atan;
90
91 ---------
92 -- Exp --
93 ---------
94
95 function Exp (X : Double) return Double is
96 Result : Double;
97 begin
98 Asm (Template =>
99 "fldl2e " & NL
100 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E)
101 & "fld %%st(0) " & NL
102 & "frndint " & NL -- Integer (X * Log2 (E))
103 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
104 & "fxch " & NL
105 & "f2xm1 " & NL -- 2**(...) - 1
106 & "fld1 " & NL
107 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
108 & "fscale " & NL -- E ** X
109 & "fstp %%st(1) ",
110 Outputs => Double'Asm_Output ("=t", Result),
111 Inputs => Double'Asm_Input ("0", X));
112 return Result;
113 end Exp;
114
115 ------------
116 -- Is_Nan --
117 ------------
118
119 function Is_Nan (X : Double) return Boolean is
120 begin
121 -- The IEEE NaN values are the only ones that do not equal themselves
122
123 return X /= X;
124 end Is_Nan;
125
126 ---------
127 -- Log --
128 ---------
129
130 function Log (X : Double) return Double is
131 Result : Double;
132
133 begin
134 Asm (Template =>
135 "fldln2 " & NL
136 & "fxch " & NL
137 & "fyl2x " & NL,
138 Outputs => Double'Asm_Output ("=t", Result),
139 Inputs => Double'Asm_Input ("0", X));
140 return Result;
141 end Log;
142
143 ------------
144 -- Reduce --
145 ------------
146
147 procedure Reduce (X : in out Double; Q : out Natural) is
148 Half_Pi : constant := Pi / 2.0;
149 Two_Over_Pi : constant := 2.0 / Pi;
150
151 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
152 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
153 P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
154 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
155 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
156 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
157 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
158 - P4, HM);
159 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
160 K : Double;
161 R : Integer;
162
163 begin
164 -- For X < 2.0**HM, all products below are computed exactly.
165 -- Due to cancellation effects all subtractions are exact as well.
166 -- As no double extended floating-point number has more than 75
167 -- zeros after the binary point, the result will be the correctly
168 -- rounded result of X - K * (Pi / 2.0).
169
170 K := X * Two_Over_Pi;
171 while abs K >= 2.0**HM loop
172 K := K * M - (K * M - K);
173 X :=
174 (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
175 K := X * Two_Over_Pi;
176 end loop;
177
178 -- If K is not a number (because X was not finite) raise exception
179
180 if Is_Nan (K) then
181 raise Constraint_Error;
182 end if;
183
184 -- Go through an integer temporary so as to use machine instructions
185
186 R := Integer (Double'Rounding (K));
187 Q := R mod 4;
188 K := Double (R);
189 X := (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
190 end Reduce;
191
192 ----------
193 -- Sqrt --
194 ----------
195
196 function Sqrt (X : Double) return Double is
197 Result : Double;
198
199 begin
200 if X < 0.0 then
201 raise Argument_Error;
202 end if;
203
204 Asm (Template => "fsqrt",
205 Outputs => Double'Asm_Output ("=t", Result),
206 Inputs => Double'Asm_Input ("0", X));
207
208 return Result;
209 end Sqrt;
210
211 --------------------------------
212 -- Other Elementary Functions --
213 --------------------------------
214
215 -- These are built using the previously implemented basic functions
216
217 ----------
218 -- Acos --
219 ----------
220
221 function Acos (X : Double) return Double is
222 Result : Double;
223
224 begin
225 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
226
227 -- The result value is NaN iff input was invalid
228
229 if Is_Nan (Result) then
230 raise Argument_Error;
231 end if;
232
233 return Result;
234 end Acos;
235
236 ----------
237 -- Asin --
238 ----------
239
240 function Asin (X : Double) return Double is
241 Result : Double;
242
243 begin
244 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
245
246 -- The result value is NaN iff input was invalid
247
248 if Is_Nan (Result) then
249 raise Argument_Error;
250 end if;
251
252 return Result;
253 end Asin;
254
255 ---------
256 -- Cos --
257 ---------
258
259 function Cos (X : Double) return Double is
260 Reduced_X : Double := abs X;
261 Result : Double;
262 Quadrant : Natural range 0 .. 3;
263
264 begin
265 if Reduced_X > Pi / 4.0 then
266 Reduce (Reduced_X, Quadrant);
267
268 case Quadrant is
269 when 0 =>
270 Asm (Template => "fcos",
271 Outputs => Double'Asm_Output ("=t", Result),
272 Inputs => Double'Asm_Input ("0", Reduced_X));
273
274 when 1 =>
275 Asm (Template => "fsin",
276 Outputs => Double'Asm_Output ("=t", Result),
277 Inputs => Double'Asm_Input ("0", -Reduced_X));
278
279 when 2 =>
280 Asm (Template => "fcos ; fchs",
281 Outputs => Double'Asm_Output ("=t", Result),
282 Inputs => Double'Asm_Input ("0", Reduced_X));
283
284 when 3 =>
285 Asm (Template => "fsin",
286 Outputs => Double'Asm_Output ("=t", Result),
287 Inputs => Double'Asm_Input ("0", Reduced_X));
288 end case;
289
290 else
291 Asm (Template => "fcos",
292 Outputs => Double'Asm_Output ("=t", Result),
293 Inputs => Double'Asm_Input ("0", Reduced_X));
294 end if;
295
296 return Result;
297 end Cos;
298
299 ---------------------
300 -- Logarithmic_Pow --
301 ---------------------
302
303 function Logarithmic_Pow (X, Y : Double) return Double is
304 Result : Double;
305 begin
306 Asm (Template => "" -- X : Y
307 & "fyl2x " & NL -- Y * Log2 (X)
308 & "fld %%st(0) " & NL -- Y * Log2 (X) : Y * Log2 (X)
309 & "frndint " & NL -- Int (...) : Y * Log2 (X)
310 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...)
311 & "fxch " & NL -- Fract (...) : Int (...)
312 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...)
313 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...)
314 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...)
315 & "fscale ", -- 2**(Fract (...) + Int (...))
316 Outputs => Double'Asm_Output ("=t", Result),
317 Inputs =>
318 (Double'Asm_Input ("0", X),
319 Double'Asm_Input ("u", Y)));
320 return Result;
321 end Logarithmic_Pow;
322
323 ---------
324 -- Pow --
325 ---------
326
327 function Pow (X, Y : Double) return Double is
328 type Mantissa_Type is mod 2**Double'Machine_Mantissa;
329 -- Modular type that can hold all bits of the mantissa of Double
330
331 -- For negative exponents, do divide at the end of the processing
332
333 Negative_Y : constant Boolean := Y < 0.0;
334 Abs_Y : constant Double := abs Y;
335
336 -- During this function the following invariant is kept:
337 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
338
339 Base : Double := X;
340
341 Exp_High : Double := Double'Floor (Abs_Y);
342 Exp_Mid : Double;
343 Exp_Low : Double;
344 Exp_Int : Mantissa_Type;
345
346 Factor : Double := 1.0;
347
348 begin
349 -- Select algorithm for calculating Pow (integer cases fall through)
350
351 if Exp_High >= 2.0**Double'Machine_Mantissa then
352
353 -- In case of Y that is IEEE infinity, just raise constraint error
354
355 if Exp_High > Double'Safe_Last then
356 raise Constraint_Error;
357 end if;
358
359 -- Large values of Y are even integers and will stay integer
360 -- after division by two.
361
362 loop
363 -- Exp_Mid and Exp_Low are zero, so
364 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
365
366 Exp_High := Exp_High / 2.0;
367 Base := Base * Base;
368 exit when Exp_High < 2.0**Double'Machine_Mantissa;
369 end loop;
370
371 elsif Exp_High /= Abs_Y then
372 Exp_Low := Abs_Y - Exp_High;
373 Factor := 1.0;
374
375 if Exp_Low /= 0.0 then
376
377 -- Exp_Low now is in interval (0.0, 1.0)
378 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
379
380 Exp_Mid := 0.0;
381 Exp_Low := Exp_Low - Exp_Mid;
382
383 if Exp_Low >= 0.5 then
384 Factor := Sqrt (X);
385 Exp_Low := Exp_Low - 0.5; -- exact
386
387 if Exp_Low >= 0.25 then
388 Factor := Factor * Sqrt (Factor);
389 Exp_Low := Exp_Low - 0.25; -- exact
390 end if;
391
392 elsif Exp_Low >= 0.25 then
393 Factor := Sqrt (Sqrt (X));
394 Exp_Low := Exp_Low - 0.25; -- exact
395 end if;
396
397 -- Exp_Low now is in interval (0.0, 0.25)
398
399 -- This means it is safe to call Logarithmic_Pow
400 -- for the remaining part.
401
402 Factor := Factor * Logarithmic_Pow (X, Exp_Low);
403 end if;
404
405 elsif X = 0.0 then
406 return 0.0;
407 end if;
408
409 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
410
411 Exp_Int := Mantissa_Type (Exp_High);
412
413 -- Standard way for processing integer powers > 0
414
415 while Exp_Int > 1 loop
416 if (Exp_Int and 1) = 1 then
417
418 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
419
420 Factor := Factor * Base;
421 end if;
422
423 -- Exp_Int is even and Exp_Int > 0, so
424 -- Base**Y = (Base**2)**(Exp_Int / 2)
425
426 Base := Base * Base;
427 Exp_Int := Exp_Int / 2;
428 end loop;
429
430 -- Exp_Int = 1 or Exp_Int = 0
431
432 if Exp_Int = 1 then
433 Factor := Base * Factor;
434 end if;
435
436 if Negative_Y then
437 Factor := 1.0 / Factor;
438 end if;
439
440 return Factor;
441 end Pow;
442
443 ---------
444 -- Sin --
445 ---------
446
447 function Sin (X : Double) return Double is
448 Reduced_X : Double := X;
449 Result : Double;
450 Quadrant : Natural range 0 .. 3;
451
452 begin
453 if abs X > Pi / 4.0 then
454 Reduce (Reduced_X, Quadrant);
455
456 case Quadrant is
457 when 0 =>
458 Asm (Template => "fsin",
459 Outputs => Double'Asm_Output ("=t", Result),
460 Inputs => Double'Asm_Input ("0", Reduced_X));
461
462 when 1 =>
463 Asm (Template => "fcos",
464 Outputs => Double'Asm_Output ("=t", Result),
465 Inputs => Double'Asm_Input ("0", Reduced_X));
466
467 when 2 =>
468 Asm (Template => "fsin",
469 Outputs => Double'Asm_Output ("=t", Result),
470 Inputs => Double'Asm_Input ("0", -Reduced_X));
471
472 when 3 =>
473 Asm (Template => "fcos ; fchs",
474 Outputs => Double'Asm_Output ("=t", Result),
475 Inputs => Double'Asm_Input ("0", Reduced_X));
476 end case;
477
478 else
479 Asm (Template => "fsin",
480 Outputs => Double'Asm_Output ("=t", Result),
481 Inputs => Double'Asm_Input ("0", Reduced_X));
482 end if;
483
484 return Result;
485 end Sin;
486
487 ---------
488 -- Tan --
489 ---------
490
491 function Tan (X : Double) return Double is
492 Reduced_X : Double := X;
493 Result : Double;
494 Quadrant : Natural range 0 .. 3;
495
496 begin
497 if abs X > Pi / 4.0 then
498 Reduce (Reduced_X, Quadrant);
499
500 if Quadrant mod 2 = 0 then
501 Asm (Template => "fptan" & NL
502 & "ffree %%st(0)" & NL
503 & "fincstp",
504 Outputs => Double'Asm_Output ("=t", Result),
505 Inputs => Double'Asm_Input ("0", Reduced_X));
506 else
507 Asm (Template => "fsincos" & NL
508 & "fdivp %%st, %%st(1)" & NL
509 & "fchs",
510 Outputs => Double'Asm_Output ("=t", Result),
511 Inputs => Double'Asm_Input ("0", Reduced_X));
512 end if;
513
514 else
515 Asm (Template =>
516 "fptan " & NL
517 & "ffree %%st(0) " & NL
518 & "fincstp ",
519 Outputs => Double'Asm_Output ("=t", Result),
520 Inputs => Double'Asm_Input ("0", Reduced_X));
521 end if;
522
523 return Result;
524 end Tan;
525
526 ----------
527 -- Sinh --
528 ----------
529
530 function Sinh (X : Double) return Double is
531 begin
532 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
533
534 if abs X < 25.0 then
535 return (Exp (X) - Exp (-X)) / 2.0;
536 else
537 return Exp (X) / 2.0;
538 end if;
539 end Sinh;
540
541 ----------
542 -- Cosh --
543 ----------
544
545 function Cosh (X : Double) return Double is
546 begin
547 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
548
549 if abs X < 22.0 then
550 return (Exp (X) + Exp (-X)) / 2.0;
551 else
552 return Exp (X) / 2.0;
553 end if;
554 end Cosh;
555
556 ----------
557 -- Tanh --
558 ----------
559
560 function Tanh (X : Double) return Double is
561 begin
562 -- Return the Hyperbolic Tangent of x
563
564 -- x -x
565 -- e - e Sinh (X)
566 -- Tanh (X) is defined to be ----------- = --------
567 -- x -x Cosh (X)
568 -- e + e
569
570 if abs X > 23.0 then
571 return Double'Copy_Sign (1.0, X);
572 end if;
573
574 return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
575 end Tanh;
576
577 end Ada.Numerics.Aux;