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