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