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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [math_lib.adb] - Blame information for rev 750

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
--                             M A T H _ L I B                              --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--          Copyright (C) 1992-2009, 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
--  A known weakness is that on the x86, all computation is done in Double,
39
--  which means that a lot of accuracy is lost for the Long_Long_Float case.
40
 
41
--  Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
42
--  sinh, cosh, tanh from C library via math.h
43
 
44
--  This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
45
--  provides a compatible body for the DEC Math_Lib package.
46
 
47
with Ada.Numerics.Aux;
48
use type Ada.Numerics.Aux.Double;
49
with Ada.Numerics; use Ada.Numerics;
50
 
51
package body Math_Lib is
52
 
53
   Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
54
 
55
   Two_Pi     : constant Real'Base := 2.0 * Pi;
56
   Half_Pi    : constant Real'Base := Pi / 2.0;
57
   Fourth_Pi  : constant Real'Base := Pi / 4.0;
58
   Epsilon    : constant Real'Base := Real'Base'Epsilon;
59
   IEpsilon   : constant Real'Base := 1.0 / Epsilon;
60
 
61
   subtype Double is Aux.Double;
62
 
63
   DEpsilon    : constant Double := Double (Epsilon);
64
   DIEpsilon   : constant Double := Double (IEpsilon);
65
 
66
   -----------------------
67
   -- Local Subprograms --
68
   -----------------------
69
 
70
   function Arctan
71
     (Y    : Real;
72
      A    : Real := 1.0)
73
      return Real;
74
 
75
   function Arctan
76
     (Y     : Real;
77
      A     : Real := 1.0;
78
      Cycle : Real)
79
      return  Real;
80
 
81
   function Exact_Remainder
82
     (A    : Real;
83
      Y    : Real)
84
      return Real;
85
   --  Computes exact remainder of A divided by Y
86
 
87
   function Half_Log_Epsilon return Real;
88
   --  Function to provide constant: 0.5 * Log (Epsilon)
89
 
90
   function Local_Atan
91
     (Y    : Real;
92
      A    : Real := 1.0)
93
      return Real;
94
   --  Common code for arc tangent after cycle reduction
95
 
96
   function Log_Inverse_Epsilon return Real;
97
   --  Function to provide constant: Log (1.0 / Epsilon)
98
 
99
   function Square_Root_Epsilon return Real;
100
   --  Function to provide constant: Sqrt (Epsilon)
101
 
102
   ----------
103
   -- "**" --
104
   ----------
105
 
106
   function "**" (A1, A2 : Real) return Real is
107
 
108
   begin
109
      if A1 = 0.0
110
        and then A2 = 0.0
111
      then
112
         raise Argument_Error;
113
 
114
      elsif A1 < 0.0 then
115
         raise Argument_Error;
116
 
117
      elsif A2 = 0.0 then
118
         return 1.0;
119
 
120
      elsif A1 = 0.0 then
121
         if A2 < 0.0 then
122
            raise Constraint_Error;
123
         else
124
            return 0.0;
125
         end if;
126
 
127
      elsif A1 = 1.0 then
128
         return 1.0;
129
 
130
      elsif A2 = 1.0 then
131
         return A1;
132
 
133
      else
134
         begin
135
            if A2 = 2.0 then
136
               return A1 * A1;
137
            else
138
               return
139
                 Real (Aux.pow (Double (A1), Double (A2)));
140
            end if;
141
 
142
         exception
143
            when others =>
144
               raise Constraint_Error;
145
         end;
146
      end if;
147
   end "**";
148
 
149
   ------------
150
   -- Arccos --
151
   ------------
152
 
153
   --  Natural cycle
154
 
155
   function Arccos (A : Real) return Real is
156
      Temp : Real'Base;
157
 
158
   begin
159
      if abs A > 1.0 then
160
         raise Argument_Error;
161
 
162
      elsif abs A < Square_Root_Epsilon then
163
         return Pi / 2.0 - A;
