1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS --
10 -- Copyright (C) 1992-2001 Free Software Foundation, Inc.
12 -- GNAT is free software; you can redistribute it and/or modify it under --
13 -- terms of the GNU General Public License as published by the Free Soft- --
14 -- ware Foundation; either version 2, or (at your option) any later ver- --
15 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
16 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
17 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
18 -- for more details. You should have received a copy of the GNU General --
19 -- Public License distributed with GNAT; see file COPYING. If not, write --
20 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
21 -- MA 02111-1307, USA. --
23 -- As a special exception, if other files instantiate generics from this --
24 -- unit, or you link this unit with other files to produce an executable, --
25 -- this unit does not by itself cause the resulting executable to be --
26 -- covered by the GNU General Public License. This exception does not --
27 -- however invalidate any other reasons why the executable file might be --
28 -- covered by the GNU Public License. --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
33 ------------------------------------------------------------------------------
35 with Ada.Numerics.Generic_Elementary_Functions;
37 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
39 package Elementary_Functions is new
40 Ada.Numerics.Generic_Elementary_Functions (Real'Base);
41 use Elementary_Functions;
43 PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
44 PI_2 : constant := PI / 2.0;
45 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
46 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
48 subtype T is Real'Base;
50 Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);
51 Square_Root_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
52 Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
53 Root_Root_Epsilon : constant T := Sqrt_Two **
54 ((1 - T'Model_Mantissa) / 2);
55 Log_Inverse_Epsilon_2 : constant T := T (T'Model_Mantissa - 1) / 2.0;
57 Complex_Zero : constant Complex := (0.0, 0.0);
58 Complex_One : constant Complex := (1.0, 0.0);
59 Complex_I : constant Complex := (0.0, 1.0);
60 Half_Pi : constant Complex := (PI_2, 0.0);
66 function "**" (Left : Complex; Right : Complex) return Complex is
69 and then Im (Right) = 0.0
70 and then Re (Left) = 0.0
71 and then Im (Left) = 0.0
76 and then Im (Left) = 0.0
77 and then Re (Right) < 0.0
79 raise Constraint_Error;
81 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
84 elsif Right = (0.0, 0.0) then
87 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
90 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
94 return Exp (Right * Log (Left));
98 function "**" (Left : Real'Base; Right : Complex) return Complex is
100 if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
101 raise Argument_Error;
103 elsif Left = 0.0 and then Re (Right) < 0.0 then
104 raise Constraint_Error;
106 elsif Left = 0.0 then
107 return Compose_From_Cartesian (Left, 0.0);
109 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
112 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
113 return Compose_From_Cartesian (Left, 0.0);
116 return Exp (Log (Left) * Right);
120 function "**" (Left : Complex; Right : Real'Base) return Complex is
123 and then Re (Left) = 0.0
124 and then Im (Left) = 0.0
126 raise Argument_Error;
128 elsif Re (Left) = 0.0
129 and then Im (Left) = 0.0
132 raise Constraint_Error;
134 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
137 elsif Right = 0.0 then
140 elsif Right = 1.0 then
144 return Exp (Right * Log (Left));
152 function Arccos (X : Complex) return Complex is
156 if X = Complex_One then
159 elsif abs Re (X) < Square_Root_Epsilon and then
160 abs Im (X) < Square_Root_Epsilon
164 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
165 abs Im (X) > Inv_Square_Root_Epsilon
167 return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
168 Complex_I * Sqrt ((1.0 - X) / 2.0));
171 Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
174 and then abs Re (X) <= 1.00
176 Set_Im (Result, Im (X));
186 function Arccosh (X : Complex) return Complex is
190 if X = Complex_One then
193 elsif abs Re (X) < Square_Root_Epsilon and then
194 abs Im (X) < Square_Root_Epsilon
196 Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
198 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
199 abs Im (X) > Inv_Square_Root_Epsilon
201 Result := Log_Two + Log (X);
204 Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
205 Sqrt ((X - 1.0) / 2.0));
208 if Re (Result) <= 0.0 then
219 function Arccot (X : Complex) return Complex is
223 if abs Re (X) < Square_Root_Epsilon and then
224 abs Im (X) < Square_Root_Epsilon
228 elsif abs Re (X) > 1.0 / Epsilon or else
229 abs Im (X) > 1.0 / Epsilon
231 Xt := Complex_One / X;
234 Set_Re (Xt, PI - Re (Xt));
241 Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
243 if Re (Xt) < 0.0 then
254 function Arccoth (X : Complex) return Complex is
258 if X = (0.0, 0.0) then
259 return Compose_From_Cartesian (0.0, PI_2);
261 elsif abs Re (X) < Square_Root_Epsilon
262 and then abs Im (X) < Square_Root_Epsilon
264 return PI_2 * Complex_I + X;
266 elsif abs Re (X) > 1.0 / Epsilon or else
267 abs Im (X) > 1.0 / Epsilon
272 return PI * Complex_I;
275 elsif Im (X) = 0.0 and then Re (X) = 1.0 then
276 raise Constraint_Error;
278 elsif Im (X) = 0.0 and then Re (X) = -1.0 then
279 raise Constraint_Error;
283 R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
286 when Constraint_Error =>
287 R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
291 Set_Im (R, PI + Im (R));
305 function Arcsin (X : Complex) return Complex is
309 if abs Re (X) < Square_Root_Epsilon and then
310 abs Im (X) < Square_Root_Epsilon
314 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
315 abs Im (X) > Inv_Square_Root_Epsilon
317 Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
319 if Im (Result) > PI_2 then
320 Set_Im (Result, PI - Im (X));
322 elsif Im (Result) < -PI_2 then
323 Set_Im (Result, -(PI + Im (X)));
327 Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
330 Set_Re (Result, Re (X));
333 and then abs Re (X) <= 1.00
335 Set_Im (Result, Im (X));
345 function Arcsinh (X : Complex) return Complex is
349 if abs Re (X) < Square_Root_Epsilon and then
350 abs Im (X) < Square_Root_Epsilon
354 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
355 abs Im (X) > Inv_Square_Root_Epsilon
357 Result := Log_Two + Log (X); -- may have wrong sign
359 if (Re (X) < 0.0 and Re (Result) > 0.0)
360 or else (Re (X) > 0.0 and Re (Result) < 0.0)
362 Set_Re (Result, -Re (Result));
368 Result := Log (X + Sqrt (1.0 + X * X));
371 Set_Re (Result, Re (X));
372 elsif Im (X) = 0.0 then
373 Set_Im (Result, Im (X));
383 function Arctan (X : Complex) return Complex is
385 if abs Re (X) < Square_Root_Epsilon and then
386 abs Im (X) < Square_Root_Epsilon
391 return -Complex_I * (Log (1.0 + Complex_I * X)
392 - Log (1.0 - Complex_I * X)) / 2.0;
400 function Arctanh (X : Complex) return Complex is
402 if abs Re (X) < Square_Root_Epsilon and then
403 abs Im (X) < Square_Root_Epsilon
407 return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
415 function Cos (X : Complex) return Complex is
418 Compose_From_Cartesian
419 (Cos (Re (X)) * Cosh (Im (X)),
420 -Sin (Re (X)) * Sinh (Im (X)));
427 function Cosh (X : Complex) return Complex is
430 Compose_From_Cartesian
431 (Cosh (Re (X)) * Cos (Im (X)),
432 Sinh (Re (X)) * Sin (Im (X)));
439 function Cot (X : Complex) return Complex is
441 if abs Re (X) < Square_Root_Epsilon and then
442 abs Im (X) < Square_Root_Epsilon
444 return Complex_One / X;
446 elsif Im (X) > Log_Inverse_Epsilon_2 then
449 elsif Im (X) < -Log_Inverse_Epsilon_2 then
453 return Cos (X) / Sin (X);
460 function Coth (X : Complex) return Complex is
462 if abs Re (X) < Square_Root_Epsilon and then
463 abs Im (X) < Square_Root_Epsilon
465 return Complex_One / X;
467 elsif Re (X) > Log_Inverse_Epsilon_2 then
470 elsif Re (X) < -Log_Inverse_Epsilon_2 then
474 return Cosh (X) / Sinh (X);
482 function Exp (X : Complex) return Complex is
483 EXP_RE_X : Real'Base := Exp (Re (X));
486 return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
487 EXP_RE_X * Sin (Im (X)));
491 function Exp (X : Imaginary) return Complex is
492 ImX : Real'Base := Im (X);
495 return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
502 function Log (X : Complex) return Complex is
508 if Re (X) = 0.0 and then Im (X) = 0.0 then
509 raise Constraint_Error;
511 elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
512 and then abs Im (X) < Root_Root_Epsilon
515 Set_Re (Z, Re (Z) - 1.0);
517 return (1.0 - (1.0 / 2.0 -
518 (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
522 ReX := Log (Modulus (X));
525 when Constraint_Error =>
526 ReX := Log (Modulus (X / 2.0)) - Log_Two;
529 ImX := Arctan (Im (X), Re (X));
532 ImX := ImX - 2.0 * PI;
535 return Compose_From_Cartesian (ReX, ImX);
542 function Sin (X : Complex) return Complex is
544 if abs Re (X) < Square_Root_Epsilon and then
545 abs Im (X) < Square_Root_Epsilon then
550 Compose_From_Cartesian
551 (Sin (Re (X)) * Cosh (Im (X)),
552 Cos (Re (X)) * Sinh (Im (X)));
559 function Sinh (X : Complex) return Complex is
561 if abs Re (X) < Square_Root_Epsilon and then
562 abs Im (X) < Square_Root_Epsilon
567 return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
568 Cosh (Re (X)) * Sin (Im (X)));
576 function Sqrt (X : Complex) return Complex is
577 ReX : constant Real'Base := Re (X);
578 ImX : constant Real'Base := Im (X);
579 XR : constant Real'Base := abs Re (X);
580 YR : constant Real'Base := abs Im (X);
586 -- Deal with pure real case, see (RM G.1.2(39))
591 Compose_From_Cartesian
599 Compose_From_Cartesian
600 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
604 R_X := Sqrt (YR / 2.0);
607 return Compose_From_Cartesian (R_X, R_X);
609 return Compose_From_Cartesian (R_X, -R_X);
613 R := Sqrt (XR ** 2 + YR ** 2);
615 -- If the square of the modulus overflows, try rescaling the
616 -- real and imaginary parts. We cannot depend on an exception
617 -- being raised on all targets.
619 if R > Real'Base'Last then
620 raise Constraint_Error;
623 -- We are solving the system
625 -- XR = R_X ** 2 - Y_R ** 2 (1)
626 -- YR = 2.0 * R_X * R_Y (2)
628 -- The symmetric solution involves square roots for both R_X and
629 -- R_Y, but it is more accurate to use the square root with the
630 -- larger argument for either R_X or R_Y, and equation (2) for the
634 R_Y := Sqrt (0.5 * (R - ReX));
635 R_X := YR / (2.0 * R_Y);
638 R_X := Sqrt (0.5 * (R + ReX));
639 R_Y := YR / (2.0 * R_X);
643 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
646 return Compose_From_Cartesian (R_X, R_Y);
649 when Constraint_Error =>
651 -- Rescale and try again.
653 R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
654 R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
655 R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
657 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
661 return Compose_From_Cartesian (R_X, R_Y);
668 function Tan (X : Complex) return Complex is
670 if abs Re (X) < Square_Root_Epsilon and then
671 abs Im (X) < Square_Root_Epsilon
675 elsif Im (X) > Log_Inverse_Epsilon_2 then
678 elsif Im (X) < -Log_Inverse_Epsilon_2 then
682 return Sin (X) / Cos (X);
690 function Tanh (X : Complex) return Complex is
692 if abs Re (X) < Square_Root_Epsilon and then
693 abs Im (X) < Square_Root_Epsilon
697 elsif Re (X) > Log_Inverse_Epsilon_2 then
700 elsif Re (X) < -Log_Inverse_Epsilon_2 then
704 return Sinh (X) / Cosh (X);
708 end Ada.Numerics.Generic_Complex_Elementary_Functions;