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