164
 
165
      elsif A = 1.0 then
166
         return 0.0;
167
 
168
      elsif A = -1.0 then
169
         return Pi;
170
      end if;
171
 
172
      Temp := Real (Aux.acos (Double (A)));
173
 
174
      if Temp < 0.0 then
175
         Temp := Pi + Temp;
176
      end if;
177
 
178
      return Temp;
179
   end Arccos;
180
 
181
   --  Arbitrary cycle
182
 
183
   function Arccos (A, Cycle : Real) return Real is
184
      Temp : Real'Base;
185
 
186
   begin
187
      if Cycle <= 0.0 then
188
         raise Argument_Error;
189
 
190
      elsif abs A > 1.0 then
191
         raise Argument_Error;
192
 
193
      elsif abs A < Square_Root_Epsilon then
194
         return Cycle / 4.0;
195
 
196
      elsif A = 1.0 then
197
         return 0.0;
198
 
199
      elsif A = -1.0 then
200
         return Cycle / 2.0;
201
      end if;
202
 
203
      Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
204
 
205
      if Temp < 0.0 then
206
         Temp := Cycle / 2.0 + Temp;
207
      end if;
208
 
209
      return Temp;
210
   end Arccos;
211
 
212
   -------------
213
   -- Arccosh --
214
   -------------
215
 
216
   function Arccosh (A : Real) return Real is
217
   begin
218
      --  Return Log (A - Sqrt (A * A - 1.0));  double valued,
219
      --    only positive value returned
220
      --  What is this comment ???
221
 
222
      if A < 1.0 then
223
         raise Argument_Error;
224
 
225
      elsif A < 1.0 + Square_Root_Epsilon then
226
         return A - 1.0;
227
 
228
      elsif abs A > 1.0 / Square_Root_Epsilon then
229
         return Log (A) + Log_Two;
230
 
231
      else
232
         return Log (A + Sqrt (A * A - 1.0));
233
      end if;
234
   end Arccosh;
235
 
236
   ------------
237
   -- Arccot --
238
   ------------
239
 
240
   --  Natural cycle
241
 
242
   function Arccot
243
     (A    : Real;
244
      Y    : Real := 1.0)
245
      return Real
246
   is
247
   begin
248
      --  Just reverse arguments
249
 
250
      return Arctan (Y, A);
251
   end Arccot;
252
 
253
   --  Arbitrary cycle
254
 
255
   function Arccot
256
     (A     : Real;
257
      Y     : Real := 1.0;
258
      Cycle : Real)
259
      return  Real
260
   is
261
   begin
262
      --  Just reverse arguments
263
 
264
      return Arctan (Y, A, Cycle);
265
   end Arccot;
266
 
267
   -------------
268
   -- Arccoth --
269
   -------------
270
 
271
   function Arccoth (A : Real) return Real is
272
   begin
273
      if abs A = 1.0 then
274
         raise Constraint_Error;
275
 
276
      elsif abs A < 1.0 then
277
         raise Argument_Error;
278
 
279
      elsif abs A > 1.0 / Epsilon then
280
         return 0.0;
281
 
282
      else
283
         return 0.5 * Log ((1.0 + A) / (A - 1.0));
284
      end if;
285
   end Arccoth;
286
 
287
   ------------
288
   -- Arcsin --
289
   ------------
290
 
291
   --  Natural cycle
292
 
293
   function Arcsin (A : Real) return Real is
294
   begin
295
      if abs A > 1.0 then
296
         raise Argument_Error;
297
 
298
      elsif abs A < Square_Root_Epsilon then
299
         return A;
300
 
301
      elsif A = 1.0 then
302
         return Pi / 2.0;
303
 
304
      elsif A = -1.0 then
305
         return -Pi / 2.0;
306
      end if;
307
 
308
      return Real (Aux.asin (Double (A)));
309
   end Arcsin;
310
 
311
   --  Arbitrary cycle
312
 
