1 // Written in the D programming language. 2 3 /** 4 This is a submodule of $(MREF std, math). 5 6 It contains classical algebraic functions like `abs`, `sqrt`, and `poly`. 7 8 Copyright: Copyright The D Language Foundation 2000 - 2011. 9 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 10 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 11 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 12 Source: $(PHOBOSSRC std/math/algebraic.d) 13 14 Macros: 15 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> 16 <caption>Special Values</caption> 17 $0</table> 18 NAN = $(RED NAN) 19 POWER = $1<sup>$2</sup> 20 SUB = $1<sub>$2</sub> 21 PLUSMN = ± 22 INFIN = ∞ 23 PLUSMNINF = ±∞ 24 LT = < 25 26 */ 27 28 module std.math.algebraic; 29 30 static import core.math; 31 static import core.stdc.math; 32 import std.traits : CommonType, isFloatingPoint, isIntegral, isSigned, Unqual; 33 34 /*********************************** 35 * Calculates the absolute value of a number. 36 * 37 * Params: 38 * Num = (template parameter) type of number 39 * x = real number value 40 * 41 * Returns: 42 * The absolute value of the number. If floating-point or integral, 43 * the return type will be the same as the input. 44 * 45 * Limitations: 46 * Does not work correctly for signed intergal types and value `Num`.min. 47 */ 48 auto abs(Num)(Num x) @nogc pure nothrow 49 if ((is(immutable Num == immutable short) || is(immutable Num == immutable byte)) || 50 (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)))) 51 { 52 static if (isFloatingPoint!(Num)) 53 return fabs(x); 54 else 55 { 56 static if (is(immutable Num == immutable short) || is(immutable Num == immutable byte)) 57 return x >= 0 ? x : cast(Num) -int(x); 58 else 59 return x >= 0 ? x : -x; 60 } 61 } 62 63 /// ditto 64 @safe pure nothrow @nogc unittest 65 { 66 import std.math.traits : isIdentical, isNaN; 67 68 assert(isIdentical(abs(-0.0L), 0.0L)); 69 assert(isNaN(abs(real.nan))); 70 assert(abs(-real.infinity) == real.infinity); 71 assert(abs(-56) == 56); 72 assert(abs(2321312L) == 2321312L); 73 } 74 75 @safe pure nothrow @nogc unittest 76 { 77 short s = -8; 78 byte b = -8; 79 assert(abs(s) == 8); 80 assert(abs(b) == 8); 81 immutable(byte) c = -8; 82 assert(abs(c) == 8); 83 } 84 85 @safe pure nothrow @nogc unittest 86 { 87 import std.meta : AliasSeq; 88 static foreach (T; AliasSeq!(float, double, real)) 89 {{ 90 T f = 3; 91 assert(abs(f) == f); 92 assert(abs(-f) == f); 93 }} 94 } 95 96 // see https://issues.dlang.org/show_bug.cgi?id=20205 97 // to avoid falling into the trap again 98 @safe pure nothrow @nogc unittest 99 { 100 assert(50 - abs(-100) == -50); 101 } 102 103 // https://issues.dlang.org/show_bug.cgi?id=19162 104 @safe unittest 105 { 106 struct Vector(T, int size) 107 { 108 T x, y, z; 109 } 110 111 static auto abs(T, int size)(auto ref const Vector!(T, size) v) 112 { 113 return v; 114 } 115 Vector!(int, 3) v; 116 assert(abs(v) == v); 117 } 118 119 /******************************* 120 * Returns |x| 121 * 122 * $(TABLE_SV 123 * $(TR $(TH x) $(TH fabs(x))) 124 * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) ) 125 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) ) 126 * ) 127 */ 128 pragma(inline, true) 129 real fabs(real x) @safe pure nothrow @nogc { return core.math.fabs(x); } 130 131 ///ditto 132 pragma(inline, true) 133 double fabs(double x) @safe pure nothrow @nogc { return core.math.fabs(x); } 134 135 ///ditto 136 pragma(inline, true) 137 float fabs(float x) @safe pure nothrow @nogc { return core.math.fabs(x); } 138 139 /// 140 @safe unittest 141 { 142 import std.math.traits : isIdentical; 143 144 assert(isIdentical(fabs(0.0f), 0.0f)); 145 assert(isIdentical(fabs(-0.0f), 0.0f)); 146 assert(fabs(-10.0f) == 10.0f); 147 148 assert(isIdentical(fabs(0.0), 0.0)); 149 assert(isIdentical(fabs(-0.0), 0.0)); 150 assert(fabs(-10.0) == 10.0); 151 152 assert(isIdentical(fabs(0.0L), 0.0L)); 153 assert(isIdentical(fabs(-0.0L), 0.0L)); 154 assert(fabs(-10.0L) == 10.0L); 155 } 156 157 @safe unittest 158 { 159 real function(real) pfabs = &fabs; 160 assert(pfabs != null); 161 } 162 163 @safe pure nothrow @nogc unittest 164 { 165 float f = fabs(-2.0f); 166 assert(f == 2); 167 168 double d = fabs(-2.0); 169 assert(d == 2); 170 171 real r = fabs(-2.0L); 172 assert(r == 2); 173 } 174 175 /*************************************** 176 * Compute square root of x. 177 * 178 * $(TABLE_SV 179 * $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?)) 180 * $(TR $(TD -0.0) $(TD -0.0) $(TD no)) 181 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes)) 182 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no)) 183 * ) 184 */ 185 pragma(inline, true) 186 float sqrt(float x) @nogc @safe pure nothrow { return core.math.sqrt(x); } 187 188 /// ditto 189 pragma(inline, true) 190 double sqrt(double x) @nogc @safe pure nothrow { return core.math.sqrt(x); } 191 192 /// ditto 193 pragma(inline, true) 194 real sqrt(real x) @nogc @safe pure nothrow { return core.math.sqrt(x); } 195 196 /// 197 @safe pure nothrow @nogc unittest 198 { 199 import std.math.operations : feqrel; 200 import std.math.traits : isNaN; 201 202 assert(sqrt(2.0).feqrel(1.4142) > 16); 203 assert(sqrt(9.0).feqrel(3.0) > 16); 204 205 assert(isNaN(sqrt(-1.0f))); 206 assert(isNaN(sqrt(-1.0))); 207 assert(isNaN(sqrt(-1.0L))); 208 } 209 210 @safe unittest 211 { 212 // https://issues.dlang.org/show_bug.cgi?id=5305 213 float function(float) psqrtf = &sqrt; 214 assert(psqrtf != null); 215 double function(double) psqrtd = &sqrt; 216 assert(psqrtd != null); 217 real function(real) psqrtr = &sqrt; 218 assert(psqrtr != null); 219 220 //ctfe 221 enum ZX80 = sqrt(7.0f); 222 enum ZX81 = sqrt(7.0); 223 enum ZX82 = sqrt(7.0L); 224 } 225 226 @safe pure nothrow @nogc unittest 227 { 228 float f = sqrt(2.0f); 229 assert(fabs(f * f - 2.0f) < .00001); 230 231 double d = sqrt(2.0); 232 assert(fabs(d * d - 2.0) < .00001); 233 234 real r = sqrt(2.0L); 235 assert(fabs(r * r - 2.0) < .00001); 236 } 237 238 /*************** 239 * Calculates the cube root of x. 240 * 241 * $(TABLE_SV 242 * $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?)) 243 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) ) 244 * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) 245 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) ) 246 * ) 247 */ 248 real cbrt(real x) @trusted nothrow @nogc 249 { 250 version (CRuntime_Microsoft) 251 { 252 import std.math.traits : copysign; 253 import std.math.exponential : exp2; 254 255 version (INLINE_YL2X) 256 return copysign(exp2(core.math.yl2x(fabs(x), 1.0L/3.0L)), x); 257 else 258 return core.stdc.math.cbrtl(x); 259 } 260 else 261 return core.stdc.math.cbrtl(x); 262 } 263 264 /// 265 @safe unittest 266 { 267 import std.math.operations : feqrel; 268 269 assert(cbrt(1.0).feqrel(1.0) > 16); 270 assert(cbrt(27.0).feqrel(3.0) > 16); 271 assert(cbrt(15.625).feqrel(2.5) > 16); 272 } 273 274 /*********************************************************************** 275 * Calculates the length of the 276 * hypotenuse of a right-angled triangle with sides of length x and y. 277 * The hypotenuse is the value of the square root of 278 * the sums of the squares of x and y: 279 * 280 * sqrt($(POWER x, 2) + $(POWER y, 2)) 281 * 282 * Note that hypot(x, y), hypot(y, x) and 283 * hypot(x, -y) are equivalent. 284 * 285 * $(TABLE_SV 286 * $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?)) 287 * $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no)) 288 * $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no)) 289 * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no)) 290 * ) 291 */ 292 T hypot(T)(const T x, const T y) @safe pure nothrow @nogc 293 if (isFloatingPoint!T) 294 { 295 // Scale x and y to avoid underflow and overflow. 296 // If one is huge and the other tiny, return the larger. 297 // If both are huge, avoid overflow by scaling by 2^^-N. 298 // If both are tiny, avoid underflow by scaling by 2^^N. 299 import core.math : fabs, sqrt; 300 import std.math : floatTraits, RealFormat; 301 302 alias F = floatTraits!T; 303 304 T u = fabs(x); 305 T v = fabs(y); 306 if (!(u >= v)) // check for NaN as well. 307 { 308 v = u; 309 u = fabs(y); 310 if (u == T.infinity) return u; // hypot(inf, nan) == inf 311 if (v == T.infinity) return v; // hypot(nan, inf) == inf 312 } 313 314 static if (F.realFormat == RealFormat.ieeeSingle) 315 { 316 enum SQRTMIN = 0x1p-60f; 317 enum SQRTMAX = 0x1p+60f; 318 enum SCALE_UNDERFLOW = 0x1p+90f; 319 enum SCALE_OVERFLOW = 0x1p-90f; 320 } 321 else static if (F.realFormat == RealFormat.ieeeDouble || 322 F.realFormat == RealFormat.ieeeExtended53 || 323 F.realFormat == RealFormat.ibmExtended) 324 { 325 enum SQRTMIN = 0x1p-450L; 326 enum SQRTMAX = 0x1p+500L; 327 enum SCALE_UNDERFLOW = 0x1p+600L; 328 enum SCALE_OVERFLOW = 0x1p-600L; 329 } 330 else static if (F.realFormat == RealFormat.ieeeExtended || 331 F.realFormat == RealFormat.ieeeQuadruple) 332 { 333 enum SQRTMIN = 0x1p-8000L; 334 enum SQRTMAX = 0x1p+8000L; 335 enum SCALE_UNDERFLOW = 0x1p+10000L; 336 enum SCALE_OVERFLOW = 0x1p-10000L; 337 } 338 else 339 assert(0, "hypot not implemented"); 340 341 // Now u >= v, or else one is NaN. 342 T ratio = 1.0; 343 if (v >= SQRTMAX) 344 { 345 // hypot(huge, huge) -- avoid overflow 346 ratio = SCALE_UNDERFLOW; 347 u *= SCALE_OVERFLOW; 348 v *= SCALE_OVERFLOW; 349 } 350 else if (u <= SQRTMIN) 351 { 352 // hypot (tiny, tiny) -- avoid underflow 353 // This is only necessary to avoid setting the underflow 354 // flag. 355 ratio = SCALE_OVERFLOW; 356 u *= SCALE_UNDERFLOW; 357 v *= SCALE_UNDERFLOW; 358 } 359 360 if (u * T.epsilon > v) 361 { 362 // hypot (huge, tiny) = huge 363 return u; 364 } 365 366 // both are in the normal range 367 return ratio * sqrt(u*u + v*v); 368 } 369 370 /// 371 @safe unittest 372 { 373 import std.math.operations : feqrel; 374 375 assert(hypot(1.0, 1.0).feqrel(1.4142) > 16); 376 assert(hypot(3.0, 4.0).feqrel(5.0) > 16); 377 assert(hypot(real.infinity, 1.0L) == real.infinity); 378 assert(hypot(real.infinity, real.nan) == real.infinity); 379 } 380 381 @safe unittest 382 { 383 import std.math.operations : feqrel; 384 385 assert(hypot(1.0f, 1.0f).feqrel(1.4142f) > 16); 386 assert(hypot(3.0f, 4.0f).feqrel(5.0f) > 16); 387 assert(hypot(float.infinity, 1.0f) == float.infinity); 388 assert(hypot(float.infinity, float.nan) == float.infinity); 389 390 assert(hypot(1.0L, 1.0L).feqrel(1.4142L) > 16); 391 assert(hypot(3.0L, 4.0L).feqrel(5.0L) > 16); 392 assert(hypot(double.infinity, 1.0) == double.infinity); 393 assert(hypot(double.infinity, double.nan) == double.infinity); 394 } 395 396 @safe unittest 397 { 398 import std.math.operations : feqrel; 399 import std.math.traits : isIdentical; 400 import std.meta : AliasSeq; 401 402 static foreach (T; AliasSeq!(float, double, real)) 403 {{ 404 static T[3][] vals = // x,y,hypot 405 [ 406 [ 0.0, 0.0, 0.0], 407 [ 0.0, -0.0, 0.0], 408 [ -0.0, -0.0, 0.0], 409 [ 3.0, 4.0, 5.0], 410 [ -300, -400, 500], 411 [0.0, 7.0, 7.0], 412 [9.0, 9*T.epsilon, 9.0], 413 [88/(64*sqrt(T.min_normal)), 105/(64*sqrt(T.min_normal)), 137/(64*sqrt(T.min_normal))], 414 [88/(128*sqrt(T.min_normal)), 105/(128*sqrt(T.min_normal)), 137/(128*sqrt(T.min_normal))], 415 [3*T.min_normal*T.epsilon, 4*T.min_normal*T.epsilon, 5*T.min_normal*T.epsilon], 416 [ T.min_normal, T.min_normal, sqrt(2.0L)*T.min_normal], 417 [ T.max/sqrt(2.0L), T.max/sqrt(2.0L), T.max], 418 [ T.infinity, T.nan, T.infinity], 419 [ T.nan, T.infinity, T.infinity], 420 [ T.nan, T.nan, T.nan], 421 [ T.nan, T.max, T.nan], 422 [ T.max, T.nan, T.nan], 423 ]; 424 for (int i = 0; i < vals.length; i++) 425 { 426 T x = vals[i][0]; 427 T y = vals[i][1]; 428 T z = vals[i][2]; 429 T h = hypot(x, y); 430 assert(isIdentical(z,h) || feqrel(z, h) >= T.mant_dig - 1); 431 } 432 }} 433 } 434 435 /*********************************************************************** 436 * Calculates the distance of the point (x, y, z) from the origin (0, 0, 0) 437 * in three-dimensional space. 438 * The distance is the value of the square root of the sums of the squares 439 * of x, y, and z: 440 * 441 * sqrt($(POWER x, 2) + $(POWER y, 2) + $(POWER z, 2)) 442 * 443 * Note that the distance between two points (x1, y1, z1) and (x2, y2, z2) 444 * in three-dimensional space can be calculated as hypot(x2-x1, y2-y1, z2-z1). 445 * 446 * Params: 447 * x = floating point value 448 * y = floating point value 449 * z = floating point value 450 * 451 * Returns: 452 * The square root of the sum of the squares of the given arguments. 453 */ 454 T hypot(T)(const T x, const T y, const T z) @safe pure nothrow @nogc 455 if (isFloatingPoint!T) 456 { 457 import core.math : fabs, sqrt; 458 import std.math.operations : fmax; 459 const absx = fabs(x); 460 const absy = fabs(y); 461 const absz = fabs(z); 462 463 // Scale all parameters to avoid overflow. 464 const ratio = fmax(absx, fmax(absy, absz)); 465 if (ratio == 0.0) 466 return ratio; 467 468 return ratio * sqrt((absx / ratio) * (absx / ratio) 469 + (absy / ratio) * (absy / ratio) 470 + (absz / ratio) * (absz / ratio)); 471 } 472 473 /// 474 @safe unittest 475 { 476 import std.math.operations : isClose; 477 478 assert(isClose(hypot(1.0, 2.0, 2.0), 3.0)); 479 assert(isClose(hypot(2.0, 3.0, 6.0), 7.0)); 480 assert(isClose(hypot(1.0, 4.0, 8.0), 9.0)); 481 } 482 483 @safe unittest 484 { 485 import std.meta : AliasSeq; 486 import std.math.traits : isIdentical; 487 import std.math.operations : isClose; 488 static foreach (T; AliasSeq!(float, double, real)) 489 {{ 490 static T[4][] vals = [ 491 [ 0.0L, 0.0L, 0.0L, 0.0L ], 492 [ 0.0L, 1.0L, 1.0L, sqrt(2.0L) ], 493 [ 1.0L, 1.0L, 1.0L, sqrt(3.0L) ], 494 [ 1.0L, 2.0L, 2.0L, 3.0L ], 495 [ 2.0L, 3.0L, 6.0L, 7.0L ], 496 [ 1.0L, 4.0L, 8.0L, 9.0L ], 497 [ 4.0L, 4.0L, 7.0L, 9.0L ], 498 [ 12.0L, 16.0L, 21.0L, 29.0L ], 499 [ 1e+8L, 1.0L, 1e-8L, 1e+8L+5e-9L ], 500 [ 1.0L, 1e+8L, 1e-8L, 1e+8L+5e-9L ], 501 [ 1e-8L, 1.0L, 1e+8L, 1e+8L+5e-9L ], 502 [ 1e-2L, 1e-4L, 1e-4L, 0.010000999950004999375L ], 503 [ 2147483647.0L, 2147483647.0L, 2147483647.0L, 3719550785.027307813987L ] 504 ]; 505 for (int i = 0; i < vals.length; i++) 506 { 507 T x = vals[i][0]; 508 T y = vals[i][1]; 509 T z = vals[i][2]; 510 T r = vals[i][3]; 511 T a = hypot(x, y, z); 512 assert(isIdentical(r, a) || isClose(r, a)); 513 } 514 }} 515 } 516 517 /*********************************** 518 * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2) + 519 * $(SUB a,3)$(POWER x,3); ... 520 * 521 * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2) + 522 * x($(SUB a, 3) + ...))) 523 * Params: 524 * x = the value to evaluate. 525 * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc. 526 */ 527 Unqual!(CommonType!(T1, T2)) poly(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc 528 if (isFloatingPoint!T1 && isFloatingPoint!T2) 529 in 530 { 531 assert(A.length > 0); 532 } 533 do 534 { 535 static if (is(immutable T2 == immutable real)) 536 { 537 return polyImpl(x, A); 538 } 539 else 540 { 541 return polyImplBase(x, A); 542 } 543 } 544 545 /// ditto 546 Unqual!(CommonType!(T1, T2)) poly(T1, T2, int N)(T1 x, ref const T2[N] A) @safe pure nothrow @nogc 547 if (isFloatingPoint!T1 && isFloatingPoint!T2 && N > 0 && N <= 10) 548 { 549 // statically unrolled version for up to 10 coefficients 550 typeof(return) r = A[N - 1]; 551 static foreach (i; 1 .. N) 552 { 553 r *= x; 554 r += A[N - 1 - i]; 555 } 556 return r; 557 } 558 559 /// 560 @safe nothrow @nogc unittest 561 { 562 real x = 3.1L; 563 static real[] pp = [56.1L, 32.7L, 6]; 564 565 assert(poly(x, pp) == (56.1L + (32.7L + 6.0L * x) * x)); 566 } 567 568 @safe nothrow @nogc unittest 569 { 570 double x = 3.1; 571 static double[] pp = [56.1, 32.7, 6]; 572 double y = x; 573 y *= 6.0; 574 y += 32.7; 575 y *= x; 576 y += 56.1; 577 assert(poly(x, pp) == y); 578 } 579 580 @safe unittest 581 { 582 static assert(poly(3.0, [1.0, 2.0, 3.0]) == 34); 583 } 584 585 private Unqual!(CommonType!(T1, T2)) polyImplBase(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc 586 if (isFloatingPoint!T1 && isFloatingPoint!T2) 587 { 588 ptrdiff_t i = A.length - 1; 589 typeof(return) r = A[i]; 590 while (--i >= 0) 591 { 592 r *= x; 593 r += A[i]; 594 } 595 return r; 596 } 597 598 version (linux) version = GenericPosixVersion; 599 else version (FreeBSD) version = GenericPosixVersion; 600 else version (OpenBSD) version = GenericPosixVersion; 601 else version (Solaris) version = GenericPosixVersion; 602 else version (DragonFlyBSD) version = GenericPosixVersion; 603 604 private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc 605 { 606 version (D_InlineAsm_X86) 607 { 608 if (__ctfe) 609 { 610 return polyImplBase(x, A); 611 } 612 version (Windows) 613 { 614 // BUG: This code assumes a frame pointer in EBP. 615 asm pure nothrow @nogc // assembler by W. Bright 616 { 617 // EDX = (A.length - 1) * real.sizeof 618 mov ECX,A[EBP] ; // ECX = A.length 619 dec ECX ; 620 lea EDX,[ECX][ECX*8] ; 621 add EDX,ECX ; 622 add EDX,A+4[EBP] ; 623 fld real ptr [EDX] ; // ST0 = coeff[ECX] 624 jecxz return_ST ; 625 fld x[EBP] ; // ST0 = x 626 fxch ST(1) ; // ST1 = x, ST0 = r 627 align 4 ; 628 L2: fmul ST,ST(1) ; // r *= x 629 fld real ptr -10[EDX] ; 630 sub EDX,10 ; // deg-- 631 faddp ST(1),ST ; 632 dec ECX ; 633 jne L2 ; 634 fxch ST(1) ; // ST1 = r, ST0 = x 635 fstp ST(0) ; // dump x 636 align 4 ; 637 return_ST: ; 638 } 639 } 640 else version (GenericPosixVersion) 641 { 642 asm pure nothrow @nogc // assembler by W. Bright 643 { 644 // EDX = (A.length - 1) * real.sizeof 645 mov ECX,A[EBP] ; // ECX = A.length 646 dec ECX ; 647 lea EDX,[ECX*8] ; 648 lea EDX,[EDX][ECX*4] ; 649 add EDX,A+4[EBP] ; 650 fld real ptr [EDX] ; // ST0 = coeff[ECX] 651 jecxz return_ST ; 652 fld x[EBP] ; // ST0 = x 653 fxch ST(1) ; // ST1 = x, ST0 = r 654 align 4 ; 655 L2: fmul ST,ST(1) ; // r *= x 656 fld real ptr -12[EDX] ; 657 sub EDX,12 ; // deg-- 658 faddp ST(1),ST ; 659 dec ECX ; 660 jne L2 ; 661 fxch ST(1) ; // ST1 = r, ST0 = x 662 fstp ST(0) ; // dump x 663 align 4 ; 664 return_ST: ; 665 } 666 } 667 else version (OSX) 668 { 669 asm pure nothrow @nogc // assembler by W. Bright 670 { 671 // EDX = (A.length - 1) * real.sizeof 672 mov ECX,A[EBP] ; // ECX = A.length 673 dec ECX ; 674 lea EDX,[ECX*8] ; 675 add EDX,EDX ; 676 add EDX,A+4[EBP] ; 677 fld real ptr [EDX] ; // ST0 = coeff[ECX] 678 jecxz return_ST ; 679 fld x[EBP] ; // ST0 = x 680 fxch ST(1) ; // ST1 = x, ST0 = r 681 align 4 ; 682 L2: fmul ST,ST(1) ; // r *= x 683 fld real ptr -16[EDX] ; 684 sub EDX,16 ; // deg-- 685 faddp ST(1),ST ; 686 dec ECX ; 687 jne L2 ; 688 fxch ST(1) ; // ST1 = r, ST0 = x 689 fstp ST(0) ; // dump x 690 align 4 ; 691 return_ST: ; 692 } 693 } 694 else 695 { 696 static assert(0); 697 } 698 } 699 else 700 { 701 return polyImplBase(x, A); 702 } 703 } 704 705 /** 706 * Gives the next power of two after `val`. `T` can be any built-in 707 * numerical type. 708 * 709 * If the operation would lead to an over/underflow, this function will 710 * return `0`. 711 * 712 * Params: 713 * val = any number 714 * 715 * Returns: 716 * the next power of two after `val` 717 */ 718 T nextPow2(T)(const T val) 719 if (isIntegral!T) 720 { 721 return powIntegralImpl!(PowType.ceil)(val); 722 } 723 724 /// ditto 725 T nextPow2(T)(const T val) 726 if (isFloatingPoint!T) 727 { 728 return powFloatingPointImpl!(PowType.ceil)(val); 729 } 730 731 /// 732 @safe @nogc pure nothrow unittest 733 { 734 assert(nextPow2(2) == 4); 735 assert(nextPow2(10) == 16); 736 assert(nextPow2(4000) == 4096); 737 738 assert(nextPow2(-2) == -4); 739 assert(nextPow2(-10) == -16); 740 741 assert(nextPow2(uint.max) == 0); 742 assert(nextPow2(uint.min) == 0); 743 assert(nextPow2(size_t.max) == 0); 744 assert(nextPow2(size_t.min) == 0); 745 746 assert(nextPow2(int.max) == 0); 747 assert(nextPow2(int.min) == 0); 748 assert(nextPow2(long.max) == 0); 749 assert(nextPow2(long.min) == 0); 750 } 751 752 /// 753 @safe @nogc pure nothrow unittest 754 { 755 assert(nextPow2(2.1) == 4.0); 756 assert(nextPow2(-2.0) == -4.0); 757 assert(nextPow2(0.25) == 0.5); 758 assert(nextPow2(-4.0) == -8.0); 759 760 assert(nextPow2(double.max) == 0.0); 761 assert(nextPow2(double.infinity) == double.infinity); 762 } 763 764 @safe @nogc pure nothrow unittest 765 { 766 assert(nextPow2(ubyte(2)) == 4); 767 assert(nextPow2(ubyte(10)) == 16); 768 769 assert(nextPow2(byte(2)) == 4); 770 assert(nextPow2(byte(10)) == 16); 771 772 assert(nextPow2(short(2)) == 4); 773 assert(nextPow2(short(10)) == 16); 774 assert(nextPow2(short(4000)) == 4096); 775 776 assert(nextPow2(ushort(2)) == 4); 777 assert(nextPow2(ushort(10)) == 16); 778 assert(nextPow2(ushort(4000)) == 4096); 779 } 780 781 @safe @nogc pure nothrow unittest 782 { 783 foreach (ulong i; 1 .. 62) 784 { 785 assert(nextPow2(1UL << i) == 2UL << i); 786 assert(nextPow2((1UL << i) - 1) == 1UL << i); 787 assert(nextPow2((1UL << i) + 1) == 2UL << i); 788 assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i); 789 } 790 } 791 792 @safe @nogc pure nothrow unittest 793 { 794 import std.math.traits : isNaN; 795 import std.meta : AliasSeq; 796 797 static foreach (T; AliasSeq!(float, double, real)) 798 {{ 799 enum T subNormal = T.min_normal / 2; 800 801 static if (subNormal) assert(nextPow2(subNormal) == T.min_normal); 802 803 assert(nextPow2(T(0.0)) == 0.0); 804 805 assert(nextPow2(T(2.0)) == 4.0); 806 assert(nextPow2(T(2.1)) == 4.0); 807 assert(nextPow2(T(3.1)) == 4.0); 808 assert(nextPow2(T(4.0)) == 8.0); 809 assert(nextPow2(T(0.25)) == 0.5); 810 811 assert(nextPow2(T(-2.0)) == -4.0); 812 assert(nextPow2(T(-2.1)) == -4.0); 813 assert(nextPow2(T(-3.1)) == -4.0); 814 assert(nextPow2(T(-4.0)) == -8.0); 815 assert(nextPow2(T(-0.25)) == -0.5); 816 817 assert(nextPow2(T.max) == 0); 818 assert(nextPow2(-T.max) == 0); 819 820 assert(nextPow2(T.infinity) == T.infinity); 821 assert(nextPow2(T.init).isNaN); 822 }} 823 } 824 825 // https://issues.dlang.org/show_bug.cgi?id=15973 826 @safe @nogc pure nothrow unittest 827 { 828 assert(nextPow2(uint.max / 2) == uint.max / 2 + 1); 829 assert(nextPow2(uint.max / 2 + 2) == 0); 830 assert(nextPow2(int.max / 2) == int.max / 2 + 1); 831 assert(nextPow2(int.max / 2 + 2) == 0); 832 assert(nextPow2(int.min + 1) == int.min); 833 } 834 835 /** 836 * Gives the last power of two before `val`. $(T) can be any built-in 837 * numerical type. 838 * 839 * Params: 840 * val = any number 841 * 842 * Returns: 843 * the last power of two before `val` 844 */ 845 T truncPow2(T)(const T val) 846 if (isIntegral!T) 847 { 848 return powIntegralImpl!(PowType.floor)(val); 849 } 850 851 /// ditto 852 T truncPow2(T)(const T val) 853 if (isFloatingPoint!T) 854 { 855 return powFloatingPointImpl!(PowType.floor)(val); 856 } 857 858 /// 859 @safe @nogc pure nothrow unittest 860 { 861 assert(truncPow2(3) == 2); 862 assert(truncPow2(4) == 4); 863 assert(truncPow2(10) == 8); 864 assert(truncPow2(4000) == 2048); 865 866 assert(truncPow2(-5) == -4); 867 assert(truncPow2(-20) == -16); 868 869 assert(truncPow2(uint.max) == int.max + 1); 870 assert(truncPow2(uint.min) == 0); 871 assert(truncPow2(ulong.max) == long.max + 1); 872 assert(truncPow2(ulong.min) == 0); 873 874 assert(truncPow2(int.max) == (int.max / 2) + 1); 875 assert(truncPow2(int.min) == int.min); 876 assert(truncPow2(long.max) == (long.max / 2) + 1); 877 assert(truncPow2(long.min) == long.min); 878 } 879 880 /// 881 @safe @nogc pure nothrow unittest 882 { 883 assert(truncPow2(2.1) == 2.0); 884 assert(truncPow2(7.0) == 4.0); 885 assert(truncPow2(-1.9) == -1.0); 886 assert(truncPow2(0.24) == 0.125); 887 assert(truncPow2(-7.0) == -4.0); 888 889 assert(truncPow2(double.infinity) == double.infinity); 890 } 891 892 @safe @nogc pure nothrow unittest 893 { 894 assert(truncPow2(ubyte(3)) == 2); 895 assert(truncPow2(ubyte(4)) == 4); 896 assert(truncPow2(ubyte(10)) == 8); 897 898 assert(truncPow2(byte(3)) == 2); 899 assert(truncPow2(byte(4)) == 4); 900 assert(truncPow2(byte(10)) == 8); 901 902 assert(truncPow2(ushort(3)) == 2); 903 assert(truncPow2(ushort(4)) == 4); 904 assert(truncPow2(ushort(10)) == 8); 905 assert(truncPow2(ushort(4000)) == 2048); 906 907 assert(truncPow2(short(3)) == 2); 908 assert(truncPow2(short(4)) == 4); 909 assert(truncPow2(short(10)) == 8); 910 assert(truncPow2(short(4000)) == 2048); 911 } 912 913 @safe @nogc pure nothrow unittest 914 { 915 foreach (ulong i; 1 .. 62) 916 { 917 assert(truncPow2(2UL << i) == 2UL << i); 918 assert(truncPow2((2UL << i) + 1) == 2UL << i); 919 assert(truncPow2((2UL << i) - 1) == 1UL << i); 920 assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i); 921 } 922 } 923 924 @safe @nogc pure nothrow unittest 925 { 926 import std.math.traits : isNaN; 927 import std.meta : AliasSeq; 928 929 static foreach (T; AliasSeq!(float, double, real)) 930 { 931 assert(truncPow2(T(0.0)) == 0.0); 932 933 assert(truncPow2(T(4.0)) == 4.0); 934 assert(truncPow2(T(2.1)) == 2.0); 935 assert(truncPow2(T(3.5)) == 2.0); 936 assert(truncPow2(T(7.0)) == 4.0); 937 assert(truncPow2(T(0.24)) == 0.125); 938 939 assert(truncPow2(T(-2.0)) == -2.0); 940 assert(truncPow2(T(-2.1)) == -2.0); 941 assert(truncPow2(T(-3.1)) == -2.0); 942 assert(truncPow2(T(-7.0)) == -4.0); 943 assert(truncPow2(T(-0.24)) == -0.125); 944 945 assert(truncPow2(T.infinity) == T.infinity); 946 assert(truncPow2(T.init).isNaN); 947 } 948 } 949 950 private enum PowType 951 { 952 floor, 953 ceil 954 } 955 956 pragma(inline, true) 957 private T powIntegralImpl(PowType type, T)(T val) 958 { 959 import core.bitop : bsr; 960 961 if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min))) 962 return 0; 963 else 964 { 965 static if (isSigned!T) 966 return cast(Unqual!T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type); 967 else 968 return cast(Unqual!T) (T(1) << bsr(val) + type); 969 } 970 } 971 972 private T powFloatingPointImpl(PowType type, T)(T x) 973 { 974 import std.math.traits : copysign, isFinite; 975 import std.math.exponential : frexp; 976 977 if (!x.isFinite) 978 return x; 979 980 if (!x) 981 return x; 982 983 int exp; 984 auto y = frexp(x, exp); 985 986 static if (type == PowType.ceil) 987 y = core.math.ldexp(cast(T) 0.5, exp + 1); 988 else 989 y = core.math.ldexp(cast(T) 0.5, exp); 990 991 if (!y.isFinite) 992 return cast(T) 0.0; 993 994 y = copysign(y, x); 995 996 return y; 997 }