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-- 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.Aux; use Ada.Numerics.Aux; 33 34package 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 681end Ada.Numerics.Generic_Complex_Types; 682