313
   function Arcsin (A, Cycle : Real) return Real is
314
   begin
315
      if Cycle <= 0.0 then
316
         raise Argument_Error;
317
 
318
      elsif abs A > 1.0 then
319
         raise Argument_Error;
320
 
321
      elsif A = 0.0 then
322
         return A;
323
 
324
      elsif A = 1.0 then
325
         return Cycle / 4.0;
326
 
327
      elsif A = -1.0 then
328
         return -Cycle / 4.0;
329
      end if;
330
 
331
      return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
332
   end Arcsin;
333
 
334
   -------------
335
   -- Arcsinh --
336
   -------------
337
 
338
   function Arcsinh (A : Real) return Real is
339
   begin
340
      if abs A < Square_Root_Epsilon then
341
         return A;
342
 
343
      elsif A > 1.0 / Square_Root_Epsilon then
344
         return Log (A) + Log_Two;
345
 
346
      elsif A < -1.0 / Square_Root_Epsilon then
347
         return -(Log (-A) + Log_Two);
348
 
349
      elsif A < 0.0 then
350
         return -Log (abs A + Sqrt (A * A + 1.0));
351
 
352
      else
353
         return Log (A + Sqrt (A * A + 1.0));
354
      end if;
355
   end Arcsinh;
356
 
357
   ------------
358
   -- Arctan --
359
   ------------
360
 
361
   --  Natural cycle
362
 
363
   function Arctan
364
     (Y    : Real;
365
      A    : Real := 1.0)
366
      return Real
367
   is
368
   begin
369
      if A = 0.0
370
        and then Y = 0.0
371
      then
372
         raise Argument_Error;
373
 
374
      elsif Y = 0.0 then
375
         if A > 0.0 then
376
            return 0.0;
377
         else -- A < 0.0
378
            return Pi;
379
         end if;
380
 
381
      elsif A = 0.0 then
382
         if Y > 0.0 then
383
            return Half_Pi;
384
         else -- Y < 0.0
385
            return -Half_Pi;
386
         end if;
387
 
388
      else
389
         return Local_Atan (Y, A);
390
      end if;
391
   end Arctan;
392
 
393
   --  Arbitrary cycle
394
 
395
   function Arctan
396
     (Y     : Real;
397
      A     : Real := 1.0;
398
      Cycle : Real)
399
      return  Real
400
   is
401
   begin
402
      if Cycle <= 0.0 then
403
         raise Argument_Error;
404
 
405
      elsif A = 0.0
406
        and then Y = 0.0
407
      then
408
         raise Argument_Error;
409
 
410
      elsif Y = 0.0 then
411
         if A > 0.0 then
412
            return 0.0;
413
         else -- A < 0.0
414
            return Cycle / 2.0;
415
         end if;
416
 
417
      elsif A = 0.0 then
418
         if Y > 0.0 then
419
            return Cycle / 4.0;
420
         else -- Y < 0.0
421
            return -Cycle / 4.0;
422
         end if;
423
 
424
      else
425
         return Local_Atan (Y, A) *  Cycle / Two_Pi;
426
      end if;
427
   end Arctan;
428
 
429
   -------------
430
   -- Arctanh --
431
   -------------
432
 
433
   function Arctanh (A : Real) return Real is
434
   begin
435
      if abs A = 1.0 then
436
         raise Constraint_Error;
437
 
438
      elsif abs A > 1.0 then
439
         raise Argument_Error;
440
 
441
      elsif abs A < Square_Root_Epsilon then
442
         return A;
443
 
444
      else
445
         return 0.5 * Log ((1.0 + A) / (1.0 - A));
446
      end if;
447
   end Arctanh;
448
 
449
   ---------
450
   -- Cos --
451
   ---------
452
 
453
   --  Natural cycle
454
 
455
   function Cos (A : Real) return Real is
456
   begin
457
      if A = 0.0 then
458
         return 1.0;
459
 
