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