|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
|
9 -- Copyright (C) 1992-2017, 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.Generic_Elementary_Functions;
|
|
|
33
|
|
|
34 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
|
|
|
35
|
|
|
36 package Elementary_Functions is new
|
|
|
37 Ada.Numerics.Generic_Elementary_Functions (Real'Base);
|
|
|
38 use Elementary_Functions;
|
|
|
39
|
|
|
40 PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
|
|
|
41 PI_2 : constant := PI / 2.0;
|
|
|
42 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
|
|
|
43 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
|
|
|
44
|
|
|
45 subtype T is Real'Base;
|
|
|
46
|
|
|
47 Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);
|
|
|
48 Square_Root_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
|
|
|
49 Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
|
|
|
50 Root_Root_Epsilon : constant T := Sqrt_Two **
|
|
|
51 ((1 - T'Model_Mantissa) / 2);
|
|
|
52 Log_Inverse_Epsilon_2 : constant T := T (T'Model_Mantissa - 1) / 2.0;
|
|
|
53
|
|
|
54 Complex_Zero : constant Complex := (0.0, 0.0);
|
|
|
55 Complex_One : constant Complex := (1.0, 0.0);
|
|
|
56 Complex_I : constant Complex := (0.0, 1.0);
|
|
|
57 Half_Pi : constant Complex := (PI_2, 0.0);
|
|
|
58
|
|
|
59 --------
|
|
|
60 -- ** --
|
|
|
61 --------
|
|
|
62
|
|
|
63 function "**" (Left : Complex; Right : Complex) return Complex is
|
|
|
64 begin
|
|
|
65 if Re (Right) = 0.0
|
|
|
66 and then Im (Right) = 0.0
|
|
|
67 and then Re (Left) = 0.0
|
|
|
68 and then Im (Left) = 0.0
|
|
|
69 then
|
|
|
70 raise Argument_Error;
|
|
|
71
|
|
|
72 elsif Re (Left) = 0.0
|
|
|
73 and then Im (Left) = 0.0
|
|
|
74 and then Re (Right) < 0.0
|
|
|
75 then
|
|
|
76 raise Constraint_Error;
|
|
|
77
|
|
|
78 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
|
|
|
79 return Left;
|
|
|
80
|
|
|
81 elsif Right = (0.0, 0.0) then
|
|
|
82 return Complex_One;
|
|
|
83
|
|
|
84 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
|
|
|
85 return 1.0 + Right;
|
|
|
86
|
|
|
87 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
|
|
|
88 return Left;
|
|
|
89
|
|
|
90 else
|
|
|
91 return Exp (Right * Log (Left));
|
|
|
92 end if;
|
|
|
93 end "**";
|
|
|
94
|
|
|
95 function "**" (Left : Real'Base; Right : Complex) return Complex is
|
|
|
96 begin
|
|
|
97 if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
|
|
|
98 raise Argument_Error;
|
|
|
99
|
|
|
100 elsif Left = 0.0 and then Re (Right) < 0.0 then
|
|
|
101 raise Constraint_Error;
|
|
|
102
|
|
|
103 elsif Left = 0.0 then
|
|
|
104 return Compose_From_Cartesian (Left, 0.0);
|
|
|
105
|
|
|
106 elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
|
|
|
107 return Complex_One;
|
|
|
108
|
|
|
109 elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
|
|
|
110 return Compose_From_Cartesian (Left, 0.0);
|
|
|
111
|
|
|
112 else
|
|
|
113 return Exp (Log (Left) * Right);
|
|
|
114 end if;
|
|
|
115 end "**";
|
|
|
116
|
|
|
117 function "**" (Left : Complex; Right : Real'Base) return Complex is
|
|
|
118 begin
|
|
|
119 if Right = 0.0
|
|
|
120 and then Re (Left) = 0.0
|
|
|
121 and then Im (Left) = 0.0
|
|
|
122 then
|
|
|
123 raise Argument_Error;
|
|
|
124
|
|
|
125 elsif Re (Left) = 0.0
|
|
|
126 and then Im (Left) = 0.0
|
|
|
127 and then Right < 0.0
|
|
|
128 then
|
|
|
129 raise Constraint_Error;
|
|
|
130
|
|
|
131 elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
|
|
|
132 return Left;
|
|
|
133
|
|
|
134 elsif Right = 0.0 then
|
|
|
135 return Complex_One;
|
|
|
136
|
|
|
137 elsif Right = 1.0 then
|
|
|
138 return Left;
|
|
|
139
|
|
|
140 else
|
|
|
141 return Exp (Right * Log (Left));
|
|
|
142 end if;
|
|
|
143 end "**";
|
|
|
144
|
|
|
145 ------------
|
|
|
146 -- Arccos --
|
|
|
147 ------------
|
|
|
148
|
|
|
149 function Arccos (X : Complex) return Complex is
|
|
|
150 Result : Complex;
|
|
|
151
|
|
|
152 begin
|
|
|
153 if X = Complex_One then
|
|
|
154 return Complex_Zero;
|
|
|
155
|
|
|
156 elsif abs Re (X) < Square_Root_Epsilon and then
|
|
|
157 abs Im (X) < Square_Root_Epsilon
|
|
|
158 then
|
|
|
159 return Half_Pi - X;
|
|
|
160
|
|
|
161 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
|
162 abs Im (X) > Inv_Square_Root_Epsilon
|
|
|
163 then
|
|
|
164 return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
|
|
|
165 Complex_I * Sqrt ((1.0 - X) / 2.0));
|
|
|
166 end if;
|
|
|
167
|
|
|
168 Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
|
|
|
169
|
|
|
170 if Im (X) = 0.0
|
|
|
171 and then abs Re (X) <= 1.00
|
|
|
172 then
|
|
|
173 Set_Im (Result, Im (X));
|
|
|
174 end if;
|
|
|
175
|
|
|
176 return Result;
|
|
|
177 end Arccos;
|
|
|
178
|
|
|
179 -------------
|
|
|
180 -- Arccosh --
|
|
|
181 -------------
|
|
|
182
|
|
|
183 function Arccosh (X : Complex) return Complex is
|
|
|
184 Result : Complex;
|
|
|
185
|
|
|
186 begin
|
|
|
187 if X = Complex_One then
|
|
|
188 return Complex_Zero;
|
|
|
189
|
|
|
190 elsif abs Re (X) < Square_Root_Epsilon and then
|
|
|
191 abs Im (X) < Square_Root_Epsilon
|
|
|
192 then
|
|
|
193 Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
|
|
|
194
|
|
|
195 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
|
196 abs Im (X) > Inv_Square_Root_Epsilon
|
|
|
197 then
|
|
|
198 Result := Log_Two + Log (X);
|
|
|
199
|
|
|
200 else
|
|
|
201 Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
|
|
|
202 Sqrt ((X - 1.0) / 2.0));
|
|
|
203 end if;
|
|
|
204
|
|
|
205 if Re (Result) <= 0.0 then
|
|
|
206 Result := -Result;
|
|
|
207 end if;
|
|
|
208
|
|
|
209 return Result;
|
|
|
210 end Arccosh;
|
|
|
211
|
|
|
212 ------------
|
|
|
213 -- Arccot --
|
|
|
214 ------------
|
|
|
215
|
|
|
216 function Arccot (X : Complex) return Complex is
|
|
|
217 Xt : Complex;
|
|
|
218
|
|
|
219 begin
|
|
|
220 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
221 abs Im (X) < Square_Root_Epsilon
|
|
|
222 then
|
|
|
223 return Half_Pi - X;
|
|
|
224
|
|
|
225 elsif abs Re (X) > 1.0 / Epsilon or else
|
|
|
226 abs Im (X) > 1.0 / Epsilon
|
|
|
227 then
|
|
|
228 Xt := Complex_One / X;
|
|
|
229
|
|
|
230 if Re (X) < 0.0 then
|
|
|
231 Set_Re (Xt, PI - Re (Xt));
|
|
|
232 return Xt;
|
|
|
233 else
|
|
|
234 return Xt;
|
|
|
235 end if;
|
|
|
236 end if;
|
|
|
237
|
|
|
238 Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
|
|
|
239
|
|
|
240 if Re (Xt) < 0.0 then
|
|
|
241 Xt := PI + Xt;
|
|
|
242 end if;
|
|
|
243
|
|
|
244 return Xt;
|
|
|
245 end Arccot;
|
|
|
246
|
|
|
247 --------------
|
|
|
248 -- Arccoth --
|
|
|
249 --------------
|
|
|
250
|
|
|
251 function Arccoth (X : Complex) return Complex is
|
|
|
252 R : Complex;
|
|
|
253
|
|
|
254 begin
|
|
|
255 if X = (0.0, 0.0) then
|
|
|
256 return Compose_From_Cartesian (0.0, PI_2);
|
|
|
257
|
|
|
258 elsif abs Re (X) < Square_Root_Epsilon
|
|
|
259 and then abs Im (X) < Square_Root_Epsilon
|
|
|
260 then
|
|
|
261 return PI_2 * Complex_I + X;
|
|
|
262
|
|
|
263 elsif abs Re (X) > 1.0 / Epsilon or else
|
|
|
264 abs Im (X) > 1.0 / Epsilon
|
|
|
265 then
|
|
|
266 if Im (X) > 0.0 then
|
|
|
267 return (0.0, 0.0);
|
|
|
268 else
|
|
|
269 return PI * Complex_I;
|
|
|
270 end if;
|
|
|
271
|
|
|
272 elsif Im (X) = 0.0 and then Re (X) = 1.0 then
|
|
|
273 raise Constraint_Error;
|
|
|
274
|
|
|
275 elsif Im (X) = 0.0 and then Re (X) = -1.0 then
|
|
|
276 raise Constraint_Error;
|
|
|
277 end if;
|
|
|
278
|
|
|
279 begin
|
|
|
280 R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
|
|
|
281
|
|
|
282 exception
|
|
|
283 when Constraint_Error =>
|
|
|
284 R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
|
|
|
285 end;
|
|
|
286
|
|
|
287 if Im (R) < 0.0 then
|
|
|
288 Set_Im (R, PI + Im (R));
|
|
|
289 end if;
|
|
|
290
|
|
|
291 if Re (X) = 0.0 then
|
|
|
292 Set_Re (R, Re (X));
|
|
|
293 end if;
|
|
|
294
|
|
|
295 return R;
|
|
|
296 end Arccoth;
|
|
|
297
|
|
|
298 ------------
|
|
|
299 -- Arcsin --
|
|
|
300 ------------
|
|
|
301
|
|
|
302 function Arcsin (X : Complex) return Complex is
|
|
|
303 Result : Complex;
|
|
|
304
|
|
|
305 begin
|
|
|
306 -- For very small argument, sin (x) = x
|
|
|
307
|
|
|
308 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
309 abs Im (X) < Square_Root_Epsilon
|
|
|
310 then
|
|
|
311 return X;
|
|
|
312
|
|
|
313 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
|
314 abs Im (X) > Inv_Square_Root_Epsilon
|
|
|
315 then
|
|
|
316 Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
|
|
|
317
|
|
|
318 if Im (Result) > PI_2 then
|
|
|
319 Set_Im (Result, PI - Im (X));
|
|
|
320
|
|
|
321 elsif Im (Result) < -PI_2 then
|
|
|
322 Set_Im (Result, -(PI + Im (X)));
|
|
|
323 end if;
|
|
|
324
|
|
|
325 return Result;
|
|
|
326 end if;
|
|
|
327
|
|
|
328 Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
|
|
|
329
|
|
|
330 if Re (X) = 0.0 then
|
|
|
331 Set_Re (Result, Re (X));
|
|
|
332
|
|
|
333 elsif Im (X) = 0.0
|
|
|
334 and then abs Re (X) <= 1.00
|
|
|
335 then
|
|
|
336 Set_Im (Result, Im (X));
|
|
|
337 end if;
|
|
|
338
|
|
|
339 return Result;
|
|
|
340 end Arcsin;
|
|
|
341
|
|
|
342 -------------
|
|
|
343 -- Arcsinh --
|
|
|
344 -------------
|
|
|
345
|
|
|
346 function Arcsinh (X : Complex) return Complex is
|
|
|
347 Result : Complex;
|
|
|
348
|
|
|
349 begin
|
|
|
350 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
351 abs Im (X) < Square_Root_Epsilon
|
|
|
352 then
|
|
|
353 return X;
|
|
|
354
|
|
|
355 elsif abs Re (X) > Inv_Square_Root_Epsilon or else
|
|
|
356 abs Im (X) > Inv_Square_Root_Epsilon
|
|
|
357 then
|
|
|
358 Result := Log_Two + Log (X); -- may have wrong sign
|
|
|
359
|
|
|
360 if (Re (X) < 0.0 and then Re (Result) > 0.0)
|
|
|
361 or else (Re (X) > 0.0 and then Re (Result) < 0.0)
|
|
|
362 then
|
|
|
363 Set_Re (Result, -Re (Result));
|
|
|
364 end if;
|
|
|
365
|
|
|
366 return Result;
|
|
|
367 end if;
|
|
|
368
|
|
|
369 Result := Log (X + Sqrt (1.0 + X * X));
|
|
|
370
|
|
|
371 if Re (X) = 0.0 then
|
|
|
372 Set_Re (Result, Re (X));
|
|
|
373 elsif Im (X) = 0.0 then
|
|
|
374 Set_Im (Result, Im (X));
|
|
|
375 end if;
|
|
|
376
|
|
|
377 return Result;
|
|
|
378 end Arcsinh;
|
|
|
379
|
|
|
380 ------------
|
|
|
381 -- Arctan --
|
|
|
382 ------------
|
|
|
383
|
|
|
384 function Arctan (X : Complex) return Complex is
|
|
|
385 begin
|
|
|
386 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
387 abs Im (X) < Square_Root_Epsilon
|
|
|
388 then
|
|
|
389 return X;
|
|
|
390
|
|
|
391 else
|
|
|
392 return -Complex_I * (Log (1.0 + Complex_I * X)
|
|
|
393 - Log (1.0 - Complex_I * X)) / 2.0;
|
|
|
394 end if;
|
|
|
395 end Arctan;
|
|
|
396
|
|
|
397 -------------
|
|
|
398 -- Arctanh --
|
|
|
399 -------------
|
|
|
400
|
|
|
401 function Arctanh (X : Complex) return Complex is
|
|
|
402 begin
|
|
|
403 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
404 abs Im (X) < Square_Root_Epsilon
|
|
|
405 then
|
|
|
406 return X;
|
|
|
407 else
|
|
|
408 return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
|
|
|
409 end if;
|
|
|
410 end Arctanh;
|
|
|
411
|
|
|
412 ---------
|
|
|
413 -- Cos --
|
|
|
414 ---------
|
|
|
415
|
|
|
416 function Cos (X : Complex) return Complex is
|
|
|
417 begin
|
|
|
418 return
|
|
|
419 Compose_From_Cartesian
|
|
|
420 (Cos (Re (X)) * Cosh (Im (X)),
|
|
|
421 -(Sin (Re (X)) * Sinh (Im (X))));
|
|
|
422 end Cos;
|
|
|
423
|
|
|
424 ----------
|
|
|
425 -- Cosh --
|
|
|
426 ----------
|
|
|
427
|
|
|
428 function Cosh (X : Complex) return Complex is
|
|
|
429 begin
|
|
|
430 return
|
|
|
431 Compose_From_Cartesian
|
|
|
432 (Cosh (Re (X)) * Cos (Im (X)),
|
|
|
433 Sinh (Re (X)) * Sin (Im (X)));
|
|
|
434 end Cosh;
|
|
|
435
|
|
|
436 ---------
|
|
|
437 -- Cot --
|
|
|
438 ---------
|
|
|
439
|
|
|
440 function Cot (X : Complex) return Complex is
|
|
|
441 begin
|
|
|
442 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
443 abs Im (X) < Square_Root_Epsilon
|
|
|
444 then
|
|
|
445 return Complex_One / X;
|
|
|
446
|
|
|
447 elsif Im (X) > Log_Inverse_Epsilon_2 then
|
|
|
448 return -Complex_I;
|
|
|
449
|
|
|
450 elsif Im (X) < -Log_Inverse_Epsilon_2 then
|
|
|
451 return Complex_I;
|
|
|
452 end if;
|
|
|
453
|
|
|
454 return Cos (X) / Sin (X);
|
|
|
455 end Cot;
|
|
|
456
|
|
|
457 ----------
|
|
|
458 -- Coth --
|
|
|
459 ----------
|
|
|
460
|
|
|
461 function Coth (X : Complex) return Complex is
|
|
|
462 begin
|
|
|
463 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
464 abs Im (X) < Square_Root_Epsilon
|
|
|
465 then
|
|
|
466 return Complex_One / X;
|
|
|
467
|
|
|
468 elsif Re (X) > Log_Inverse_Epsilon_2 then
|
|
|
469 return Complex_One;
|
|
|
470
|
|
|
471 elsif Re (X) < -Log_Inverse_Epsilon_2 then
|
|
|
472 return -Complex_One;
|
|
|
473
|
|
|
474 else
|
|
|
475 return Cosh (X) / Sinh (X);
|
|
|
476 end if;
|
|
|
477 end Coth;
|
|
|
478
|
|
|
479 ---------
|
|
|
480 -- Exp --
|
|
|
481 ---------
|
|
|
482
|
|
|
483 function Exp (X : Complex) return Complex is
|
|
|
484 EXP_RE_X : constant Real'Base := Exp (Re (X));
|
|
|
485
|
|
|
486 begin
|
|
|
487 return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
|
|
|
488 EXP_RE_X * Sin (Im (X)));
|
|
|
489 end Exp;
|
|
|
490
|
|
|
491 function Exp (X : Imaginary) return Complex is
|
|
|
492 ImX : constant Real'Base := Im (X);
|
|
|
493
|
|
|
494 begin
|
|
|
495 return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
|
|
|
496 end Exp;
|
|
|
497
|
|
|
498 ---------
|
|
|
499 -- Log --
|
|
|
500 ---------
|
|
|
501
|
|
|
502 function Log (X : Complex) return Complex is
|
|
|
503 ReX : Real'Base;
|
|
|
504 ImX : Real'Base;
|
|
|
505 Z : Complex;
|
|
|
506
|
|
|
507 begin
|
|
|
508 if Re (X) = 0.0 and then Im (X) = 0.0 then
|
|
|
509 raise Constraint_Error;
|
|
|
510
|
|
|
511 elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
|
|
|
512 and then abs Im (X) < Root_Root_Epsilon
|
|
|
513 then
|
|
|
514 Z := X;
|
|
|
515 Set_Re (Z, Re (Z) - 1.0);
|
|
|
516
|
|
|
517 return (1.0 - (1.0 / 2.0 -
|
|
|
518 (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
|
|
|
519 end if;
|
|
|
520
|
|
|
521 begin
|
|
|
522 ReX := Log (Modulus (X));
|
|
|
523
|
|
|
524 exception
|
|
|
525 when Constraint_Error =>
|
|
|
526 ReX := Log (Modulus (X / 2.0)) - Log_Two;
|
|
|
527 end;
|
|
|
528
|
|
|
529 ImX := Arctan (Im (X), Re (X));
|
|
|
530
|
|
|
531 if ImX > PI then
|
|
|
532 ImX := ImX - 2.0 * PI;
|
|
|
533 end if;
|
|
|
534
|
|
|
535 return Compose_From_Cartesian (ReX, ImX);
|
|
|
536 end Log;
|
|
|
537
|
|
|
538 ---------
|
|
|
539 -- Sin --
|
|
|
540 ---------
|
|
|
541
|
|
|
542 function Sin (X : Complex) return Complex is
|
|
|
543 begin
|
|
|
544 if abs Re (X) < Square_Root_Epsilon
|
|
|
545 and then
|
|
|
546 abs Im (X) < Square_Root_Epsilon
|
|
|
547 then
|
|
|
548 return X;
|
|
|
549 end if;
|
|
|
550
|
|
|
551 return
|
|
|
552 Compose_From_Cartesian
|
|
|
553 (Sin (Re (X)) * Cosh (Im (X)),
|
|
|
554 Cos (Re (X)) * Sinh (Im (X)));
|
|
|
555 end Sin;
|
|
|
556
|
|
|
557 ----------
|
|
|
558 -- Sinh --
|
|
|
559 ----------
|
|
|
560
|
|
|
561 function Sinh (X : Complex) return Complex is
|
|
|
562 begin
|
|
|
563 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
564 abs Im (X) < Square_Root_Epsilon
|
|
|
565 then
|
|
|
566 return X;
|
|
|
567
|
|
|
568 else
|
|
|
569 return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
|
|
|
570 Cosh (Re (X)) * Sin (Im (X)));
|
|
|
571 end if;
|
|
|
572 end Sinh;
|
|
|
573
|
|
|
574 ----------
|
|
|
575 -- Sqrt --
|
|
|
576 ----------
|
|
|
577
|
|
|
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);
|
|
|
583 R : Real'Base;
|
|
|
584 R_X : Real'Base;
|
|
|
585 R_Y : Real'Base;
|
|
|
586
|
|
|
587 begin
|
|
|
588 -- Deal with pure real case, see (RM G.1.2(39))
|
|
|
589
|
|
|
590 if ImX = 0.0 then
|
|
|
591 if ReX > 0.0 then
|
|
|
592 return
|
|
|
593 Compose_From_Cartesian
|
|
|
594 (Sqrt (ReX), 0.0);
|
|
|
595
|
|
|
596 elsif ReX = 0.0 then
|
|
|
597 return X;
|
|
|
598
|
|
|
599 else
|
|
|
600 return
|
|
|
601 Compose_From_Cartesian
|
|
|
602 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
|
|
|
603 end if;
|
|
|
604
|
|
|
605 elsif ReX = 0.0 then
|
|
|
606 R_X := Sqrt (YR / 2.0);
|
|
|
607
|
|
|
608 if ImX > 0.0 then
|
|
|
609 return Compose_From_Cartesian (R_X, R_X);
|
|
|
610 else
|
|
|
611 return Compose_From_Cartesian (R_X, -R_X);
|
|
|
612 end if;
|
|
|
613
|
|
|
614 else
|
|
|
615 R := Sqrt (XR ** 2 + YR ** 2);
|
|
|
616
|
|
|
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.
|
|
|
620
|
|
|
621 if R > Real'Base'Last then
|
|
|
622 raise Constraint_Error;
|
|
|
623 end if;
|
|
|
624
|
|
|
625 -- We are solving the system
|
|
|
626
|
|
|
627 -- XR = R_X ** 2 - Y_R ** 2 (1)
|
|
|
628 -- YR = 2.0 * R_X * R_Y (2)
|
|
|
629 --
|
|
|
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
|
|
|
633 -- other.
|
|
|
634
|
|
|
635 if ReX < 0.0 then
|
|
|
636 R_Y := Sqrt (0.5 * (R - ReX));
|
|
|
637 R_X := YR / (2.0 * R_Y);
|
|
|
638
|
|
|
639 else
|
|
|
640 R_X := Sqrt (0.5 * (R + ReX));
|
|
|
641 R_Y := YR / (2.0 * R_X);
|
|
|
642 end if;
|
|
|
643 end if;
|
|
|
644
|
|
|
645 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
|
|
|
646 R_Y := -R_Y;
|
|
|
647 end if;
|
|
|
648 return Compose_From_Cartesian (R_X, R_Y);
|
|
|
649
|
|
|
650 exception
|
|
|
651 when Constraint_Error =>
|
|
|
652
|
|
|
653 -- Rescale and try again
|
|
|
654
|
|
|
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));
|
|
|
658
|
|
|
659 if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
|
|
|
660 R_Y := -R_Y;
|
|
|
661 end if;
|
|
|
662
|
|
|
663 return Compose_From_Cartesian (R_X, R_Y);
|
|
|
664 end Sqrt;
|
|
|
665
|
|
|
666 ---------
|
|
|
667 -- Tan --
|
|
|
668 ---------
|
|
|
669
|
|
|
670 function Tan (X : Complex) return Complex is
|
|
|
671 begin
|
|
|
672 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
673 abs Im (X) < Square_Root_Epsilon
|
|
|
674 then
|
|
|
675 return X;
|
|
|
676
|
|
|
677 elsif Im (X) > Log_Inverse_Epsilon_2 then
|
|
|
678 return Complex_I;
|
|
|
679
|
|
|
680 elsif Im (X) < -Log_Inverse_Epsilon_2 then
|
|
|
681 return -Complex_I;
|
|
|
682
|
|
|
683 else
|
|
|
684 return Sin (X) / Cos (X);
|
|
|
685 end if;
|
|
|
686 end Tan;
|
|
|
687
|
|
|
688 ----------
|
|
|
689 -- Tanh --
|
|
|
690 ----------
|
|
|
691
|
|
|
692 function Tanh (X : Complex) return Complex is
|
|
|
693 begin
|
|
|
694 if abs Re (X) < Square_Root_Epsilon and then
|
|
|
695 abs Im (X) < Square_Root_Epsilon
|
|
|
696 then
|
|
|
697 return X;
|
|
|
698
|
|
|
699 elsif Re (X) > Log_Inverse_Epsilon_2 then
|
|
|
700 return Complex_One;
|
|
|
701
|
|
|
702 elsif Re (X) < -Log_Inverse_Epsilon_2 then
|
|
|
703 return -Complex_One;
|
|
|
704
|
|
|
705 else
|
|
|
706 return Sinh (X) / Cosh (X);
|
|
|
707 end if;
|
|
|
708 end Tanh;
|
|
|
709
|
|
|
710 end Ada.Numerics.Generic_Complex_Elementary_Functions;
|