460
      elsif abs A < Square_Root_Epsilon then
461
         return 1.0;
462
 
463
      end if;
464
 
465
      return Real (Aux.Cos (Double (A)));
466
   end Cos;
467
 
468
   --  Arbitrary cycle
469
 
470
   function Cos (A, Cycle : Real) return Real is
471
      T : Real'Base;
472
 
473
   begin
474
      if Cycle <= 0.0 then
475
         raise Argument_Error;
476
 
477
      elsif A = 0.0 then
478
         return 1.0;
479
      end if;
480
 
481
      T := Exact_Remainder (abs (A), Cycle) / Cycle;
482
 
483
      if T = 0.25
484
        or else T = 0.75
485
        or else T = -0.25
486
        or else T = -0.75
487
      then
488
         return 0.0;
489
 
490
      elsif T = 0.5 or T = -0.5 then
491
         return -1.0;
492
      end if;
493
 
494
      return Real (Aux.Cos (Double (T * Two_Pi)));
495
   end Cos;
496
 
497
   ----------
498
   -- Cosh --
499
   ----------
500
 
501
   function Cosh (A : Real) return Real is
502
   begin
503
      if abs A < Square_Root_Epsilon then
504
         return 1.0;
505
 
506
      elsif abs A > Log_Inverse_Epsilon then
507
         return Exp ((abs A) - Log_Two);
508
      end if;
509
 
510
      return Real (Aux.cosh (Double (A)));
511
 
512
   exception
513
      when others =>
514
         raise Constraint_Error;
515
   end Cosh;
516
 
517
   ---------
518
   -- Cot --
519
   ---------
520
 
521
   --  Natural cycle
522
 
523
   function Cot (A : Real) return Real is
524
   begin
525
      if A = 0.0 then
526
         raise Constraint_Error;
527
 
528
      elsif abs A < Square_Root_Epsilon then
529
         return 1.0 / A;
530
      end if;
531
 
