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 : 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 pragma(inline, true) // LDC 616 private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc 617 { 618 version (LDC) 619 { 620 return polyImplBase(x, A); 621 } 622 else version (D_InlineAsm_X86) 623 { 624 if (__ctfe) 625 { 626 return polyImplBase(x, A); 627 } 628 version (Windows) 629 { 630 // BUG: This code assumes a frame pointer in EBP. 631 asm pure nothrow @nogc // assembler by W. Bright 632 { 633 // EDX = (A.length - 1) * real.sizeof 634 mov ECX,A[EBP] ; // ECX = A.length 635 dec ECX ; 636 lea EDX,[ECX][ECX*8] ; 637 add EDX,ECX ; 638 add EDX,A+4[EBP] ; 639 fld real ptr [EDX] ; // ST0 = coeff[ECX] 640 jecxz return_ST ; 641 fld x[EBP] ; // ST0 = x 642 fxch ST(1) ; // ST1 = x, ST0 = r 643 align 4 ; 644 L2: fmul ST,ST(1) ; // r *= x 645 fld real ptr -10[EDX] ; 646 sub EDX,10 ; // deg-- 647 faddp ST(1),ST ; 648 dec ECX ; 649 jne L2 ; 650 fxch ST(1) ; // ST1 = r, ST0 = x 651 fstp ST(0) ; // dump x 652 align 4 ; 653 return_ST: ; 654 } 655 } 656 else version (GenericPosixVersion) 657 { 658 asm pure nothrow @nogc // assembler by W. Bright 659 { 660 // EDX = (A.length - 1) * real.sizeof 661 mov ECX,A[EBP] ; // ECX = A.length 662 dec ECX ; 663 lea EDX,[ECX*8] ; 664 lea EDX,[EDX][ECX*4] ; 665 add EDX,A+4[EBP] ; 666 fld real ptr [EDX] ; // ST0 = coeff[ECX] 667 jecxz return_ST ; 668 fld x[EBP] ; // ST0 = x 669 fxch ST(1) ; // ST1 = x, ST0 = r 670 align 4 ; 671 L2: fmul ST,ST(1) ; // r *= x 672 fld real ptr -12[EDX] ; 673 sub EDX,12 ; // deg-- 674 faddp ST(1),ST ; 675 dec ECX ; 676 jne L2 ; 677 fxch ST(1) ; // ST1 = r, ST0 = x 678 fstp ST(0) ; // dump x 679 align 4 ; 680 return_ST: ; 681 } 682 } 683 else version (OSX) 684 { 685 asm pure nothrow @nogc // assembler by W. Bright 686 { 687 // EDX = (A.length - 1) * real.sizeof 688 mov ECX,A[EBP] ; // ECX = A.length 689 dec ECX ; 690 lea EDX,[ECX*8] ; 691 add EDX,EDX ; 692 add EDX,A+4[EBP] ; 693 fld real ptr [EDX] ; // ST0 = coeff[ECX] 694 jecxz return_ST ; 695 fld x[EBP] ; // ST0 = x 696 fxch ST(1) ; // ST1 = x, ST0 = r 697 align 4 ; 698 L2: fmul ST,ST(1) ; // r *= x 699 fld real ptr -16[EDX] ; 700 sub EDX,16 ; // deg-- 701 faddp ST(1),ST ; 702 dec ECX ; 703 jne L2 ; 704 fxch ST(1) ; // ST1 = r, ST0 = x 705 fstp ST(0) ; // dump x 706 align 4 ; 707 return_ST: ; 708 } 709 } 710 else 711 { 712 static assert(0); 713 } 714 } 715 else 716 { 717 return polyImplBase(x, A); 718 } 719 } 720 721 /** 722 * Gives the next power of two after `val`. `T` can be any built-in 723 * numerical type. 724 * 725 * If the operation would lead to an over/underflow, this function will 726 * return `0`. 727 * 728 * Params: 729 * val = any number 730 * 731 * Returns: 732 * the next power of two after `val` 733 */ 734 T nextPow2(T)(const T val) 735 if (isIntegral!T) 736 { 737 return powIntegralImpl!(PowType.ceil)(val); 738 } 739 740 /// ditto 741 T nextPow2(T)(const T val) 742 if (isFloatingPoint!T) 743 { 744 return powFloatingPointImpl!(PowType.ceil)(val); 745 } 746 747 /// 748 @safe @nogc pure nothrow unittest 749 { 750 assert(nextPow2(2) == 4); 751 assert(nextPow2(10) == 16); 752 assert(nextPow2(4000) == 4096); 753 754 assert(nextPow2(-2) == -4); 755 assert(nextPow2(-10) == -16); 756 757 assert(nextPow2(uint.max) == 0); 758 assert(nextPow2(uint.min) == 0); 759 assert(nextPow2(size_t.max) == 0); 760 assert(nextPow2(size_t.min) == 0); 761 762 assert(nextPow2(int.max) == 0); 763 assert(nextPow2(int.min) == 0); 764 assert(nextPow2(long.max) == 0); 765 assert(nextPow2(long.min) == 0); 766 } 767 768 /// 769 @safe @nogc pure nothrow unittest 770 { 771 assert(nextPow2(2.1) == 4.0); 772 assert(nextPow2(-2.0) == -4.0); 773 assert(nextPow2(0.25) == 0.5); 774 assert(nextPow2(-4.0) == -8.0); 775 776 assert(nextPow2(double.max) == 0.0); 777 assert(nextPow2(double.infinity) == double.infinity); 778 } 779 780 @safe @nogc pure nothrow unittest 781 { 782 assert(nextPow2(ubyte(2)) == 4); 783 assert(nextPow2(ubyte(10)) == 16); 784 785 assert(nextPow2(byte(2)) == 4); 786 assert(nextPow2(byte(10)) == 16); 787 788 assert(nextPow2(short(2)) == 4); 789 assert(nextPow2(short(10)) == 16); 790 assert(nextPow2(short(4000)) == 4096); 791 792 assert(nextPow2(ushort(2)) == 4); 793 assert(nextPow2(ushort(10)) == 16); 794 assert(nextPow2(ushort(4000)) == 4096); 795 } 796 797 @safe @nogc pure nothrow unittest 798 { 799 foreach (ulong i; 1 .. 62) 800 { 801 assert(nextPow2(1UL << i) == 2UL << i); 802 assert(nextPow2((1UL << i) - 1) == 1UL << i); 803 assert(nextPow2((1UL << i) + 1) == 2UL << i); 804 assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i); 805 } 806 } 807 808 @safe @nogc pure nothrow unittest 809 { 810 import std.math.traits : isNaN; 811 import std.meta : AliasSeq; 812 813 static foreach (T; AliasSeq!(float, double, real)) 814 {{ 815 enum T subNormal = T.min_normal / 2; 816 817 static if (subNormal) assert(nextPow2(subNormal) == T.min_normal); 818 819 assert(nextPow2(T(0.0)) == 0.0); 820 821 assert(nextPow2(T(2.0)) == 4.0); 822 assert(nextPow2(T(2.1)) == 4.0); 823 assert(nextPow2(T(3.1)) == 4.0); 824 assert(nextPow2(T(4.0)) == 8.0); 825 assert(nextPow2(T(0.25)) == 0.5); 826 827 assert(nextPow2(T(-2.0)) == -4.0); 828 assert(nextPow2(T(-2.1)) == -4.0); 829 assert(nextPow2(T(-3.1)) == -4.0); 830 assert(nextPow2(T(-4.0)) == -8.0); 831 assert(nextPow2(T(-0.25)) == -0.5); 832 833 assert(nextPow2(T.max) == 0); 834 assert(nextPow2(-T.max) == 0); 835 836 assert(nextPow2(T.infinity) == T.infinity); 837 assert(nextPow2(T.init).isNaN); 838 }} 839 } 840 841 // https://issues.dlang.org/show_bug.cgi?id=15973 842 @safe @nogc pure nothrow unittest 843 { 844 assert(nextPow2(uint.max / 2) == uint.max / 2 + 1); 845 assert(nextPow2(uint.max / 2 + 2) == 0); 846 assert(nextPow2(int.max / 2) == int.max / 2 + 1); 847 assert(nextPow2(int.max / 2 + 2) == 0); 848 assert(nextPow2(int.min + 1) == int.min); 849 } 850 851 /** 852 * Gives the last power of two before `val`. $(T) can be any built-in 853 * numerical type. 854 * 855 * Params: 856 * val = any number 857 * 858 * Returns: 859 * the last power of two before `val` 860 */ 861 T truncPow2(T)(const T val) 862 if (isIntegral!T) 863 { 864 return powIntegralImpl!(PowType.floor)(val); 865 } 866 867 /// ditto 868 T truncPow2(T)(const T val) 869 if (isFloatingPoint!T) 870 { 871 return powFloatingPointImpl!(PowType.floor)(val); 872 } 873 874 /// 875 @safe @nogc pure nothrow unittest 876 { 877 assert(truncPow2(3) == 2); 878 assert(truncPow2(4) == 4); 879 assert(truncPow2(10) == 8); 880 assert(truncPow2(4000) == 2048); 881 882 assert(truncPow2(-5) == -4); 883 assert(truncPow2(-20) == -16); 884 885 assert(truncPow2(uint.max) == int.max + 1); 886 assert(truncPow2(uint.min) == 0); 887 assert(truncPow2(ulong.max) == long.max + 1); 888 assert(truncPow2(ulong.min) == 0); 889 890 assert(truncPow2(int.max) == (int.max / 2) + 1); 891 version (LDC) 892 { 893 // this test relies on undefined behaviour, i.e. (1 << 63) == int.min 894 // that fails for LDC with optimizations enabled 895 } 896 else 897 { 898 assert(truncPow2(int.min) == int.min); 899 } 900 assert(truncPow2(long.max) == (long.max / 2) + 1); 901 assert(truncPow2(long.min) == long.min); 902 } 903 904 /// 905 @safe @nogc pure nothrow unittest 906 { 907 assert(truncPow2(2.1) == 2.0); 908 assert(truncPow2(7.0) == 4.0); 909 assert(truncPow2(-1.9) == -1.0); 910 assert(truncPow2(0.24) == 0.125); 911 assert(truncPow2(-7.0) == -4.0); 912 913 assert(truncPow2(double.infinity) == double.infinity); 914 } 915 916 @safe @nogc pure nothrow unittest 917 { 918 assert(truncPow2(ubyte(3)) == 2); 919 assert(truncPow2(ubyte(4)) == 4); 920 assert(truncPow2(ubyte(10)) == 8); 921 922 assert(truncPow2(byte(3)) == 2); 923 assert(truncPow2(byte(4)) == 4); 924 assert(truncPow2(byte(10)) == 8); 925 926 assert(truncPow2(ushort(3)) == 2); 927 assert(truncPow2(ushort(4)) == 4); 928 assert(truncPow2(ushort(10)) == 8); 929 assert(truncPow2(ushort(4000)) == 2048); 930 931 assert(truncPow2(short(3)) == 2); 932 assert(truncPow2(short(4)) == 4); 933 assert(truncPow2(short(10)) == 8); 934 assert(truncPow2(short(4000)) == 2048); 935 } 936 937 @safe @nogc pure nothrow unittest 938 { 939 foreach (ulong i; 1 .. 62) 940 { 941 assert(truncPow2(2UL << i) == 2UL << i); 942 assert(truncPow2((2UL << i) + 1) == 2UL << i); 943 assert(truncPow2((2UL << i) - 1) == 1UL << i); 944 assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i); 945 } 946 } 947 948 @safe @nogc pure nothrow unittest 949 { 950 import std.math.traits : isNaN; 951 import std.meta : AliasSeq; 952 953 static foreach (T; AliasSeq!(float, double, real)) 954 { 955 assert(truncPow2(T(0.0)) == 0.0); 956 957 assert(truncPow2(T(4.0)) == 4.0); 958 assert(truncPow2(T(2.1)) == 2.0); 959 assert(truncPow2(T(3.5)) == 2.0); 960 assert(truncPow2(T(7.0)) == 4.0); 961 assert(truncPow2(T(0.24)) == 0.125); 962 963 assert(truncPow2(T(-2.0)) == -2.0); 964 assert(truncPow2(T(-2.1)) == -2.0); 965 assert(truncPow2(T(-3.1)) == -2.0); 966 assert(truncPow2(T(-7.0)) == -4.0); 967 assert(truncPow2(T(-0.24)) == -0.125); 968 969 assert(truncPow2(T.infinity) == T.infinity); 970 assert(truncPow2(T.init).isNaN); 971 } 972 } 973 974 private enum PowType 975 { 976 floor, 977 ceil 978 } 979 980 pragma(inline, true) 981 private T powIntegralImpl(PowType type, T)(T val) 982 { 983 import core.bitop : bsr; 984 985 if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min))) 986 return 0; 987 else 988 { 989 static if (isSigned!T) 990 return cast(Unqual!T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type); 991 else 992 return cast(Unqual!T) (T(1) << bsr(val) + type); 993 } 994 } 995 996 private T powFloatingPointImpl(PowType type, T)(T x) 997 { 998 import std.math.traits : copysign, isFinite; 999 import std.math.exponential : frexp; 1000 1001 if (!x.isFinite) 1002 return x; 1003 1004 if (!x) 1005 return x; 1006 1007 int exp; 1008 auto y = frexp(x, exp); 1009 1010 static if (type == PowType.ceil) 1011 y = core.math.ldexp(cast(T) 0.5, exp + 1); 1012 else 1013 y = core.math.ldexp(cast(T) 0.5, exp); 1014 1015 if (!y.isFinite) 1016 return cast(T) 0.0; 1017 1018 y = copysign(y, x); 1019 1020 return y; 1021 }