1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS --
9 -- Copyright (C) 1992-2003 Free Software Foundation, Inc. --
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 2, 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. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
20 -- MA 02111-1307, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 with Ada.Numerics.Generic_Elementary_Functions;
36 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
38 package Elementary_Functions is new
39 Ada.Numerics.Generic_Elementary_Functions (Real'Base);
40 use Elementary_Functions;
42 PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
43 PI_2 : constant := PI / 2.0;
44 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
45 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
47 subtype T is Real'Base;
49 Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);
50 Square_Root_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
51 Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
52 Root_Root_Epsilon : constant T := Sqrt_Two **
53 ((1 - T'Model_Mantissa) / 2);
54 Log_Inverse_Epsilon_2 : constant T := T (T'Model_Mantissa - 1) / 2.0;
56 Complex_Zero : constant Complex := (0.0, 0.0);
57 Complex_One : constant Complex := (1.0, 0.0);
58 Complex_I : constant Complex := (0.0, 1.0);
59 Half_Pi : constant Complex := (PI_2, 0.0);
65 function "**" (Left : Complex; Right : Complex) return Complex is
68 and then Im (Right) = 0.0
69 and then Re (Left) = 0.0
70 and then Im (Left) = 0.0
75 and then Im (Left) = 0.0
76 and then Re (Right) < 0.0
78 raise Constraint_Error;
80 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
83 elsif Right = (0.0, 0.0) then
86 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
89 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
93 return Exp (Right * Log (Left));
97 function "**" (Left : Real'Base; Right : Complex) return Complex is
99 if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
100 raise Argument_Error;
102 elsif Left = 0.0 and then Re (Right) < 0.0 then
103 raise Constraint_Error;
105 elsif Left = 0.0 then
106 return Compose_From_Cartesian (Left, 0.0);
108 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
111 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
112 return Compose_From_Cartesian (Left, 0.0);
115 return Exp (Log (Left) * Right);
119 function "**" (Left : Complex; Right : Real'Base) return Complex is
122 and then Re (Left) = 0.0
123 and then Im (Left) = 0.0
125 raise Argument_Error;
127 elsif Re (Left) = 0.0
128 and then Im (Left) = 0.0
131 raise Constraint_Error;
133 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
136 elsif Right = 0.0 then
139 elsif Right = 1.0 then
143 return Exp (Right * Log (Left));
151 function Arccos (X : Complex) return Complex is
155 if X = Complex_One then
158 elsif abs Re (X) < Square_Root_Epsilon and then
159 abs Im (X) < Square_Root_Epsilon
163 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
164 abs Im (X) > Inv_Square_Root_Epsilon
166 return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
167 Complex_I * Sqrt ((1.0 - X) / 2.0));
170 Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
173 and then abs Re (X) <= 1.00
175 Set_Im (Result, Im (X));
185 function Arccosh (X : Complex) return Complex is
189 if X = Complex_One then
192 elsif abs Re (X) < Square_Root_Epsilon and then
193 abs Im (X) < Square_Root_Epsilon
195 Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
197 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
198 abs Im (X) > Inv_Square_Root_Epsilon
200 Result := Log_Two + Log (X);
203 Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
204 Sqrt ((X - 1.0) / 2.0));
207 if Re (Result) <= 0.0 then
218 function Arccot (X : Complex) return Complex is
222 if abs Re (X) < Square_Root_Epsilon and then
223 abs Im (X) < Square_Root_Epsilon
227 elsif abs Re (X) > 1.0 / Epsilon or else
228 abs Im (X) > 1.0 / Epsilon
230 Xt := Complex_One / X;
233 Set_Re (Xt, PI - Re (Xt));
240 Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
242 if Re (Xt) < 0.0 then
253 function Arccoth (X : Complex) return Complex is
257 if X = (0.0, 0.0) then
258 return Compose_From_Cartesian (0.0, PI_2);
260 elsif abs Re (X) < Square_Root_Epsilon
261 and then abs Im (X) < Square_Root_Epsilon
263 return PI_2 * Complex_I + X;
265 elsif abs Re (X) > 1.0 / Epsilon or else
266 abs Im (X) > 1.0 / Epsilon
271 return PI * Complex_I;
274 elsif Im (X) = 0.0 and then Re (X) = 1.0 then
275 raise Constraint_Error;
277 elsif Im (X) = 0.0 and then Re (X) = -1.0 then
278 raise Constraint_Error;
282 R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
285 when Constraint_Error =>
286 R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
290 Set_Im (R, PI + Im (R));
304 function Arcsin (X : Complex) return Complex is
308 -- For very small argument, sin (x) = x.
310 if abs Re (X) < Square_Root_Epsilon and then
311 abs Im (X) < Square_Root_Epsilon
315 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
316 abs Im (X) > Inv_Square_Root_Epsilon
318 Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
320 if Im (Result) > PI_2 then
321 Set_Im (Result, PI - Im (X));
323 elsif Im (Result) < -PI_2 then
324 Set_Im (Result, -(PI + Im (X)));
330 Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
333 Set_Re (Result, Re (X));
336 and then abs Re (X) <= 1.00
338 Set_Im (Result, Im (X));
348 function Arcsinh (X : Complex) return Complex is
352 if abs Re (X) < Square_Root_Epsilon and then
353 abs Im (X) < Square_Root_Epsilon
357 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
358 abs Im (X) > Inv_Square_Root_Epsilon
360 Result := Log_Two + Log (X); -- may have wrong sign
362 if (Re (X) < 0.0 and Re (Result) > 0.0)
363 or else (Re (X) > 0.0 and Re (Result) < 0.0)
365 Set_Re (Result, -Re (Result));
371 Result := Log (X + Sqrt (1.0 + X * X));
374 Set_Re (Result, Re (X));
375 elsif Im (X) = 0.0 then
376 Set_Im (Result, Im (X));
386 function Arctan (X : Complex) return Complex is
388 if abs Re (X) < Square_Root_Epsilon and then
389 abs Im (X) < Square_Root_Epsilon
394 return -Complex_I * (Log (1.0 + Complex_I * X)
395 - Log (1.0 - Complex_I * X)) / 2.0;
403 function Arctanh (X : Complex) return Complex is
405 if abs Re (X) < Square_Root_Epsilon and then
406 abs Im (X) < Square_Root_Epsilon
410 return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
418 function Cos (X : Complex) return Complex is
421 Compose_From_Cartesian
422 (Cos (Re (X)) * Cosh (Im (X)),
423 -Sin (Re (X)) * Sinh (Im (X)));
430 function Cosh (X : Complex) return Complex is
433 Compose_From_Cartesian
434 (Cosh (Re (X)) * Cos (Im (X)),
435 Sinh (Re (X)) * Sin (Im (X)));
442 function Cot (X : Complex) return Complex is
444 if abs Re (X) < Square_Root_Epsilon and then
445 abs Im (X) < Square_Root_Epsilon
447 return Complex_One / X;
449 elsif Im (X) > Log_Inverse_Epsilon_2 then
452 elsif Im (X) < -Log_Inverse_Epsilon_2 then
456 return Cos (X) / Sin (X);
463 function Coth (X : Complex) return Complex is
465 if abs Re (X) < Square_Root_Epsilon and then
466 abs Im (X) < Square_Root_Epsilon
468 return Complex_One / X;
470 elsif Re (X) > Log_Inverse_Epsilon_2 then
473 elsif Re (X) < -Log_Inverse_Epsilon_2 then
477 return Cosh (X) / Sinh (X);
485 function Exp (X : Complex) return Complex is
486 EXP_RE_X : constant Real'Base := Exp (Re (X));
489 return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
490 EXP_RE_X * Sin (Im (X)));
493 function Exp (X : Imaginary) return Complex is
494 ImX : constant Real'Base := Im (X);
497 return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
504 function Log (X : Complex) return Complex is
510 if Re (X) = 0.0 and then Im (X) = 0.0 then
511 raise Constraint_Error;
513 elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
514 and then abs Im (X) < Root_Root_Epsilon
517 Set_Re (Z, Re (Z) - 1.0);
519 return (1.0 - (1.0 / 2.0 -
520 (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
524 ReX := Log (Modulus (X));
527 when Constraint_Error =>
528 ReX := Log (Modulus (X / 2.0)) - Log_Two;
531 ImX := Arctan (Im (X), Re (X));
534 ImX := ImX - 2.0 * PI;
537 return Compose_From_Cartesian (ReX, ImX);
544 function Sin (X : Complex) return Complex is
546 if abs Re (X) < Square_Root_Epsilon and then
547 abs Im (X) < Square_Root_Epsilon then
552 Compose_From_Cartesian
553 (Sin (Re (X)) * Cosh (Im (X)),
554 Cos (Re (X)) * Sinh (Im (X)));
561 function Sinh (X : Complex) return Complex is
563 if abs Re (X) < Square_Root_Epsilon and then
564 abs Im (X) < Square_Root_Epsilon
569 return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
570 Cosh (Re (X)) * Sin (Im (X)));
578 function Sqrt (X : Complex) return Complex is
579 ReX : constant Real'Base := Re (X);
580 ImX : constant Real'Base := Im (X);
581 XR : constant Real'Base := abs Re (X);
582 YR : constant Real'Base := abs Im (X);
588 -- Deal with pure real case, see (RM G.1.2(39))
593 Compose_From_Cartesian
601 Compose_From_Cartesian
602 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
606 R_X := Sqrt (YR / 2.0);
609 return Compose_From_Cartesian (R_X, R_X);
611 return Compose_From_Cartesian (R_X, -R_X);
615 R := Sqrt (XR ** 2 + YR ** 2);
617 -- If the square of the modulus overflows, try rescaling the
618 -- real and imaginary parts. We cannot depend on an exception
619 -- being raised on all targets.
621 if R > Real'Base'Last then
622 raise Constraint_Error;
625 -- We are solving the system
627 -- XR = R_X ** 2 - Y_R ** 2 (1)
628 -- YR = 2.0 * R_X * R_Y (2)
630 -- The symmetric solution involves square roots for both R_X and
631 -- R_Y, but it is more accurate to use the square root with the
632 -- larger argument for either R_X or R_Y, and equation (2) for the
636 R_Y := Sqrt (0.5 * (R - ReX));
637 R_X := YR / (2.0 * R_Y);
640 R_X := Sqrt (0.5 * (R + ReX));
641 R_Y := YR / (2.0 * R_X);
645 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
648 return Compose_From_Cartesian (R_X, R_Y);
651 when Constraint_Error =>
653 -- Rescale and try again.
655 R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
656 R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
657 R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
659 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
663 return Compose_From_Cartesian (R_X, R_Y);
670 function Tan (X : Complex) return Complex is
672 if abs Re (X) < Square_Root_Epsilon and then
673 abs Im (X) < Square_Root_Epsilon
677 elsif Im (X) > Log_Inverse_Epsilon_2 then
680 elsif Im (X) < -Log_Inverse_Epsilon_2 then
684 return Sin (X) / Cos (X);
692 function Tanh (X : Complex) return Complex is
694 if abs Re (X) < Square_Root_Epsilon and then
695 abs Im (X) < Square_Root_Epsilon
699 elsif Re (X) > Log_Inverse_Epsilon_2 then
702 elsif Re (X) < -Log_Inverse_Epsilon_2 then
706 return Sinh (X) / Cosh (X);
710 end Ada.Numerics.Generic_Complex_Elementary_Functions;