OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [a-ngelfu.adb] - Blame information for rev 717

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 706 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT RUN-TIME COMPONENTS                         --
4
--                                                                          --
5
--                ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS                 --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--          Copyright (C) 1992-2011, Free Software Foundation, Inc.         --
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. This is not a "strict" implementation, but takes full
35
--  advantage of the C functions, e.g. in providing interface to hardware
36
--  provided versions of the elementary functions.
37
 
38
--  Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, sinh,
39
--  cosh, tanh from C library via math.h
40
 
41
with Ada.Numerics.Aux;
42
 
43
package body Ada.Numerics.Generic_Elementary_Functions 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 X = 0.0 then
515
         return 1.0;
516
 
517
      elsif abs X < Sqrt_Epsilon then
518
         return 1.0;
519
 
520
      end if;
521
 
522
      return Float_Type'Base (Aux.Cos (Double (X)));
523
   end Cos;
524
 
525
   --  Arbitrary cycle
526
 
527
   function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
528
   begin
529
      --  Just reuse the code for Sin. The potential small loss of speed is
530
      --  negligible with proper (front-end) inlining.
531
 
532
      return -Sin (abs X - Cycle * 0.25, Cycle);
533
   end Cos;
534
 
535
   ----------
536
   -- Cosh --
537
   ----------
538
 
