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-2013, 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 32with Ada.Numerics.Generic_Elementary_Functions; 33 34package 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 710end Ada.Numerics.Generic_Complex_Elementary_Functions; 711