1 /** 2 * Elementary Mathematical Functions 3 * 4 * Copyright: 5 * Portions Copyright (C) 2001-2005 Digital Mars. 6 * Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH. 7 * All rights reserved. 8 * 9 * License: 10 * Tango Dual License: 3-Clause BSD License / Academic Free License v3.0. 11 * See LICENSE_TANGO.txt for details. 12 * 13 * Authors: Walter Bright, Don Clugston, Sean Kelly 14 * 15 */ 16 17 /** 18 * Macros: 19 * NAN = $(RED NAN) 20 * TEXTNAN = $(RED NAN:$1 ) 21 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span> 22 * GAMMA = Γ 23 * INTEGRAL = ∫ 24 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) 25 * POWER = $1<sup>$2</sup> 26 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>) 27 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG )) 28 * PLUSMN = ± 29 * INFIN = ∞ 30 * PLUSMNINF = ±∞ 31 * PI = π 32 * LT = < 33 * GT = > 34 * SQRT = &radix; 35 * HALF = ½ 36 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> 37 * <caption>Special Values</caption> 38 * $0</table> 39 * SVH = $(TR $(TH $1) $(TH $2)) 40 * SV = $(TR $(TD $1) $(TD $2)) 41 * TABLE_DOMRG = <table border=1 cellpadding=4 cellspacing=0>$0</table> 42 * DOMAIN = $(TR $(TD Domain) $(TD $0)) 43 * RANGE = $(TR $(TD Range) $(TD $0)) 44 */ 45 46 module ocean.math.Math; 47 48 import ocean.meta.types.Qualifiers; 49 50 static import core.stdc.math; 51 import ocean.math.IEEE; 52 import ocean.core.Verify; 53 54 version (unittest) import ocean.core.Test; 55 56 version(TangoNoAsm) { 57 58 } else version(D_InlineAsm_X86) { 59 version = Naked_D_InlineAsm_X86; 60 } else version(D_InlineAsm_X86_64) { 61 version = Naked_D_InlineAsm_X86_64; 62 } 63 64 /* 65 * Constants 66 */ 67 68 static immutable real E = 2.7182818284590452354L; /** e */ // 3.32193 fldl2t 0x1.5BF0A8B1_45769535_5FF5p+1L 69 static immutable real LOG2T = 0x1.a934f0979a3715fcp+1; /** $(SUB log, 2)10 */ // 1.4427 fldl2e 70 static immutable real LOG2E = 0x1.71547652b82fe178p+0; /** $(SUB log, 2)e */ // 0.30103 fldlg2 71 static immutable real LOG2 = 0x1.34413509f79fef32p-2; /** $(SUB log, 10)2 */ 72 static immutable real LOG10E = 0.43429448190325182765; /** $(SUB log, 10)e */ 73 static immutable real LN2 = 0x1.62e42fefa39ef358p-1; /** ln 2 */ // 0.693147 fldln2 74 static immutable real LN10 = 2.30258509299404568402; /** ln 10 */ 75 static immutable real PI = 0x1.921fb54442d1846ap+1; /** $(_PI) */ // 3.14159 fldpi 76 static immutable real PI_2 = 1.57079632679489661923; /** $(PI) / 2 */ 77 static immutable real PI_4 = 0.78539816339744830962; /** $(PI) / 4 */ 78 static immutable real M_1_PI = 0.31830988618379067154; /** 1 / $(PI) */ 79 static immutable real M_2_PI = 0.63661977236758134308; /** 2 / $(PI) */ 80 static immutable real M_2_SQRTPI = 1.12837916709551257390; /** 2 / $(SQRT)$(PI) */ 81 static immutable real SQRT2 = 1.41421356237309504880; /** $(SQRT)2 */ 82 static immutable real SQRT1_2 = 0.70710678118654752440; /** $(SQRT)$(HALF) */ 83 84 //const real SQRTPI = 1.77245385090551602729816748334114518279754945612238L; /** √π */ 85 //const real SQRT2PI = 2.50662827463100050242E0L; /** √(2 π) */ 86 //const real SQRTE = 1.64872127070012814684865078781416357L; /** √(e) */ 87 88 static immutable real MAXLOG = 0x1.62e42fefa39ef358p+13L; /** log(real.max) */ 89 static immutable real MINLOG = -0x1.6436716d5406e6d8p+13L; /** log(real.min*real.epsilon) */ 90 static immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */ 91 92 /* 93 * Primitives 94 */ 95 96 /** 97 * Calculates the absolute value 98 * 99 * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) ) 100 * = hypot(z.re, z.im). 101 */ 102 real abs(real x) 103 { 104 return ocean.math.IEEE.fabs(x); 105 } 106 107 /** ditto */ 108 long abs(long x) 109 { 110 return x>=0 ? x : -x; 111 } 112 113 /** ditto */ 114 int abs(int x) 115 { 116 return x>=0 ? x : -x; 117 } 118 119 /** ditto */ 120 deprecated("Use `std.complex.Complex` instead") 121 real abs(creal z) 122 { 123 return hypot(z.re, z.im); 124 } 125 126 /** ditto */ 127 deprecated("Use `real` instead") 128 real abs(ireal y) 129 { 130 return ocean.math.IEEE.fabs(y.im); 131 } 132 133 unittest 134 { 135 test(isIdentical(0.0L,abs(-0.0L))); 136 test(isNaN(abs(real.nan))); 137 test(abs(-real.infinity) == real.infinity); 138 test(abs(-56) == 56); 139 test(abs(2321312L) == 2321312L); 140 } 141 142 deprecated unittest 143 { 144 test(abs(-3.2Li) == 3.2L); 145 test(abs(71.6Li) == 71.6L); 146 test(abs(-1.0L+1.0Li) == sqrt(2.0L)); 147 } 148 149 /** 150 * Complex conjugate 151 * 152 * conj(x + iy) = x - iy 153 * 154 * Note that z * conj(z) = $(POWER z.re, 2) + $(POWER z.im, 2) 155 * is always a real number 156 */ 157 deprecated("Use `std.complex.Complex` instead") 158 creal conj(creal z) 159 { 160 return z.re - z.im*1i; 161 } 162 163 /** ditto */ 164 deprecated("Use `real` instead") 165 ireal conj(ireal y) 166 { 167 return -y; 168 } 169 170 deprecated unittest 171 { 172 test(conj(7 + 3i) == 7-3i); 173 ireal z = -3.2Li; 174 test(conj(z) == -z); 175 } 176 177 private { 178 // Return the type which would be returned by a max or min operation 179 template minmaxtype(T...){ 180 static if(T.length == 1) alias T[0] minmaxtype; 181 else static if(T.length > 2) 182 alias minmaxtype!(minmaxtype!(T[0..2]), T[2..$]) minmaxtype; 183 else alias typeof (T[1].init > T[0].init ? T[1].init : T[0].init) minmaxtype; 184 } 185 } 186 187 unittest 188 { 189 static assert (is(minmaxtype!(int, long) == long)); 190 } 191 192 /** Return the minimum of the supplied arguments. 193 * 194 * Note: If the arguments are floating-point numbers, and at least one is a NaN, 195 * the result is undefined. 196 */ 197 minmaxtype!(T) min(T...)(T arg){ 198 static if(arg.length == 1) return arg[0]; 199 else static if(arg.length == 2) return arg[1] < arg[0] ? arg[1] : arg[0]; 200 static if(arg.length > 2) return min(arg[1] < arg[0] ? arg[1] : arg[0], arg[2..$]); 201 } 202 203 /** Return the maximum of the supplied arguments. 204 * 205 * Note: If the arguments are floating-point numbers, and at least one is a NaN, 206 * the result is undefined. 207 */ 208 minmaxtype!(T) max(T...)(T arg){ 209 static if(arg.length == 1) return arg[0]; 210 else static if(arg.length == 2) return arg[1] > arg[0] ? arg[1] : arg[0]; 211 static if(arg.length > 2) return max(arg[1] > arg[0] ? arg[1] : arg[0], arg[2..$]); 212 } 213 unittest 214 { 215 test(max('e', 'f')=='f'); 216 test(min(3.5, 3.8)==3.5); 217 // check implicit conversion to integer. 218 test(min(3.5, 18)==3.5); 219 220 } 221 222 /** Returns the minimum number of x and y, favouring numbers over NaNs. 223 * 224 * If both x and y are numbers, the minimum is returned. 225 * If both parameters are NaN, either will be returned. 226 * If one parameter is a NaN and the other is a number, the number is 227 * returned (this behaviour is mandated by IEEE 754R, and is useful 228 * for determining the range of a function). 229 */ 230 real minNum(real x, real y) { 231 if (x<=y || isNaN(y)) return x; else return y; 232 } 233 234 /** Returns the maximum number of x and y, favouring numbers over NaNs. 235 * 236 * If both x and y are numbers, the maximum is returned. 237 * If both parameters are NaN, either will be returned. 238 * If one parameter is a NaN and the other is a number, the number is 239 * returned (this behaviour is mandated by IEEE 754-2008, and is useful 240 * for determining the range of a function). 241 */ 242 real maxNum(real x, real y) { 243 if (x>=y || isNaN(y)) return x; else return y; 244 } 245 246 /** Returns the minimum of x and y, favouring NaNs over numbers 247 * 248 * If both x and y are numbers, the minimum is returned. 249 * If both parameters are NaN, either will be returned. 250 * If one parameter is a NaN and the other is a number, the NaN is returned. 251 */ 252 real minNaN(real x, real y) { 253 return (x<=y || isNaN(x))? x : y; 254 } 255 256 /** Returns the maximum of x and y, favouring NaNs over numbers 257 * 258 * If both x and y are numbers, the maximum is returned. 259 * If both parameters are NaN, either will be returned. 260 * If one parameter is a NaN and the other is a number, the NaN is returned. 261 */ 262 real maxNaN(real x, real y) { 263 return (x>=y || isNaN(x))? x : y; 264 } 265 266 unittest 267 { 268 test(maxNum(NaN(0xABC), 56.1L)== 56.1L); 269 test(isIdentical(maxNaN(NaN(1389), 56.1L), NaN(1389))); 270 test(maxNum(28.0, NaN(0xABC))== 28.0); 271 test(minNum(1e12, NaN(0xABC))== 1e12); 272 test(isIdentical(minNaN(1e12, NaN(23454)), NaN(23454))); 273 test(isIdentical(minNum(NaN(489), NaN(23)), NaN(489))); 274 } 275 276 /* 277 * Trig Functions 278 */ 279 280 /*********************************** 281 * Returns cosine of x. x is in radians. 282 * 283 * $(TABLE_SV 284 * $(TR $(TH x) $(TH cos(x)) $(TH invalid?)) 285 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) 286 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes) ) 287 * ) 288 * Bugs: 289 * Results are undefined if |x| >= $(POWER 2,64). 290 */ 291 292 real cos(real x) /* intrinsic */ 293 { 294 version(D_InlineAsm_X86) 295 { 296 asm 297 { 298 fld x; 299 fcos; 300 } 301 } 302 else 303 { 304 return core.stdc.math.cosl(x); 305 } 306 } 307 308 unittest { 309 // NaN payloads 310 test(isIdentical(cos(NaN(314)), NaN(314))); 311 } 312 313 /*********************************** 314 * Returns sine of x. x is in radians. 315 * 316 * $(TABLE_SV 317 * $(TR $(TH x) $(TH sin(x)) $(TH invalid?)) 318 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) 319 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 320 * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes)) 321 * ) 322 * Bugs: 323 * Results are undefined if |x| >= $(POWER 2,64). 324 */ 325 real sin(real x) /* intrinsic */ 326 { 327 version(D_InlineAsm_X86) 328 { 329 asm 330 { 331 fld x; 332 fsin; 333 } 334 } 335 else 336 { 337 return core.stdc.math.sinl(x); 338 } 339 } 340 341 unittest { 342 // NaN payloads 343 test(isIdentical(sin(NaN(314)), NaN(314))); 344 } 345 346 /** 347 * Returns tangent of x. x is in radians. 348 * 349 * $(TABLE_SV 350 * $(TR $(TH x) $(TH tan(x)) $(TH invalid?)) 351 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) 352 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 353 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)) 354 * ) 355 */ 356 real tan(real x) 357 { 358 asm 359 { 360 fld x[EBP] ; // load theta 361 fxam ; // test for oddball values 362 fstsw AX ; 363 sahf ; 364 jc trigerr ; // x is NAN, infinity, or empty 365 // 387's can handle denormals 366 SC18: fptan ; 367 fstp ST(0) ; // dump X, which is always 1 368 fstsw AX ; 369 sahf ; 370 jnp Lret ; // C2 = 1 (x is out of range) 371 372 // Do argument reduction to bring x into range 373 fldpi ; 374 fxch ; 375 SC17: fprem1 ; 376 fstsw AX ; 377 sahf ; 378 jp SC17 ; 379 fstp ST(1) ; // remove pi from stack 380 jmp SC18 ; 381 382 trigerr: 383 jnp Lret ; // if x is NaN, return x. 384 fstp ST(0) ; // dump x, which will be infinity 385 } 386 return NaN(TANGO_NAN.TAN_DOMAIN); 387 Lret: 388 ; 389 } 390 391 unittest 392 { 393 static real[2][] vals = // angle,tan 394 [ 395 [ 0, 0], 396 [ .5, .5463024898], 397 [ 1, 1.557407725], 398 [ 1.5, 14.10141995], 399 [ 2, -2.185039863], 400 [ 2.5,-.7470222972], 401 [ 3, -.1425465431], 402 [ 3.5, .3745856402], 403 [ 4, 1.157821282], 404 [ 4.5, 4.637332055], 405 [ 5, -3.380515006], 406 [ 5.5,-.9955840522], 407 [ 6, -.2910061914], 408 [ 6.5, .2202772003], 409 [ 10, .6483608275], 410 411 // special angles 412 [ PI_4, 1], 413 //[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2. 414 [ 3*PI_4, -1], 415 [ PI, 0], 416 [ 5*PI_4, 1], 417 //[ 3*PI_2, -real.infinity], 418 [ 7*PI_4, -1], 419 [ 2*PI, 0], 420 ]; 421 int i; 422 423 for (i = 0; i < vals.length; i++) 424 { 425 real x = vals[i][0]; 426 real r = vals[i][1]; 427 real t = tan(x); 428 429 //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); 430 if (!isIdentical(r, t)) test(fabs(r-t) <= .0000001); 431 432 x = -x; 433 r = -r; 434 t = tan(x); 435 //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); 436 if (!isIdentical(r, t) && !(isNaN(r) && isNaN(t))) 437 test(fabs(r-t) <= .0000001); 438 } 439 // overflow 440 test(isNaN(tan(real.infinity))); 441 test(isNaN(tan(-real.infinity))); 442 // NaN propagation 443 test(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) )); 444 } 445 446 /***************************************** 447 * Sine, cosine, and arctangent of multiple of π 448 * 449 * Accuracy is preserved for large values of x. 450 */ 451 real cosPi(real x) 452 { 453 return cos((x%2.0)*PI); 454 } 455 456 /** ditto */ 457 real sinPi(real x) 458 { 459 return sin((x%2.0)*PI); 460 } 461 462 /** ditto */ 463 real atanPi(real x) 464 { 465 return PI * atan(x); // BUG: Fix this. 466 } 467 468 unittest { 469 test(isIdentical(sinPi(0.0), 0.0)); 470 test(isIdentical(sinPi(-0.0), -0.0)); 471 test(isIdentical(atanPi(0.0), 0.0)); 472 test(isIdentical(atanPi(-0.0), -0.0)); 473 } 474 475 /*********************************** 476 * sine, complex and imaginary 477 * 478 * sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i 479 * 480 * If both sin(θ) and cos(θ) are required, 481 * it is most efficient to use expi(&theta). 482 */ 483 deprecated("Use `std.complex.Complex` instead") 484 creal sin(creal z) 485 { 486 creal cs = expi(z.re); 487 return cs.im * cosh(z.im) + cs.re * sinh(z.im) * 1i; 488 } 489 490 /** ditto */ 491 deprecated("Use `ireal` instead") 492 ireal sin(ireal y) 493 { 494 return cosh(y.im)*1i; 495 } 496 497 deprecated unittest 498 { 499 test(sin(0.0+0.0i) == 0.0); 500 test(sin(2.0+0.0i) == sin(2.0L) ); 501 } 502 503 /*********************************** 504 * cosine, complex and imaginary 505 * 506 * cos(z) = cos(z.re)*cosh(z.im) + sin(z.re)*sinh(z.im)i 507 */ 508 deprecated("Use `std.complex.Complex` instead") 509 creal cos(creal z) 510 { 511 creal cs = expi(z.re); 512 return cs.re * cosh(z.im) - cs.im * sinh(z.im) * 1i; 513 } 514 515 /** ditto */ 516 deprecated("Use `real` instead") 517 real cos(ireal y) 518 { 519 return cosh(y.im); 520 } 521 522 deprecated unittest 523 { 524 test(cos(0.0+0.0i)==1.0); 525 test(cos(1.3L+0.0i)==cos(1.3L)); 526 test(cos(5.2Li)== cosh(5.2L)); 527 } 528 529 /*************** 530 * Calculates the arc cosine of x, 531 * returning a value ranging from 0 to $(PI). 532 * 533 * $(TABLE_SV 534 * $(TR $(TH x) $(TH acos(x)) $(TH invalid?)) 535 * $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes)) 536 * $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes)) 537 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) 538 * ) 539 */ 540 real acos(real x) 541 { 542 return core.stdc.math.acosl(x); 543 } 544 545 unittest { 546 // NaN payloads 547 version(darwin){} 548 else { 549 test(isIdentical(acos(NaN(254)), NaN(254))); 550 } 551 } 552 553 /*************** 554 * Calculates the arc sine of x, 555 * returning a value ranging from -$(PI)/2 to $(PI)/2. 556 * 557 * $(TABLE_SV 558 * $(TR $(TH x) $(TH asin(x)) $(TH invalid?)) 559 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 560 * $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes)) 561 * $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes)) 562 * ) 563 */ 564 real asin(real x) 565 { 566 return core.stdc.math.asinl(x); 567 } 568 569 unittest { 570 // NaN payloads 571 version(darwin){} 572 else{ 573 test(isIdentical(asin(NaN(7249)), NaN(7249))); 574 } 575 } 576 577 /*************** 578 * Calculates the arc tangent of x, 579 * returning a value ranging from -$(PI)/2 to $(PI)/2. 580 * 581 * $(TABLE_SV 582 * $(TR $(TH x) $(TH atan(x)) $(TH invalid?)) 583 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 584 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)) 585 * ) 586 */ 587 real atan(real x) 588 { 589 return core.stdc.math.atanl(x); 590 } 591 592 unittest 593 { 594 // NaN payloads 595 test(isIdentical(atan(NaN(9876)), NaN(9876))); 596 } 597 598 /*************** 599 * Calculates the arc tangent of y / x, 600 * returning a value ranging from -$(PI) to $(PI). 601 * 602 * $(TABLE_SV 603 * $(TR $(TH y) $(TH x) $(TH atan(y, x))) 604 * $(TR $(TD $(NAN)) $(TD anything) $(TD $(NAN)) ) 605 * $(TR $(TD anything) $(TD $(NAN)) $(TD $(NAN)) ) 606 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) ) 607 * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) $(TD $(PLUSMN)0.0) ) 608 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT)0.0) $(TD $(PLUSMN)$(PI))) 609 * $(TR $(TD $(PLUSMN)0.0) $(TD -0.0) $(TD $(PLUSMN)$(PI))) 610 * $(TR $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) $(TD $(PI)/2) ) 611 * $(TR $(TD $(LT)0.0) $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) ) 612 * $(TR $(TD $(GT)0.0) $(TD $(INFIN)) $(TD $(PLUSMN)0.0) ) 613 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything) $(TD $(PLUSMN)$(PI)/2)) 614 * $(TR $(TD $(GT)0.0) $(TD -$(INFIN)) $(TD $(PLUSMN)$(PI)) ) 615 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN)) $(TD $(PLUSMN)$(PI)/4)) 616 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN)) $(TD $(PLUSMN)3$(PI)/4)) 617 * ) 618 */ 619 real atan2(real y, real x) 620 { 621 return core.stdc.math.atan2l(y,x); 622 } 623 624 unittest 625 { 626 // NaN payloads 627 test(isIdentical(atan2(5.3, NaN(9876)), NaN(9876))); 628 test(isIdentical(atan2(NaN(9876), 2.18), NaN(9876))); 629 } 630 631 /*********************************** 632 * Complex inverse sine 633 * 634 * asin(z) = -i log( sqrt(1-$(POWER z, 2)) + iz) 635 * where both log and sqrt are complex. 636 */ 637 deprecated("Use `std.complex.Complex` instead") 638 creal asin(creal z) 639 { 640 return -log(sqrt(1-z*z) + z*1i)*1i; 641 } 642 643 deprecated unittest 644 { 645 test(asin(sin(0+0i)) == 0 + 0i); 646 } 647 648 /*********************************** 649 * Complex inverse cosine 650 * 651 * acos(z) = $(PI)/2 - asin(z) 652 */ 653 deprecated("Use `std.complex.Complex` instead") 654 creal acos(creal z) 655 { 656 return PI_2 - asin(z); 657 } 658 659 660 /*********************************** 661 * Calculates the hyperbolic cosine of x. 662 * 663 * $(TABLE_SV 664 * $(TR $(TH x) $(TH cosh(x)) $(TH invalid?)) 665 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) ) 666 * ) 667 */ 668 real cosh(real x) 669 { 670 // cosh = (exp(x)+exp(-x))/2. 671 // The naive implementation works correctly. 672 real y = exp(x); 673 return (y + 1.0/y) * 0.5; 674 } 675 676 unittest 677 { 678 // NaN payloads 679 test(isIdentical(cosh(NaN(432)), NaN(432))); 680 } 681 682 /*********************************** 683 * Calculates the hyperbolic sine of x. 684 * 685 * $(TABLE_SV 686 * $(TR $(TH x) $(TH sinh(x)) $(TH invalid?)) 687 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) 688 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no)) 689 * ) 690 */ 691 real sinh(real x) 692 { 693 // sinh(x) = (exp(x)-exp(-x))/2; 694 // Very large arguments could cause an overflow, but 695 // the maximum value of x for which exp(x) + exp(-x)) != exp(x) 696 // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80. 697 if (fabs(x) > real.mant_dig * LN2) { 698 return copysign(0.5*exp(fabs(x)), x); 699 } 700 real y = expm1(x); 701 return 0.5 * y / (y+1) * (y+2); 702 } 703 704 unittest 705 { 706 // NaN payloads 707 test(isIdentical(sinh(NaN(0xABC)), NaN(0xABC))); 708 } 709 710 /*********************************** 711 * Calculates the hyperbolic tangent of x. 712 * 713 * $(TABLE_SV 714 * $(TR $(TH x) $(TH tanh(x)) $(TH invalid?)) 715 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) ) 716 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no)) 717 * ) 718 */ 719 real tanh(real x) 720 { 721 // tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x)) 722 if (fabs(x)> real.mant_dig * LN2){ 723 return copysign(1, x); 724 } 725 real y = expm1(2*x); 726 return y/(y + 2); 727 } 728 729 unittest 730 { 731 // NaN payloads 732 test(isIdentical(tanh(NaN(0xABC)), NaN(0xABC))); 733 } 734 735 /*********************************** 736 * hyperbolic sine, complex and imaginary 737 * 738 * sinh(z) = cos(z.im)*sinh(z.re) + sin(z.im)*cosh(z.re)i 739 */ 740 deprecated("Use `std.complex.Complex` instead") 741 creal sinh(creal z) 742 { 743 creal cs = expi(z.im); 744 return cs.re * sinh(z.re) + cs.im * cosh(z.re) * 1i; 745 } 746 747 /** ditto */ 748 deprecated("Use `real` instead") 749 ireal sinh(ireal y) 750 { 751 return sin(y.im)*1i; 752 } 753 754 deprecated unittest 755 { 756 test(sinh(4.2L + 0i)==sinh(4.2L)); 757 } 758 759 /*********************************** 760 * hyperbolic cosine, complex and imaginary 761 * 762 * cosh(z) = cos(z.im)*cosh(z.re) + sin(z.im)*sinh(z.re)i 763 */ 764 deprecated("Use `std.complex.Complex` instead") 765 creal cosh(creal z) 766 { 767 creal cs = expi(z.im); 768 return cs.re * cosh(z.re) + cs.im * sinh(z.re) * 1i; 769 } 770 771 /** ditto */ 772 deprecated("Use `ireal` instead") 773 real cosh(ireal y) 774 { 775 return cos(y.im); 776 } 777 778 deprecated unittest 779 { 780 test(cosh(8.3L + 0i)==cosh(8.3L)); 781 } 782 783 784 /*********************************** 785 * Calculates the inverse hyperbolic cosine of x. 786 * 787 * Mathematically, acosh(x) = log(x + sqrt( x*x - 1)) 788 * 789 * $(TABLE_SV 790 * $(SVH x, acosh(x) ) 791 * $(SV $(NAN), $(NAN) ) 792 * $(SV $(LT)1, $(NAN) ) 793 * $(SV 1, 0 ) 794 * $(SV +$(INFIN),+$(INFIN)) 795 * ) 796 */ 797 real acosh(real x) 798 { 799 if (x > 1/real.epsilon) 800 return LN2 + log(x); 801 else 802 return log(x + sqrt(x*x - 1)); 803 } 804 805 unittest 806 { 807 test(isNaN(acosh(0.9))); 808 test(isNaN(acosh(real.nan))); 809 test(acosh(1)==0.0); 810 test(acosh(real.infinity) == real.infinity); 811 // NaN payloads 812 test(isIdentical(acosh(NaN(0xABC)), NaN(0xABC))); 813 } 814 815 /*********************************** 816 * Calculates the inverse hyperbolic sine of x. 817 * 818 * Mathematically, 819 * --------------- 820 * asinh(x) = log( x + sqrt( x*x + 1 )) // if x >= +0 821 * asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0 822 * ------------- 823 * 824 * $(TABLE_SV 825 * $(SVH x, asinh(x) ) 826 * $(SV $(NAN), $(NAN) ) 827 * $(SV $(PLUSMN)0, $(PLUSMN)0 ) 828 * $(SV $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN)) 829 * ) 830 */ 831 real asinh(real x) 832 { 833 if (ocean.math.IEEE.fabs(x) > 1 / real.epsilon) // beyond this point, x*x + 1 == x*x 834 return ocean.math.IEEE.copysign(LN2 + log(ocean.math.IEEE.fabs(x)), x); 835 else 836 { 837 // sqrt(x*x + 1) == 1 + x * x / ( 1 + sqrt(x*x + 1) ) 838 return ocean.math.IEEE.copysign(log1p(ocean.math.IEEE.fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x); 839 } 840 } 841 842 unittest 843 { 844 test(isIdentical(0.0L,asinh(0.0))); 845 test(isIdentical(-0.0L,asinh(-0.0))); 846 test(asinh(real.infinity) == real.infinity); 847 test(asinh(-real.infinity) == -real.infinity); 848 test(isNaN(asinh(real.nan))); 849 // NaN payloads 850 test(isIdentical(asinh(NaN(0xABC)), NaN(0xABC))); 851 } 852 853 /*********************************** 854 * Calculates the inverse hyperbolic tangent of x, 855 * returning a value from ranging from -1 to 1. 856 * 857 * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2 858 * 859 * 860 * $(TABLE_SV 861 * $(SVH x, acosh(x) ) 862 * $(SV $(NAN), $(NAN) ) 863 * $(SV $(PLUSMN)0, $(PLUSMN)0) 864 * $(SV -$(INFIN), -0) 865 * ) 866 */ 867 real atanh(real x) 868 { 869 // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) ) 870 return 0.5 * log1p( 2 * x / (1 - x) ); 871 } 872 873 unittest 874 { 875 test(isIdentical(0.0L, atanh(0.0))); 876 test(isIdentical(-0.0L,atanh(-0.0))); 877 test(isIdentical(atanh(-1),-real.infinity)); 878 // Fails with -O because of DMD : https://issues.dlang.org/show_bug.cgi?id=13743 879 // Nothing can be done about it without feedback from DMD upstream 880 //assert(isIdentical(atanh(1),real.infinity)); 881 test(isNaN(atanh(-real.infinity))); 882 // NaN payloads 883 test(isIdentical(atanh(NaN(0xABC)), NaN(0xABC))); 884 } 885 886 /** ditto */ 887 deprecated("Use `std.complex.Complex` / `real` instead") 888 creal atanh(ireal y) 889 { 890 // Not optimised for accuracy or speed 891 return 0.5*(log(1+y) - log(1-y)); 892 } 893 894 /** ditto */ 895 deprecated("Use `std.complex.Complex` instead") 896 creal atanh(creal z) 897 { 898 // Not optimised for accuracy or speed 899 return 0.5 * (log(1 + z) - log(1-z)); 900 } 901 902 /* 903 * Powers and Roots 904 */ 905 906 /*************************************** 907 * Compute square root of x. 908 * 909 * $(TABLE_SV 910 * $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?)) 911 * $(TR $(TD -0.0) $(TD -0.0) $(TD no)) 912 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes)) 913 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no)) 914 * ) 915 */ 916 float sqrt(float x) /* intrinsic */ 917 { 918 version(D_InlineAsm_X86) 919 { 920 asm 921 { 922 fld x; 923 fsqrt; 924 } 925 } 926 else 927 { 928 return core.stdc.math.sqrtf(x); 929 } 930 } 931 932 double sqrt(double x) /* intrinsic */ /// ditto 933 { 934 version(D_InlineAsm_X86) 935 { 936 asm 937 { 938 fld x; 939 fsqrt; 940 } 941 } 942 else 943 { 944 return core.stdc.math.sqrt(x); 945 } 946 } 947 948 real sqrt(real x) /* intrinsic */ /// ditto 949 { 950 version(D_InlineAsm_X86) 951 { 952 asm 953 { 954 fld x; 955 fsqrt; 956 } 957 } 958 else 959 { 960 return core.stdc.math.sqrtl(x); 961 } 962 } 963 964 /** ditto */ 965 deprecated("Use `std.complex.Complex` instead") 966 creal sqrt(creal z) 967 { 968 969 if (z == 0.0) return z; 970 real x,y,w,r; 971 creal c; 972 973 x = ocean.math.IEEE.fabs(z.re); 974 y = ocean.math.IEEE.fabs(z.im); 975 if (x >= y) { 976 r = y / x; 977 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); 978 } else { 979 r = x / y; 980 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); 981 } 982 983 if (z.re >= 0) { 984 c = w + (z.im / (w + w)) * 1.0i; 985 } else { 986 if (z.im < 0) w = -w; 987 c = z.im / (w + w) + w * 1.0i; 988 } 989 return c; 990 } 991 992 unittest { 993 // NaN payloads 994 test(isIdentical(sqrt(NaN(0xABC)), NaN(0xABC))); 995 } 996 997 deprecated unittest { 998 test(sqrt(-1+0i) == 1i); 999 test(isIdentical(sqrt(0-0i), 0-0i)); 1000 test(cfeqrel(sqrt(4+16i)*sqrt(4+16i), 4+16i)>=real.mant_dig-2); 1001 } 1002 1003 /*************** 1004 * Calculates the cube root of x. 1005 * 1006 * $(TABLE_SV 1007 * $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?)) 1008 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) ) 1009 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) 1010 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) ) 1011 * ) 1012 */ 1013 real cbrt(real x) 1014 { 1015 return core.stdc.math.cbrtl(x); 1016 } 1017 1018 1019 unittest { 1020 // NaN payloads 1021 test(isIdentical(cbrt(NaN(0xABC)), NaN(0xABC))); 1022 } 1023 1024 public: 1025 1026 /** 1027 * Calculates e$(SUP x). 1028 * 1029 * $(TABLE_SV 1030 * $(TR $(TH x) $(TH e$(SUP x)) ) 1031 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1032 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 1033 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1034 * ) 1035 */ 1036 real exp(real x) { 1037 version(Naked_D_InlineAsm_X86) { 1038 // e^x = 2^(LOG2E*x) 1039 // (This is valid because the overflow & underflow limits for exp 1040 // and exp2 are so similar). 1041 return exp2(LOG2E*x); 1042 } else { 1043 return core.stdc.math.expl(x); 1044 } 1045 } 1046 1047 /** 1048 * Calculates the value of the natural logarithm base (e) 1049 * raised to the power of x, minus 1. 1050 * 1051 * For very small x, expm1(x) is more accurate 1052 * than exp(x)-1. 1053 * 1054 * $(TABLE_SV 1055 * $(TR $(TH x) $(TH e$(SUP x)-1) ) 1056 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) 1057 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1058 * $(TR $(TD -$(INFIN)) $(TD -1.0) ) 1059 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1060 * ) 1061 */ 1062 real expm1(real x) 1063 { 1064 version(Naked_D_InlineAsm_X86) { 1065 enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4 1066 asm { 1067 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1068 * Author: Don Clugston. 1069 * 1070 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. 1071 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) 1072 * and 2ym1 = (2^(y-rndint(y))-1). 1073 * If 2rndy < 0.5*real.epsilon, result is -1. 1074 * Implementation is otherwise the same as for exp2() 1075 */ 1076 naked; 1077 fld real ptr [ESP+4] ; // x 1078 mov AX, [ESP+4+8]; // AX = exponent and sign 1079 sub ESP, 12+8; // Create scratch space on the stack 1080 // [ESP,ESP+2] = scratchint 1081 // [ESP+4..+6, +8..+10, +10] = scratchreal 1082 // set scratchreal mantissa = 1.0 1083 mov dword ptr [ESP+8], 0; 1084 mov dword ptr [ESP+8+4], 0x80000000; 1085 and AX, 0x7FFF; // drop sign bit 1086 cmp AX, 0x401D; // avoid InvalidException in fist 1087 jae L_extreme; 1088 fldl2e; 1089 fmul ; // y = x*log2(e) 1090 fist dword ptr [ESP]; // scratchint = rndint(y) 1091 fisub dword ptr [ESP]; // y - rndint(y) 1092 // and now set scratchreal exponent 1093 mov EAX, [ESP]; 1094 add EAX, 0x3fff; 1095 jle short L_largenegative; 1096 cmp EAX,0x8000; 1097 jge short L_largepositive; 1098 mov [ESP+8+8],AX; 1099 f2xm1; // 2^(y-rndint(y)) -1 1100 fld real ptr [ESP+8] ; // 2^rndint(y) 1101 fmul ST(1), ST; 1102 fld1; 1103 fsubp ST(1), ST; 1104 fadd; 1105 add ESP,12+8; 1106 ret PARAMSIZE; 1107 1108 L_extreme: // Extreme exponent. X is very large positive, very 1109 // large negative, infinity, or NaN. 1110 fxam; 1111 fstsw AX; 1112 test AX, 0x0400; // NaN_or_zero, but we already know x!=0 1113 jz L_was_nan; // if x is NaN, returns x 1114 test AX, 0x0200; 1115 jnz L_largenegative; 1116 L_largepositive: 1117 // Set scratchreal = real.max. 1118 // squaring it will create infinity, and set overflow flag. 1119 mov word ptr [ESP+8+8], 0x7FFE; 1120 fstp ST(0); 1121 fld real ptr [ESP+8]; // load scratchreal 1122 fmul ST(0), ST; // square it, to create havoc! 1123 L_was_nan: 1124 add ESP,12+8; 1125 ret PARAMSIZE; 1126 L_largenegative: 1127 fstp ST(0); 1128 fld1; 1129 fchs; // return -1. Underflow flag is not set. 1130 add ESP,12+8; 1131 ret PARAMSIZE; 1132 } 1133 } else version(D_InlineAsm_X86_64) { 1134 asm 1135 { 1136 naked; 1137 } 1138 asm 1139 { 1140 fld real ptr [RSP+8]; // x 1141 mov AX,[RSP+8+8]; // AX = exponent and sign 1142 } 1143 asm 1144 { 1145 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. 1146 * Author: Don Clugston. 1147 * 1148 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. 1149 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) 1150 * and 2ym1 = (2^(y-rndint(y))-1). 1151 * If 2rndy < 0.5*real.epsilon, result is -1. 1152 * Implementation is otherwise the same as for exp2() 1153 */ 1154 sub RSP, 24; // Create scratch space on the stack 1155 // [RSP,RSP+2] = scratchint 1156 // [RSP+4..+6, +8..+10, +10] = scratchreal 1157 // set scratchreal mantissa = 1.0 1158 mov dword ptr [RSP+8], 0; 1159 mov dword ptr [RSP+8+4], 0x80000000; 1160 and AX, 0x7FFF; // drop sign bit 1161 cmp AX, 0x401D; // avoid InvalidException in fist 1162 jae L_extreme; 1163 fldl2e; 1164 fmul ; // y = x*log2(e) 1165 fist dword ptr [RSP]; // scratchint = rndint(y) 1166 fisub dword ptr [RSP]; // y - rndint(y) 1167 // and now set scratchreal exponent 1168 mov EAX, [RSP]; 1169 add EAX, 0x3fff; 1170 jle short L_largenegative; 1171 cmp EAX,0x8000; 1172 jge short L_largepositive; 1173 mov [RSP+8+8],AX; 1174 f2xm1; // 2^(y-rndint(y)) -1 1175 fld real ptr [RSP+8] ; // 2^rndint(y) 1176 fmul ST(1), ST; 1177 fld1; 1178 fsubp ST(1), ST; 1179 fadd; 1180 add RSP,24; 1181 ret; 1182 1183 L_extreme: // Extreme exponent. X is very large positive, very 1184 // large negative, infinity, or NaN. 1185 fxam; 1186 fstsw AX; 1187 test AX, 0x0400; // NaN_or_zero, but we already know x!=0 1188 jz L_was_nan; // if x is NaN, returns x 1189 test AX, 0x0200; 1190 jnz L_largenegative; 1191 L_largepositive: 1192 // Set scratchreal = real.max. 1193 // squaring it will create infinity, and set overflow flag. 1194 mov word ptr [RSP+8+8], 0x7FFE; 1195 fstp ST(0); 1196 fld real ptr [RSP+8]; // load scratchreal 1197 fmul ST(0), ST; // square it, to create havoc! 1198 L_was_nan: 1199 add RSP,24; 1200 ret; 1201 1202 L_largenegative: 1203 fstp ST(0); 1204 fld1; 1205 fchs; // return -1. Underflow flag is not set. 1206 add RSP,24; 1207 ret; 1208 } 1209 } else { 1210 return core.stdc.math.expm1l(x); 1211 } 1212 } 1213 1214 /** 1215 * Calculates 2$(SUP x). 1216 * 1217 * $(TABLE_SV 1218 * $(TR $(TH x) $(TH exp2(x)) ) 1219 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) 1220 * $(TR $(TD -$(INFIN)) $(TD +0.0) ) 1221 * $(TR $(TD $(NAN)) $(TD $(NAN)) ) 1222 * ) 1223 */ 1224 real exp2(real x) 1225 { 1226 version(Naked_D_InlineAsm_X86) { 1227 enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4 1228 asm { 1229 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1230 * Author: Don Clugston. 1231 * 1232 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) 1233 * The trick for high performance is to avoid the fscale(28cycles on core2), 1234 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1235 * 1236 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1237 * because it will set the Invalid Operation flag is overflow or NaN occurs. 1238 * Fortunately, whenever this happens the result would be zero or infinity. 1239 * 1240 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1241 * work for the (very rare) cases where the result is subnormal. So we fall back 1242 * to the slow method in that case. 1243 */ 1244 naked; 1245 fld real ptr [ESP+4] ; // x 1246 mov AX, [ESP+4+8]; // AX = exponent and sign 1247 sub ESP, 12+8; // Create scratch space on the stack 1248 // [ESP,ESP+2] = scratchint 1249 // [ESP+4..+6, +8..+10, +10] = scratchreal 1250 // set scratchreal mantissa = 1.0 1251 mov dword ptr [ESP+8], 0; 1252 mov dword ptr [ESP+8+4], 0x80000000; 1253 and AX, 0x7FFF; // drop sign bit 1254 cmp AX, 0x401D; // avoid InvalidException in fist 1255 jae L_extreme; 1256 fist dword ptr [ESP]; // scratchint = rndint(x) 1257 fisub dword ptr [ESP]; // x - rndint(x) 1258 // and now set scratchreal exponent 1259 mov EAX, [ESP]; 1260 add EAX, 0x3fff; 1261 jle short L_subnormal; 1262 cmp EAX,0x8000; 1263 jge short L_overflow; 1264 mov [ESP+8+8],AX; 1265 L_normal: 1266 f2xm1; 1267 fld1; 1268 fadd; // 2^(x-rndint(x)) 1269 fld real ptr [ESP+8] ; // 2^rndint(x) 1270 add ESP,12+8; 1271 fmulp ST(1), ST; 1272 ret PARAMSIZE; 1273 1274 L_subnormal: 1275 // Result will be subnormal. 1276 // In this rare case, the simple poking method doesn't work. 1277 // The speed doesn't matter, so use the slow fscale method. 1278 fild dword ptr [ESP]; // scratchint 1279 fld1; 1280 fscale; 1281 fstp real ptr [ESP+8]; // scratchreal = 2^scratchint 1282 fstp ST(0); // drop scratchint 1283 jmp L_normal; 1284 1285 L_extreme: // Extreme exponent. X is very large positive, very 1286 // large negative, infinity, or NaN. 1287 fxam; 1288 fstsw AX; 1289 test AX, 0x0400; // NaN_or_zero, but we already know x!=0 1290 jz L_was_nan; // if x is NaN, returns x 1291 // set scratchreal = real.min 1292 // squaring it will return 0, setting underflow flag 1293 mov word ptr [ESP+8+8], 1; 1294 test AX, 0x0200; 1295 jnz L_waslargenegative; 1296 L_overflow: 1297 // Set scratchreal = real.max. 1298 // squaring it will create infinity, and set overflow flag. 1299 mov word ptr [ESP+8+8], 0x7FFE; 1300 L_waslargenegative: 1301 fstp ST(0); 1302 fld real ptr [ESP+8]; // load scratchreal 1303 fmul ST(0), ST; // square it, to create havoc! 1304 L_was_nan: 1305 add ESP,12+8; 1306 ret PARAMSIZE; 1307 } 1308 } else version(D_InlineAsm_X86_64) { 1309 asm 1310 { 1311 naked; 1312 } 1313 asm 1314 { 1315 fld real ptr [RSP+8]; // x 1316 mov AX,[RSP+8+8]; // AX = exponent and sign 1317 } 1318 asm 1319 { 1320 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. 1321 * Author: Don Clugston. 1322 * 1323 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) 1324 * The trick for high performance is to avoid the fscale(28cycles on core2), 1325 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. 1326 * 1327 * We can do frndint by using fist. BUT we can't use it for huge numbers, 1328 * because it will set the Invalid Operation flag is overflow or NaN occurs. 1329 * Fortunately, whenever this happens the result would be zero or infinity. 1330 * 1331 * We can perform fscale by directly poking into the exponent. BUT this doesn't 1332 * work for the (very rare) cases where the result is subnormal. So we fall back 1333 * to the slow method in that case. 1334 */ 1335 sub RSP, 24; // Create scratch space on the stack 1336 // [RSP,RSP+2] = scratchint 1337 // [RSP+4..+6, +8..+10, +10] = scratchreal 1338 // set scratchreal mantissa = 1.0 1339 mov dword ptr [RSP+8], 0; 1340 mov dword ptr [RSP+8+4], 0x80000000; 1341 and AX, 0x7FFF; // drop sign bit 1342 cmp AX, 0x401D; // avoid InvalidException in fist 1343 jae L_extreme; 1344 fist dword ptr [RSP]; // scratchint = rndint(x) 1345 fisub dword ptr [RSP]; // x - rndint(x) 1346 // and now set scratchreal exponent 1347 mov EAX, [RSP]; 1348 add EAX, 0x3fff; 1349 jle short L_subnormal; 1350 cmp EAX,0x8000; 1351 jge short L_overflow; 1352 mov [RSP+8+8],AX; 1353 L_normal: 1354 f2xm1; 1355 fld1; 1356 fadd; // 2^(x-rndint(x)) 1357 fld real ptr [RSP+8] ; // 2^rndint(x) 1358 add RSP,24; 1359 fmulp ST(1), ST; 1360 ret; 1361 1362 L_subnormal: 1363 // Result will be subnormal. 1364 // In this rare case, the simple poking method doesn't work. 1365 // The speed doesn't matter, so use the slow fscale method. 1366 fild dword ptr [RSP]; // scratchint 1367 fld1; 1368 fscale; 1369 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint 1370 fstp ST(0); // drop scratchint 1371 jmp L_normal; 1372 1373 L_extreme: // Extreme exponent. X is very large positive, very 1374 // large negative, infinity, or NaN. 1375 fxam; 1376 fstsw AX; 1377 test AX, 0x0400; // NaN_or_zero, but we already know x!=0 1378 jz L_was_nan; // if x is NaN, returns x 1379 // set scratchreal = real.min 1380 // squaring it will return 0, setting underflow flag 1381 mov word ptr [RSP+8+8], 1; 1382 test AX, 0x0200; 1383 jnz L_waslargenegative; 1384 L_overflow: 1385 // Set scratchreal = real.max. 1386 // squaring it will create infinity, and set overflow flag. 1387 mov word ptr [RSP+8+8], 0x7FFE; 1388 L_waslargenegative: 1389 fstp ST(0); 1390 fld real ptr [RSP+8]; // load scratchreal 1391 fmul ST(0), ST; // square it, to create havoc! 1392 L_was_nan: 1393 add RSP,24; 1394 ret; 1395 } 1396 } else { 1397 return core.stdc.math.exp2l(x); 1398 } 1399 } 1400 1401 unittest { 1402 // NaN payloads 1403 test(isIdentical(exp(NaN(0xABC)), NaN(0xABC))); 1404 } 1405 1406 unittest { 1407 // NaN payloads 1408 test(isIdentical(expm1(NaN(0xABC)), NaN(0xABC))); 1409 } 1410 1411 unittest { 1412 // NaN payloads 1413 test(isIdentical(exp2(NaN(0xABC)), NaN(0xABC))); 1414 } 1415 1416 /* 1417 * Powers and Roots 1418 */ 1419 1420 /************************************** 1421 * Calculate the natural logarithm of x. 1422 * 1423 * $(TABLE_SV 1424 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?)) 1425 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 1426 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 1427 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 1428 * ) 1429 */ 1430 real log(real x) 1431 { 1432 return core.stdc.math.logl(x); 1433 } 1434 1435 unittest { 1436 // NaN payloads 1437 test(isIdentical(log(NaN(0xABC)), NaN(0xABC))); 1438 } 1439 1440 /****************************************** 1441 * Calculates the natural logarithm of 1 + x. 1442 * 1443 * For very small x, log1p(x) will be more accurate than 1444 * log(1 + x). 1445 * 1446 * $(TABLE_SV 1447 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?)) 1448 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no)) 1449 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 1450 * $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD no) $(TD yes)) 1451 * $(TR $(TD +$(INFIN)) $(TD -$(INFIN)) $(TD no) $(TD no)) 1452 * ) 1453 */ 1454 real log1p(real x) 1455 { 1456 return core.stdc.math.log1pl(x); 1457 } 1458 1459 unittest { 1460 // NaN payloads 1461 test(isIdentical(log1p(NaN(0xABC)), NaN(0xABC))); 1462 } 1463 1464 /*************************************** 1465 * Calculates the base-2 logarithm of x: 1466 * $(SUB log, 2)x 1467 * 1468 * $(TABLE_SV 1469 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?)) 1470 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) ) 1471 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) ) 1472 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) ) 1473 * ) 1474 */ 1475 real log2(real x) 1476 { 1477 return core.stdc.math.log2l(x); 1478 } 1479 1480 unittest { 1481 // NaN payloads 1482 test(isIdentical(log2(NaN(0xABC)), NaN(0xABC))); 1483 } 1484 1485 /************************************** 1486 * Calculate the base-10 logarithm of x. 1487 * 1488 * $(TABLE_SV 1489 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?)) 1490 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) 1491 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) 1492 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) 1493 * ) 1494 */ 1495 real log10(real x) 1496 { 1497 return core.stdc.math.log10l(x); 1498 } 1499 1500 unittest { 1501 // NaN payloads 1502 test(isIdentical(log10(NaN(0xABC)), NaN(0xABC))); 1503 } 1504 1505 /*********************************** 1506 * Exponential, complex and imaginary 1507 * 1508 * For complex numbers, the exponential function is defined as 1509 * 1510 * exp(z) = exp(z.re)cos(z.im) + exp(z.re)sin(z.im)i. 1511 * 1512 * For a pure imaginary argument, 1513 * exp(θi) = cos(θ) + sin(θ)i. 1514 * 1515 */ 1516 deprecated("Use `std.complex.Complex` / `real` instead") 1517 creal exp(ireal y) 1518 { 1519 return expi(y.im); 1520 } 1521 1522 /** ditto */ 1523 deprecated("Use `std.complex.Complex` instead") 1524 creal exp(creal z) 1525 { 1526 return expi(z.im) * exp(z.re); 1527 } 1528 1529 deprecated unittest { 1530 test(exp(1.3e5Li)==cos(1.3e5L)+sin(1.3e5L)*1i); 1531 test(exp(0.0Li)==1L+0.0Li); 1532 test(exp(7.2 + 0.0i) == exp(7.2L)); 1533 creal c = exp(ireal.nan); 1534 test(isNaN(c.re) && isNaN(c.im)); 1535 c = exp(ireal.infinity); 1536 test(isNaN(c.re) && isNaN(c.im)); 1537 } 1538 1539 /*********************************** 1540 * Natural logarithm, complex 1541 * 1542 * Returns complex logarithm to the base e (2.718...) of 1543 * the complex argument x. 1544 * 1545 * If z = x + iy, then 1546 * log(z) = log(abs(z)) + i arctan(y/x). 1547 * 1548 * The arctangent ranges from -PI to +PI. 1549 * There are branch cuts along both the negative real and negative 1550 * imaginary axes. For pure imaginary arguments, use one of the 1551 * following forms, depending on which branch is required. 1552 * ------------ 1553 * log( 0.0 + yi) = log(-y) + PI_2i // y<=-0.0 1554 * log(-0.0 + yi) = log(-y) - PI_2i // y<=-0.0 1555 * ------------ 1556 */ 1557 deprecated("Use `std.complex.Complex` instead") 1558 creal log(creal z) 1559 { 1560 return log(abs(z)) + atan2(z.im, z.re)*1i; 1561 } 1562 1563 /* 1564 * feqrel for complex numbers. Returns the worst relative 1565 * equality of the two components. 1566 */ 1567 deprecated private int cfeqrel(creal a, creal b) 1568 { 1569 int intmin(int a, int b) { return a<b? a: b; } 1570 return intmin(feqrel(a.re, b.re), feqrel(a.im, b.im)); 1571 } 1572 1573 deprecated unittest { 1574 1575 test(log(3.0L +0i) == log(3.0L)+0i); 1576 test(cfeqrel(log(0.0L-2i),( log(2.0L)-PI_2*1i)) >= real.mant_dig-10); 1577 test(cfeqrel(log(0.0L+2i),( log(2.0L)+PI_2*1i)) >= real.mant_dig-10); 1578 } 1579 1580 /** 1581 * Fast integral powers. 1582 */ 1583 real pow(real x, uint n) 1584 { 1585 real p; 1586 1587 switch (n) 1588 { 1589 case 0: 1590 p = 1.0; 1591 break; 1592 1593 case 1: 1594 p = x; 1595 break; 1596 1597 case 2: 1598 p = x * x; 1599 break; 1600 1601 default: 1602 p = 1.0; 1603 while (1){ 1604 if (n & 1) 1605 p *= x; 1606 n >>= 1; 1607 if (!n) 1608 break; 1609 x *= x; 1610 } 1611 break; 1612 } 1613 return p; 1614 } 1615 1616 /** ditto */ 1617 real pow(real x, int n) 1618 { 1619 if (n < 0) return pow(x, cast(real)n); 1620 else return pow(x, cast(uint)n); 1621 } 1622 1623 /********************************************* 1624 * Calculates x$(SUP y). 1625 * 1626 * $(TABLE_SV 1627 * $(TR $(TH x) $(TH y) $(TH pow(x, y)) 1628 * $(TH div 0) $(TH invalid?)) 1629 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0) 1630 * $(TD no) $(TD no) ) 1631 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN)) 1632 * $(TD no) $(TD no) ) 1633 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0) 1634 * $(TD no) $(TD no) ) 1635 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0) 1636 * $(TD no) $(TD no) ) 1637 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN)) 1638 * $(TD no) $(TD no) ) 1639 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN)) 1640 * $(TD no) $(TD no) ) 1641 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0) 1642 * $(TD no) $(TD no) ) 1643 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN)) 1644 * $(TD no) $(TD no) ) 1645 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN)) 1646 * $(TD no) $(TD no)) 1647 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0) 1648 * $(TD no) $(TD no) ) 1649 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0) 1650 * $(TD no) $(TD no) ) 1651 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) 1652 * $(TD no) $(TD yes) ) 1653 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN)) 1654 * $(TD no) $(TD yes)) 1655 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF)) 1656 * $(TD yes) $(TD no) ) 1657 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN)) 1658 * $(TD yes) $(TD no)) 1659 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0) 1660 * $(TD no) $(TD no) ) 1661 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0) 1662 * $(TD no) $(TD no) ) 1663 * ) 1664 */ 1665 real pow(real x, real y) 1666 { 1667 if (isNaN(y)) 1668 return y; 1669 1670 if (y == 0) 1671 return 1; // even if x is $(NAN) 1672 if (isNaN(x) && y != 0) 1673 return x; 1674 if (isInfinity(y)) 1675 { 1676 if (ocean.math.IEEE.fabs(x) > 1) 1677 { 1678 if (signbit(y)) 1679 return +0.0; 1680 else 1681 return real.infinity; 1682 } 1683 else if (ocean.math.IEEE.fabs(x) == 1) 1684 { 1685 return NaN(TANGO_NAN.POW_DOMAIN); 1686 } 1687 else // < 1 1688 { 1689 if (signbit(y)) 1690 return real.infinity; 1691 else 1692 return +0.0; 1693 } 1694 } 1695 if (isInfinity(x)) 1696 { 1697 if (signbit(x)) 1698 { 1699 long i; 1700 i = cast(long)y; 1701 if (y > 0) 1702 { 1703 if (i == y && i & 1) 1704 return -real.infinity; 1705 else 1706 return real.infinity; 1707 } 1708 else if (y < 0) 1709 { 1710 if (i == y && i & 1) 1711 return -0.0; 1712 else 1713 return +0.0; 1714 } 1715 } 1716 else 1717 { 1718 if (y > 0) 1719 return real.infinity; 1720 else if (y < 0) 1721 return +0.0; 1722 } 1723 } 1724 1725 if (x == 0.0) 1726 { 1727 if (signbit(x)) 1728 { 1729 long i; 1730 1731 i = cast(long)y; 1732 if (y > 0) 1733 { 1734 if (i == y && i & 1) 1735 return -0.0; 1736 else 1737 return +0.0; 1738 } 1739 else if (y < 0) 1740 { 1741 if (i == y && i & 1) 1742 return -real.infinity; 1743 else 1744 return real.infinity; 1745 } 1746 } 1747 else 1748 { 1749 if (y > 0) 1750 return +0.0; 1751 else if (y < 0) 1752 return real.infinity; 1753 } 1754 } 1755 1756 return core.stdc.math.powl(x, y); 1757 } 1758 1759 unittest 1760 { 1761 real x = 46; 1762 1763 test(pow(x,0) == 1.0); 1764 test(pow(x,1) == x); 1765 test(pow(x,2) == x * x); 1766 test(pow(x,3) == x * x * x); 1767 test(pow(x,8) == (x * x) * (x * x) * (x * x) * (x * x)); 1768 // NaN payloads 1769 test(isIdentical(pow(NaN(0xABC), 19), NaN(0xABC))); 1770 } 1771 1772 /*********************************************************************** 1773 * Calculates the length of the 1774 * hypotenuse of a right-angled triangle with sides of length x and y. 1775 * The hypotenuse is the value of the square root of 1776 * the sums of the squares of x and y: 1777 * 1778 * sqrt($(POW x, 2) + $(POW y, 2)) 1779 * 1780 * Note that hypot(x, y), hypot(y, x) and 1781 * hypot(x, -y) are equivalent. 1782 * 1783 * $(TABLE_SV 1784 * $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?)) 1785 * $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no)) 1786 * $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no)) 1787 * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no)) 1788 * ) 1789 */ 1790 real hypot(real x, real y) 1791 { 1792 // Scale x and y to avoid underflow and overflow. 1793 // If one is huge and the other tiny, return the larger. 1794 // If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2). 1795 // If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon). 1796 1797 static immutable real SQRTMIN = 0x8.0p-8195L; // 0.5 * sqrt(real.min_normal); // This is a power of 2. 1798 static immutable real SQRTMAX = 1.0L / SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max)) 1799 1800 static assert(2 * (SQRTMAX / 2) * (SQRTMAX / 2) <= real.max); 1801 1802 // Proves that sqrt(real.max) ~~ 0.5/sqrt(real.min_normal) 1803 static assert(real.min_normal * real.max > 2 && real.min_normal * real.max <= 4); 1804 1805 real u = fabs(x); 1806 real v = fabs(y); 1807 if (!(u >= v)) // check for NaN as well. 1808 { 1809 v = u; 1810 u = fabs(y); 1811 if (u == real.infinity) return u; // hypot(inf, nan) == inf 1812 if (v == real.infinity) return v; // hypot(nan, inf) == inf 1813 } 1814 1815 // Now u >= v, or else one is NaN. 1816 if (v >= SQRTMAX * 0.5) 1817 { 1818 // hypot(huge, huge) -- avoid overflow 1819 u *= SQRTMIN * 0.5; 1820 v *= SQRTMIN * 0.5; 1821 return sqrt(u * u + v * v) * SQRTMAX * 2.0; 1822 } 1823 1824 if (u <= SQRTMIN) 1825 { 1826 // hypot (tiny, tiny) -- avoid underflow 1827 // This is only necessary to avoid setting the underflow 1828 // flag. 1829 u *= SQRTMAX / real.epsilon; 1830 v *= SQRTMAX / real.epsilon; 1831 return sqrt(u * u + v * v) * SQRTMIN * real.epsilon; 1832 } 1833 1834 if (u * real.epsilon > v) 1835 { 1836 // hypot (huge, tiny) = huge 1837 return u; 1838 } 1839 1840 // both are in the normal range 1841 return sqrt(u * u + v * v); 1842 } 1843 1844 unittest 1845 { 1846 static real[3][] vals = // x,y,hypot 1847 [ 1848 [ 0.0, 0.0, 0.0], 1849 [ 0.0, -0.0, 0.0], 1850 [-0.0, -0.0, 0.0], 1851 [ 3.0, 4.0, 5.0], 1852 [-300, -400, 500], 1853 [ 0.0, 7.0, 7.0], 1854 [ 9.0, 9.0 * real.epsilon, 9.0], 1855 [ 0xb.0p+8188L /+88 / (64 * sqrt(real.min_normal))+/, 0xd.2p+8188L /+105 / (64 * sqrt(real.min_normal))+/, 0x8.9p+8189L /+137 / (64 * sqrt(real.min_normal))+/], 1856 [ 0xb.0p+8187L /+88 / (128 * sqrt(real.min_normal))+/, 0xd.2p+8187L /+105 / (128 * sqrt(real.min_normal))+/, 0x8.9p+8188L /+137 / (128 * sqrt(real.min_normal))+/], 1857 [ 3 * real.min_normal * real.epsilon, 4 * real.min_normal * real.epsilon, 5 * real.min_normal * real.epsilon], 1858 [ real.min_normal, real.min_normal, 0x1.6a09e667f3bcc908p-16382L /+sqrt(2.0L) * real,min_normal+/], 1859 [ real.max / 2.0, real.max / 2.0, 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/], 1860 [ 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/, 0x1.6a09e667f3bcc908p+16383L /+real.max / sqrt(2.0L)+/, real.max], 1861 [ real.max, 1.0, real.max], 1862 [ real.infinity, real.nan, real.infinity], 1863 [ real.nan, real.infinity, real.infinity], 1864 [ real.nan, real.nan, real.nan], 1865 [ real.nan, real.max, real.nan], 1866 [ real.max, real.nan, real.nan] 1867 ]; 1868 1869 for (int i = 0; i < vals.length; i++) 1870 { 1871 real x = vals[i][0]; 1872 real y = vals[i][1]; 1873 real z = vals[i][2]; 1874 real h = hypot(x, y); 1875 1876 test(isIdentical(z, h)); 1877 } 1878 // NaN payloads 1879 test(isIdentical(hypot(NaN(0xABC), 3.14), NaN(0xABC))); 1880 test(isIdentical(hypot(7.6e39, NaN(0xABC)), NaN(0xABC))); 1881 } 1882 1883 /*********************************** 1884 * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2) 1885 * + $(SUB a,3)$(POWER x,3); ... 1886 * 1887 * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2) 1888 * + x($(SUB a, 3) + ...))) 1889 * Params: 1890 * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc. 1891 * x = point in which to evaluate polynomial 1892 */ 1893 T poly(T)(T x, T[] A) 1894 { 1895 verify(A.length > 0); 1896 1897 version (Naked_D_InlineAsm_X86) { 1898 static immutable bool Use_D_InlineAsm_X86 = true; 1899 } else static immutable bool Use_D_InlineAsm_X86 = false; 1900 1901 // BUG (Inherited from Phobos): This code assumes a frame pointer in EBP. 1902 // This is not in the spec. 1903 static if (Use_D_InlineAsm_X86 && is(T==real) && T.sizeof == 10) { 1904 asm // assembler by W. Bright 1905 { 1906 // EDX = (A.length - 1) * real.sizeof 1907 mov ECX,A[EBP] ; // ECX = A.length 1908 dec ECX ; 1909 lea EDX,[ECX][ECX*8] ; 1910 add EDX,ECX ; 1911 add EDX,A+4[EBP] ; 1912 fld real ptr [EDX] ; // ST0 = coeff[ECX] 1913 jecxz return_ST ; 1914 fld x[EBP] ; // ST0 = x 1915 fxch ST(1) ; // ST1 = x, ST0 = r 1916 align 4 ; 1917 L2: fmul ST,ST(1) ; // r *= x 1918 fld real ptr -10[EDX] ; 1919 sub EDX,10 ; // deg-- 1920 faddp ST(1),ST ; 1921 dec ECX ; 1922 jne L2 ; 1923 fxch ST(1) ; // ST1 = r, ST0 = x 1924 fstp ST(0) ; // dump x 1925 align 4 ; 1926 return_ST: ; 1927 ; 1928 } 1929 } else static if ( Use_D_InlineAsm_X86 && is(T==real) && T.sizeof==12){ 1930 asm // assembler by W. Bright 1931 { 1932 // EDX = (A.length - 1) * real.sizeof 1933 mov ECX,A[EBP] ; // ECX = A.length 1934 dec ECX ; 1935 lea EDX,[ECX*8] ; 1936 lea EDX,[EDX][ECX*4] ; 1937 add EDX,A+4[EBP] ; 1938 fld real ptr [EDX] ; // ST0 = coeff[ECX] 1939 jecxz return_ST ; 1940 fld x ; // ST0 = x 1941 fxch ST(1) ; // ST1 = x, ST0 = r 1942 align 4 ; 1943 L2: fmul ST,ST(1) ; // r *= x 1944 fld real ptr -12[EDX] ; 1945 sub EDX,12 ; // deg-- 1946 faddp ST(1),ST ; 1947 dec ECX ; 1948 jne L2 ; 1949 fxch ST(1) ; // ST1 = r, ST0 = x 1950 fstp ST(0) ; // dump x 1951 align 4 ; 1952 return_ST: ; 1953 ; 1954 } 1955 } else { 1956 ptrdiff_t i = A.length - 1; 1957 real r = A[i]; 1958 while (--i >= 0) 1959 { 1960 r *= x; 1961 r += A[i]; 1962 } 1963 return r; 1964 } 1965 } 1966 1967 unittest 1968 { 1969 real x = 3.1; 1970 static immutable real[] pp = [56.1L, 32.7L, 6L]; 1971 1972 test( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) ); 1973 1974 test(isIdentical(poly(NaN(0xABC), pp), NaN(0xABC))); 1975 } 1976 1977 package { 1978 T rationalPoly(T)(T x, T [] numerator, T [] denominator) 1979 { 1980 return poly(x, numerator)/poly(x, denominator); 1981 } 1982 } 1983 1984 /* 1985 * Rounding (returning real) 1986 */ 1987 1988 /** 1989 * Returns the value of x rounded downward to the next integer 1990 * (toward negative infinity). 1991 */ 1992 real floor(real x) 1993 { 1994 return core.stdc.math.floorl(x); 1995 } 1996 1997 unittest { 1998 test(isIdentical(floor(NaN(0xABC)), NaN(0xABC))); 1999 } 2000 2001 /** 2002 * Returns the value of x rounded upward to the next integer 2003 * (toward positive infinity). 2004 */ 2005 real ceil(real x) 2006 { 2007 return core.stdc.math.ceill(x); 2008 } 2009 2010 unittest { 2011 test(isIdentical(ceil(NaN(0xABC)), NaN(0xABC))); 2012 } 2013 2014 /** 2015 * Return the value of x rounded to the nearest integer. 2016 * If the fractional part of x is exactly 0.5, the return value is rounded to 2017 * the even integer. 2018 */ 2019 real round(real x) 2020 { 2021 return core.stdc.math.roundl(x); 2022 } 2023 2024 unittest { 2025 test(isIdentical(round(NaN(0xABC)), NaN(0xABC))); 2026 } 2027 2028 /** 2029 * Returns the integer portion of x, dropping the fractional portion. 2030 * 2031 * This is also known as "chop" rounding. 2032 */ 2033 real trunc(real x) 2034 { 2035 return core.stdc.math.truncl(x); 2036 } 2037 2038 unittest { 2039 test(isIdentical(trunc(NaN(0xABC)), NaN(0xABC))); 2040 } 2041 2042 /** 2043 * Rounds x to the nearest int or long. 2044 * 2045 * This is generally the fastest method to convert a floating-point number 2046 * to an integer. Note that the results from this function 2047 * depend on the rounding mode, if the fractional part of x is exactly 0.5. 2048 * If using the default rounding mode (ties round to even integers) 2049 * rndint(4.5) == 4, rndint(5.5)==6. 2050 */ 2051 int rndint(real x) 2052 { 2053 version(Naked_D_InlineAsm_X86) 2054 { 2055 int n; 2056 asm 2057 { 2058 fld x; 2059 fistp n; 2060 } 2061 return n; 2062 } 2063 else 2064 { 2065 return cast(int) core.stdc.math.lrintl(x); 2066 } 2067 } 2068 2069 /** ditto */ 2070 long rndlong(real x) 2071 { 2072 version(Naked_D_InlineAsm_X86) 2073 { 2074 long n; 2075 asm 2076 { 2077 fld x; 2078 fistp n; 2079 } 2080 return n; 2081 } 2082 else 2083 { 2084 return core.stdc.math.llrintl(x); 2085 } 2086 } 2087 2088 version(D_InlineAsm_X86) 2089 { 2090 // Won't work for anything else yet 2091 2092 unittest 2093 { 2094 int r = getIeeeRounding; 2095 test(r==RoundingMode.ROUNDTONEAREST); 2096 real b = 5.5; 2097 int cnear = ocean.math.Math.rndint(b); 2098 test(cnear == 6); 2099 auto oldrounding = setIeeeRounding(RoundingMode.ROUNDDOWN); 2100 scope (exit) setIeeeRounding(oldrounding); 2101 2102 test(getIeeeRounding==RoundingMode.ROUNDDOWN); 2103 2104 int cdown = ocean.math.Math.rndint(b); 2105 test(cdown==5); 2106 } 2107 2108 unittest 2109 { 2110 // Check that the previous test correctly restored the rounding mode 2111 test(getIeeeRounding==RoundingMode.ROUNDTONEAREST); 2112 } 2113 } 2114 2115 /*************************************************************************** 2116 2117 Integer pow function. Returns the power'th power of base 2118 2119 Params: 2120 base = base number 2121 power = power 2122 2123 Returns: 2124 the power'th power of base 2125 2126 ***************************************************************************/ 2127 2128 public ulong pow ( ulong base, ulong power ) 2129 { 2130 ulong res = void; 2131 2132 switch (power) 2133 { 2134 case 0: 2135 res = 1; 2136 break; 2137 case 1: 2138 res = base; 2139 break; 2140 case 2: 2141 res = base * base; 2142 break; 2143 2144 default: 2145 res = 1; 2146 2147 while (1) 2148 { 2149 if (power & 1) res *= base; 2150 power >>= 1; 2151 if (!power) break; 2152 base *= base; 2153 } 2154 break; 2155 } 2156 2157 return res; 2158 } 2159 2160 unittest 2161 { 2162 ulong x = 46; 2163 2164 test(pow(x,0UL) == 1); 2165 test(pow(x,1UL) == x); 2166 test(pow(x,2UL) == x * x); 2167 test(pow(x,3UL) == x * x * x); 2168 test(pow(x,8UL) == (x * x) * (x * x) * (x * x) * (x * x)); 2169 } 2170 2171 /******************************************************************************* 2172 2173 Does an integer division, rounding towards the nearest integer. 2174 Rounds to the even one if both integers are equal near. 2175 2176 Params: 2177 a = number to divide 2178 b = number to divide by 2179 2180 Returns: 2181 number divided according to given description 2182 2183 *******************************************************************************/ 2184 2185 T divRoundEven(T)(T a, T b) 2186 { 2187 // both integers equal near? 2188 if (b % 2 == 0 && (a % b == b / 2 || a % b == -b / 2)) 2189 { 2190 auto div_rounded_down = a / b; 2191 2192 auto add = div_rounded_down < 0 ? -1 : 1; 2193 2194 return div_rounded_down % 2 == 0 ? 2195 div_rounded_down : div_rounded_down + add; 2196 } 2197 2198 if ( (a >= 0) != (b >= 0) ) 2199 { 2200 return (a - b / 2) / b; 2201 } 2202 else 2203 { 2204 return (a + b / 2) / b; 2205 } 2206 } 2207 2208 version (unittest) 2209 { 2210 import ocean.math.Math : rndlong; 2211 } 2212 2213 unittest 2214 { 2215 long roundDivCheat ( long a, long b ) 2216 { 2217 real x = cast(real)a / cast(real)b; 2218 return rndlong(x); 2219 } 2220 2221 test(divRoundEven(-3, 2) == -2); 2222 test(divRoundEven(3, 2) == 2); 2223 test(divRoundEven(-3, -2) == 2); 2224 test(divRoundEven(3, -2) == -2); 2225 2226 test(divRoundEven(7, 11) == 1); 2227 test(divRoundEven(11, 11) == 1); 2228 test(divRoundEven(16, 11) == 1); 2229 test(divRoundEven(-17, 11) == -2); 2230 test(divRoundEven(-17, 11) == -2); 2231 test(divRoundEven(-16, 11) == -1); 2232 2233 test(divRoundEven(17, -11) == -2); 2234 test(divRoundEven(16, -11) == -1); 2235 test(divRoundEven(-17, -11) == 2); 2236 test(divRoundEven(-16, -11) == 1); 2237 2238 for (int i = -100; i <= 100; ++i) for (int j = -100; j <= 100; ++j) 2239 { 2240 if (j != 0) 2241 { 2242 test (divRoundEven(i,j) == roundDivCheat(i,j)); 2243 } 2244 } 2245 }