539
   function Cosh (X : Float_Type'Base) return Float_Type'Base is
540
      Lnv      : constant Float_Type'Base := 8#0.542714#;
541
      V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
542
      Y        : constant Float_Type'Base := abs X;
543
      Z        : Float_Type'Base;
544
 
545
   begin
546
      if Y < Sqrt_Epsilon then
547
         return 1.0;
548
 
549
      elsif  Y > Log_Inverse_Epsilon then
550
         Z := Exp_Strict (Y - Lnv);
551
         return (Z + V2minus1 * Z);
552
 
553
      else
554
         Z := Exp_Strict (Y);
555
         return 0.5 * (Z + 1.0 / Z);
556
      end if;
557
 
558
   end Cosh;
559
 
560
   ---------
561
   -- Cot --
562
   ---------
563
 
564
   --  Natural cycle
565
 
566
   function Cot (X : Float_Type'Base) return Float_Type'Base is
567
   begin
568
      if X = 0.0 then
569
         raise Constraint_Error;
570
 
571
      elsif abs X < Sqrt_Epsilon then
572
         return 1.0 / X;
573
      end if;
574
 
575
      return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
576
   end Cot;
577
 
578
   --  Arbitrary cycle
579
 
580
   function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
581
      T : Float_Type'Base;
582
 
583
   begin
584
      if Cycle <= 0.0 then
585
         raise Argument_Error;
586
      end if;
587
 
588
      T := Float_Type'Base'Remainder (X, Cycle);
589
 
590
      if T = 0.0 or else abs T = 0.5 * Cycle then
591
         raise Constraint_Error;
592
 
593
      elsif abs T < Sqrt_Epsilon then
594
         return 1.0 / T;
595
 
596
      elsif abs T = 0.25 * Cycle then
597
         return 0.0;
598
 
599
      else
600
         T := T / Cycle * Two_Pi;
601
         return Cos (T) / Sin (T);
602
      end if;
603
   end Cot;
604
 
605
   ----------
606
   -- Coth --
607
   ----------
608
 
609
   function Coth (X : Float_Type'Base) return Float_Type'Base is
610
   begin
611
      if X = 0.0 then
612
         raise Constraint_Error;
613
 
614
      elsif X < Half_Log_Epsilon then
615
         return -1.0;
616
 
617
      elsif X > -Half_Log_Epsilon then
618
         return 1.0;
619
 
620
      elsif abs X < Sqrt_Epsilon then
621
         return 1.0 / X;
622
      end if;
623
 
624
      return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
625
   end Coth;
626
 
627
   ---------
628
   -- Exp --
629
   ---------
630
 
631
   function Exp (X : Float_Type'Base) return Float_Type'Base is
632
      Result : Float_Type'Base;
633
 
634
   begin
635
      if X = 0.0 then
636
         return 1.0;
637
      end if;
638
 
639
      Result := Float_Type'Base (Aux.Exp (Double (X)));
640
 
641
      --  Deal with case of Exp returning IEEE infinity. If Machine_Overflows
642
      --  is False, then we can just leave it as an infinity (and indeed we
643
      --  prefer to do so). But if Machine_Overflows is True, then we have
644
      --  to raise a Constraint_Error exception as required by the RM.
645
 
646
      if Float_Type'Machine_Overflows and then not Result'Valid then
647
         raise Constraint_Error;
648
      end if;
649
 
650
      return Result;
651
   end Exp;
652
 
653
   ----------------
654
   -- Exp_Strict --
655
   ----------------
656
 
657
   function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
658
      G : Float_Type'Base;
659
      Z : Float_Type'Base;
660
 
661
      P0 : constant := 0.25000_00000_00000_00000;
662
      P1 : constant := 0.75753_18015_94227_76666E-2;
663
      P2 : constant := 0.31555_19276_56846_46356E-4;
664
 
665
      Q0 : constant := 0.5;
666
      Q1 : constant := 0.56817_30269_85512_21787E-1;
667
      Q2 : constant := 0.63121_89437_43985_02557E-3;
668
      Q3 : constant := 0.75104_02839_98700_46114E-6;
669
 
670
      C1 : constant := 8#0.543#;
671
      C2 : constant := -2.1219_44400_54690_58277E-4;
672
      Le : constant := 1.4426_95040_88896_34074;
673
 
674
      XN : Float_Type'Base;
675
      P, Q, R : Float_Type'Base;
676
 
677
   begin
678
      if X = 0.0 then
679
         return 1.0;
680
      end if;
681
 
682
      XN := Float_Type'Base'Rounding (X * Le);
683
      G := (X - XN * C1) - XN * C2;
684
      Z := G * G;
685
      P := G * ((P2 * Z + P1) * Z + P0);
686
      Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
687
      R := 0.5 + P / (Q - P);
688
 
689
      R := Float_Type'Base'Scaling (R, Integer (XN) + 1);
690
 
691
      --  Deal with case of Exp returning IEEE infinity. If Machine_Overflows
692
      --  is False, then we can just leave it as an infinity (and indeed we
693
      --  prefer to do so). But if Machine_Overflows is True, then we have to
694
      --  raise a Constraint_Error exception as required by the RM.
695
 
696
      if Float_Type'Machine_Overflows and then not R'Valid then
697
         raise Constraint_Error;
698
      else
699
         return R;
700
      end if;
701
 
702
   end Exp_Strict;
703
 
704
   ----------------
705
   -- Local_Atan --
706
   ----------------
707
 
708
   function Local_Atan
709
     (Y : Float_Type'Base;
710
      X : Float_Type'Base := 1.0) return Float_Type'Base
711
   is
712
      Z        : Float_Type'Base;
713
      Raw_Atan : Float_Type'Base;
714
 
715
   begin
716
      Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));
717
 
718
      Raw_Atan :=
719
        (if Z < Sqrt_Epsilon then Z
720
         elsif Z = 1.0 then Pi / 4.0
721
         else Float_Type'Base (Aux.Atan (Double (Z))));
722
 
723
      if abs Y > abs X then
724
         Raw_Atan := Half_Pi - Raw_Atan;
725
      end if;
726
 
727
      if X > 0.0 then
728
         return Float_Type'Copy_Sign (Raw_Atan, Y);
729
      else
730
         return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
731
      end if;
732
   end Local_Atan;
733
 
734
   ---------
735
   -- Log --
736
   ---------
737
 
738
   --  Natural base
739
 
740
   function Log (X : Float_Type'Base) return Float_Type'Base is
741
   begin
742
      if X < 0.0 then
743
         raise Argument_Error;
744
 
745
      elsif X = 0.0 then
746
         raise Constraint_Error;
747
 
748
      elsif X = 1.0 then
749
         return 0.0;
750
      end if;
751
 
752
      return Float_Type'Base (Aux.Log (Double (X)));
753
   end Log;
754
 
755
   --  Arbitrary base
756
 
757
   function Log (X, Base : Float_Type'Base) return Float_Type'Base is
758
   begin
759
      if X < 0.0 then
760
         raise Argument_Error;
761
 
762
      elsif Base <= 0.0 or else Base = 1.0 then
763
         raise Argument_Error;
764
 
765
      elsif X = 0.0 then
766
         raise Constraint_Error;
767
 
768
      elsif X = 1.0 then
769
         return 0.0;
770
      end if;
771
 
772
      return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
773
   end Log;
774
 
775
   ---------
776
   -- Sin --
777
   ---------
778
 
779
   --  Natural cycle
780
 
781
   function Sin (X : Float_Type'Base) return Float_Type'Base is
782
   begin
783
      if abs X < Sqrt_Epsilon then
784
         return X;
785
      end if;
786
 
787
      return Float_Type'Base (Aux.Sin (Double (X)));
788
   end Sin;
789
 
790
   --  Arbitrary cycle
791
 
792
   function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
793
      T : Float_Type'Base;
794
 
795
   begin
796
      if Cycle <= 0.0 then
797
         raise Argument_Error;
798
 
799
      --  If X is zero, return it as the result, preserving the argument sign.
800
      --  Is this test really needed on any machine ???
801
 
802
      elsif X = 0.0 then
803
         return X;
804
      end if;
805
 
806
      T := Float_Type'Base'Remainder (X, Cycle);
807
 
808
      --  The following two reductions reduce the argument to the interval
809
      --  [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed
810
      --  to prevent inaccuracy that may result if the sine function uses a
811
      --  different (more accurate) value of Pi in its reduction than is used
812
      --  in the multiplication with Two_Pi.
813
 
814
      if abs T > 0.25 * Cycle then
815
         T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
816
      end if;
817
 
818
      --  Could test for 12.0 * abs T = Cycle, and return an exact value in
819
      --  those cases. It is not clear this is worth the extra test though.
820
 
821
      return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
822
   end Sin;
823
 
824
   ----------
825
   -- Sinh --
826
   ----------
827
 
828
   function Sinh (X : Float_Type'Base) return Float_Type'Base is
829
      Lnv      : constant Float_Type'Base := 8#0.542714#;
830
      V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
831
      Y        : constant Float_Type'Base := abs X;
832
      F        : constant Float_Type'Base := Y * Y;
833
      Z        : Float_Type'Base;
834
 
835
      Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;
836
 
837
   begin
838
      if Y < Sqrt_Epsilon then
839
         return X;
840
 
841
      elsif  Y > Log_Inverse_Epsilon then
842
         Z := Exp_Strict (Y - Lnv);
843
         Z := Z + V2minus1 * Z;
844
 
845
      elsif Y < 1.0 then
846
 
847
         if Float_Digits_1_6 then
848
 
849
            --  Use expansion provided by Cody and Waite, p. 226. Note that
850
            --  leading term of the polynomial in Q is exactly 1.0.
851
 
852
            declare
853
               P0 : constant := -0.71379_3159E+1;
854
               P1 : constant := -0.19033_3399E+0;
855
               Q0 : constant := -0.42827_7109E+2;
856
 
857
            begin
858
               Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
859
            end;
860
 
861
         else
862
            declare
863
               P0 : constant := -0.35181_28343_01771_17881E+6;
864
               P1 : constant := -0.11563_52119_68517_68270E+5;
865
               P2 : constant := -0.16375_79820_26307_51372E+3;
866
               P3 : constant := -0.78966_12741_73570_99479E+0;
867
               Q0 : constant := -0.21108_77005_81062_71242E+7;
868
               Q1 : constant :=  0.36162_72310_94218_36460E+5;
869
               Q2 : constant := -0.27773_52311_96507_01667E+3;
870
 
871
            begin
872
               Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
873
                              / (((F + Q2) * F + Q1) * F + Q0);
874
            end;
875
         end if;
876
 
877
      else
878
         Z := Exp_Strict (Y);
879
         Z := 0.5 * (Z - 1.0 / Z);
880
      end if;
881
 
882
      if X > 0.0 then
883
         return Z;
884
      else
885
         return -Z;
886
      end if;
887
   end Sinh;
888
 
889
   ----------
890
   -- Sqrt --
891
   ----------
892
 
893
   function Sqrt (X : Float_Type'Base) return Float_Type'Base is
894
   begin
895
      if X < 0.0 then
896
         raise Argument_Error;
897
 
898
      --  Special case Sqrt (0.0) to preserve possible minus sign per IEEE
899
 
900
      elsif X = 0.0 then
901
         return X;
902
      end if;
903
 
904
      return Float_Type'Base (Aux.Sqrt (Double (X)));
905
   end Sqrt;
906
 
907
   ---------
908
   -- Tan --
909
   ---------
910
 
911
   --  Natural cycle
912
 
913
   function Tan (X : Float_Type'Base) return Float_Type'Base is
914
   begin
915
      if abs X < Sqrt_Epsilon then
916
         return X;
917
      end if;
918
 
919
      --  Note: if X is exactly pi/2, then we should raise an exception, since
920
      --  the result would overflow. But for all floating-point formats we deal
921
      --  with, it is impossible for X to be exactly pi/2, and the result is
922
      --  always in range.
923
 
924
      return Float_Type'Base (Aux.Tan (Double (X)));
925
   end Tan;
926
 
927
   --  Arbitrary cycle
928
 
929
   function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
930
      T : Float_Type'Base;
931
 
932
   begin
933
      if Cycle <= 0.0 then
934
         raise Argument_Error;
935
 
936
      elsif X = 0.0 then
937
         return X;
938
      end if;
939
 
940
      T := Float_Type'Base'Remainder (X, Cycle);
941
 
942
      if abs T = 0.25 * Cycle then
943
         raise Constraint_Error;
944
 
945
      elsif abs T = 0.5 * Cycle then
946
         return 0.0;
947
 
948
      else
949
         T := T / Cycle * Two_Pi;
950
         return Sin (T) / Cos (T);
951
      end if;
952
 
953
   end Tan;
954
 
955
   ----------
956
   -- Tanh --
957
   ----------
958
 
959
   function Tanh (X : Float_Type'Base) return Float_Type'Base is
960
      P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
961
      P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
962
      P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;
963
 
964
      Q0 : constant Float_Type'Base :=  0.48402_35707_19886_88686E+4;
965
      Q1 : constant Float_Type'Base :=  0.22337_72071_89623_12926E+4;
966
      Q2 : constant Float_Type'Base :=  0.11274_47438_05349_49335E+3;
967
      Q3 : constant Float_Type'Base :=  0.10000_00000_00000_00000E+1;
968
 
969
      Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
970
 
971
      P, Q, R : Float_Type'Base;
972
      Y : constant Float_Type'Base := abs X;
973
      G : constant Float_Type'Base := Y * Y;
974
 
975
      Float_Type_Digits_15_Or_More : constant Boolean :=
976
                                       Float_Type'Digits > 14;
977
 
978
   begin
979
      if X < Half_Log_Epsilon then
980
         return -1.0;
981
 
982
      elsif X > -Half_Log_Epsilon then
983
         return 1.0;
984
 
985
      elsif Y < Sqrt_Epsilon then
986
         return X;
987
 
988
      elsif Y < Half_Ln3
989
        and then Float_Type_Digits_15_Or_More
990
      then
991
         P := (P2 * G + P1) * G + P0;
992
         Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
993
         R := G * (P / Q);
994
         return X + X * R;
995
 
996
      else
997
         return Float_Type'Base (Aux.Tanh (Double (X)));
998
      end if;
999
   end Tanh;
1000
 
1001
end Ada.Numerics.Generic_Elementary_Functions;

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.