1------------------------------------------------------------------------------ 2-- -- 3-- GNAT RUNTIME 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-2001 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 2, 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. See the GNU General Public License -- 18-- for more details. You should have received a copy of the GNU General -- 19-- Public License distributed with GNAT; see file COPYING. If not, write -- 20-- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, -- 21-- MA 02111-1307, USA. -- 22-- -- 23-- As a special exception, if other files instantiate generics from this -- 24-- unit, or you link this unit with other files to produce an executable, -- 25-- this unit does not by itself cause the resulting executable to be -- 26-- covered by the GNU General Public License. This exception does not -- 27-- however invalidate any other reasons why the executable file might be -- 28-- covered by the GNU Public License. -- 29-- -- 30-- GNAT was originally developed by the GNAT team at New York University. -- 31-- Extensive contributions were provided by Ada Core Technologies Inc. -- 32-- -- 33------------------------------------------------------------------------------ 34 35-- File a-numaux.adb <- 86numaux.adb 36 37-- This version of Numerics.Aux is for the IEEE Double Extended floating 38-- point format on x86. 39 40with System.Machine_Code; use System.Machine_Code; 41 42package body Ada.Numerics.Aux is 43 44 NL : constant String := ASCII.LF & ASCII.HT; 45 46 type FPU_Stack_Pointer is range 0 .. 7; 47 for FPU_Stack_Pointer'Size use 3; 48 49 type FPU_Status_Word is record 50 B : Boolean; -- FPU Busy (for 8087 compatibility only) 51 ES : Boolean; -- Error Summary Status 52 SF : Boolean; -- Stack Fault 53 54 Top : FPU_Stack_Pointer; 55 56 -- Condition Code Flags 57 58 -- C2 is set by FPREM and FPREM1 to indicate incomplete reduction. 59 -- In case of successfull recorction, C0, C3 and C1 are set to the 60 -- three least significant bits of the result (resp. Q2, Q1 and Q0). 61 62 -- C2 is used by FPTAN, FSIN, FCOS, and FSINCOS to indicate that 63 -- that source operand is beyond the allowable range of 64 -- -2.0**63 .. 2.0**63. 65 66 C3 : Boolean; 67 C2 : Boolean; 68 C1 : Boolean; 69 C0 : Boolean; 70 71 -- Exception Flags 72 73 PE : Boolean; -- Precision 74 UE : Boolean; -- Underflow 75 OE : Boolean; -- Overflow 76 ZE : Boolean; -- Zero Divide 77 DE : Boolean; -- Denormalized Operand 78 IE : Boolean; -- Invalid Operation 79 end record; 80 81 for FPU_Status_Word use record 82 B at 0 range 15 .. 15; 83 C3 at 0 range 14 .. 14; 84 Top at 0 range 11 .. 13; 85 C2 at 0 range 10 .. 10; 86 C1 at 0 range 9 .. 9; 87 C0 at 0 range 8 .. 8; 88 ES at 0 range 7 .. 7; 89 SF at 0 range 6 .. 6; 90 PE at 0 range 5 .. 5; 91 UE at 0 range 4 .. 4; 92 OE at 0 range 3 .. 3; 93 ZE at 0 range 2 .. 2; 94 DE at 0 range 1 .. 1; 95 IE at 0 range 0 .. 0; 96 end record; 97 98 for FPU_Status_Word'Size use 16; 99 100 ----------------------- 101 -- Local subprograms -- 102 ----------------------- 103 104 function Is_Nan (X : Double) return Boolean; 105 -- Return True iff X is a IEEE NaN value 106 107 function Logarithmic_Pow (X, Y : Double) return Double; 108 -- Implementation of X**Y using Exp and Log functions (binary base) 109 -- to calculate the exponentiation. This is used by Pow for values 110 -- for values of Y in the open interval (-0.25, 0.25) 111 112 function Reduce (X : Double) return Double; 113 -- Implement partial reduction of X by Pi in the x86. 114 115 -- Note that for the Sin, Cos and Tan functions completely accurate 116 -- reduction of the argument is done for arguments in the range of 117 -- -2.0**63 .. 2.0**63, using a 66-bit approximation of Pi. 118 119 pragma Inline (Is_Nan); 120 pragma Inline (Reduce); 121 122 --------------------------------- 123 -- Basic Elementary Functions -- 124 --------------------------------- 125 126 -- This section implements a few elementary functions that are 127 -- used to build the more complex ones. This ordering enables 128 -- better inlining. 129 130 ---------- 131 -- Atan -- 132 ---------- 133 134 function Atan (X : Double) return Double is 135 Result : Double; 136 137 begin 138 Asm (Template => 139 "fld1" & NL 140 & "fpatan", 141 Outputs => Double'Asm_Output ("=t", Result), 142 Inputs => Double'Asm_Input ("0", X)); 143 144 -- The result value is NaN iff input was invalid 145 146 if not (Result = Result) then 147 raise Argument_Error; 148 end if; 149 150 return Result; 151 end Atan; 152 153 --------- 154 -- Exp -- 155 --------- 156 157 function Exp (X : Double) return Double is 158 Result : Double; 159 begin 160 Asm (Template => 161 "fldl2e " & NL 162 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E) 163 & "fld %%st(0) " & NL 164 & "frndint " & NL -- Integer (X * Log2 (E)) 165 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E)) 166 & "fxch " & NL 167 & "f2xm1 " & NL -- 2**(...) - 1 168 & "fld1 " & NL 169 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E))) 170 & "fscale " & NL -- E ** X 171 & "fstp %%st(1) ", 172 Outputs => Double'Asm_Output ("=t", Result), 173 Inputs => Double'Asm_Input ("0", X)); 174 return Result; 175 end Exp; 176 177 ------------ 178 -- Is_Nan -- 179 ------------ 180 181 function Is_Nan (X : Double) return Boolean is 182 begin 183 -- The IEEE NaN values are the only ones that do not equal themselves 184 185 return not (X = X); 186 end Is_Nan; 187 188 --------- 189 -- Log -- 190 --------- 191 192 function Log (X : Double) return Double is 193 Result : Double; 194 195 begin 196 Asm (Template => 197 "fldln2 " & NL 198 & "fxch " & NL 199 & "fyl2x " & NL, 200 Outputs => Double'Asm_Output ("=t", Result), 201 Inputs => Double'Asm_Input ("0", X)); 202 return Result; 203 end Log; 204 205 ------------ 206 -- Reduce -- 207 ------------ 208 209 function Reduce (X : Double) return Double is 210 Result : Double; 211 begin 212 Asm 213 (Template => 214 -- Partial argument reduction 215 "fldpi " & NL 216 & "fadd %%st(0), %%st" & NL 217 & "fxch %%st(1) " & NL 218 & "fprem1 " & NL 219 & "fstp %%st(1) ", 220 Outputs => Double'Asm_Output ("=t", Result), 221 Inputs => Double'Asm_Input ("0", X)); 222 return Result; 223 end Reduce; 224 225 ---------- 226 -- Sqrt -- 227 ---------- 228 229 function Sqrt (X : Double) return Double is 230 Result : Double; 231 232 begin 233 if X < 0.0 then 234 raise Argument_Error; 235 end if; 236 237 Asm (Template => "fsqrt", 238 Outputs => Double'Asm_Output ("=t", Result), 239 Inputs => Double'Asm_Input ("0", X)); 240 241 return Result; 242 end Sqrt; 243 244 --------------------------------- 245 -- Other Elementary Functions -- 246 --------------------------------- 247 248 -- These are built using the previously implemented basic functions 249 250 ---------- 251 -- Acos -- 252 ---------- 253 254 function Acos (X : Double) return Double is 255 Result : Double; 256 begin 257 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X))); 258 259 -- The result value is NaN iff input was invalid 260 261 if Is_Nan (Result) then 262 raise Argument_Error; 263 end if; 264 265 return Result; 266 end Acos; 267 268 ---------- 269 -- Asin -- 270 ---------- 271 272 function Asin (X : Double) return Double is 273 Result : Double; 274 begin 275 276 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X))); 277 278 -- The result value is NaN iff input was invalid 279 280 if Is_Nan (Result) then 281 raise Argument_Error; 282 end if; 283 284 return Result; 285 end Asin; 286 287 --------- 288 -- Cos -- 289 --------- 290 291 function Cos (X : Double) return Double is 292 Reduced_X : Double := X; 293 Result : Double; 294 Status : FPU_Status_Word; 295 296 begin 297 298 loop 299 Asm 300 (Template => 301 "fcos " & NL 302 & "xorl %%eax, %%eax " & NL 303 & "fnstsw %%ax ", 304 Outputs => (Double'Asm_Output ("=t", Result), 305 FPU_Status_Word'Asm_Output ("=a", Status)), 306 Inputs => Double'Asm_Input ("0", Reduced_X)); 307 308 exit when not Status.C2; 309 310 -- Original argument was not in range and the result 311 -- is the unmodified argument. 312 313 Reduced_X := Reduce (Result); 314 end loop; 315 316 return Result; 317 end Cos; 318 319 --------------------- 320 -- Logarithmic_Pow -- 321 --------------------- 322 323 function Logarithmic_Pow (X, Y : Double) return Double is 324 Result : Double; 325 326 begin 327 Asm (Template => "" -- X : Y 328 & "fyl2x " & NL -- Y * Log2 (X) 329 & "fst %%st(1) " & NL -- Y * Log2 (X) : Y * Log2 (X) 330 & "frndint " & NL -- Int (...) : Y * Log2 (X) 331 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...) 332 & "fxch " & NL -- Fract (...) : Int (...) 333 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...) 334 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...) 335 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...) 336 & "fscale " & NL -- 2**(Fract (...) + Int (...)) 337 & "fstp %%st(1) ", 338 Outputs => Double'Asm_Output ("=t", Result), 339 Inputs => 340 (Double'Asm_Input ("0", X), 341 Double'Asm_Input ("u", Y))); 342 343 return Result; 344 end Logarithmic_Pow; 345 346 --------- 347 -- Pow -- 348 --------- 349 350 function Pow (X, Y : Double) return Double is 351 type Mantissa_Type is mod 2**Double'Machine_Mantissa; 352 -- Modular type that can hold all bits of the mantissa of Double 353 354 -- For negative exponents, a division is done 355 -- at the end of the processing. 356 357 Negative_Y : constant Boolean := Y < 0.0; 358 Abs_Y : constant Double := abs Y; 359 360 -- During this function the following invariant is kept: 361 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor 362 363 Base : Double := X; 364 365 Exp_High : Double := Double'Floor (Abs_Y); 366 Exp_Mid : Double; 367 Exp_Low : Double; 368 Exp_Int : Mantissa_Type; 369 370 Factor : Double := 1.0; 371 372 begin 373 -- Select algorithm for calculating Pow: 374 -- integer cases fall through 375 376 if Exp_High >= 2.0**Double'Machine_Mantissa then 377 378 -- In case of Y that is IEEE infinity, just raise constraint error 379 380 if Exp_High > Double'Safe_Last then 381 raise Constraint_Error; 382 end if; 383 384 -- Large values of Y are even integers and will stay integer 385 -- after division by two. 386 387 loop 388 -- Exp_Mid and Exp_Low are zero, so 389 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2) 390 391 Exp_High := Exp_High / 2.0; 392 Base := Base * Base; 393 exit when Exp_High < 2.0**Double'Machine_Mantissa; 394 end loop; 395 396 elsif Exp_High /= Abs_Y then 397 Exp_Low := Abs_Y - Exp_High; 398 399 Factor := 1.0; 400 401 if Exp_Low /= 0.0 then 402 403 -- Exp_Low now is in interval (0.0, 1.0) 404 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0; 405 406 Exp_Mid := 0.0; 407 Exp_Low := Exp_Low - Exp_Mid; 408 409 if Exp_Low >= 0.5 then 410 Factor := Sqrt (X); 411 Exp_Low := Exp_Low - 0.5; -- exact 412 413 if Exp_Low >= 0.25 then 414 Factor := Factor * Sqrt (Factor); 415 Exp_Low := Exp_Low - 0.25; -- exact 416 end if; 417 418 elsif Exp_Low >= 0.25 then 419 Factor := Sqrt (Sqrt (X)); 420 Exp_Low := Exp_Low - 0.25; -- exact 421 end if; 422 423 -- Exp_Low now is in interval (0.0, 0.25) 424 425 -- This means it is safe to call Logarithmic_Pow 426 -- for the remaining part. 427 428 Factor := Factor * Logarithmic_Pow (X, Exp_Low); 429 end if; 430 431 elsif X = 0.0 then 432 return 0.0; 433 end if; 434 435 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa 436 437 Exp_Int := Mantissa_Type (Exp_High); 438 439 -- Standard way for processing integer powers > 0 440 441 while Exp_Int > 1 loop 442 if (Exp_Int and 1) = 1 then 443 444 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0 445 446 Factor := Factor * Base; 447 end if; 448 449 -- Exp_Int is even and Exp_Int > 0, so 450 -- Base**Y = (Base**2)**(Exp_Int / 2) 451 452 Base := Base * Base; 453 Exp_Int := Exp_Int / 2; 454 end loop; 455 456 -- Exp_Int = 1 or Exp_Int = 0 457 458 if Exp_Int = 1 then 459 Factor := Base * Factor; 460 end if; 461 462 if Negative_Y then 463 Factor := 1.0 / Factor; 464 end if; 465 466 return Factor; 467 end Pow; 468 469 --------- 470 -- Sin -- 471 --------- 472 473 function Sin (X : Double) return Double is 474 Reduced_X : Double := X; 475 Result : Double; 476 Status : FPU_Status_Word; 477 478 begin 479 480 loop 481 Asm 482 (Template => 483 "fsin " & NL 484 & "xorl %%eax, %%eax " & NL 485 & "fnstsw %%ax ", 486 Outputs => (Double'Asm_Output ("=t", Result), 487 FPU_Status_Word'Asm_Output ("=a", Status)), 488 Inputs => Double'Asm_Input ("0", Reduced_X)); 489 490 exit when not Status.C2; 491 492 -- Original argument was not in range and the result 493 -- is the unmodified argument. 494 495 Reduced_X := Reduce (Result); 496 end loop; 497 498 return Result; 499 end Sin; 500 501 --------- 502 -- Tan -- 503 --------- 504 505 function Tan (X : Double) return Double is 506 Reduced_X : Double := X; 507 Result : Double; 508 Status : FPU_Status_Word; 509 510 begin 511 512 loop 513 Asm 514 (Template => 515 "fptan " & NL 516 & "xorl %%eax, %%eax " & NL 517 & "fnstsw %%ax " & NL 518 & "ffree %%st(0) " & NL 519 & "fincstp ", 520 521 Outputs => (Double'Asm_Output ("=t", Result), 522 FPU_Status_Word'Asm_Output ("=a", Status)), 523 Inputs => Double'Asm_Input ("0", Reduced_X)); 524 525 exit when not Status.C2; 526 527 -- Original argument was not in range and the result 528 -- is the unmodified argument. 529 530 Reduced_X := Reduce (Result); 531 end loop; 532 533 return Result; 534 end Tan; 535 536 ---------- 537 -- Sinh -- 538 ---------- 539 540 function Sinh (X : Double) return Double is 541 begin 542 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0 543 544 if abs X < 25.0 then 545 return (Exp (X) - Exp (-X)) / 2.0; 546 547 else 548 return Exp (X) / 2.0; 549 end if; 550 551 end Sinh; 552 553 ---------- 554 -- Cosh -- 555 ---------- 556 557 function Cosh (X : Double) return Double is 558 begin 559 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0 560 561 if abs X < 22.0 then 562 return (Exp (X) + Exp (-X)) / 2.0; 563 564 else 565 return Exp (X) / 2.0; 566 end if; 567 568 end Cosh; 569 570 ---------- 571 -- Tanh -- 572 ---------- 573 574 function Tanh (X : Double) return Double is 575 begin 576 -- Return the Hyperbolic Tangent of x 577 -- 578 -- x -x 579 -- e - e Sinh (X) 580 -- Tanh (X) is defined to be ----------- = -------- 581 -- x -x Cosh (X) 582 -- e + e 583 584 if abs X > 23.0 then 585 return Double'Copy_Sign (1.0, X); 586 end if; 587 588 return 1.0 / (1.0 + Exp (-2.0 * X)) - 1.0 / (1.0 + Exp (2.0 * X)); 589 590 end Tanh; 591 592end Ada.Numerics.Aux; 593