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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [gcc-4.5.1/] [gcc/] [ada/] [a-ngcoty.adb] - Blame information for rev 414

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

Line No. Rev Author Line
1 281 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT RUN-TIME COMPONENTS                         --
4
--                                                                          --
5
--   A D A . N U M E R I C S . G E N E R I C _ C O M P L E X _ T Y P E S    --
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
with Ada.Numerics.Aux; use Ada.Numerics.Aux;
33
 
34
package body Ada.Numerics.Generic_Complex_Types is
35
 
36
   subtype R is Real'Base;
37
 
38
   Two_Pi  : constant R := R (2.0) * Pi;
39
   Half_Pi : constant R := Pi / R (2.0);
40
 
41
   ---------
42
   -- "*" --
43
   ---------
44
 
45
   function "*" (Left, Right : Complex) return Complex is
46
      X : R;
47
      Y : R;
48
 
49
   begin
50
      X := Left.Re * Right.Re - Left.Im * Right.Im;
51
      Y := Left.Re * Right.Im + Left.Im * Right.Re;
52
 
53
      --  If either component overflows, try to scale (skip in fast math mode)
54
 
55
      if not Standard'Fast_Math then
56
         if abs (X) > R'Last then
57
            X := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Re / 2.0)
58
                            - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
59
         end if;
60
 
61
         if abs (Y) > R'Last then
62
            Y := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Im / 2.0)
63
                            - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
64
         end if;
65
      end if;
66
 
67
      return (X, Y);
68
   end "*";
69
 
70
   function "*" (Left, Right : Imaginary) return Real'Base is
71
   begin
72
      return -(R (Left) * R (Right));
73
   end "*";
74
 
