1------------------------------------------------------------------------------ 2-- -- 3-- GNAT RUN-TIME COMPONENTS -- 4-- -- 5-- A D A . N U M E R I C S . A U X -- 6-- -- 7-- B o d y -- 8-- (Machine Version for x86) -- 9-- -- 10-- Copyright (C) 1998-2009, Free Software Foundation, Inc. -- 11-- -- 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 3, 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. -- 18-- -- 19-- As a special exception under Section 7 of GPL version 3, you are granted -- 20-- additional permissions described in the GCC Runtime Library Exception, -- 21-- version 3.1, as published by the Free Software Foundation. -- 22-- -- 23-- You should have received a copy of the GNU General Public License and -- 24-- a copy of the GCC Runtime Library Exception along with this program; -- 25-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- 26-- <http://www.gnu.org/licenses/>. -- 27-- -- 28-- GNAT was originally developed by the GNAT team at New York University. -- 29-- Extensive contributions were provided by Ada Core Technologies Inc. -- 30-- -- 31------------------------------------------------------------------------------ 32 33-- File a-numaux.adb <- 86numaux.adb 34 35-- This version of Numerics.Aux is for the IEEE Double Extended floating 36-- point format on x86. 37 38with System.Machine_Code; use System.Machine_Code; 39 40package body Ada.Numerics.Aux is 41 42 NL : constant String := ASCII.LF & ASCII.HT; 43 44 ----------------------- 45 -- Local subprograms -- 46 ----------------------- 47 48 function Is_Nan (X : Double) return Boolean; 49 -- Return True iff X is a IEEE NaN value 50 51 function Logarithmic_Pow (X, Y : Double) return Double; 52 -- Implementation of X**Y using Exp and Log functions (binary base) 53 -- to calculate the exponentiation. This is used by Pow for values 54 -- for values of Y in the open interval (-0.25, 0.25) 55 56 procedure Reduce (X : in out Double; Q : out Natural); 57 -- Implements reduction of X by Pi/2. Q is the quadrant of the final 58 -- result in the range 0 .. 3. The absolute value of X is at most Pi. 59 60 pragma Inline (Is_Nan); 61 pragma Inline (Reduce); 62 63 -------------------------------- 64 -- Basic Elementary Functions -- 65 -------------------------------- 66 67 -- This section implements a few elementary functions that are used to 68 -- build the more complex ones. This ordering enables better inlining. 69 70 ---------- 71 -- Atan -- 72 ---------- 73 74 function Atan (X : Double) return Double is 75 Result : Double; 76 77 begin 78 Asm (Template => 79 "fld1" & NL 80 & "fpatan", 81 Outputs => Double'Asm_Output ("=t", Result), 82 Inputs => Double'Asm_Input ("0", X)); 83 84 -- The result value is NaN iff input was invalid 85 86 if not (Result = Result) then 87 raise Argument_Error; 88 end if; 89 90 return Result; 91 end Atan; 92 93 --------- 94 -- Exp -- 95 --------- 96 97 function Exp (X : Double) return Double is 98 Result : Double; 99 begin 100 Asm (Template => 101 "fldl2e " & NL 102 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E) 103 & "fld %%st(0) " & NL 104 & "frndint " & NL -- Integer (X * Log2 (E)) 105 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E)) 106 & "fxch " & NL 107 & "f2xm1 " & NL -- 2**(...) - 1 108 & "fld1 " & NL 109 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E))) 110 & "fscale " & NL -- E ** X 111 & "fstp %%st(1) ", 112 Outputs => Double'Asm_Output ("=t", Result), 113 Inputs => Double'Asm_Input ("0", X)); 114 return Result; 115 end Exp; 116 117 ------------ 118 -- Is_Nan -- 119 ------------ 120 121 function Is_Nan (X : Double) return Boolean is 122 begin 123 -- The IEEE NaN values are the only ones that do not equal themselves 124 125 return not (X = X); 126 end Is_Nan; 127 128 --------- 129 -- Log -- 130 --------- 131 132 function Log (X : Double) return Double is 133 Result : Double; 134 135 begin 136 Asm (Template => 137 "fldln2 " & NL 138 & "fxch " & NL 139 & "fyl2x " & NL, 140 Outputs => Double'Asm_Output ("=t", Result), 141 Inputs => Double'Asm_Input ("0", X)); 142 return Result; 143 end Log; 144 145 ------------ 146 -- Reduce -- 147 ------------ 148 149 procedure Reduce (X : in out Double; Q : out Natural) is 150 Half_Pi : constant := Pi / 2.0; 151 Two_Over_Pi : constant := 2.0 / Pi; 152 153 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size); 154 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant 155 P1 : constant Double := Double'Leading_Part (Half_Pi, HM); 156 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM); 157 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM); 158 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM); 159 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3 160 - P4, HM); 161 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5); 162 K : Double := X * Two_Over_Pi; 163 begin 164 -- For X < 2.0**32, all products below are computed exactly. 165 -- Due to cancellation effects all subtractions are exact as well. 166 -- As no double extended floating-point number has more than 75 167 -- zeros after the binary point, the result will be the correctly 168 -- rounded result of X - K * (Pi / 2.0). 169 170 while abs K >= 2.0**HM loop 171 K := K * M - (K * M - K); 172 X := (((((X - K * P1) - K * P2) - K * P3) 173 - K * P4) - K * P5) - K * P6; 174 K := X * Two_Over_Pi; 175 end loop; 176 177 if K /= K then 178 179 -- K is not a number, because X was not finite 180 181 raise Constraint_Error; 182 end if; 183 184 K := Double'Rounding (K); 185 Q := Integer (K) mod 4; 186 X := (((((X - K * P1) - K * P2) - K * P3) 187 - K * P4) - K * P5) - K * P6; 188 end Reduce; 189 190 ---------- 191 -- Sqrt -- 192 ---------- 193 194 function Sqrt (X : Double) return Double is 195 Result : Double; 196 197 begin 198 if X < 0.0 then 199 raise Argument_Error; 200 end if; 201 202 Asm (Template => "fsqrt", 203 Outputs => Double'Asm_Output ("=t", Result), 204 Inputs => Double'Asm_Input ("0", X)); 205 206 return Result; 207 end Sqrt; 208 209 -------------------------------- 210 -- Other Elementary Functions -- 211 -------------------------------- 212 213 -- These are built using the previously implemented basic functions 214 215 ---------- 216 -- Acos -- 217 ---------- 218 219 function Acos (X : Double) return Double is 220 Result : Double; 221 222 begin 223 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X))); 224 225 -- The result value is NaN iff input was invalid 226 227 if Is_Nan (Result) then 228 raise Argument_Error; 229 end if; 230 231 return Result; 232 end Acos; 233 234 ---------- 235 -- Asin -- 236 ---------- 237 238 function Asin (X : Double) return Double is 239 Result : Double; 240 241 begin 242 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X))); 243 244 -- The result value is NaN iff input was invalid 245 246 if Is_Nan (Result) then 247 raise Argument_Error; 248 end if; 249 250 return Result; 251 end Asin; 252 253 --------- 254 -- Cos -- 255 --------- 256 257 function Cos (X : Double) return Double is 258 Reduced_X : Double := abs X; 259 Result : Double; 260 Quadrant : Natural range 0 .. 3; 261 262 begin 263 if Reduced_X > Pi / 4.0 then 264 Reduce (Reduced_X, Quadrant); 265 266 case Quadrant is 267 when 0 => 268 Asm (Template => "fcos", 269 Outputs => Double'Asm_Output ("=t", Result), 270 Inputs => Double'Asm_Input ("0", Reduced_X)); 271 when 1 => 272 Asm (Template => "fsin", 273 Outputs => Double'Asm_Output ("=t", Result), 274 Inputs => Double'Asm_Input ("0", -Reduced_X)); 275 when 2 => 276 Asm (Template => "fcos ; fchs", 277 Outputs => Double'Asm_Output ("=t", Result), 278 Inputs => Double'Asm_Input ("0", Reduced_X)); 279 when 3 => 280 Asm (Template => "fsin", 281 Outputs => Double'Asm_Output ("=t", Result), 282 Inputs => Double'Asm_Input ("0", Reduced_X)); 283 end case; 284 285 else 286 Asm (Template => "fcos", 287 Outputs => Double'Asm_Output ("=t", Result), 288 Inputs => Double'Asm_Input ("0", Reduced_X)); 289 end if; 290 291 return Result; 292 end Cos; 293 294 --------------------- 295 -- Logarithmic_Pow -- 296 --------------------- 297 298 function Logarithmic_Pow (X, Y : Double) return Double is 299 Result : Double; 300 begin 301 Asm (Template => "" -- X : Y 302 & "fyl2x " & NL -- Y * Log2 (X) 303 & "fld %%st(0) " & NL -- Y * Log2 (X) : Y * Log2 (X) 304 & "frndint " & NL -- Int (...) : Y * Log2 (X) 305 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...) 306 & "fxch " & NL -- Fract (...) : Int (...) 307 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...) 308 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...) 309 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...) 310 & "fscale ", -- 2**(Fract (...) + Int (...)) 311 Outputs => Double'Asm_Output ("=t", Result), 312 Inputs => 313 (Double'Asm_Input ("0", X), 314 Double'Asm_Input ("u", Y))); 315 return Result; 316 end Logarithmic_Pow; 317 318 --------- 319 -- Pow -- 320 --------- 321 322 function Pow (X, Y : Double) return Double is 323 type Mantissa_Type is mod 2**Double'Machine_Mantissa; 324 -- Modular type that can hold all bits of the mantissa of Double 325 326 -- For negative exponents, do divide at the end of the processing 327 328 Negative_Y : constant Boolean := Y < 0.0; 329 Abs_Y : constant Double := abs Y; 330 331 -- During this function the following invariant is kept: 332 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor 333 334 Base : Double := X; 335 336 Exp_High : Double := Double'Floor (Abs_Y); 337 Exp_Mid : Double; 338 Exp_Low : Double; 339 Exp_Int : Mantissa_Type; 340 341 Factor : Double := 1.0; 342 343 begin 344 -- Select algorithm for calculating Pow (integer cases fall through) 345 346 if Exp_High >= 2.0**Double'Machine_Mantissa then 347 348 -- In case of Y that is IEEE infinity, just raise constraint error 349 350 if Exp_High > Double'Safe_Last then 351 raise Constraint_Error; 352 end if; 353 354 -- Large values of Y are even integers and will stay integer 355 -- after division by two. 356 357 loop 358 -- Exp_Mid and Exp_Low are zero, so 359 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2) 360 361 Exp_High := Exp_High / 2.0; 362 Base := Base * Base; 363 exit when Exp_High < 2.0**Double'Machine_Mantissa; 364 end loop; 365 366 elsif Exp_High /= Abs_Y then 367 Exp_Low := Abs_Y - Exp_High; 368 Factor := 1.0; 369 370 if Exp_Low /= 0.0 then 371 372 -- Exp_Low now is in interval (0.0, 1.0) 373 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0; 374 375 Exp_Mid := 0.0; 376 Exp_Low := Exp_Low - Exp_Mid; 377 378 if Exp_Low >= 0.5 then 379 Factor := Sqrt (X); 380 Exp_Low := Exp_Low - 0.5; -- exact 381 382 if Exp_Low >= 0.25 then 383 Factor := Factor * Sqrt (Factor); 384 Exp_Low := Exp_Low - 0.25; -- exact 385 end if; 386 387 elsif Exp_Low >= 0.25 then 388 Factor := Sqrt (Sqrt (X)); 389 Exp_Low := Exp_Low - 0.25; -- exact 390 end if; 391 392 -- Exp_Low now is in interval (0.0, 0.25) 393 394 -- This means it is safe to call Logarithmic_Pow 395 -- for the remaining part. 396 397 Factor := Factor * Logarithmic_Pow (X, Exp_Low); 398 end if; 399 400 elsif X = 0.0 then 401 return 0.0; 402 end if; 403 404 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa 405 406 Exp_Int := Mantissa_Type (Exp_High); 407 408 -- Standard way for processing integer powers > 0 409 410 while Exp_Int > 1 loop 411 if (Exp_Int and 1) = 1 then 412 413 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0 414 415 Factor := Factor * Base; 416 end if; 417 418 -- Exp_Int is even and Exp_Int > 0, so 419 -- Base**Y = (Base**2)**(Exp_Int / 2) 420 421 Base := Base * Base; 422 Exp_Int := Exp_Int / 2; 423 end loop; 424 425 -- Exp_Int = 1 or Exp_Int = 0 426 427 if Exp_Int = 1 then 428 Factor := Base * Factor; 429 end if; 430 431 if Negative_Y then 432 Factor := 1.0 / Factor; 433 end if; 434 435 return Factor; 436 end Pow; 437 438 --------- 439 -- Sin -- 440 --------- 441 442 function Sin (X : Double) return Double is 443 Reduced_X : Double := X; 444 Result : Double; 445 Quadrant : Natural range 0 .. 3; 446 447 begin 448 if abs X > Pi / 4.0 then 449 Reduce (Reduced_X, Quadrant); 450 451 case Quadrant is 452 when 0 => 453 Asm (Template => "fsin", 454 Outputs => Double'Asm_Output ("=t", Result), 455 Inputs => Double'Asm_Input ("0", Reduced_X)); 456 when 1 => 457 Asm (Template => "fcos", 458 Outputs => Double'Asm_Output ("=t", Result), 459 Inputs => Double'Asm_Input ("0", Reduced_X)); 460 when 2 => 461 Asm (Template => "fsin", 462 Outputs => Double'Asm_Output ("=t", Result), 463 Inputs => Double'Asm_Input ("0", -Reduced_X)); 464 when 3 => 465 Asm (Template => "fcos ; fchs", 466 Outputs => Double'Asm_Output ("=t", Result), 467 Inputs => Double'Asm_Input ("0", Reduced_X)); 468 end case; 469 470 else 471 Asm (Template => "fsin", 472 Outputs => Double'Asm_Output ("=t", Result), 473 Inputs => Double'Asm_Input ("0", Reduced_X)); 474 end if; 475 476 return Result; 477 end Sin; 478 479 --------- 480 -- Tan -- 481 --------- 482 483 function Tan (X : Double) return Double is 484 Reduced_X : Double := X; 485 Result : Double; 486 Quadrant : Natural range 0 .. 3; 487 488 begin 489 if abs X > Pi / 4.0 then 490 Reduce (Reduced_X, Quadrant); 491 492 if Quadrant mod 2 = 0 then 493 Asm (Template => "fptan" & NL 494 & "ffree %%st(0)" & NL 495 & "fincstp", 496 Outputs => Double'Asm_Output ("=t", Result), 497 Inputs => Double'Asm_Input ("0", Reduced_X)); 498 else 499 Asm (Template => "fsincos" & NL 500 & "fdivp %%st, %%st(1)" & NL 501 & "fchs", 502 Outputs => Double'Asm_Output ("=t", Result), 503 Inputs => Double'Asm_Input ("0", Reduced_X)); 504 end if; 505 506 else 507 Asm (Template => 508 "fptan " & NL 509 & "ffree %%st(0) " & NL 510 & "fincstp ", 511 Outputs => Double'Asm_Output ("=t", Result), 512 Inputs => Double'Asm_Input ("0", Reduced_X)); 513 end if; 514 515 return Result; 516 end Tan; 517 518 ---------- 519 -- Sinh -- 520 ---------- 521 522 function Sinh (X : Double) return Double is 523 begin 524 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0 525 526 if abs X < 25.0 then 527 return (Exp (X) - Exp (-X)) / 2.0; 528 else 529 return Exp (X) / 2.0; 530 end if; 531 end Sinh; 532 533 ---------- 534 -- Cosh -- 535 ---------- 536 537 function Cosh (X : Double) return Double is 538 begin 539 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0 540 541 if abs X < 22.0 then 542 return (Exp (X) + Exp (-X)) / 2.0; 543 else 544 return Exp (X) / 2.0; 545 end if; 546 end Cosh; 547 548 ---------- 549 -- Tanh -- 550 ---------- 551 552 function Tanh (X : Double) return Double is 553 begin 554 -- Return the Hyperbolic Tangent of x 555 556 -- x -x 557 -- e - e Sinh (X) 558 -- Tanh (X) is defined to be ----------- = -------- 559 -- x -x Cosh (X) 560 -- e + e 561 562 if abs X > 23.0 then 563 return Double'Copy_Sign (1.0, X); 564 end if; 565 566 return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X)); 567 end Tanh; 568 569end Ada.Numerics.Aux; 570