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