75
   function "*" (Left : Complex; Right : Real'Base) return Complex is
76
   begin
77
      return Complex'(Left.Re * Right, Left.Im * Right);
78
   end "*";
79
 
80
   function "*" (Left : Real'Base; Right : Complex) return Complex is
81
   begin
82
      return (Left * Right.Re, Left * Right.Im);
83
   end "*";
84
 
85
   function "*" (Left : Complex; Right : Imaginary) return Complex is
86
   begin
87
      return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
88
   end "*";
89
 
90
   function "*" (Left : Imaginary; Right : Complex) return Complex is
91
   begin
92
      return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
93
   end "*";
94
 
95
   function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
96
   begin
97
      return Left * Imaginary (Right);
98
   end "*";
99
 
100
   function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
101
   begin
102
      return Imaginary (Left * R (Right));
103
   end "*";
104
 
105
   ----------
106
   -- "**" --
107
   ----------
108
 
109
   function "**" (Left : Complex; Right : Integer) return Complex is
110
      Result : Complex := (1.0, 0.0);
111
      Factor : Complex := Left;
112
      Exp    : Integer := Right;
113
 
114
   begin
115
      --  We use the standard logarithmic approach, Exp gets shifted right
116
      --  testing successive low order bits and Factor is the value of the
117
      --  base raised to the next power of 2. For positive exponents we
118
      --  multiply the result by this factor, for negative exponents, we
119
      --  divide by this factor.
120
 
121
      if Exp >= 0 then
122
 
123
         --  For a positive exponent, if we get a constraint error during
124
         --  this loop, it is an overflow, and the constraint error will
125
         --  simply be passed on to the caller.
126
 
127
         while Exp /= 0 loop
128
            if Exp rem 2 /= 0 then
129
               Result := Result * Factor;
130
            end if;
131
 
132
            Factor := Factor * Factor;
133
            Exp := Exp / 2;
134
         end loop;
135
 
136
         return Result;
137
 
138
      else -- Exp < 0 then
139
 
140
         --  For the negative exponent case, a constraint error during this
141
         --  calculation happens if Factor gets too large, and the proper
142
         --  response is to return 0.0, since what we essentially have is
143
         --  1.0 / infinity, and the closest model number will be zero.
144
 
145
         begin
146
            while Exp /= 0 loop
147
               if Exp rem 2 /= 0 then
148
                  Result := Result * Factor;
149
               end if;
150
 
151
               Factor := Factor * Factor;
152
               Exp := Exp / 2;
153
            end loop;
154
 
155
            return R'(1.0) / Result;
156
 
157
         exception
158
            when Constraint_Error =>
159
               return (0.0, 0.0);
160
         end;
161
      end if;
162
   end "**";
163
 
164
   function "**" (Left : Imaginary; Right : Integer) return Complex is
165
      M : constant R := R (Left) ** Right;
166
   begin
167
      case Right mod 4 is
168
         when 0 => return (M,   0.0);
169
         when 1 => return (0.0, M);
170
         when 2 => return (-M,  0.0);
171
         when 3 => return (0.0, -M);
172
         when others => raise Program_Error;
173
      end case;
174
   end "**";
175
 
176
   ---------
177
   -- "+" --
178
   ---------
179
 
180
   function "+" (Right : Complex) return Complex is
181
   begin
182
      return Right;
183
   end "+";
184
 
185
   function "+" (Left, Right : Complex) return Complex is
186
   begin
187
      return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
188
   end "+";
189
 
190
   function "+" (Right : Imaginary) return Imaginary is
191
   begin
192
      return Right;
193
   end "+";
194
 
195
   function "+" (Left, Right : Imaginary) return Imaginary is
196
   begin
197
      return Imaginary (R (Left) + R (Right));
198
   end "+";
199
 
200
   function "+" (Left : Complex; Right : Real'Base) return Complex is
201
   begin
202
      return Complex'(Left.Re + Right, Left.Im);
203
   end "+";
204
 
205
   function "+" (Left : Real'Base; Right : Complex) return Complex is
206
   begin
207
      return Complex'(Left + Right.Re, Right.Im);
208
   end "+";
209
 
210
   function "+" (Left : Complex; Right : Imaginary) return Complex is
211
   begin
212
      return Complex'(Left.Re, Left.Im + R (Right));
213
   end "+";
214
 
215
   function "+" (Left : Imaginary; Right : Complex) return Complex is
216
   begin
217
      return Complex'(Right.Re, R (Left) + Right.Im);
218
   end "+";
219
 
220
   function "+" (Left : Imaginary; Right : Real'Base) return Complex is
221
   begin
222
      return Complex'(Right, R (Left));
223
   end "+";
224
 
225
   function "+" (Left : Real'Base; Right : Imaginary) return Complex is
226
   begin
227
      return Complex'(Left, R (Right));
228
   end "+";
229
 
230
   ---------
231
   -- "-" --
232
   ---------
233
 
234
   function "-" (Right : Complex) return Complex is
235
   begin
236
      return (-Right.Re, -Right.Im);
237
   end "-";
238
 
239
   function "-" (Left, Right : Complex) return Complex is
240
   begin
241
      return (Left.Re - Right.Re, Left.Im - Right.Im);
242
   end "-";
243
 
244
   function "-" (Right : Imaginary) return Imaginary is
245
   begin
246
      return Imaginary (-R (Right));
247
   end "-";
248
 
249
   function "-" (Left, Right : Imaginary) return Imaginary is
250
   begin
251
      return Imaginary (R (Left) - R (Right));
252
   end "-";
253
 
254
   function "-" (Left : Complex; Right : Real'Base) return Complex is
255
   begin
256
      return Complex'(Left.Re - Right, Left.Im);
257
   end "-";
258
 
259
   function "-" (Left : Real'Base; Right : Complex) return Complex is
260
   begin
261
      return Complex'(Left - Right.Re, -Right.Im);
262
   end "-";
263
 
264
   function "-" (Left : Complex; Right : Imaginary) return Complex is
265
   begin
266
      return Complex'(Left.Re, Left.Im - R (Right));
267
   end "-";
268
 
269
   function "-" (Left : Imaginary; Right : Complex) return Complex is
270
   begin
271
      return Complex'(-Right.Re, R (Left) - Right.Im);
272
   end "-";
273
 
274
   function "-" (Left : Imaginary; Right : Real'Base) return Complex is
275
   begin
276
      return Complex'(-Right, R (Left));
277
   end "-";
278
 
279
   function "-" (Left : Real'Base; Right : Imaginary) return Complex is
280
   begin
281
      return Complex'(Left, -R (Right));
282
   end "-";
283
 
284
   ---------
285
   -- "/" --
286
   ---------
287
 
288
   function "/" (Left, Right : Complex) return Complex is
289
      a : constant R := Left.Re;
290
      b : constant R := Left.Im;
291
      c : constant R := Right.Re;
292
      d : constant R := Right.Im;
293
 
294
   begin
295
      if c = 0.0 and then d = 0.0 then
296
         raise Constraint_Error;
297
      else
298
         return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
299
                         Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
300
      end if;
301
   end "/";
302
 
303
   function "/" (Left, Right : Imaginary) return Real'Base is
304
   begin
305
      return R (Left) / R (Right);
306
   end "/";
307
 
308
   function "/" (Left : Complex; Right : Real'Base) return Complex is
309
   begin
310
      return Complex'(Left.Re / Right, Left.Im / Right);
311
   end "/";
312
 
313
   function "/" (Left : Real'Base; Right : Complex) return Complex is
314
      a : constant R := Left;
315
      c : constant R := Right.Re;
316
      d : constant R := Right.Im;
317
   begin
318
      return Complex'(Re =>   (a * c) / (c ** 2 + d ** 2),
319
                      Im => -((a * d) / (c ** 2 + d ** 2)));
320
   end "/";
321
 
322
   function "/" (Left : Complex; Right : Imaginary) return Complex is
323
      a : constant R := Left.Re;
324
      b : constant R := Left.Im;
325
      d : constant R := R (Right);
326
 
327
   begin
328
      return (b / d,  -(a / d));
329
   end "/";
330
 
331
   function "/" (Left : Imaginary; Right : Complex) return Complex is
332
      b : constant R := R (Left);
333
      c : constant R := Right.Re;
334
      d : constant R := Right.Im;
335
 
336
   begin
337
      return (Re => b * d / (c ** 2 + d ** 2),
338
              Im => b * c / (c ** 2 + d ** 2));
339
   end "/";
340
 
341
   function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
342
   begin
343
      return Imaginary (R (Left) / Right);
344
   end "/";
345
 
346
   function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
347
   begin
348
      return Imaginary (-(Left / R (Right)));
349
   end "/";
350
 
351
   ---------
352
   -- "<" --
353
   ---------
354
 
355
   function "<" (Left, Right : Imaginary) return Boolean is
356
   begin
357
      return R (Left) < R (Right);
358
   end "<";
359
 
360
   ----------
361
   -- "<=" --
362
   ----------
363
 
364
   function "<=" (Left, Right : Imaginary) return Boolean is
365
   begin
366
      return R (Left) <= R (Right);
367
   end "<=";
368
 
369
   ---------
370
   -- ">" --
371
   ---------
372
 
373
   function ">" (Left, Right : Imaginary) return Boolean is
374
   begin
375
      return R (Left) > R (Right);
376
   end ">";
377
 
378
   ----------
379
   -- ">=" --
380
   ----------
381
 
382
   function ">=" (Left, Right : Imaginary) return Boolean is
383
   begin
384
      return R (Left) >= R (Right);
385
   end ">=";
386
 
387
   -----------
388
   -- "abs" --
389
   -----------
390
 
391
   function "abs" (Right : Imaginary) return Real'Base is
392
   begin
393
      return abs R (Right);
394
   end "abs";
395
 
396
   --------------
397
   -- Argument --
398
   --------------
399
 
400
   function Argument (X : Complex) return Real'Base is
401
      a   : constant R := X.Re;
402
      b   : constant R := X.Im;
403
      arg : R;
404
 
405
   begin
406
      if b = 0.0 then
407
 
408
         if a >= 0.0 then
409
            return 0.0;
410
         else
411
            return R'Copy_Sign (Pi, b);
412
         end if;
413
 
414
      elsif a = 0.0 then
415
 
416
         if b >= 0.0 then
417
            return Half_Pi;
418
         else
419
            return -Half_Pi;
420
         end if;
421
 
422
      else
423
         arg := R (Atan (Double (abs (b / a))));
424
 
425
         if a > 0.0 then
426
            if b > 0.0 then
427
               return arg;
428
            else                  --  b < 0.0
429
               return -arg;
430
            end if;
431
 
432
         else                     --  a < 0.0
433
            if b >= 0.0 then
434
               return Pi - arg;
435
            else                  --  b < 0.0
436
               return -(Pi - arg);
437
            end if;
438
         end if;
439
      end if;
440
 
441
   exception
442
      when Constraint_Error =>
443
         if b > 0.0 then
444
            return Half_Pi;
445
         else
446
            return -Half_Pi;
447
         end if;
448
   end Argument;
449
 
450
   function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
451
   begin
452
      if Cycle > 0.0 then
453
         return Argument (X) * Cycle / Two_Pi;
454
      else
455
         raise Argument_Error;
456
      end if;
457
   end Argument;
458
 
459
   ----------------------------
460
   -- Compose_From_Cartesian --
461
   ----------------------------
462
 
463
   function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
464
   begin
465
      return (Re, Im);
466
   end Compose_From_Cartesian;
467
 
468
   function Compose_From_Cartesian (Re : Real'Base) return Complex is
469
   begin
470
      return (Re, 0.0);
471
   end Compose_From_Cartesian;
472
 
473
   function Compose_From_Cartesian (Im : Imaginary) return Complex is
474
   begin
475
      return (0.0, R (Im));
476
   end Compose_From_Cartesian;
477
 
478
   ------------------------
479
   -- Compose_From_Polar --
480
   ------------------------
481
 
482
   function Compose_From_Polar (
483
     Modulus, Argument : Real'Base)
484
     return Complex
485
   is
486
   begin
487
      if Modulus = 0.0 then
488
         return (0.0, 0.0);
489
      else
490
         return (Modulus * R (Cos (Double (Argument))),
491
                 Modulus * R (Sin (Double (Argument))));
492
      end if;
493
   end Compose_From_Polar;
494
 
495
   function Compose_From_Polar (
496
     Modulus, Argument, Cycle : Real'Base)
497
     return Complex
498
   is
499
      Arg : Real'Base;
500
 
501
   begin
502
      if Modulus = 0.0 then
503
         return (0.0, 0.0);
504
 
505
      elsif Cycle > 0.0 then
506
         if Argument = 0.0 then
507
            return (Modulus, 0.0);
508
 
509
         elsif Argument = Cycle / 4.0 then
510
            return (0.0, Modulus);
511
 
512
         elsif Argument = Cycle / 2.0 then
513
            return (-Modulus, 0.0);
514
 
515
         elsif Argument = 3.0 * Cycle / R (4.0) then
516
            return (0.0, -Modulus);
517
         else
518
            Arg := Two_Pi * Argument / Cycle;
519
            return (Modulus * R (Cos (Double (Arg))),
520
                    Modulus * R (Sin (Double (Arg))));
521
         end if;
522
      else
523
         raise Argument_Error;
524
      end if;
525
   end Compose_From_Polar;
526
 
527
   ---------------
528
   -- Conjugate --
529
   ---------------
530
 
531
   function Conjugate (X : Complex) return Complex is
532
   begin
533
      return Complex'(X.Re, -X.Im);
534
   end Conjugate;
535
 
536
   --------
537
   -- Im --
538
   --------
539
 
540
   function Im (X : Complex) return Real'Base is
541
   begin
542
      return X.Im;
543
   end Im;
544
 
545
   function Im (X : Imaginary) return Real'Base is
546
   begin
547
      return R (X);
548
   end Im;
549
 
550
   -------------
551
   -- Modulus --
552
   -------------
553
 
554
   function Modulus (X : Complex) return Real'Base is
555
      Re2, Im2 : R;
556
 
557
   begin
558
 
559
      begin
560
         Re2 := X.Re ** 2;
561
 
562
         --  To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
563
         --  compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
564
         --  squaring does not raise constraint_error but generates infinity,
565
         --  we can use an explicit comparison to determine whether to use
566
         --  the scaling expression.
567
 
568
         --  The scaling expression is computed in double format throughout
569
         --  in order to prevent inaccuracies on machines where not all
570
         --  immediate expressions are rounded, such as PowerPC.
571
 
572
         if Re2 > R'Last then
573
            raise Constraint_Error;
574
         end if;
575
 
576
      exception
577
         when Constraint_Error =>
578
            return R (Double (abs (X.Re))
579
              * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
580
      end;
581
 
582
      begin
583
         Im2 := X.Im ** 2;
584
 
585
         if Im2 > R'Last then
586
            raise Constraint_Error;
587
         end if;
588
 
589
      exception
590
         when Constraint_Error =>
591
            return R (Double (abs (X.Im))
592
              * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
593
      end;
594
 
595
      --  Now deal with cases of underflow. If only one of the squares
596
      --  underflows, return the modulus of the other component. If both
597
      --  squares underflow, use scaling as above.
598
 
599
      if Re2 = 0.0 then
600
 
601
         if X.Re = 0.0 then
602
            return abs (X.Im);
603
 
604
         elsif Im2 = 0.0 then
605
 
606
            if X.Im = 0.0 then
607
               return abs (X.Re);
608
 
609
            else
610
               if abs (X.Re) > abs (X.Im) then
611
                  return
612
                    R (Double (abs (X.Re))
613
                      * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
614
               else
615
                  return
616
                    R (Double (abs (X.Im))
617
                      * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
618
               end if;
619
            end if;
620
 
621
         else
622
            return abs (X.Im);
623
         end if;
624
 
625
      elsif Im2 = 0.0 then
626
         return abs (X.Re);
627
 
628
      --  In all other cases, the naive computation will do
629
 
630
      else
631
         return R (Sqrt (Double (Re2 + Im2)));
632
      end if;
633
   end Modulus;
634
 
635
   --------
636
   -- Re --
637
   --------
638
 
639
   function Re (X : Complex) return Real'Base is
640
   begin
641
      return X.Re;
642
   end Re;
643
 
644
   ------------
645
   -- Set_Im --
646
   ------------
647
 
648
   procedure Set_Im (X : in out Complex; Im : Real'Base) is
649
   begin
650
      X.Im := Im;
651
   end Set_Im;
652
 
653
   procedure Set_Im (X : out Imaginary; Im : Real'Base) is
654
   begin
655
      X := Imaginary (Im);
656
   end Set_Im;
657
 
658
   ------------
659
   -- Set_Re --
660
   ------------
661
 
662
   procedure Set_Re (X : in out Complex; Re : Real'Base) is
663
   begin
664
      X.Re := Re;
665
   end Set_Re;
666
 
667
end Ada.Numerics.Generic_Complex_Types;

powered by: WebSVN 2.1.0

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