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 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 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 * ) 302 */ 303 T hypot(T)(const T x, const T y) @safe pure nothrow @nogc 304 if (isFloatingPoint!T) 305 { 306 // Scale x and y to avoid underflow and overflow. 307 // If one is huge and the other tiny, return the larger. 308 // If both are huge, avoid overflow by scaling by 2^^-N. 309 // If both are tiny, avoid underflow by scaling by 2^^N. 310 import core.math : fabs, sqrt; 311 import std.math.traits : floatTraits, RealFormat; 312 313 alias F = floatTraits!T; 314 315 T u = fabs(x); 316 T v = fabs(y); 317 if (!(u >= v)) // check for NaN as well. 318 { 319 v = u; 320 u = fabs(y); 321 if (u == T.infinity) return u; // hypot(inf, nan) == inf 322 if (v == T.infinity) return v; // hypot(nan, inf) == inf 323 } 324 325 static if (F.realFormat == RealFormat.ieeeSingle) 326 { 327 enum SQRTMIN = 0x1p-60f; 328 enum SQRTMAX = 0x1p+60f; 329 enum SCALE_UNDERFLOW = 0x1p+90f; 330 enum SCALE_OVERFLOW = 0x1p-90f; 331 } 332 else static if (F.realFormat == RealFormat.ieeeDouble || 333 F.realFormat == RealFormat.ieeeExtended53 || 334 F.realFormat == RealFormat.ibmExtended) 335 { 336 enum SQRTMIN = 0x1p-450L; 337 enum SQRTMAX = 0x1p+500L; 338 enum SCALE_UNDERFLOW = 0x1p+600L; 339 enum SCALE_OVERFLOW = 0x1p-600L; 340 } 341 else static if (F.realFormat == RealFormat.ieeeExtended || 342 F.realFormat == RealFormat.ieeeQuadruple) 343 { 344 enum SQRTMIN = 0x1p-8000L; 345 enum SQRTMAX = 0x1p+8000L; 346 enum SCALE_UNDERFLOW = 0x1p+10000L; 347 enum SCALE_OVERFLOW = 0x1p-10000L; 348 } 349 else 350 assert(0, "hypot not implemented"); 351 352 // Now u >= v, or else one is NaN. 353 T ratio = 1.0; 354 if (v >= SQRTMAX) 355 { 356 // hypot(huge, huge) -- avoid overflow 357 ratio = SCALE_UNDERFLOW; 358 u *= SCALE_OVERFLOW; 359 v *= SCALE_OVERFLOW; 360 } 361 else if (u <= SQRTMIN) 362 { 363 // hypot (tiny, tiny) -- avoid underflow 364 // This is only necessary to avoid setting the underflow 365 // flag. 366 ratio = SCALE_OVERFLOW; 367 u *= SCALE_UNDERFLOW; 368 v *= SCALE_UNDERFLOW; 369 } 370 371 if (u * T.epsilon > v) 372 { 373 // hypot (huge, tiny) = huge 374 return u; 375 } 376 377 // both are in the normal range 378 return ratio * sqrt(u*u + v*v); 379 } 380 381 /// 382 @safe unittest 383 { 384 import std.math.operations : feqrel; 385 386 assert(hypot(1.0, 1.0).feqrel(1.4142) > 16); 387 assert(hypot(3.0, 4.0).feqrel(5.0) > 16); 388 assert(hypot(real.infinity, 1.0L) == real.infinity); 389 assert(hypot(real.infinity, real.nan) == real.infinity); 390 } 391 392 @safe unittest 393 { 394 import std.math.operations : feqrel; 395 396 assert(hypot(1.0f, 1.0f).feqrel(1.4142f) > 16); 397 assert(hypot(3.0f, 4.0f).feqrel(5.0f) > 16); 398 assert(hypot(float.infinity, 1.0f) == float.infinity); 399 assert(hypot(float.infinity, float.nan) == float.infinity); 400 401 assert(hypot(1.0L, 1.0L).feqrel(1.4142L) > 16); 402 assert(hypot(3.0L, 4.0L).feqrel(5.0L) > 16); 403 assert(hypot(double.infinity, 1.0) == double.infinity); 404 assert(hypot(double.infinity, double.nan) == double.infinity); 405 } 406 407 @safe unittest 408 { 409 import std.math.operations : feqrel; 410 import std.math.traits : isIdentical; 411 import std.meta : AliasSeq; 412 413 static foreach (T; AliasSeq!(float, double, real)) 414 {{ 415 static T[3][] vals = // x,y,hypot 416 [ 417 [ 0.0, 0.0, 0.0], 418 [ 0.0, -0.0, 0.0], 419 [ -0.0, -0.0, 0.0], 420 [ 3.0, 4.0, 5.0], 421 [ -300, -400, 500], 422 [0.0, 7.0, 7.0], 423 [9.0, 9*T.epsilon, 9.0], 424 [88/(64*sqrt(T.min_normal)), 105/(64*sqrt(T.min_normal)), 137/(64*sqrt(T.min_normal))], 425 [88/(128*sqrt(T.min_normal)), 105/(128*sqrt(T.min_normal)), 137/(128*sqrt(T.min_normal))], 426 [3*T.min_normal*T.epsilon, 4*T.min_normal*T.epsilon, 5*T.min_normal*T.epsilon], 427 [ T.min_normal, T.min_normal, sqrt(2.0L)*T.min_normal], 428 [ T.max/sqrt(2.0L), T.max/sqrt(2.0L), T.max], 429 [ T.infinity, T.nan, T.infinity], 430 [ T.nan, T.infinity, T.infinity], 431 [ T.nan, T.nan, T.nan], 432 [ T.nan, T.max, T.nan], 433 [ T.max, T.nan, T.nan], 434 ]; 435 for (int i = 0; i < vals.length; i++) 436 { 437 T x = vals[i][0]; 438 T y = vals[i][1]; 439 T z = vals[i][2]; 440 T h = hypot(x, y); 441 assert(isIdentical(z,h) || feqrel(z, h) >= T.mant_dig - 1); 442 } 443 }} 444 } 445 446 /*********************************************************************** 447 * Calculates the distance of the point (x, y, z) from the origin (0, 0, 0) 448 * in three-dimensional space. 449 * The distance is the value of the square root of the sums of the squares 450 * of x, y, and z: 451 * 452 * sqrt($(POWER x, 2) + $(POWER y, 2) + $(POWER z, 2)) 453 * 454 * Note that the distance between two points (x1, y1, z1) and (x2, y2, z2) 455 * in three-dimensional space can be calculated as hypot(x2-x1, y2-y1, z2-z1). 456 * 457 * Params: 458 * x = floating point value 459 * y = floating point value 460 * z = floating point value 461 * 462 * Returns: 463 * The square root of the sum of the squares of the given arguments. 464 */ 465 T hypot(T)(const T x, const T y, const T z) @safe pure nothrow @nogc 466 if (isFloatingPoint!T) 467 { 468 import core.math : fabs, sqrt; 469 import std.math.operations : fmax; 470 const absx = fabs(x); 471 const absy = fabs(y); 472 const absz = fabs(z); 473 474 // Scale all parameters to avoid overflow. 475 const ratio = fmax(absx, fmax(absy, absz)); 476 if (ratio == 0.0) 477 return ratio; 478 479 return ratio * sqrt((absx / ratio) * (absx / ratio) 480 + (absy / ratio) * (absy / ratio) 481 + (absz / ratio) * (absz / ratio)); 482 } 483 484 /// 485 @safe unittest 486 { 487 import std.math.operations : isClose; 488 489 assert(isClose(hypot(1.0, 2.0, 2.0), 3.0)); 490 assert(isClose(hypot(2.0, 3.0, 6.0), 7.0)); 491 assert(isClose(hypot(1.0, 4.0, 8.0), 9.0)); 492 } 493 494 @safe unittest 495 { 496 import std.meta : AliasSeq; 497 import std.math.traits : isIdentical; 498 import std.math.operations : isClose; 499 static foreach (T; AliasSeq!(float, double, real)) 500 {{ 501 static T[4][] vals = [ 502 [ 0.0L, 0.0L, 0.0L, 0.0L ], 503 [ 0.0L, 1.0L, 1.0L, sqrt(2.0L) ], 504 [ 1.0L, 1.0L, 1.0L, sqrt(3.0L) ], 505 [ 1.0L, 2.0L, 2.0L, 3.0L ], 506 [ 2.0L, 3.0L, 6.0L, 7.0L ], 507 [ 1.0L, 4.0L, 8.0L, 9.0L ], 508 [ 4.0L, 4.0L, 7.0L, 9.0L ], 509 [ 12.0L, 16.0L, 21.0L, 29.0L ], 510 [ 1e+8L, 1.0L, 1e-8L, 1e+8L+5e-9L ], 511 [ 1.0L, 1e+8L, 1e-8L, 1e+8L+5e-9L ], 512 [ 1e-8L, 1.0L, 1e+8L, 1e+8L+5e-9L ], 513 [ 1e-2L, 1e-4L, 1e-4L, 0.010000999950004999375L ], 514 [ 2147483647.0L, 2147483647.0L, 2147483647.0L, 3719550785.027307813987L ] 515 ]; 516 for (int i = 0; i < vals.length; i++) 517 { 518 T x = vals[i][0]; 519 T y = vals[i][1]; 520 T z = vals[i][2]; 521 T r = vals[i][3]; 522 T a = hypot(x, y, z); 523 assert(isIdentical(r, a) || isClose(r, a)); 524 } 525 }} 526 } 527 528 /*********************************** 529 * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2) + 530 * $(SUB a,3)$(POWER x,3); ... 531 * 532 * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2) + 533 * x($(SUB a, 3) + ...))) 534 * Params: 535 * x = the value to evaluate. 536 * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc. 537 */ 538 Unqual!(CommonType!(T1, T2)) poly(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc 539 if (isFloatingPoint!T1 && isFloatingPoint!T2) 540 in 541 { 542 assert(A.length > 0); 543 } 544 do 545 { 546 static if (is(immutable T2 == immutable real)) 547 { 548 return polyImpl(x, A); 549 } 550 else 551 { 552 return polyImplBase(x, A); 553 } 554 } 555 556 /// ditto 557 Unqual!(CommonType!(T1, T2)) poly(T1, T2, int N)(T1 x, ref const T2[N] A) @safe pure nothrow @nogc 558 if (isFloatingPoint!T1 && isFloatingPoint!T2 && N > 0 && N <= 10) 559 { 560 // statically unrolled version for up to 10 coefficients 561 typeof(return) r = A[N - 1]; 562 static foreach (i; 1 .. N) 563 { 564 r *= x; 565 r += A[N - 1 - i]; 566 } 567 return r; 568 } 569 570 /// 571 @safe nothrow @nogc unittest 572 { 573 real x = 3.1L; 574 static real[] pp = [56.1L, 32.7L, 6]; 575 576 assert(poly(x, pp) == (56.1L + (32.7L + 6.0L * x) * x)); 577 } 578 579 @safe nothrow @nogc unittest 580 { 581 double x = 3.1; 582 static double[] pp = [56.1, 32.7, 6]; 583 double y = x; 584 y *= 6.0; 585 y += 32.7; 586 y *= x; 587 y += 56.1; 588 assert(poly(x, pp) == y); 589 } 590 591 @safe unittest 592 { 593 static assert(poly(3.0, [1.0, 2.0, 3.0]) == 34); 594 } 595 596 private Unqual!(CommonType!(T1, T2)) polyImplBase(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc 597 if (isFloatingPoint!T1 && isFloatingPoint!T2) 598 { 599 ptrdiff_t i = A.length - 1; 600 typeof(return) r = A[i]; 601 while (--i >= 0) 602 { 603 r *= x; 604 r += A[i]; 605 } 606 return r; 607 } 608 609 version (linux) version = GenericPosixVersion; 610 else version (FreeBSD) version = GenericPosixVersion; 611 else version (OpenBSD) version = GenericPosixVersion; 612 else version (Solaris) version = GenericPosixVersion; 613 else version (DragonFlyBSD) version = GenericPosixVersion; 614 615 private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc 616 { 617 version (D_InlineAsm_X86) 618 { 619 if (__ctfe) 620 { 621 return polyImplBase(x, A); 622 } 623 version (Windows) 624 { 625 // BUG: This code assumes a frame pointer in EBP. 626 asm pure nothrow @nogc // assembler by W. Bright 627 { 628 // EDX = (A.length - 1) * real.sizeof 629 mov ECX,A[EBP] ; // ECX = A.length 630 dec ECX ; 631 lea EDX,[ECX][ECX*8] ; 632 add EDX,ECX ; 633 add EDX,A+4[EBP] ; 634 fld real ptr [EDX] ; // ST0 = coeff[ECX] 635 jecxz return_ST ; 636 fld x[EBP] ; // ST0 = x 637 fxch ST(1) ; // ST1 = x, ST0 = r 638 align 4 ; 639 L2: fmul ST,ST(1) ; // r *= x 640 fld real ptr -10[EDX] ; 641 sub EDX,10 ; // deg-- 642 faddp ST(1),ST ; 643 dec ECX ; 644 jne L2 ; 645 fxch ST(1) ; // ST1 = r, ST0 = x 646 fstp ST(0) ; // dump x 647 align 4 ; 648 return_ST: ; 649 } 650 } 651 else version (GenericPosixVersion) 652 { 653 asm pure nothrow @nogc // assembler by W. Bright 654 { 655 // EDX = (A.length - 1) * real.sizeof 656 mov ECX,A[EBP] ; // ECX = A.length 657 dec ECX ; 658 lea EDX,[ECX*8] ; 659 lea EDX,[EDX][ECX*4] ; 660 add EDX,A+4[EBP] ; 661 fld real ptr [EDX] ; // ST0 = coeff[ECX] 662 jecxz return_ST ; 663 fld x[EBP] ; // ST0 = x 664 fxch ST(1) ; // ST1 = x, ST0 = r 665 align 4 ; 666 L2: fmul ST,ST(1) ; // r *= x 667 fld real ptr -12[EDX] ; 668 sub EDX,12 ; // deg-- 669 faddp ST(1),ST ; 670 dec ECX ; 671 jne L2 ; 672 fxch ST(1) ; // ST1 = r, ST0 = x 673 fstp ST(0) ; // dump x 674 align 4 ; 675 return_ST: ; 676 } 677 } 678 else version (OSX) 679 { 680 asm pure nothrow @nogc // assembler by W. Bright 681 { 682 // EDX = (A.length - 1) * real.sizeof 683 mov ECX,A[EBP] ; // ECX = A.length 684 dec ECX ; 685 lea EDX,[ECX*8] ; 686 add EDX,EDX ; 687 add EDX,A+4[EBP] ; 688 fld real ptr [EDX] ; // ST0 = coeff[ECX] 689 jecxz return_ST ; 690 fld x[EBP] ; // ST0 = x 691 fxch ST(1) ; // ST1 = x, ST0 = r 692 align 4 ; 693 L2: fmul ST,ST(1) ; // r *= x 694 fld real ptr -16[EDX] ; 695 sub EDX,16 ; // deg-- 696 faddp ST(1),ST ; 697 dec ECX ; 698 jne L2 ; 699 fxch ST(1) ; // ST1 = r, ST0 = x 700 fstp ST(0) ; // dump x 701 align 4 ; 702 return_ST: ; 703 } 704 } 705 else 706 { 707 static assert(0); 708 } 709 } 710 else 711 { 712 return polyImplBase(x, A); 713 } 714 } 715 716 /** 717 * Gives the next power of two after `val`. `T` can be any built-in 718 * numerical type. 719 * 720 * If the operation would lead to an over/underflow, this function will 721 * return `0`. 722 * 723 * Params: 724 * val = any number 725 * 726 * Returns: 727 * the next power of two after `val` 728 */ 729 T nextPow2(T)(const T val) 730 if (isIntegral!T) 731 { 732 return powIntegralImpl!(PowType.ceil)(val); 733 } 734 735 /// ditto 736 T nextPow2(T)(const T val) 737 if (isFloatingPoint!T) 738 { 739 return powFloatingPointImpl!(PowType.ceil)(val); 740 } 741 742 /// 743 @safe @nogc pure nothrow unittest 744 { 745 assert(nextPow2(2) == 4); 746 assert(nextPow2(10) == 16); 747 assert(nextPow2(4000) == 4096); 748 749 assert(nextPow2(-2) == -4); 750 assert(nextPow2(-10) == -16); 751 752 assert(nextPow2(uint.max) == 0); 753 assert(nextPow2(uint.min) == 0); 754 assert(nextPow2(size_t.max) == 0); 755 assert(nextPow2(size_t.min) == 0); 756 757 assert(nextPow2(int.max) == 0); 758 assert(nextPow2(int.min) == 0); 759 assert(nextPow2(long.max) == 0); 760 assert(nextPow2(long.min) == 0); 761 } 762 763 /// 764 @safe @nogc pure nothrow unittest 765 { 766 assert(nextPow2(2.1) == 4.0); 767 assert(nextPow2(-2.0) == -4.0); 768 assert(nextPow2(0.25) == 0.5); 769 assert(nextPow2(-4.0) == -8.0); 770 771 assert(nextPow2(double.max) == 0.0); 772 assert(nextPow2(double.infinity) == double.infinity); 773 } 774 775 @safe @nogc pure nothrow unittest 776 { 777 assert(nextPow2(ubyte(2)) == 4); 778 assert(nextPow2(ubyte(10)) == 16); 779 780 assert(nextPow2(byte(2)) == 4); 781 assert(nextPow2(byte(10)) == 16); 782 783 assert(nextPow2(short(2)) == 4); 784 assert(nextPow2(short(10)) == 16); 785 assert(nextPow2(short(4000)) == 4096); 786 787 assert(nextPow2(ushort(2)) == 4); 788 assert(nextPow2(ushort(10)) == 16); 789 assert(nextPow2(ushort(4000)) == 4096); 790 } 791 792 @safe @nogc pure nothrow unittest 793 { 794 foreach (ulong i; 1 .. 62) 795 { 796 assert(nextPow2(1UL << i) == 2UL << i); 797 assert(nextPow2((1UL << i) - 1) == 1UL << i); 798 assert(nextPow2((1UL << i) + 1) == 2UL << i); 799 assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i); 800 } 801 } 802 803 @safe @nogc pure nothrow unittest 804 { 805 import std.math.traits : isNaN; 806 import std.meta : AliasSeq; 807 808 static foreach (T; AliasSeq!(float, double, real)) 809 {{ 810 enum T subNormal = T.min_normal / 2; 811 812 static if (subNormal) assert(nextPow2(subNormal) == T.min_normal); 813 814 assert(nextPow2(T(0.0)) == 0.0); 815 816 assert(nextPow2(T(2.0)) == 4.0); 817 assert(nextPow2(T(2.1)) == 4.0); 818 assert(nextPow2(T(3.1)) == 4.0); 819 assert(nextPow2(T(4.0)) == 8.0); 820 assert(nextPow2(T(0.25)) == 0.5); 821 822 assert(nextPow2(T(-2.0)) == -4.0); 823 assert(nextPow2(T(-2.1)) == -4.0); 824 assert(nextPow2(T(-3.1)) == -4.0); 825 assert(nextPow2(T(-4.0)) == -8.0); 826 assert(nextPow2(T(-0.25)) == -0.5); 827 828 assert(nextPow2(T.max) == 0); 829 assert(nextPow2(-T.max) == 0); 830 831 assert(nextPow2(T.infinity) == T.infinity); 832 assert(nextPow2(T.init).isNaN); 833 }} 834 } 835 836 // https://issues.dlang.org/show_bug.cgi?id=15973 837 @safe @nogc pure nothrow unittest 838 { 839 assert(nextPow2(uint.max / 2) == uint.max / 2 + 1); 840 assert(nextPow2(uint.max / 2 + 2) == 0); 841 assert(nextPow2(int.max / 2) == int.max / 2 + 1); 842 assert(nextPow2(int.max / 2 + 2) == 0); 843 assert(nextPow2(int.min + 1) == int.min); 844 } 845 846 /** 847 * Gives the last power of two before `val`. $(T) can be any built-in 848 * numerical type. 849 * 850 * Params: 851 * val = any number 852 * 853 * Returns: 854 * the last power of two before `val` 855 */ 856 T truncPow2(T)(const T val) 857 if (isIntegral!T) 858 { 859 return powIntegralImpl!(PowType.floor)(val); 860 } 861 862 /// ditto 863 T truncPow2(T)(const T val) 864 if (isFloatingPoint!T) 865 { 866 return powFloatingPointImpl!(PowType.floor)(val); 867 } 868 869 /// 870 @safe @nogc pure nothrow unittest 871 { 872 assert(truncPow2(3) == 2); 873 assert(truncPow2(4) == 4); 874 assert(truncPow2(10) == 8); 875 assert(truncPow2(4000) == 2048); 876 877 assert(truncPow2(-5) == -4); 878 assert(truncPow2(-20) == -16); 879 880 assert(truncPow2(uint.max) == int.max + 1); 881 assert(truncPow2(uint.min) == 0); 882 assert(truncPow2(ulong.max) == long.max + 1); 883 assert(truncPow2(ulong.min) == 0); 884 885 assert(truncPow2(int.max) == (int.max / 2) + 1); 886 assert(truncPow2(int.min) == int.min); 887 assert(truncPow2(long.max) == (long.max / 2) + 1); 888 assert(truncPow2(long.min) == long.min); 889 } 890 891 /// 892 @safe @nogc pure nothrow unittest 893 { 894 assert(truncPow2(2.1) == 2.0); 895 assert(truncPow2(7.0) == 4.0); 896 assert(truncPow2(-1.9) == -1.0); 897 assert(truncPow2(0.24) == 0.125); 898 assert(truncPow2(-7.0) == -4.0); 899 900 assert(truncPow2(double.infinity) == double.infinity); 901 } 902 903 @safe @nogc pure nothrow unittest 904 { 905 assert(truncPow2(ubyte(3)) == 2); 906 assert(truncPow2(ubyte(4)) == 4); 907 assert(truncPow2(ubyte(10)) == 8); 908 909 assert(truncPow2(byte(3)) == 2); 910 assert(truncPow2(byte(4)) == 4); 911 assert(truncPow2(byte(10)) == 8); 912 913 assert(truncPow2(ushort(3)) == 2); 914 assert(truncPow2(ushort(4)) == 4); 915 assert(truncPow2(ushort(10)) == 8); 916 assert(truncPow2(ushort(4000)) == 2048); 917 918 assert(truncPow2(short(3)) == 2); 919 assert(truncPow2(short(4)) == 4); 920 assert(truncPow2(short(10)) == 8); 921 assert(truncPow2(short(4000)) == 2048); 922 } 923 924 @safe @nogc pure nothrow unittest 925 { 926 foreach (ulong i; 1 .. 62) 927 { 928 assert(truncPow2(2UL << i) == 2UL << i); 929 assert(truncPow2((2UL << i) + 1) == 2UL << i); 930 assert(truncPow2((2UL << i) - 1) == 1UL << i); 931 assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i); 932 } 933 } 934 935 @safe @nogc pure nothrow unittest 936 { 937 import std.math.traits : isNaN; 938 import std.meta : AliasSeq; 939 940 static foreach (T; AliasSeq!(float, double, real)) 941 { 942 assert(truncPow2(T(0.0)) == 0.0); 943 944 assert(truncPow2(T(4.0)) == 4.0); 945 assert(truncPow2(T(2.1)) == 2.0); 946 assert(truncPow2(T(3.5)) == 2.0); 947 assert(truncPow2(T(7.0)) == 4.0); 948 assert(truncPow2(T(0.24)) == 0.125); 949 950 assert(truncPow2(T(-2.0)) == -2.0); 951 assert(truncPow2(T(-2.1)) == -2.0); 952 assert(truncPow2(T(-3.1)) == -2.0); 953 assert(truncPow2(T(-7.0)) == -4.0); 954 assert(truncPow2(T(-0.24)) == -0.125); 955 956 assert(truncPow2(T.infinity) == T.infinity); 957 assert(truncPow2(T.init).isNaN); 958 } 959 } 960 961 private enum PowType 962 { 963 floor, 964 ceil 965 } 966 967 pragma(inline, true) 968 private T powIntegralImpl(PowType type, T)(T val) 969 { 970 import core.bitop : bsr; 971 972 if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min))) 973 return 0; 974 else 975 { 976 static if (isSigned!T) 977 return cast() cast(T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type); 978 else 979 return cast() cast(T) (T(1) << bsr(val) + type); 980 } 981 } 982 983 private T powFloatingPointImpl(PowType type, T)(T x) 984 { 985 import std.math.traits : copysign, isFinite; 986 import std.math.exponential : frexp; 987 988 if (!x.isFinite) 989 return x; 990 991 if (!x) 992 return x; 993 994 int exp; 995 auto y = frexp(x, exp); 996 997 static if (type == PowType.ceil) 998 y = core.math.ldexp(cast(T) 0.5, exp + 1); 999 else 1000 y = core.math.ldexp(cast(T) 0.5, exp); 1001 1002 if (!y.isFinite) 1003 return cast(T) 0.0; 1004 1005 y = copysign(y, x); 1006 1007 return y; 1008 }