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