532
      return Real (1.0 / Real'Base (Aux.tan (Double (A))));
533
   end Cot;
534
 
535
   --  Arbitrary cycle
536
 
537
   function Cot (A, Cycle : Real) return Real is
538
      T : Real'Base;
539
 
540
   begin
541
      if Cycle <= 0.0 then
542
         raise Argument_Error;
543
 
544
      elsif A = 0.0 then
545
         raise Constraint_Error;
546
 
547
      elsif abs A < Square_Root_Epsilon then
548
         return 1.0 / A;
549
      end if;
550
 
551
      T := Exact_Remainder (A, Cycle) / Cycle;
552
 
553
      if T = 0.0 or T = 0.5 or T = -0.5 then
554
         raise Constraint_Error;
555
      else
556
         return  Cos (T * Two_Pi) / Sin (T * Two_Pi);
557
      end if;
558
   end Cot;
559
 
560
   ----------
561
   -- Coth --
562
   ----------
563
 
564
   function Coth (A : Real) return Real is
565
   begin
566
      if A = 0.0 then
567
         raise Constraint_Error;
568
 
569
      elsif A < Half_Log_Epsilon then
570
         return -1.0;
571
 
572
      elsif A > -Half_Log_Epsilon then
573
         return 1.0;
574
 
575
      elsif abs A < Square_Root_Epsilon then
576
         return 1.0 / A;
577
      end if;
578
 
579
      return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
580
   end Coth;
581
 
582
   ---------------------
583
   -- Exact_Remainder --
584
   ---------------------
585
 
586
   function Exact_Remainder
587
     (A    : Real;
588
      Y    : Real)
589
      return Real
590
   is
591
      Denominator : Real'Base := abs A;
592
      Divisor     : Real'Base := abs Y;
593
      Reducer     : Real'Base;
594
      Sign        : Real'Base := 1.0;
595
 
596
   begin
597
      if Y = 0.0 then
598
         raise Constraint_Error;
599
 
600
      elsif A = 0.0 then
601
         return 0.0;
602
 
603
      elsif A = Y then
604
         return 0.0;
605
 
606
      elsif Denominator < Divisor then
607
         return A;
608
      end if;
609
 
610
      while Denominator >= Divisor loop
611
 
612
         --  Put divisors mantissa with denominators exponent to make reducer
613
 
614
         Reducer := Divisor;
615
 
616
         begin
617
            while Reducer * 1_048_576.0 < Denominator loop
618
               Reducer := Reducer * 1_048_576.0;
619
            end loop;
620
 
621
         exception
622
            when others => null;
623
         end;
624
 
625
         begin
626
            while Reducer * 1_024.0 < Denominator loop
627
               Reducer := Reducer * 1_024.0;
628
            end loop;
629
 
630
         exception
631
            when others => null;
632
         end;
633
 
634
         begin
635
            while Reducer * 2.0 < Denominator loop
636
               Reducer := Reducer * 2.0;
637
            end loop;
638
 
639
         exception
640
            when others => null;
641
         end;
642
 
643
         Denominator := Denominator - Reducer;
644
      end loop;
645
 
646
      if A < 0.0 then
647
         return -Denominator;
648
      else
649
         return Denominator;
650
      end if;
651
   end Exact_Remainder;
652
 
653
   ---------
654
   -- Exp --
655
   ---------
656
 
657
   function Exp (A : Real) return Real is
658
      Result : Real'Base;
659
 
660
   begin
661
      if A = 0.0 then
662
         return 1.0;
663
 
664
      else
665
         Result := Real (Aux.Exp (Double (A)));
666
 
667
         --  The check here catches the case of Exp returning IEEE infinity
668
 
669
         if Result > Real'Last then
670
            raise Constraint_Error;
671
         else
672
            return Result;
673
         end if;
674
      end if;
675
   end Exp;
676
 
677
   ----------------------
678
   -- Half_Log_Epsilon --
679
   ----------------------
680
 
681
   --  Cannot precompute this constant, because this is required to be a
682
   --  pure package, which allows no state. A pity, but no way around it!
683
 
684
   function Half_Log_Epsilon return Real is
685
   begin
686
      return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
687
   end Half_Log_Epsilon;
688
 
689
   ----------------
690
   -- Local_Atan --
691
   ----------------
692
 
693
   function Local_Atan
694
     (Y    : Real;
695
      A    : Real := 1.0)
696
      return Real
697
   is
698
      Z        : Real'Base;
699
      Raw_Atan : Real'Base;
700
 
701
   begin
702
      if abs Y > abs A then
703
         Z := abs (A / Y);
704
      else
705
         Z := abs (Y / A);
706
      end if;
707
 
708
      if Z < Square_Root_Epsilon then
709
         Raw_Atan := Z;
710
 
711
      elsif Z = 1.0 then
712
         Raw_Atan := Pi / 4.0;
713
 
714
      elsif Z < Square_Root_Epsilon then
715
         Raw_Atan := Z;
716
 
717
      else
718
         Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
719
      end if;
720
 
721
      if abs Y > abs A then
722
         Raw_Atan := Half_Pi - Raw_Atan;
723
      end if;
724
 
725
      if A > 0.0 then
726
         if Y > 0.0 then
727
            return Raw_Atan;
728
         else                 --  Y < 0.0
729
            return -Raw_Atan;
730
         end if;
731
 
732
      else                    --  A < 0.0
733
         if Y > 0.0 then
734
            return Pi - Raw_Atan;
735
         else                  --  Y < 0.0
736
            return -(Pi - Raw_Atan);
737
         end if;
738
      end if;
739
   end Local_Atan;
740
 
741
   ---------
742
   -- Log --
743
   ---------
744
 
745
   --  Natural base
746
 
747
   function Log (A : Real) return Real is
748
   begin
749
      if A < 0.0 then
750
         raise Argument_Error;
751
 
752
      elsif A = 0.0 then
753
         raise Constraint_Error;
754
 
755
      elsif A = 1.0 then
756
         return 0.0;
757
      end if;
758
 
759
      return Real (Aux.Log (Double (A)));
760
   end Log;
761
 
762
   --  Arbitrary base
763
 
764
   function Log (A, Base : Real) return Real is
765
   begin
766
      if A < 0.0 then
767
         raise Argument_Error;
768
 
769
      elsif Base <= 0.0 or else Base = 1.0 then
770
         raise Argument_Error;
771
 
772
      elsif A = 0.0 then
773
         raise Constraint_Error;
774
 
775
      elsif A = 1.0 then
776
         return 0.0;
777
      end if;
778
 
779
      return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
780
   end Log;
781
 
782
   -------------------------
783
   -- Log_Inverse_Epsilon --
784
   -------------------------
785
 
786
   --  Cannot precompute this constant, because this is required to be a
787
   --  pure package, which allows no state. A pity, but no way around it!
788
 
789
   function Log_Inverse_Epsilon return Real is
790
   begin
791
      return Real (Aux.Log (DIEpsilon));
792
   end Log_Inverse_Epsilon;
793
 
794
   ---------
795
   -- Sin --
796
   ---------
797
 
798
   --  Natural cycle
799
 
800
   function Sin (A : Real) return Real is
801
   begin
802
      if abs A < Square_Root_Epsilon then
803
         return A;
804
      end if;
805
 
806
      return Real (Aux.Sin (Double (A)));
807
   end Sin;
808
 
809
   --  Arbitrary cycle
810
 
811
   function Sin (A, Cycle : Real) return Real is
812
      T : Real'Base;
813
 
814
   begin
815
      if Cycle <= 0.0 then
816
         raise Argument_Error;
817
 
818
      elsif A = 0.0 then
819
         return A;
820
      end if;
821
 
822
      T := Exact_Remainder (A, Cycle) / Cycle;
823
 
824
      if T = 0.0 or T = 0.5 or T = -0.5 then
825
         return 0.0;
826
 
827
      elsif T = 0.25 or T = -0.75 then
828
         return 1.0;
829
 
830
      elsif T = -0.25 or T = 0.75 then
831
         return -1.0;
832
 
833
      end if;
834
 
835
      return Real (Aux.Sin (Double (T * Two_Pi)));
836
   end Sin;
837
 
838
   ----------
839
   -- Sinh --
840
   ----------
841
 
842
   function Sinh (A : Real) return Real is
843
   begin
844
      if abs A < Square_Root_Epsilon then
845
         return A;
846
 
847
      elsif  A > Log_Inverse_Epsilon then
848
         return Exp (A - Log_Two);
849
 
850
      elsif A < -Log_Inverse_Epsilon then
851
         return -Exp ((-A) - Log_Two);
852
      end if;
853
 
854
      return Real (Aux.Sinh (Double (A)));
855
 
856
   exception
857
      when others =>
858
         raise Constraint_Error;
859
   end Sinh;
860
 
861
   -------------------------
862
   -- Square_Root_Epsilon --
863
   -------------------------
864
 
865
   --  Cannot precompute this constant, because this is required to be a
866
   --  pure package, which allows no state. A pity, but no way around it!
867
 
868
   function Square_Root_Epsilon return Real is
869
   begin
870
      return Real (Aux.Sqrt (DEpsilon));
871
   end Square_Root_Epsilon;
872
 
873
   ----------
874
   -- Sqrt --
875
   ----------
876
 
877
   function Sqrt (A : Real) return Real is
878
   begin
879
      if A < 0.0 then
880
         raise Argument_Error;
881
 
882
      --  Special case Sqrt (0.0) to preserve possible minus sign per IEEE
883
 
884
      elsif A = 0.0 then
885
         return A;
886
 
887
      --  Sqrt (1.0) must be exact for good complex accuracy
888
 
889
      elsif A = 1.0 then
890
         return 1.0;
891
 
892
      end if;
893
 
894
      return Real (Aux.Sqrt (Double (A)));
895
   end Sqrt;
896
 
897
   ---------
898
   -- Tan --
899
   ---------
900
 
901
   --  Natural cycle
902
 
903
   function Tan (A : Real) return Real is
904
   begin
905
      if abs A < Square_Root_Epsilon then
906
         return A;
907
 
908
      elsif abs A = Pi / 2.0 then
909
         raise Constraint_Error;
910
      end if;
911
 
912
      return Real (Aux.tan (Double (A)));
913
   end Tan;
914
 
915
   --  Arbitrary cycle
916
 
917
   function Tan (A, Cycle : Real) return Real is
918
      T : Real'Base;
919
 
920
   begin
921
      if Cycle <= 0.0 then
922
         raise Argument_Error;
923
 
924
      elsif A = 0.0 then
925
         return A;
926
      end if;
927
 
928
      T := Exact_Remainder (A, Cycle) / Cycle;
929
 
930
      if T = 0.25
931
        or else T = 0.75
932
        or else T = -0.25
933
        or else T = -0.75
934
      then
935
         raise Constraint_Error;
936
 
937
      else
938
         return  Sin (T * Two_Pi) / Cos (T * Two_Pi);
939
      end if;
940
   end Tan;
941
 
942
   ----------
943
   -- Tanh --
944
   ----------
945
 
946
   function Tanh (A : Real) return Real is
947
   begin
948
      if A < Half_Log_Epsilon then
949
         return -1.0;
950
 
951
      elsif A > -Half_Log_Epsilon then
952
         return 1.0;
953
 
954
      elsif abs A < Square_Root_Epsilon then
955
         return A;
956
      end if;
957
 
958
      return Real (Aux.tanh (Double (A)));
959
   end Tanh;
960
 
961
   ----------------------------
962
   -- DEC-Specific functions --
963
   ----------------------------
964
 
965
   function LOG10  (A : REAL) return REAL is
966
   begin
967
      return Log (A, 10.0);
968
   end LOG10;
969
 
970
   function LOG2   (A : REAL) return REAL is
971
   begin
972
      return Log (A, 2.0);
973
   end LOG2;
974
 
975
   function ASIN (A : REAL) return REAL renames Arcsin;
976
   function ACOS (A : REAL) return REAL renames Arccos;
977
 
978
   function ATAN (A : REAL) return REAL is
979
   begin
980
      return Arctan (A, 1.0);
981
   end ATAN;
982
 
983
   function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
984
 
985
   function SIND   (A : REAL) return REAL is
986
   begin
987
      return Sin (A, 360.0);
988
   end SIND;
989
 
990
   function COSD   (A : REAL) return REAL is
991
   begin
992
      return  Cos (A, 360.0);
993
   end COSD;
994
 
995
   function TAND   (A : REAL) return REAL is
996
   begin
997
      return  Tan (A, 360.0);
998
   end TAND;
999
 
1000
   function ASIND  (A : REAL) return REAL is
1001
   begin
1002
      return  Arcsin (A, 360.0);
1003
   end ASIND;
1004
 
1005
   function ACOSD  (A : REAL) return REAL is
1006
   begin
1007
      return  Arccos (A, 360.0);
1008
   end ACOSD;
1009
 
1010
   function Arctan  (A : REAL) return REAL is
1011
   begin
1012
      return  Arctan (A, 1.0, 360.0);
1013
   end Arctan;
1014
 
1015
   function ATAND (A : REAL) return REAL is
1016
   begin
1017
      return Arctan (A, 1.0, 360.0);
1018
   end ATAND;
1019
 
1020
   function ATAN2D (A1, A2 : REAL) return REAL is
1021
   begin
1022
      return Arctan (A1, A2, 360.0);
1023
   end ATAN2D;
1024
 
1025
end Math_Lib;

powered by: WebSVN 2.1.0

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