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