1 // Written in the D programming language.
2
3 /**
4 This is a submodule of $(MREF std, math).
5
6 It contains several exponential and logarithm functions.
7
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9 D implementations of exp, expm1, exp2, log, log10, log1p, and log2
10 functions are based on the CEPHES math library, which is Copyright
11 (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) and are
12 incorporated herein by permission of the author. The author reserves
13 the right to distribute this material elsewhere under different
14 copying permissions. These modifications are distributed here under
15 the following terms:
16 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
17 Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston,
18 Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
19 Source: $(PHOBOSSRC std/math/exponential.d)
20
21 Macros:
22 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
23 <caption>Special Values</caption>
24 $0</table>
25 NAN = $(RED NAN)
26 PLUSMN = ±
27 INFIN = ∞
28 PLUSMNINF = ±∞
29 LT = <
30 GT = >
31 */
32
33 module std.math.exponential;
34
35 import std.traits : isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual;
36
37 static import core.math;
38 static import core.stdc.math;
39
version(DigitalMars)40 version (DigitalMars)
41 {
42 version = INLINE_YL2X; // x87 has opcodes for these
43 }
44
45 version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
46 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
47
48 version (InlineAsm_X86_Any) version = InlineAsm_X87;
version(InlineAsm_X87)49 version (InlineAsm_X87)
50 {
51 static assert(real.mant_dig == 64);
52 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
53 }
54
version(D_HardFloat)55 version (D_HardFloat)
56 {
57 // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport
58 version (IeeeFlagsSupport) version = FloatingPointControlSupport;
59 }
60
61 /**
62 * Compute the value of x $(SUPERSCRIPT n), where n is an integer
63 */
64 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow
65 if (isFloatingPoint!(F) && isIntegral!(G))
66 {
67 import std.traits : Unsigned;
68
69 real p = 1.0, v = void;
70 Unsigned!(Unqual!G) m = n;
71
72 if (n < 0)
73 {
74 if (n == -1) return 1 / x;
75
76 m = cast(typeof(m))(0 - n);
77 v = p / x;
78 }
79 else
80 {
81 switch (n)
82 {
83 case 0:
84 return 1.0;
85 case 1:
86 return x;
87 case 2:
88 return x * x;
89 default:
90 }
91
92 v = x;
93 }
94
95 while (1)
96 {
97 if (m & 1)
98 p *= v;
99 m >>= 1;
100 if (!m)
101 break;
102 v *= v;
103 }
104 return p;
105 }
106
107 ///
108 @safe pure nothrow @nogc unittest
109 {
110 import std.math.operations : feqrel;
111
112 assert(pow(2.0, 5) == 32.0);
113 assert(pow(1.5, 9).feqrel(38.4433) > 16);
114 assert(pow(real.nan, 2) is real.nan);
115 assert(pow(real.infinity, 2) == real.infinity);
116 }
117
118 @safe pure nothrow @nogc unittest
119 {
120 import std.math.operations : isClose, feqrel;
121
122 // Make sure it instantiates and works properly on immutable values and
123 // with various integer and float types.
124 immutable real x = 46;
125 immutable float xf = x;
126 immutable double xd = x;
127 immutable uint one = 1;
128 immutable ushort two = 2;
129 immutable ubyte three = 3;
130 immutable ulong eight = 8;
131
132 immutable int neg1 = -1;
133 immutable short neg2 = -2;
134 immutable byte neg3 = -3;
135 immutable long neg8 = -8;
136
137
138 assert(pow(x,0) == 1.0);
139 assert(pow(xd,one) == x);
140 assert(pow(xf,two) == x * x);
141 assert(pow(x,three) == x * x * x);
142 assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
143
144 assert(pow(x, neg1) == 1 / x);
145
146 assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15));
147 assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15));
148
149 assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1);
150 }
151
152 @safe @nogc nothrow unittest
153 {
154 import std.math.operations : isClose;
155
156 assert(isClose(pow(2.0L, 10L), 1024, 1e-18));
157 }
158
159 // https://issues.dlang.org/show_bug.cgi?id=21601
160 @safe @nogc nothrow pure unittest
161 {
162 // When reals are large enough the results of pow(b, e) can be
163 // calculated correctly, if b is of type float or double and e is
164 // not too large.
165 static if (real.mant_dig >= 64)
166 {
167 // expected result: 3.790e-42
168 assert(pow(-513645318757045764096.0f, -2) > 0.0);
169
170 // expected result: 3.763915357831797e-309
171 assert(pow(-1.6299717435255677e+154, -2) > 0.0);
172 }
173 }
174
175 @safe @nogc nothrow unittest
176 {
177 import std.math.operations : isClose;
178 import std.math.traits : isInfinity;
179
180 static float f1 = 19100.0f;
181 static float f2 = 0.000012f;
182
183 assert(isClose(pow(f1,9), 3.3829868e+38f));
184 assert(isInfinity(pow(f1,10)));
185 assert(pow(f2,9) > 0.0f);
186 assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
187
188 static double d1 = 21800.0;
189 static double d2 = 0.000012;
190
191 assert(isClose(pow(d1,71), 1.0725339442974e+308));
192 assert(isInfinity(pow(d1,72)));
193 assert(pow(d2,65) > 0.0f);
194 assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
195
196 static if (real.mant_dig == 64) // x87
197 {
198 static real r1 = 21950.0L;
199 static real r2 = 0.000011L;
200
201 assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
202 assert(isInfinity(pow(r1,1137)));
203 assert(pow(r2,998) > 0.0L);
204 assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
205 }
206 }
207
208 @safe @nogc nothrow pure unittest
209 {
210 import std.math.operations : isClose;
211
212 enum f1 = 19100.0f;
213 enum f2 = 0.000012f;
214
215 static assert(isClose(pow(f1,9), 3.3829868e+38f));
216 static assert(pow(f1,10) > float.max);
217 static assert(pow(f2,9) > 0.0f);
218 static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
219
220 enum d1 = 21800.0;
221 enum d2 = 0.000012;
222
223 static assert(isClose(pow(d1,71), 1.0725339442974e+308));
224 static assert(pow(d1,72) > double.max);
225 static assert(pow(d2,65) > 0.0f);
226 static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
227
228 static if (real.mant_dig == 64) // x87
229 {
230 enum r1 = 21950.0L;
231 enum r2 = 0.000011L;
232
233 static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
234 static assert(pow(r1,1137) > real.max);
235 static assert(pow(r2,998) > 0.0L);
236 static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
237 }
238 }
239
240 /**
241 * Compute the power of two integral numbers.
242 *
243 * Params:
244 * x = base
245 * n = exponent
246 *
247 * Returns:
248 * x raised to the power of n. If n is negative the result is 1 / pow(x, -n),
249 * which is calculated as integer division with remainder. This may result in
250 * a division by zero error.
251 *
252 * If both x and n are 0, the result is 1.
253 *
254 * Throws:
255 * If x is 0 and n is negative, the result is the same as the result of a
256 * division by zero.
257 */
258 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow
259 if (isIntegral!(F) && isIntegral!(G))
260 {
261 import std.traits : isSigned;
262
263 typeof(return) p, v = void;
264 Unqual!G m = n;
265
266 static if (isSigned!(F))
267 {
268 if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1);
269 }
270 static if (isSigned!(G))
271 {
272 if (x == 0 && m <= -1) return x / 0;
273 }
274 if (x == 1) return 1;
275 static if (isSigned!(G))
276 {
277 if (m < 0) return 0;
278 }
279
280 switch (m)
281 {
282 case 0:
283 p = 1;
284 break;
285
286 case 1:
287 p = x;
288 break;
289
290 case 2:
291 p = x * x;
292 break;
293
294 default:
295 v = x;
296 p = 1;
297 while (1)
298 {
299 if (m & 1)
300 p *= v;
301 m >>= 1;
302 if (!m)
303 break;
304 v *= v;
305 }
306 break;
307 }
308 return p;
309 }
310
311 ///
312 @safe pure nothrow @nogc unittest
313 {
314 assert(pow(2, 3) == 8);
315 assert(pow(3, 2) == 9);
316
317 assert(pow(2, 10) == 1_024);
318 assert(pow(2, 20) == 1_048_576);
319 assert(pow(2, 30) == 1_073_741_824);
320
321 assert(pow(0, 0) == 1);
322
323 assert(pow(1, -5) == 1);
324 assert(pow(1, -6) == 1);
325 assert(pow(-1, -5) == -1);
326 assert(pow(-1, -6) == 1);
327
328 assert(pow(-2, 5) == -32);
329 assert(pow(-2, -5) == 0);
330 assert(pow(cast(double) -2, -5) == -0.03125);
331 }
332
333 @safe pure nothrow @nogc unittest
334 {
335 immutable int one = 1;
336 immutable byte two = 2;
337 immutable ubyte three = 3;
338 immutable short four = 4;
339 immutable long ten = 10;
340
341 assert(pow(two, three) == 8);
342 assert(pow(two, ten) == 1024);
343 assert(pow(one, ten) == 1);
344 assert(pow(ten, four) == 10_000);
345 assert(pow(four, 10) == 1_048_576);
346 assert(pow(three, four) == 81);
347 }
348
349 // https://issues.dlang.org/show_bug.cgi?id=7006
350 @safe pure nothrow @nogc unittest
351 {
352 assert(pow(5, -1) == 0);
353 assert(pow(-5, -1) == 0);
354 assert(pow(5, -2) == 0);
355 assert(pow(-5, -2) == 0);
356 assert(pow(-1, int.min) == 1);
357 assert(pow(-2, int.min) == 0);
358
359 assert(pow(4294967290UL,2) == 18446744022169944100UL);
360 assert(pow(0,uint.max) == 0);
361 }
362
363 /**Computes integer to floating point powers.*/
364 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow
365 if (isIntegral!I && isFloatingPoint!F)
366 {
367 return pow(cast(real) x, cast(Unqual!F) y);
368 }
369
370 ///
371 @safe pure nothrow @nogc unittest
372 {
373 assert(pow(2, 5.0) == 32.0);
374 assert(pow(7, 3.0) == 343.0);
375 assert(pow(2, real.nan) is real.nan);
376 assert(pow(2, real.infinity) == real.infinity);
377 }
378
379 /**
380 * Calculates x$(SUPERSCRIPT y).
381 *
382 * $(TABLE_SV
383 * $(TR $(TH x) $(TH y) $(TH pow(x, y))
384 * $(TH div 0) $(TH invalid?))
385 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0)
386 * $(TD no) $(TD no) )
387 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN))
388 * $(TD no) $(TD no) )
389 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0)
390 * $(TD no) $(TD no) )
391 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0)
392 * $(TD no) $(TD no) )
393 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN))
394 * $(TD no) $(TD no) )
395 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN))
396 * $(TD no) $(TD no) )
397 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0)
398 * $(TD no) $(TD no) )
399 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN))
400 * $(TD no) $(TD no) )
401 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
402 * $(TD no) $(TD no))
403 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0)
404 * $(TD no) $(TD no) )
405 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
406 * $(TD no) $(TD no) )
407 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN))
408 * $(TD no) $(TD yes) )
409 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN))
410 * $(TD no) $(TD yes))
411 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF))
412 * $(TD yes) $(TD no) )
413 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
414 * $(TD yes) $(TD no))
415 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0)
416 * $(TD no) $(TD no) )
417 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
418 * $(TD no) $(TD no) )
419 * )
420 */
421 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow
422 if (isFloatingPoint!(F) && isFloatingPoint!(G))
423 {
424 import core.math : fabs, sqrt;
425 import std.math.traits : isInfinity, isNaN, signbit;
426
427 alias Float = typeof(return);
428
impl(real x,real y)429 static real impl(real x, real y) @nogc pure nothrow
430 {
431 // Special cases.
432 if (isNaN(y))
433 return y;
434 if (isNaN(x) && y != 0.0)
435 return x;
436
437 // Even if x is NaN.
438 if (y == 0.0)
439 return 1.0;
440 if (y == 1.0)
441 return x;
442
443 if (isInfinity(y))
444 {
445 if (isInfinity(x))
446 {
447 if (!signbit(y) && !signbit(x))
448 return F.infinity;
449 else
450 return F.nan;
451 }
452 else if (fabs(x) > 1)
453 {
454 if (signbit(y))
455 return +0.0;
456 else
457 return F.infinity;
458 }
459 else if (fabs(x) == 1)
460 {
461 return F.nan;
462 }
463 else // < 1
464 {
465 if (signbit(y))
466 return F.infinity;
467 else
468 return +0.0;
469 }
470 }
471 if (isInfinity(x))
472 {
473 if (signbit(x))
474 {
475 long i = cast(long) y;
476 if (y > 0.0)
477 {
478 if (i == y && i & 1)
479 return -F.infinity;
480 else if (i == y)
481 return F.infinity;
482 else
483 return -F.nan;
484 }
485 else if (y < 0.0)
486 {
487 if (i == y && i & 1)
488 return -0.0;
489 else if (i == y)
490 return +0.0;
491 else
492 return F.nan;
493 }
494 }
495 else
496 {
497 if (y > 0.0)
498 return F.infinity;
499 else if (y < 0.0)
500 return +0.0;
501 }
502 }
503
504 if (x == 0.0)
505 {
506 if (signbit(x))
507 {
508 long i = cast(long) y;
509 if (y > 0.0)
510 {
511 if (i == y && i & 1)
512 return -0.0;
513 else
514 return +0.0;
515 }
516 else if (y < 0.0)
517 {
518 if (i == y && i & 1)
519 return -F.infinity;
520 else
521 return F.infinity;
522 }
523 }
524 else
525 {
526 if (y > 0.0)
527 return +0.0;
528 else if (y < 0.0)
529 return F.infinity;
530 }
531 }
532 if (x == 1.0)
533 return 1.0;
534
535 if (y >= F.max)
536 {
537 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
538 return 0.0;
539 if (x > 1.0 || x < -1.0)
540 return F.infinity;
541 }
542 if (y <= -F.max)
543 {
544 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
545 return F.infinity;
546 if (x > 1.0 || x < -1.0)
547 return 0.0;
548 }
549
550 if (x >= F.max)
551 {
552 if (y > 0.0)
553 return F.infinity;
554 else
555 return 0.0;
556 }
557 if (x <= -F.max)
558 {
559 long i = cast(long) y;
560 if (y > 0.0)
561 {
562 if (i == y && i & 1)
563 return -F.infinity;
564 else
565 return F.infinity;
566 }
567 else if (y < 0.0)
568 {
569 if (i == y && i & 1)
570 return -0.0;
571 else
572 return +0.0;
573 }
574 }
575
576 // Integer power of x.
577 long iy = cast(long) y;
578 if (iy == y && fabs(y) < 32_768.0)
579 return pow(x, iy);
580
581 real sign = 1.0;
582 if (x < 0)
583 {
584 // Result is real only if y is an integer
585 // Check for a non-zero fractional part
586 enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
587 static if (maxOdd > ulong.max)
588 {
589 // Generic method, for any FP type
590 import std.math.rounding : floor;
591 if (floor(y) != y)
592 return sqrt(x); // Complex result -- create a NaN
593
594 const hy = 0.5 * y;
595 if (floor(hy) != hy)
596 sign = -1.0;
597 }
598 else
599 {
600 // Much faster, if ulong has enough precision
601 const absY = fabs(y);
602 if (absY <= maxOdd)
603 {
604 const uy = cast(ulong) absY;
605 if (uy != absY)
606 return sqrt(x); // Complex result -- create a NaN
607
608 if (uy & 1)
609 sign = -1.0;
610 }
611 }
612 x = -x;
613 }
614 version (INLINE_YL2X)
615 {
616 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
617 // TODO: This is not accurate in practice. A fast and accurate
618 // (though complicated) method is described in:
619 // "An efficient rounding boundary test for pow(x, y)
620 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
621 return sign * exp2( core.math.yl2x(x, y) );
622 }
623 else
624 {
625 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
626 // TODO: This is not accurate in practice. A fast and accurate
627 // (though complicated) method is described in:
628 // "An efficient rounding boundary test for pow(x, y)
629 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
630 Float w = exp2(y * log2(x));
631 return sign * w;
632 }
633 }
634 return impl(x, y);
635 }
636
637 ///
638 @safe pure nothrow @nogc unittest
639 {
640 import std.math.operations : isClose;
641
642 assert(isClose(pow(2.0, 3.0), 8.0));
643 assert(isClose(pow(1.5, 10.0), 57.6650390625));
644
645 // square root of 9
646 assert(isClose(pow(9.0, 0.5), 3.0));
647 // 10th root of 1024
648 assert(isClose(pow(1024.0, 0.1), 2.0));
649
650 assert(isClose(pow(-4.0, 3.0), -64.0));
651
652 // reciprocal of 4 ^^ 2
653 assert(isClose(pow(4.0, -2.0), 0.0625));
654 // reciprocal of (-2) ^^ 3
655 assert(isClose(pow(-2.0, -3.0), -0.125));
656
657 assert(isClose(pow(-2.5, 3.0), -15.625));
658 // reciprocal of 2.5 ^^ 3
659 assert(isClose(pow(2.5, -3.0), 0.064));
660 // reciprocal of (-2.5) ^^ 3
661 assert(isClose(pow(-2.5, -3.0), -0.064));
662
663 // reciprocal of square root of 4
664 assert(isClose(pow(4.0, -0.5), 0.5));
665
666 // per definition
667 assert(isClose(pow(0.0, 0.0), 1.0));
668 }
669
670 ///
671 @safe pure nothrow @nogc unittest
672 {
673 import std.math.operations : isClose;
674
675 // the result is a complex number
676 // which cannot be represented as floating point number
677 import std.math.traits : isNaN;
678 assert(isNaN(pow(-2.5, -1.5)));
679
680 // use the ^^-operator of std.complex instead
681 import std.complex : complex;
682 auto c1 = complex(-2.5, 0.0);
683 auto c2 = complex(-1.5, 0.0);
684 auto result = c1 ^^ c2;
685 // exact result apparently depends on `real` precision => increased tolerance
686 assert(isClose(result.re, -4.64705438e-17, 2e-4));
687 assert(isClose(result.im, 2.52982e-1, 2e-4));
688 }
689
690 @safe pure nothrow @nogc unittest
691 {
692 import std.math.traits : isNaN;
693
694 assert(pow(1.5, real.infinity) == real.infinity);
695 assert(pow(0.5, real.infinity) == 0.0);
696 assert(pow(1.5, -real.infinity) == 0.0);
697 assert(pow(0.5, -real.infinity) == real.infinity);
698 assert(pow(real.infinity, 1.0) == real.infinity);
699 assert(pow(real.infinity, -1.0) == 0.0);
700 assert(pow(real.infinity, real.infinity) == real.infinity);
701 assert(pow(-real.infinity, 1.0) == -real.infinity);
702 assert(pow(-real.infinity, 2.0) == real.infinity);
703 assert(pow(-real.infinity, -1.0) == -0.0);
704 assert(pow(-real.infinity, -2.0) == 0.0);
705 assert(isNaN(pow(1.0, real.infinity)));
706 assert(pow(0.0, -1.0) == real.infinity);
707 assert(pow(real.nan, 0.0) == 1.0);
708 assert(isNaN(pow(real.nan, 3.0)));
709 assert(isNaN(pow(3.0, real.nan)));
710 }
711
712 @safe @nogc nothrow unittest
713 {
714 import std.math.operations : isClose;
715
716 assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18));
717 }
718
719 @safe pure nothrow @nogc unittest
720 {
721 import std.math.operations : isClose;
722 import std.math.traits : isIdentical, isNaN;
723 import std.math.constants : PI;
724
725 // Test all the special values. These unittests can be run on Windows
726 // by temporarily changing the version (linux) to version (all).
727 immutable float zero = 0;
728 immutable real one = 1;
729 immutable double two = 2;
730 immutable float three = 3;
731 immutable float fnan = float.nan;
732 immutable double dnan = double.nan;
733 immutable real rnan = real.nan;
734 immutable dinf = double.infinity;
735 immutable rninf = -real.infinity;
736
737 assert(pow(fnan, zero) == 1);
738 assert(pow(dnan, zero) == 1);
739 assert(pow(rnan, zero) == 1);
740
741 assert(pow(two, dinf) == double.infinity);
742 assert(isIdentical(pow(0.2f, dinf), +0.0));
743 assert(pow(0.99999999L, rninf) == real.infinity);
744 assert(isIdentical(pow(1.000000001, rninf), +0.0));
745 assert(pow(dinf, 0.001) == dinf);
746 assert(isIdentical(pow(dinf, -0.001), +0.0));
747 assert(pow(rninf, 3.0L) == rninf);
748 assert(pow(rninf, 2.0L) == real.infinity);
749 assert(isIdentical(pow(rninf, -3.0), -0.0));
750 assert(isIdentical(pow(rninf, -2.0), +0.0));
751
752 // @@@BUG@@@ somewhere
version(OSX)753 version (OSX) {} else assert(isNaN(pow(one, dinf)));
version(OSX)754 version (OSX) {} else assert(isNaN(pow(-one, dinf)));
755 assert(isNaN(pow(-0.2, PI)));
756 // boundary cases. Note that epsilon == 2^^-n for some n,
757 // so 1/epsilon == 2^^n is always even.
758 assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
759 assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
760 assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
761 assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
762
763 assert(pow(0.0, -3.0) == double.infinity);
764 assert(pow(-0.0, -3.0) == -double.infinity);
765 assert(pow(0.0, -PI) == double.infinity);
766 assert(pow(-0.0, -PI) == double.infinity);
767 assert(isIdentical(pow(0.0, 5.0), 0.0));
768 assert(isIdentical(pow(-0.0, 5.0), -0.0));
769 assert(isIdentical(pow(0.0, 6.0), 0.0));
770 assert(isIdentical(pow(-0.0, 6.0), 0.0));
771
772 // https://issues.dlang.org/show_bug.cgi?id=14786 fixed
773 immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
774 assert(pow(-1.0L, maxOdd) == -1.0L);
775 assert(pow(-1.0L, -maxOdd) == -1.0L);
776 assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L);
777 assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L);
778 assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L);
779 assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L);
780
781 // Now, actual numbers.
782 assert(isClose(pow(two, three), 8.0));
783 assert(isClose(pow(two, -2.5), 0.1767766953));
784
785 // Test integer to float power.
786 immutable uint twoI = 2;
787 assert(isClose(pow(twoI, three), 8.0));
788 }
789
790 // https://issues.dlang.org/show_bug.cgi?id=20508
791 @safe pure nothrow @nogc unittest
792 {
793 import std.math.traits : isNaN;
794
795 assert(isNaN(pow(-double.infinity, 0.5)));
796
797 assert(isNaN(pow(-real.infinity, real.infinity)));
798 assert(isNaN(pow(-real.infinity, -real.infinity)));
799 assert(isNaN(pow(-real.infinity, 1.234)));
800 assert(isNaN(pow(-real.infinity, -0.751)));
801 assert(pow(-real.infinity, 0.0) == 1.0);
802 }
803
804 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`.
805 *
806 * Params:
807 * x = base
808 * n = exponent
809 * m = modulus
810 *
811 * Returns:
812 * `x` to the power `n`, modulo `m`.
813 * The return type is the largest of `x`'s and `m`'s type.
814 *
815 * The function requires that all values have unsigned types.
816 */
817 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m)
818 if (isUnsigned!F && isUnsigned!G && isUnsigned!H)
819 {
820 import std.meta : AliasSeq;
821
822 alias T = Unqual!(Largest!(F, H));
823 static if (T.sizeof <= 4)
824 {
825 alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof];
826 }
827
mulmod(T a,T b,T c)828 static T mulmod(T a, T b, T c)
829 {
830 static if (T.sizeof == 8)
831 {
832 static T addmod(T a, T b, T c)
833 {
834 b = c - b;
835 if (a >= b)
836 return a - b;
837 else
838 return c - b + a;
839 }
840
841 T result = 0, tmp;
842
843 b %= c;
844 while (a > 0)
845 {
846 if (a & 1)
847 result = addmod(result, b, c);
848
849 a >>= 1;
850 b = addmod(b, b, c);
851 }
852
853 return result;
854 }
855 else
856 {
857 DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b);
858 return result % c;
859 }
860 }
861
862 T base = x, result = 1, modulus = m;
863 Unqual!G exponent = n;
864
865 while (exponent > 0)
866 {
867 if (exponent & 1)
868 result = mulmod(result, base, modulus);
869
870 base = mulmod(base, base, modulus);
871 exponent >>= 1;
872 }
873
874 return result;
875 }
876
877 ///
878 @safe pure nothrow @nogc unittest
879 {
880 assert(powmod(1U, 10U, 3U) == 1);
881 assert(powmod(3U, 2U, 6U) == 3);
882 assert(powmod(5U, 5U, 15U) == 5);
883 assert(powmod(2U, 3U, 5U) == 3);
884 assert(powmod(2U, 4U, 5U) == 1);
885 assert(powmod(2U, 5U, 5U) == 2);
886 }
887
888 @safe pure nothrow @nogc unittest
889 {
890 ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u;
891 assert(powmod(a, b, c) == 95367431640625u);
892 a = 100; b = 7919; c = 18446744073709551557u;
893 assert(powmod(a, b, c) == 18223853583554725198u);
894 a = 117; b = 7919; c = 18446744073709551557u;
895 assert(powmod(a, b, c) == 11493139548346411394u);
896 a = 134; b = 7919; c = 18446744073709551557u;
897 assert(powmod(a, b, c) == 10979163786734356774u);
898 a = 151; b = 7919; c = 18446744073709551557u;
899 assert(powmod(a, b, c) == 7023018419737782840u);
900 a = 168; b = 7919; c = 18446744073709551557u;
901 assert(powmod(a, b, c) == 58082701842386811u);
902 a = 185; b = 7919; c = 18446744073709551557u;
903 assert(powmod(a, b, c) == 17423478386299876798u);
904 a = 202; b = 7919; c = 18446744073709551557u;
905 assert(powmod(a, b, c) == 5522733478579799075u);
906 a = 219; b = 7919; c = 18446744073709551557u;
907 assert(powmod(a, b, c) == 15230218982491623487u);
908 a = 236; b = 7919; c = 18446744073709551557u;
909 assert(powmod(a, b, c) == 5198328724976436000u);
910
911 a = 0; b = 7919; c = 18446744073709551557u;
912 assert(powmod(a, b, c) == 0);
913 a = 123; b = 0; c = 18446744073709551557u;
914 assert(powmod(a, b, c) == 1);
915
916 immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u;
917 assert(powmod(a1, b1, c1) == 3883707345459248860u);
918
919 uint x = 100 ,y = 7919, z = 1844674407u;
920 assert(powmod(x, y, z) == 1613100340u);
921 x = 134; y = 7919; z = 1844674407u;
922 assert(powmod(x, y, z) == 734956622u);
923 x = 151; y = 7919; z = 1844674407u;
924 assert(powmod(x, y, z) == 1738696945u);
925 x = 168; y = 7919; z = 1844674407u;
926 assert(powmod(x, y, z) == 1247580927u);
927 x = 185; y = 7919; z = 1844674407u;
928 assert(powmod(x, y, z) == 1293855176u);
929 x = 202; y = 7919; z = 1844674407u;
930 assert(powmod(x, y, z) == 1566963682u);
931 x = 219; y = 7919; z = 1844674407u;
932 assert(powmod(x, y, z) == 181227807u);
933 x = 236; y = 7919; z = 1844674407u;
934 assert(powmod(x, y, z) == 217988321u);
935 x = 253; y = 7919; z = 1844674407u;
936 assert(powmod(x, y, z) == 1588843243u);
937
938 x = 0; y = 7919; z = 184467u;
939 assert(powmod(x, y, z) == 0);
940 x = 123; y = 0; z = 1844674u;
941 assert(powmod(x, y, z) == 1);
942
943 immutable ubyte x1 = 117;
944 immutable uint y1 = 7919;
945 immutable uint z1 = 1844674407u;
946 auto res = powmod(x1, y1, z1);
947 assert(is(typeof(res) == uint));
948 assert(res == 9479781u);
949
950 immutable ushort x2 = 123;
951 immutable uint y2 = 203;
952 immutable ubyte z2 = 113;
953 auto res2 = powmod(x2, y2, z2);
954 assert(is(typeof(res2) == ushort));
955 assert(res2 == 42u);
956 }
957
958 /**
959 * Calculates e$(SUPERSCRIPT x).
960 *
961 * $(TABLE_SV
962 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) )
963 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
964 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
965 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
966 * )
967 */
pragma(inline,true)968 pragma(inline, true)
969 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe
970 {
971 import std.math.constants : LOG2E;
972
973 version (InlineAsm_X87)
974 {
975 // e^^x = 2^^(LOG2E*x)
976 // (This is valid because the overflow & underflow limits for exp
977 // and exp2 are so similar).
978 if (!__ctfe)
979 return exp2Asm(LOG2E*x);
980 }
981 return expImpl(x);
982 }
983
984 /// ditto
pragma(inline,true)985 pragma(inline, true)
986 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); }
987
988 /// ditto
pragma(inline,true)989 pragma(inline, true)
990 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); }
991
992 ///
993 @safe unittest
994 {
995 import std.math.operations : feqrel;
996 import std.math.constants : E;
997
998 assert(exp(0.0) == 1.0);
999 assert(exp(3.0).feqrel(E * E * E) > 16);
1000 }
1001
expImpl(T)1002 private T expImpl(T)(T x) @safe pure nothrow @nogc
1003 {
1004 import std.math : floatTraits, RealFormat;
1005 import std.math.traits : isNaN;
1006 import std.math.rounding : floor;
1007 import std.math.algebraic : poly;
1008 import std.math.constants : LOG2E;
1009
1010 alias F = floatTraits!T;
1011 static if (F.realFormat == RealFormat.ieeeSingle)
1012 {
1013 static immutable T[6] P = [
1014 5.0000001201E-1,
1015 1.6666665459E-1,
1016 4.1665795894E-2,
1017 8.3334519073E-3,
1018 1.3981999507E-3,
1019 1.9875691500E-4,
1020 ];
1021
1022 enum T C1 = 0.693359375;
1023 enum T C2 = -2.12194440e-4;
1024
1025 // Overflow and Underflow limits.
1026 enum T OF = 88.72283905206835;
1027 enum T UF = -103.278929903431851103; // ln(2^-149)
1028 }
1029 else static if (F.realFormat == RealFormat.ieeeDouble)
1030 {
1031 // Coefficients for exp(x)
1032 static immutable T[3] P = [
1033 9.99999999999999999910E-1L,
1034 3.02994407707441961300E-2L,
1035 1.26177193074810590878E-4L,
1036 ];
1037 static immutable T[4] Q = [
1038 2.00000000000000000009E0L,
1039 2.27265548208155028766E-1L,
1040 2.52448340349684104192E-3L,
1041 3.00198505138664455042E-6L,
1042 ];
1043
1044 // C1 + C2 = LN2.
1045 enum T C1 = 6.93145751953125E-1;
1046 enum T C2 = 1.42860682030941723212E-6;
1047
1048 // Overflow and Underflow limits.
1049 enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024)
1050 enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
1051 }
1052 else static if (F.realFormat == RealFormat.ieeeExtended ||
1053 F.realFormat == RealFormat.ieeeExtended53)
1054 {
1055 // Coefficients for exp(x)
1056 static immutable T[3] P = [
1057 9.9999999999999999991025E-1L,
1058 3.0299440770744196129956E-2L,
1059 1.2617719307481059087798E-4L,
1060 ];
1061 static immutable T[4] Q = [
1062 2.0000000000000000000897E0L,
1063 2.2726554820815502876593E-1L,
1064 2.5244834034968410419224E-3L,
1065 3.0019850513866445504159E-6L,
1066 ];
1067
1068 // C1 + C2 = LN2.
1069 enum T C1 = 6.9314575195312500000000E-1L;
1070 enum T C2 = 1.4286068203094172321215E-6L;
1071
1072 // Overflow and Underflow limits.
1073 enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384)
1074 enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446)
1075 }
1076 else static if (F.realFormat == RealFormat.ieeeQuadruple)
1077 {
1078 // Coefficients for exp(x) - 1
1079 static immutable T[5] P = [
1080 9.999999999999999999999999999999999998502E-1L,
1081 3.508710990737834361215404761139478627390E-2L,
1082 2.708775201978218837374512615596512792224E-4L,
1083 6.141506007208645008909088812338454698548E-7L,
1084 3.279723985560247033712687707263393506266E-10L
1085 ];
1086 static immutable T[6] Q = [
1087 2.000000000000000000000000000000000000150E0,
1088 2.368408864814233538909747618894558968880E-1L,
1089 3.611828913847589925056132680618007270344E-3L,
1090 1.504792651814944826817779302637284053660E-5L,
1091 1.771372078166251484503904874657985291164E-8L,
1092 2.980756652081995192255342779918052538681E-12L
1093 ];
1094
1095 // C1 + C2 = LN2.
1096 enum T C1 = 6.93145751953125E-1L;
1097 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1098
1099 // Overflow and Underflow limits.
1100 enum T OF = 1.135583025911358400418251384584930671458833e4L;
1101 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1102 }
1103 else
1104 static assert(0, "Not implemented for this architecture");
1105
1106 // Special cases.
1107 if (isNaN(x))
1108 return x;
1109 if (x > OF)
1110 return real.infinity;
1111 if (x < UF)
1112 return 0.0;
1113
1114 // Express: e^^x = e^^g * 2^^n
1115 // = e^^g * e^^(n * LOG2E)
1116 // = e^^(g + n * LOG2E)
1117 T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5);
1118 const int n = cast(int) xx;
1119 x -= xx * C1;
1120 x -= xx * C2;
1121
1122 static if (F.realFormat == RealFormat.ieeeSingle)
1123 {
1124 xx = x * x;
1125 x = poly(x, P) * xx + x + 1.0f;
1126 }
1127 else
1128 {
1129 // Rational approximation for exponential of the fractional part:
1130 // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
1131 xx = x * x;
1132 const T px = x * poly(xx, P);
1133 x = px / (poly(xx, Q) - px);
1134 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
1135 }
1136
1137 // Scale by power of 2.
1138 x = core.math.ldexp(x, n);
1139
1140 return x;
1141 }
1142
1143 @safe @nogc nothrow unittest
1144 {
1145 import std.math : floatTraits, RealFormat;
1146 import std.math.operations : NaN, feqrel, isClose;
1147 import std.math.constants : E;
1148 import std.math.traits : isIdentical;
1149 import std.math.algebraic : abs;
1150
1151 version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags;
version(FloatingPointControlSupport)1152 version (FloatingPointControlSupport)
1153 {
1154 import std.math.hardware : FloatingPointControl;
1155
1156 FloatingPointControl ctrl;
1157 if (FloatingPointControl.hasExceptionTraps)
1158 ctrl.disableExceptions(FloatingPointControl.allExceptions);
1159 ctrl.rounding = FloatingPointControl.roundToNearest;
1160 }
1161
testExp(T)1162 static void testExp(T)()
1163 {
1164 enum realFormat = floatTraits!T.realFormat;
1165 static if (realFormat == RealFormat.ieeeQuadruple)
1166 {
1167 static immutable T[2][] exptestpoints =
1168 [ // x exp(x)
1169 [ 1.0L, E ],
1170 [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ],
1171 [ 3.0L, E*E*E ],
1172 [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
1173 [ 0x1.7p+13L, T.infinity ], // close overflow
1174 [ 0x1p+80L, T.infinity ], // far overflow
1175 [ T.infinity, T.infinity ],
1176 [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
1177 [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
1178 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
1179 [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto
1180 [-0x1.655p+13L, 0 ], // close underflow
1181 [-0x1p+30L, 0 ], // far underflow
1182 ];
1183 }
1184 else static if (realFormat == RealFormat.ieeeExtended ||
1185 realFormat == RealFormat.ieeeExtended53)
1186 {
1187 static immutable T[2][] exptestpoints =
1188 [ // x exp(x)
1189 [ 1.0L, E ],
1190 [ 0.5L, 0x1.a61298e1e069bc97p+0L ],
1191 [ 3.0L, E*E*E ],
1192 [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow
1193 [ 0x1.7p+13L, T.infinity ], // close overflow
1194 [ 0x1p+80L, T.infinity ], // far overflow
1195 [ T.infinity, T.infinity ],
1196 [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow
1197 [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto
1198 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal
1199 [-0x1.643p+13L, 0x1p-16444L ], // ditto
1200 [-0x1.645p+13L, 0 ], // close underflow
1201 [-0x1p+30L, 0 ], // far underflow
1202 ];
1203 }
1204 else static if (realFormat == RealFormat.ieeeDouble)
1205 {
1206 static immutable T[2][] exptestpoints =
1207 [ // x, exp(x)
1208 [ 1.0L, E ],
1209 [ 0.5L, 0x1.a61298e1e069cp+0L ],
1210 [ 3.0L, E*E*E ],
1211 [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow
1212 [ 0x1.7p+9L, T.infinity ], // close overflow
1213 [ 0x1p+80L, T.infinity ], // far overflow
1214 [ T.infinity, T.infinity ],
1215 [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow
1216 [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal
1217 [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto
1218 [-0x1.8p+9L, 0 ], // close underflow
1219 [-0x1p+30L, 0 ], // far underflow
1220 ];
1221 }
1222 else static if (realFormat == RealFormat.ieeeSingle)
1223 {
1224 static immutable T[2][] exptestpoints =
1225 [ // x, exp(x)
1226 [ 1.0L, E ],
1227 [ 0.5L, 0x1.a61299p+0L ],
1228 [ 3.0L, E*E*E ],
1229 [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow
1230 [ 0x1.7p+6L, T.infinity ], // close overflow
1231 [ 0x1p+80L, T.infinity ], // far overflow
1232 [ T.infinity, T.infinity ],
1233 [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow
1234 [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal
1235 [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto
1236 [-0x1.ap+6L, 0 ], // close underflow
1237 [-0x1p+30L, 0 ], // far underflow
1238 ];
1239 }
1240 else
1241 static assert(0, "No exp() tests for real type!");
1242
1243 const minEqualMantissaBits = T.mant_dig - 13;
1244 T x;
1245 version (IeeeFlagsSupport) IeeeFlags f;
1246 foreach (ref pair; exptestpoints)
1247 {
1248 version (IeeeFlagsSupport) resetIeeeFlags();
1249 x = exp(pair[0]);
1250 //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]);
1251 assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
1252 }
1253
1254 // Ideally, exp(0) would not set the inexact flag.
1255 // Unfortunately, fldl2e sets it!
1256 // So it's not realistic to avoid setting it.
1257 assert(exp(cast(T) 0.0) == 1.0);
1258
1259 // NaN propagation. Doesn't set flags, bcos was already NaN.
1260 version (IeeeFlagsSupport)
1261 {
1262 resetIeeeFlags();
1263 x = exp(T.nan);
1264 f = ieeeFlags;
1265 assert(isIdentical(abs(x), T.nan));
1266 assert(f.flags == 0);
1267
1268 resetIeeeFlags();
1269 x = exp(-T.nan);
1270 f = ieeeFlags;
1271 assert(isIdentical(abs(x), T.nan));
1272 assert(f.flags == 0);
1273 }
1274 else
1275 {
1276 x = exp(T.nan);
1277 assert(isIdentical(abs(x), T.nan));
1278
1279 x = exp(-T.nan);
1280 assert(isIdentical(abs(x), T.nan));
1281 }
1282
1283 x = exp(NaN(0x123));
1284 assert(isIdentical(x, NaN(0x123)));
1285 }
1286
1287 import std.meta : AliasSeq;
1288 foreach (T; AliasSeq!(real, double, float))
1289 testExp!T();
1290
1291 // High resolution test (verified against GNU MPFR/Mathematica).
1292 assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L);
1293
1294 assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1295 }
1296
1297 /**
1298 * Calculates the value of the natural logarithm base (e)
1299 * raised to the power of x, minus 1.
1300 *
1301 * For very small x, expm1(x) is more accurate
1302 * than exp(x)-1.
1303 *
1304 * $(TABLE_SV
1305 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) )
1306 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
1307 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1308 * $(TR $(TD -$(INFIN)) $(TD -1.0) )
1309 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1310 * )
1311 */
pragma(inline,true)1312 pragma(inline, true)
1313 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe
1314 {
1315 version (InlineAsm_X87)
1316 {
1317 if (!__ctfe)
1318 return expm1Asm(x);
1319 }
1320 return expm1Impl(x);
1321 }
1322
1323 /// ditto
pragma(inline,true)1324 pragma(inline, true)
1325 double expm1(double x) @safe pure nothrow @nogc
1326 {
1327 return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x);
1328 }
1329
1330 /// ditto
pragma(inline,true)1331 pragma(inline, true)
1332 float expm1(float x) @safe pure nothrow @nogc
1333 {
1334 // no single-precision version in Cephes => use double precision
1335 return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x);
1336 }
1337
1338 ///
1339 @safe unittest
1340 {
1341 import std.math.traits : isIdentical;
1342 import std.math.operations : feqrel;
1343
1344 assert(isIdentical(expm1(0.0), 0.0));
1345 assert(expm1(1.0).feqrel(1.71828) > 16);
1346 assert(expm1(2.0).feqrel(6.3890) > 16);
1347 }
1348
version(InlineAsm_X87)1349 version (InlineAsm_X87)
1350 private real expm1Asm(real x) @trusted pure nothrow @nogc
1351 {
1352 version (X86)
1353 {
1354 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1355 asm pure nothrow @nogc
1356 {
1357 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1358 * Author: Don Clugston.
1359 *
1360 * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
1361 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
1362 * and 2ym1 = (2^^(y-rndint(y))-1).
1363 * If 2rndy < 0.5*real.epsilon, result is -1.
1364 * Implementation is otherwise the same as for exp2()
1365 */
1366 naked;
1367 fld real ptr [ESP+4] ; // x
1368 mov AX, [ESP+4+8]; // AX = exponent and sign
1369 sub ESP, 12+8; // Create scratch space on the stack
1370 // [ESP,ESP+2] = scratchint
1371 // [ESP+4..+6, +8..+10, +10] = scratchreal
1372 // set scratchreal mantissa = 1.0
1373 mov dword ptr [ESP+8], 0;
1374 mov dword ptr [ESP+8+4], 0x80000000;
1375 and AX, 0x7FFF; // drop sign bit
1376 cmp AX, 0x401D; // avoid InvalidException in fist
1377 jae L_extreme;
1378 fldl2e;
1379 fmulp ST(1), ST; // y = x*log2(e)
1380 fist dword ptr [ESP]; // scratchint = rndint(y)
1381 fisub dword ptr [ESP]; // y - rndint(y)
1382 // and now set scratchreal exponent
1383 mov EAX, [ESP];
1384 add EAX, 0x3fff;
1385 jle short L_largenegative;
1386 cmp EAX,0x8000;
1387 jge short L_largepositive;
1388 mov [ESP+8+8],AX;
1389 f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1390 fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1391 fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1
1392 fld1;
1393 fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1394 faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1395 add ESP,12+8;
1396 ret PARAMSIZE;
1397
1398 L_extreme: // Extreme exponent. X is very large positive, very
1399 // large negative, infinity, or NaN.
1400 fxam;
1401 fstsw AX;
1402 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1403 jz L_was_nan; // if x is NaN, returns x
1404 test AX, 0x0200;
1405 jnz L_largenegative;
1406 L_largepositive:
1407 // Set scratchreal = real.max.
1408 // squaring it will create infinity, and set overflow flag.
1409 mov word ptr [ESP+8+8], 0x7FFE;
1410 fstp ST(0);
1411 fld real ptr [ESP+8]; // load scratchreal
1412 fmul ST(0), ST; // square it, to create havoc!
1413 L_was_nan:
1414 add ESP,12+8;
1415 ret PARAMSIZE;
1416 L_largenegative:
1417 fstp ST(0);
1418 fld1;
1419 fchs; // return -1. Underflow flag is not set.
1420 add ESP,12+8;
1421 ret PARAMSIZE;
1422 }
1423 }
1424 else version (X86_64)
1425 {
1426 asm pure nothrow @nogc
1427 {
1428 naked;
1429 }
1430 version (Win64)
1431 {
1432 asm pure nothrow @nogc
1433 {
1434 fld real ptr [RCX]; // x
1435 mov AX,[RCX+8]; // AX = exponent and sign
1436 }
1437 }
1438 else
1439 {
1440 asm pure nothrow @nogc
1441 {
1442 fld real ptr [RSP+8]; // x
1443 mov AX,[RSP+8+8]; // AX = exponent and sign
1444 }
1445 }
1446 asm pure nothrow @nogc
1447 {
1448 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1449 * Author: Don Clugston.
1450 *
1451 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1452 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1453 * and 2ym1 = (2^(y-rndint(y))-1).
1454 * If 2rndy < 0.5*real.epsilon, result is -1.
1455 * Implementation is otherwise the same as for exp2()
1456 */
1457 sub RSP, 24; // Create scratch space on the stack
1458 // [RSP,RSP+2] = scratchint
1459 // [RSP+4..+6, +8..+10, +10] = scratchreal
1460 // set scratchreal mantissa = 1.0
1461 mov dword ptr [RSP+8], 0;
1462 mov dword ptr [RSP+8+4], 0x80000000;
1463 and AX, 0x7FFF; // drop sign bit
1464 cmp AX, 0x401D; // avoid InvalidException in fist
1465 jae L_extreme;
1466 fldl2e;
1467 fmul ; // y = x*log2(e)
1468 fist dword ptr [RSP]; // scratchint = rndint(y)
1469 fisub dword ptr [RSP]; // y - rndint(y)
1470 // and now set scratchreal exponent
1471 mov EAX, [RSP];
1472 add EAX, 0x3fff;
1473 jle short L_largenegative;
1474 cmp EAX,0x8000;
1475 jge short L_largepositive;
1476 mov [RSP+8+8],AX;
1477 f2xm1; // 2^(y-rndint(y)) -1
1478 fld real ptr [RSP+8] ; // 2^rndint(y)
1479 fmul ST(1), ST;
1480 fld1;
1481 fsubp ST(1), ST;
1482 fadd;
1483 add RSP,24;
1484 ret;
1485
1486 L_extreme: // Extreme exponent. X is very large positive, very
1487 // large negative, infinity, or NaN.
1488 fxam;
1489 fstsw AX;
1490 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1491 jz L_was_nan; // if x is NaN, returns x
1492 test AX, 0x0200;
1493 jnz L_largenegative;
1494 L_largepositive:
1495 // Set scratchreal = real.max.
1496 // squaring it will create infinity, and set overflow flag.
1497 mov word ptr [RSP+8+8], 0x7FFE;
1498 fstp ST(0);
1499 fld real ptr [RSP+8]; // load scratchreal
1500 fmul ST(0), ST; // square it, to create havoc!
1501 L_was_nan:
1502 add RSP,24;
1503 ret;
1504
1505 L_largenegative:
1506 fstp ST(0);
1507 fld1;
1508 fchs; // return -1. Underflow flag is not set.
1509 add RSP,24;
1510 ret;
1511 }
1512 }
1513 else
1514 static assert(0);
1515 }
1516
expm1Impl(T)1517 private T expm1Impl(T)(T x) @safe pure nothrow @nogc
1518 {
1519 import std.math : floatTraits, RealFormat;
1520 import std.math.rounding : floor;
1521 import std.math.algebraic : poly;
1522 import std.math.constants : LN2;
1523
1524 // Coefficients for exp(x) - 1 and overflow/underflow limits.
1525 enum realFormat = floatTraits!T.realFormat;
1526 static if (realFormat == RealFormat.ieeeQuadruple)
1527 {
1528 static immutable T[8] P = [
1529 2.943520915569954073888921213330863757240E8L,
1530 -5.722847283900608941516165725053359168840E7L,
1531 8.944630806357575461578107295909719817253E6L,
1532 -7.212432713558031519943281748462837065308E5L,
1533 4.578962475841642634225390068461943438441E4L,
1534 -1.716772506388927649032068540558788106762E3L,
1535 4.401308817383362136048032038528753151144E1L,
1536 -4.888737542888633647784737721812546636240E-1L
1537 ];
1538
1539 static immutable T[9] Q = [
1540 1.766112549341972444333352727998584753865E9L,
1541 -7.848989743695296475743081255027098295771E8L,
1542 1.615869009634292424463780387327037251069E8L,
1543 -2.019684072836541751428967854947019415698E7L,
1544 1.682912729190313538934190635536631941751E6L,
1545 -9.615511549171441430850103489315371768998E4L,
1546 3.697714952261803935521187272204485251835E3L,
1547 -8.802340681794263968892934703309274564037E1L,
1548 1.0
1549 ];
1550
1551 enum T OF = 1.1356523406294143949491931077970764891253E4L;
1552 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1553 }
1554 else static if (realFormat == RealFormat.ieeeExtended)
1555 {
1556 static immutable T[5] P = [
1557 -1.586135578666346600772998894928250240826E4L,
1558 2.642771505685952966904660652518429479531E3L,
1559 -3.423199068835684263987132888286791620673E2L,
1560 1.800826371455042224581246202420972737840E1L,
1561 -5.238523121205561042771939008061958820811E-1L,
1562 ];
1563 static immutable T[6] Q = [
1564 -9.516813471998079611319047060563358064497E4L,
1565 3.964866271411091674556850458227710004570E4L,
1566 -7.207678383830091850230366618190187434796E3L,
1567 7.206038318724600171970199625081491823079E2L,
1568 -4.002027679107076077238836622982900945173E1L,
1569 1.0
1570 ];
1571
1572 enum T OF = 1.1356523406294143949492E4L;
1573 enum T UF = -4.5054566736396445112120088E1L;
1574 }
1575 else static if (realFormat == RealFormat.ieeeDouble)
1576 {
1577 static immutable T[3] P = [
1578 9.9999999999999999991025E-1,
1579 3.0299440770744196129956E-2,
1580 1.2617719307481059087798E-4,
1581 ];
1582 static immutable T[4] Q = [
1583 2.0000000000000000000897E0,
1584 2.2726554820815502876593E-1,
1585 2.5244834034968410419224E-3,
1586 3.0019850513866445504159E-6,
1587 ];
1588 }
1589 else
1590 static assert(0, "no coefficients for expm1()");
1591
1592 static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
1593 {
1594 if (x < -0.5 || x > 0.5)
1595 return exp(x) - 1.0;
1596 if (x == 0.0)
1597 return x;
1598
1599 const T xx = x * x;
1600 x = x * poly(xx, P);
1601 x = x / (poly(xx, Q) - x);
1602 return x + x;
1603 }
1604 else
1605 {
1606 // C1 + C2 = LN2.
1607 enum T C1 = 6.9314575195312500000000E-1L;
1608 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1609
1610 // Special cases.
1611 if (x > OF)
1612 return real.infinity;
1613 if (x == cast(T) 0.0)
1614 return x;
1615 if (x < UF)
1616 return -1.0;
1617
1618 // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
1619 int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2);
1620 x -= n * C1;
1621 x -= n * C2;
1622
1623 // Rational approximation:
1624 // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
1625 T px = x * poly(x, P);
1626 T qx = poly(x, Q);
1627 const T xx = x * x;
1628 qx = x + ((cast(T) 0.5) * xx + xx * px / qx);
1629
1630 // We have qx = exp(remainder LN2) - 1, so:
1631 // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
1632 px = core.math.ldexp(cast(T) 1.0, n);
1633 x = px * qx + (px - cast(T) 1.0);
1634
1635 return x;
1636 }
1637 }
1638
1639 @safe @nogc nothrow unittest
1640 {
1641 import std.math.traits : isNaN;
1642 import std.math.operations : isClose, CommonDefaultFor;
1643
testExpm1(T)1644 static void testExpm1(T)()
1645 {
1646 // NaN
1647 assert(isNaN(expm1(cast(T) T.nan)));
1648
1649 static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ];
1650 foreach (x; xs)
1651 {
1652 const T e = expm1(x);
1653 const T r = exp(x) - 1;
1654
1655 //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
1656 assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
1657 }
1658 }
1659
1660 import std.meta : AliasSeq;
1661 foreach (T; AliasSeq!(real, double))
1662 testExpm1!T();
1663 }
1664
1665 /**
1666 * Calculates 2$(SUPERSCRIPT x).
1667 *
1668 * $(TABLE_SV
1669 * $(TR $(TH x) $(TH exp2(x)) )
1670 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1671 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
1672 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1673 * )
1674 */
pragma(inline,true)1675 pragma(inline, true)
1676 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe
1677 {
1678 version (InlineAsm_X87)
1679 {
1680 if (!__ctfe)
1681 return exp2Asm(x);
1682 }
1683 return exp2Impl(x);
1684 }
1685
1686 /// ditto
pragma(inline,true)1687 pragma(inline, true)
1688 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); }
1689
1690 /// ditto
pragma(inline,true)1691 pragma(inline, true)
1692 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); }
1693
1694 ///
1695 @safe unittest
1696 {
1697 import std.math.traits : isIdentical;
1698 import std.math.operations : feqrel;
1699
1700 assert(isIdentical(exp2(0.0), 1.0));
1701 assert(exp2(2.0).feqrel(4.0) > 16);
1702 assert(exp2(8.0).feqrel(256.0) > 16);
1703 }
1704
1705 @safe unittest
1706 {
version(CRuntime_Microsoft)1707 version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
1708 {
1709 assert( core.stdc.math.exp2f(0.0f) == 1 );
1710 assert( core.stdc.math.exp2 (0.0) == 1 );
1711 assert( core.stdc.math.exp2l(0.0L) == 1 );
1712 }
1713 }
1714
version(InlineAsm_X87)1715 version (InlineAsm_X87)
1716 private real exp2Asm(real x) @nogc @trusted pure nothrow
1717 {
1718 version (X86)
1719 {
1720 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1721
1722 asm pure nothrow @nogc
1723 {
1724 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1725 * Author: Don Clugston.
1726 *
1727 * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1728 * The trick for high performance is to avoid the fscale(28cycles on core2),
1729 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1730 *
1731 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1732 * because it will set the Invalid Operation flag if overflow or NaN occurs.
1733 * Fortunately, whenever this happens the result would be zero or infinity.
1734 *
1735 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1736 * work for the (very rare) cases where the result is subnormal. So we fall back
1737 * to the slow method in that case.
1738 */
1739 naked;
1740 fld real ptr [ESP+4] ; // x
1741 mov AX, [ESP+4+8]; // AX = exponent and sign
1742 sub ESP, 12+8; // Create scratch space on the stack
1743 // [ESP,ESP+2] = scratchint
1744 // [ESP+4..+6, +8..+10, +10] = scratchreal
1745 // set scratchreal mantissa = 1.0
1746 mov dword ptr [ESP+8], 0;
1747 mov dword ptr [ESP+8+4], 0x80000000;
1748 and AX, 0x7FFF; // drop sign bit
1749 cmp AX, 0x401D; // avoid InvalidException in fist
1750 jae L_extreme;
1751 fist dword ptr [ESP]; // scratchint = rndint(x)
1752 fisub dword ptr [ESP]; // x - rndint(x)
1753 // and now set scratchreal exponent
1754 mov EAX, [ESP];
1755 add EAX, 0x3fff;
1756 jle short L_subnormal;
1757 cmp EAX,0x8000;
1758 jge short L_overflow;
1759 mov [ESP+8+8],AX;
1760 L_normal:
1761 f2xm1;
1762 fld1;
1763 faddp ST(1), ST; // 2^^(x-rndint(x))
1764 fld real ptr [ESP+8] ; // 2^^rndint(x)
1765 add ESP,12+8;
1766 fmulp ST(1), ST;
1767 ret PARAMSIZE;
1768
1769 L_subnormal:
1770 // Result will be subnormal.
1771 // In this rare case, the simple poking method doesn't work.
1772 // The speed doesn't matter, so use the slow fscale method.
1773 fild dword ptr [ESP]; // scratchint
1774 fld1;
1775 fscale;
1776 fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1777 fstp ST(0); // drop scratchint
1778 jmp L_normal;
1779
1780 L_extreme: // Extreme exponent. X is very large positive, very
1781 // large negative, infinity, or NaN.
1782 fxam;
1783 fstsw AX;
1784 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1785 jz L_was_nan; // if x is NaN, returns x
1786 // set scratchreal = real.min_normal
1787 // squaring it will return 0, setting underflow flag
1788 mov word ptr [ESP+8+8], 1;
1789 test AX, 0x0200;
1790 jnz L_waslargenegative;
1791 L_overflow:
1792 // Set scratchreal = real.max.
1793 // squaring it will create infinity, and set overflow flag.
1794 mov word ptr [ESP+8+8], 0x7FFE;
1795 L_waslargenegative:
1796 fstp ST(0);
1797 fld real ptr [ESP+8]; // load scratchreal
1798 fmul ST(0), ST; // square it, to create havoc!
1799 L_was_nan:
1800 add ESP,12+8;
1801 ret PARAMSIZE;
1802 }
1803 }
1804 else version (X86_64)
1805 {
1806 asm pure nothrow @nogc
1807 {
1808 naked;
1809 }
1810 version (Win64)
1811 {
1812 asm pure nothrow @nogc
1813 {
1814 fld real ptr [RCX]; // x
1815 mov AX,[RCX+8]; // AX = exponent and sign
1816 }
1817 }
1818 else
1819 {
1820 asm pure nothrow @nogc
1821 {
1822 fld real ptr [RSP+8]; // x
1823 mov AX,[RSP+8+8]; // AX = exponent and sign
1824 }
1825 }
1826 asm pure nothrow @nogc
1827 {
1828 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1829 * Author: Don Clugston.
1830 *
1831 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1832 * The trick for high performance is to avoid the fscale(28cycles on core2),
1833 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1834 *
1835 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1836 * because it will set the Invalid Operation flag is overflow or NaN occurs.
1837 * Fortunately, whenever this happens the result would be zero or infinity.
1838 *
1839 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1840 * work for the (very rare) cases where the result is subnormal. So we fall back
1841 * to the slow method in that case.
1842 */
1843 sub RSP, 24; // Create scratch space on the stack
1844 // [RSP,RSP+2] = scratchint
1845 // [RSP+4..+6, +8..+10, +10] = scratchreal
1846 // set scratchreal mantissa = 1.0
1847 mov dword ptr [RSP+8], 0;
1848 mov dword ptr [RSP+8+4], 0x80000000;
1849 and AX, 0x7FFF; // drop sign bit
1850 cmp AX, 0x401D; // avoid InvalidException in fist
1851 jae L_extreme;
1852 fist dword ptr [RSP]; // scratchint = rndint(x)
1853 fisub dword ptr [RSP]; // x - rndint(x)
1854 // and now set scratchreal exponent
1855 mov EAX, [RSP];
1856 add EAX, 0x3fff;
1857 jle short L_subnormal;
1858 cmp EAX,0x8000;
1859 jge short L_overflow;
1860 mov [RSP+8+8],AX;
1861 L_normal:
1862 f2xm1;
1863 fld1;
1864 fadd; // 2^(x-rndint(x))
1865 fld real ptr [RSP+8] ; // 2^rndint(x)
1866 add RSP,24;
1867 fmulp ST(1), ST;
1868 ret;
1869
1870 L_subnormal:
1871 // Result will be subnormal.
1872 // In this rare case, the simple poking method doesn't work.
1873 // The speed doesn't matter, so use the slow fscale method.
1874 fild dword ptr [RSP]; // scratchint
1875 fld1;
1876 fscale;
1877 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1878 fstp ST(0); // drop scratchint
1879 jmp L_normal;
1880
1881 L_extreme: // Extreme exponent. X is very large positive, very
1882 // large negative, infinity, or NaN.
1883 fxam;
1884 fstsw AX;
1885 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1886 jz L_was_nan; // if x is NaN, returns x
1887 // set scratchreal = real.min
1888 // squaring it will return 0, setting underflow flag
1889 mov word ptr [RSP+8+8], 1;
1890 test AX, 0x0200;
1891 jnz L_waslargenegative;
1892 L_overflow:
1893 // Set scratchreal = real.max.
1894 // squaring it will create infinity, and set overflow flag.
1895 mov word ptr [RSP+8+8], 0x7FFE;
1896 L_waslargenegative:
1897 fstp ST(0);
1898 fld real ptr [RSP+8]; // load scratchreal
1899 fmul ST(0), ST; // square it, to create havoc!
1900 L_was_nan:
1901 add RSP,24;
1902 ret;
1903 }
1904 }
1905 else
1906 static assert(0);
1907 }
1908
exp2Impl(T)1909 private T exp2Impl(T)(T x) @nogc @safe pure nothrow
1910 {
1911 import std.math : floatTraits, RealFormat;
1912 import std.math.traits : isNaN;
1913 import std.math.rounding : floor;
1914 import std.math.algebraic : poly;
1915
1916 // Coefficients for exp2(x)
1917 enum realFormat = floatTraits!T.realFormat;
1918 static if (realFormat == RealFormat.ieeeQuadruple)
1919 {
1920 static immutable T[5] P = [
1921 9.079594442980146270952372234833529694788E12L,
1922 1.530625323728429161131811299626419117557E11L,
1923 5.677513871931844661829755443994214173883E8L,
1924 6.185032670011643762127954396427045467506E5L,
1925 1.587171580015525194694938306936721666031E2L
1926 ];
1927
1928 static immutable T[6] Q = [
1929 2.619817175234089411411070339065679229869E13L,
1930 1.490560994263653042761789432690793026977E12L,
1931 1.092141473886177435056423606755843616331E10L,
1932 2.186249607051644894762167991800811827835E7L,
1933 1.236602014442099053716561665053645270207E4L,
1934 1.0
1935 ];
1936 }
1937 else static if (realFormat == RealFormat.ieeeExtended)
1938 {
1939 static immutable T[3] P = [
1940 2.0803843631901852422887E6L,
1941 3.0286971917562792508623E4L,
1942 6.0614853552242266094567E1L,
1943 ];
1944 static immutable T[4] Q = [
1945 6.0027204078348487957118E6L,
1946 3.2772515434906797273099E5L,
1947 1.7492876999891839021063E3L,
1948 1.0000000000000000000000E0L,
1949 ];
1950 }
1951 else static if (realFormat == RealFormat.ieeeDouble)
1952 {
1953 static immutable T[3] P = [
1954 1.51390680115615096133E3L,
1955 2.02020656693165307700E1L,
1956 2.30933477057345225087E-2L,
1957 ];
1958 static immutable T[3] Q = [
1959 4.36821166879210612817E3L,
1960 2.33184211722314911771E2L,
1961 1.00000000000000000000E0L,
1962 ];
1963 }
1964 else static if (realFormat == RealFormat.ieeeSingle)
1965 {
1966 static immutable T[6] P = [
1967 6.931472028550421E-001L,
1968 2.402264791363012E-001L,
1969 5.550332471162809E-002L,
1970 9.618437357674640E-003L,
1971 1.339887440266574E-003L,
1972 1.535336188319500E-004L,
1973 ];
1974 }
1975 else
1976 static assert(0, "no coefficients for exp2()");
1977
1978 // Overflow and Underflow limits.
1979 enum T OF = T.max_exp;
1980 enum T UF = T.min_exp - 1;
1981
1982 // Special cases.
1983 if (isNaN(x))
1984 return x;
1985 if (x > OF)
1986 return real.infinity;
1987 if (x < UF)
1988 return 0.0;
1989
1990 static if (realFormat == RealFormat.ieeeSingle) // special case for single precision
1991 {
1992 // The following is necessary because range reduction blows up.
1993 if (x == 0.0f)
1994 return 1.0f;
1995
1996 // Separate into integer and fractional parts.
1997 const T i = floor(x);
1998 int n = cast(int) i;
1999 x -= i;
2000 if (x > 0.5f)
2001 {
2002 n += 1;
2003 x -= 1.0f;
2004 }
2005
2006 // Rational approximation:
2007 // exp2(x) = 1.0 + x P(x)
2008 x = 1.0f + x * poly(x, P);
2009 }
2010 else
2011 {
2012 // Separate into integer and fractional parts.
2013 const T i = floor(x + cast(T) 0.5);
2014 int n = cast(int) i;
2015 x -= i;
2016
2017 // Rational approximation:
2018 // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
2019 const T xx = x * x;
2020 const T px = x * poly(xx, P);
2021 x = px / (poly(xx, Q) - px);
2022 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
2023 }
2024
2025 // Scale by power of 2.
2026 x = core.math.ldexp(x, n);
2027
2028 return x;
2029 }
2030
2031 @safe @nogc nothrow unittest
2032 {
2033 import std.math.operations : feqrel, NaN, isClose;
2034 import std.math.traits : isIdentical;
2035 import std.math.constants : SQRT2;
2036
2037 assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
2038 assert(exp2(8.0L) == 256.0);
2039 assert(exp2(-9.0L)== 1.0L/512.0);
2040
testExp2(T)2041 static void testExp2(T)()
2042 {
2043 // NaN
2044 const T specialNaN = NaN(0x0123L);
2045 assert(isIdentical(exp2(specialNaN), specialNaN));
2046
2047 // over-/underflow
2048 enum T OF = T.max_exp;
2049 enum T UF = T.min_exp - T.mant_dig;
2050 assert(isIdentical(exp2(OF + 1), cast(T) T.infinity));
2051 assert(isIdentical(exp2(UF - 1), cast(T) 0.0));
2052
2053 static immutable T[2][] vals =
2054 [
2055 // x, exp2(x)
2056 [ 0.0, 1.0 ],
2057 [ -0.0, 1.0 ],
2058 [ 0.5, SQRT2 ],
2059 [ 8.0, 256.0 ],
2060 [ -9.0, 1.0 / 512 ],
2061 ];
2062
2063 foreach (ref val; vals)
2064 {
2065 const T x = val[0];
2066 const T r = val[1];
2067 const T e = exp2(x);
2068
2069 //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
2070 assert(isClose(r, e));
2071 }
2072 }
2073
2074 import std.meta : AliasSeq;
2075 foreach (T; AliasSeq!(real, double, float))
2076 testExp2!T();
2077 }
2078
2079 /*********************************************************************
2080 * Separate floating point value into significand and exponent.
2081 *
2082 * Returns:
2083 * Calculate and return $(I x) and $(I exp) such that
2084 * value =$(I x)*2$(SUPERSCRIPT exp) and
2085 * .5 $(LT)= |$(I x)| $(LT) 1.0
2086 *
2087 * $(I x) has same sign as value.
2088 *
2089 * $(TABLE_SV
2090 * $(TR $(TH value) $(TH returns) $(TH exp))
2091 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0))
2092 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max))
2093 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min))
2094 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
2095 * )
2096 */
2097 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
2098 if (isFloatingPoint!T)
2099 {
2100 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2101 import std.math.traits : isSubnormal;
2102
2103 if (__ctfe)
2104 {
2105 // Handle special cases.
2106 if (value == 0) { exp = 0; return value; }
2107 if (value == T.infinity) { exp = int.max; return value; }
2108 if (value == -T.infinity || value != value) { exp = int.min; return value; }
2109 // Handle ordinary cases.
2110 // In CTFE there is no performance advantage for having separate
2111 // paths for different floating point types.
2112 T absValue = value < 0 ? -value : value;
2113 int expCount;
2114 static if (T.mant_dig > double.mant_dig)
2115 {
2116 for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L)
2117 expCount += 1024;
2118 for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L)
2119 expCount -= 1021;
2120 }
2121 const double dval = cast(double) absValue;
2122 int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2;
2123 dexp++;
2124 expCount += dexp;
2125 absValue *= 2.0 ^^ -dexp;
2126 // If the original value was subnormal or if it was a real
2127 // then absValue can still be outside the [0.5, 1.0) range.
2128 if (absValue < 0.5)
2129 {
2130 assert(T.mant_dig > double.mant_dig || isSubnormal(value));
2131 do
2132 {
2133 absValue += absValue;
2134 expCount--;
2135 } while (absValue < 0.5);
2136 }
2137 else
2138 {
2139 assert(absValue < 1 || T.mant_dig > double.mant_dig);
2140 for (; absValue >= 1; absValue *= T(0.5))
2141 expCount++;
2142 }
2143 exp = expCount;
2144 return value < 0 ? -absValue : absValue;
2145 }
2146
2147 Unqual!T vf = value;
2148 ushort* vu = cast(ushort*)&vf;
2149 static if (is(immutable T == immutable float))
2150 int* vi = cast(int*)&vf;
2151 else
2152 long* vl = cast(long*)&vf;
2153 int ex;
2154 alias F = floatTraits!T;
2155
2156 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2157 static if (F.realFormat == RealFormat.ieeeExtended ||
2158 F.realFormat == RealFormat.ieeeExtended53)
2159 {
2160 if (ex)
2161 { // If exponent is non-zero
2162 if (ex == F.EXPMASK) // infinity or NaN
2163 {
2164 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2165 {
2166 *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ
2167 exp = int.min;
2168 }
2169 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2170 exp = int.min;
2171 else // positive infinity
2172 exp = int.max;
2173
2174 }
2175 else
2176 {
2177 exp = ex - F.EXPBIAS;
2178 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2179 }
2180 }
2181 else if (!*vl)
2182 {
2183 // vf is +-0.0
2184 exp = 0;
2185 }
2186 else
2187 {
2188 // subnormal
2189
2190 vf *= F.RECIP_EPSILON;
2191 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2192 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2193 vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2194 }
2195 return vf;
2196 }
2197 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2198 {
2199 if (ex) // If exponent is non-zero
2200 {
2201 if (ex == F.EXPMASK)
2202 {
2203 // infinity or NaN
2204 if (vl[MANTISSA_LSB] |
2205 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2206 {
2207 // convert NaNS to NaNQ
2208 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
2209 exp = int.min;
2210 }
2211 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2212 exp = int.min;
2213 else // positive infinity
2214 exp = int.max;
2215 }
2216 else
2217 {
2218 exp = ex - F.EXPBIAS;
2219 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2220 }
2221 }
2222 else if ((vl[MANTISSA_LSB] |
2223 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2224 {
2225 // vf is +-0.0
2226 exp = 0;
2227 }
2228 else
2229 {
2230 // subnormal
2231 vf *= F.RECIP_EPSILON;
2232 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2233 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2234 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2235 }
2236 return vf;
2237 }
2238 else static if (F.realFormat == RealFormat.ieeeDouble)
2239 {
2240 if (ex) // If exponent is non-zero
2241 {
2242 if (ex == F.EXPMASK) // infinity or NaN
2243 {
2244 if (*vl == 0x7FF0_0000_0000_0000) // positive infinity
2245 {
2246 exp = int.max;
2247 }
2248 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
2249 exp = int.min;
2250 else
2251 { // NaN
2252 *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ
2253 exp = int.min;
2254 }
2255 }
2256 else
2257 {
2258 exp = (ex - F.EXPBIAS) >> 4;
2259 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2260 }
2261 }
2262 else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
2263 {
2264 // vf is +-0.0
2265 exp = 0;
2266 }
2267 else
2268 {
2269 // subnormal
2270 vf *= F.RECIP_EPSILON;
2271 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2272 exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
2273 vu[F.EXPPOS_SHORT] =
2274 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2275 }
2276 return vf;
2277 }
2278 else static if (F.realFormat == RealFormat.ieeeSingle)
2279 {
2280 if (ex) // If exponent is non-zero
2281 {
2282 if (ex == F.EXPMASK) // infinity or NaN
2283 {
2284 if (*vi == 0x7F80_0000) // positive infinity
2285 {
2286 exp = int.max;
2287 }
2288 else if (*vi == 0xFF80_0000) // negative infinity
2289 exp = int.min;
2290 else
2291 { // NaN
2292 *vi |= 0x0040_0000; // convert NaNS to NaNQ
2293 exp = int.min;
2294 }
2295 }
2296 else
2297 {
2298 exp = (ex - F.EXPBIAS) >> 7;
2299 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
2300 }
2301 }
2302 else if (!(*vi & 0x7FFF_FFFF))
2303 {
2304 // vf is +-0.0
2305 exp = 0;
2306 }
2307 else
2308 {
2309 // subnormal
2310 vf *= F.RECIP_EPSILON;
2311 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2312 exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
2313 vu[F.EXPPOS_SHORT] =
2314 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
2315 }
2316 return vf;
2317 }
2318 else // static if (F.realFormat == RealFormat.ibmExtended)
2319 {
2320 assert(0, "frexp not implemented");
2321 }
2322 }
2323
2324 ///
2325 @safe unittest
2326 {
2327 import std.math.operations : isClose;
2328
2329 int exp;
2330 real mantissa = frexp(123.456L, exp);
2331
2332 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L));
2333
2334 assert(frexp(-real.nan, exp) && exp == int.min);
2335 assert(frexp(real.nan, exp) && exp == int.min);
2336 assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
2337 assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
2338 assert(frexp(-0.0, exp) == -0.0 && exp == 0);
2339 assert(frexp(0.0, exp) == 0.0 && exp == 0);
2340 }
2341
2342 @safe @nogc nothrow unittest
2343 {
2344 import std.math.operations : isClose;
2345
2346 int exp;
2347 real mantissa = frexp(123.456L, exp);
2348
2349 // check if values are equal to 19 decimal digits of precision
2350 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18));
2351 }
2352
2353 @safe unittest
2354 {
2355 import std.math : floatTraits, RealFormat;
2356 import std.math.traits : isIdentical;
2357 import std.meta : AliasSeq;
2358 import std.typecons : tuple, Tuple;
2359
2360 static foreach (T; AliasSeq!(real, double, float))
2361 {{
2362 Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2363 tuple(T(0.0), T( 0.0 ), 0),
2364 tuple(T(-0.0), T( -0.0), 0),
2365 tuple(T(1.0), T( .5 ), 1),
2366 tuple(T(-1.0), T( -.5 ), 1),
2367 tuple(T(2.0), T( .5 ), 2),
2368 tuple(T(float.min_normal/2.0f), T(.5), -126),
2369 tuple(T.infinity, T.infinity, int.max),
2370 tuple(-T.infinity, -T.infinity, int.min),
2371 tuple(T.nan, T.nan, int.min),
2372 tuple(-T.nan, -T.nan, int.min),
2373
2374 // https://issues.dlang.org/show_bug.cgi?id=16026:
2375 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2376 ];
2377
foreach(elem;vals)2378 foreach (elem; vals)
2379 {
2380 T x = elem[0];
2381 T e = elem[1];
2382 int exp = elem[2];
2383 int eptr;
2384 T v = frexp(x, eptr);
2385 assert(isIdentical(e, v));
2386 assert(exp == eptr);
2387 }
2388
2389 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2390 {
2391 static T[3][] extendedvals = [ // x,frexp,exp
2392 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2393 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2394 [T.min_normal, .5, -16381],
2395 [T.min_normal/2.0L, .5, -16382] // subnormal
2396 ];
foreach(elem;extendedvals)2397 foreach (elem; extendedvals)
2398 {
2399 T x = elem[0];
2400 T e = elem[1];
2401 int exp = cast(int) elem[2];
2402 int eptr;
2403 T v = frexp(x, eptr);
2404 assert(isIdentical(e, v));
2405 assert(exp == eptr);
2406 }
2407 }
2408 }}
2409
2410 // CTFE
2411 alias CtfeFrexpResult= Tuple!(real, int);
ctfeFrexp(T)2412 static CtfeFrexpResult ctfeFrexp(T)(const T value)
2413 {
2414 int exp;
2415 auto significand = frexp(value, exp);
2416 return CtfeFrexpResult(significand, exp);
2417 }
2418 static foreach (T; AliasSeq!(real, double, float))
2419 {{
2420 enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2421 tuple(T(0.0), T( 0.0 ), 0),
2422 tuple(T(-0.0), T( -0.0), 0),
2423 tuple(T(1.0), T( .5 ), 1),
2424 tuple(T(-1.0), T( -.5 ), 1),
2425 tuple(T(2.0), T( .5 ), 2),
2426 tuple(T(float.min_normal/2.0f), T(.5), -126),
2427 tuple(T.infinity, T.infinity, int.max),
2428 tuple(-T.infinity, -T.infinity, int.min),
2429 tuple(T.nan, T.nan, int.min),
2430 tuple(-T.nan, -T.nan, int.min),
2431
2432 // https://issues.dlang.org/show_bug.cgi?id=16026:
2433 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2434 ];
2435
foreach(elem;vals)2436 static foreach (elem; vals)
2437 {
2438 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2]));
2439 }
2440
2441 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2442 {
2443 enum T[3][] extendedvals = [ // x,frexp,exp
2444 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2445 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2446 [T.min_normal, .5, -16381],
2447 [T.min_normal/2.0L, .5, -16382] // subnormal
2448 ];
foreach(elem;extendedvals)2449 static foreach (elem; extendedvals)
2450 {
2451 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2]));
2452 }
2453 }
2454 }}
2455 }
2456
2457 @safe unittest
2458 {
2459 import std.meta : AliasSeq;
foo()2460 void foo() {
2461 static foreach (T; AliasSeq!(real, double, float))
2462 {{
2463 int exp;
2464 const T a = 1;
2465 immutable T b = 2;
2466 auto c = frexp(a, exp);
2467 auto d = frexp(b, exp);
2468 }}
2469 }
2470 }
2471
2472 /******************************************
2473 * Extracts the exponent of x as a signed integral value.
2474 *
2475 * If x is not a special value, the result is the same as
2476 * $(D cast(int) logb(x)).
2477 *
2478 * $(TABLE_SV
2479 * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?))
2480 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes))
2481 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no))
2482 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no))
2483 * )
2484 */
2485 int ilogb(T)(const T x) @trusted pure nothrow @nogc
2486 if (isFloatingPoint!T)
2487 {
2488 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2489
2490 import core.bitop : bsr;
2491 alias F = floatTraits!T;
2492
2493 union floatBits
2494 {
2495 T rv;
2496 ushort[T.sizeof/2] vu;
2497 uint[T.sizeof/4] vui;
2498 static if (T.sizeof >= 8)
2499 ulong[T.sizeof/8] vul;
2500 }
2501 floatBits y = void;
2502 y.rv = x;
2503
2504 int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
2505 static if (F.realFormat == RealFormat.ieeeExtended ||
2506 F.realFormat == RealFormat.ieeeExtended53)
2507 {
2508 if (ex)
2509 {
2510 // If exponent is non-zero
2511 if (ex == F.EXPMASK) // infinity or NaN
2512 {
2513 if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2514 return FP_ILOGBNAN;
2515 else // +-infinity
2516 return int.max;
2517 }
2518 else
2519 {
2520 return ex - F.EXPBIAS - 1;
2521 }
2522 }
2523 else if (!y.vul[0])
2524 {
2525 // vf is +-0.0
2526 return FP_ILOGB0;
2527 }
2528 else
2529 {
2530 // subnormal
2531 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
2532 }
2533 }
2534 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2535 {
2536 if (ex) // If exponent is non-zero
2537 {
2538 if (ex == F.EXPMASK)
2539 {
2540 // infinity or NaN
2541 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2542 return FP_ILOGBNAN;
2543 else // +- infinity
2544 return int.max;
2545 }
2546 else
2547 {
2548 return ex - F.EXPBIAS - 1;
2549 }
2550 }
2551 else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2552 {
2553 // vf is +-0.0
2554 return FP_ILOGB0;
2555 }
2556 else
2557 {
2558 // subnormal
2559 const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
2560 const ulong lsb = y.vul[MANTISSA_LSB];
2561 if (msb)
2562 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
2563 else
2564 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
2565 }
2566 }
2567 else static if (F.realFormat == RealFormat.ieeeDouble)
2568 {
2569 if (ex) // If exponent is non-zero
2570 {
2571 if (ex == F.EXPMASK) // infinity or NaN
2572 {
2573 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity
2574 return int.max;
2575 else // NaN
2576 return FP_ILOGBNAN;
2577 }
2578 else
2579 {
2580 return ((ex - F.EXPBIAS) >> 4) - 1;
2581 }
2582 }
2583 else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
2584 {
2585 // vf is +-0.0
2586 return FP_ILOGB0;
2587 }
2588 else
2589 {
2590 // subnormal
2591 enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
2592 return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
2593 }
2594 }
2595 else static if (F.realFormat == RealFormat.ieeeSingle)
2596 {
2597 if (ex) // If exponent is non-zero
2598 {
2599 if (ex == F.EXPMASK) // infinity or NaN
2600 {
2601 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity
2602 return int.max;
2603 else // NaN
2604 return FP_ILOGBNAN;
2605 }
2606 else
2607 {
2608 return ((ex - F.EXPBIAS) >> 7) - 1;
2609 }
2610 }
2611 else if (!(y.vui[0] & 0x7FFF_FFFF))
2612 {
2613 // vf is +-0.0
2614 return FP_ILOGB0;
2615 }
2616 else
2617 {
2618 // subnormal
2619 const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
2620 return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
2621 }
2622 }
2623 else // static if (F.realFormat == RealFormat.ibmExtended)
2624 {
2625 assert(0, "ilogb not implemented");
2626 }
2627 }
2628 /// ditto
2629 int ilogb(T)(const T x) @safe pure nothrow @nogc
2630 if (isIntegral!T && isUnsigned!T)
2631 {
2632 import core.bitop : bsr;
2633 if (x == 0)
2634 return FP_ILOGB0;
2635 else
2636 {
2637 static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
2638 return bsr(x);
2639 }
2640 }
2641 /// ditto
2642 int ilogb(T)(const T x) @safe pure nothrow @nogc
2643 if (isIntegral!T && isSigned!T)
2644 {
2645 import std.traits : Unsigned;
2646 // Note: abs(x) can not be used because the return type is not Unsigned and
2647 // the return value would be wrong for x == int.min
2648 Unsigned!T absx = x >= 0 ? x : -x;
2649 return ilogb(absx);
2650 }
2651
2652 ///
2653 @safe pure unittest
2654 {
2655 assert(ilogb(1) == 0);
2656 assert(ilogb(3) == 1);
2657 assert(ilogb(3.0) == 1);
2658 assert(ilogb(100_000_000) == 26);
2659
2660 assert(ilogb(0) == FP_ILOGB0);
2661 assert(ilogb(0.0) == FP_ILOGB0);
2662 assert(ilogb(double.nan) == FP_ILOGBNAN);
2663 assert(ilogb(double.infinity) == int.max);
2664 }
2665
2666 /**
2667 Special return values of $(LREF ilogb).
2668 */
2669 alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0;
2670 /// ditto
2671 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
2672
2673 ///
2674 @safe pure unittest
2675 {
2676 assert(ilogb(0) == FP_ILOGB0);
2677 assert(ilogb(0.0) == FP_ILOGB0);
2678 assert(ilogb(double.nan) == FP_ILOGBNAN);
2679 }
2680
2681 @safe nothrow @nogc unittest
2682 {
2683 import std.math : floatTraits, RealFormat;
2684 import std.math.operations : nextUp;
2685 import std.meta : AliasSeq;
2686 import std.typecons : Tuple;
2687 static foreach (F; AliasSeq!(float, double, real))
2688 {{
2689 alias T = Tuple!(F, int);
2690 T[13] vals = // x, ilogb(x)
2691 [
2692 T( F.nan , FP_ILOGBNAN ),
2693 T( -F.nan , FP_ILOGBNAN ),
2694 T( F.infinity, int.max ),
2695 T( -F.infinity, int.max ),
2696 T( 0.0 , FP_ILOGB0 ),
2697 T( -0.0 , FP_ILOGB0 ),
2698 T( 2.0 , 1 ),
2699 T( 2.0001 , 1 ),
2700 T( 1.9999 , 0 ),
2701 T( 0.5 , -1 ),
2702 T( 123.123 , 6 ),
2703 T( -123.123 , 6 ),
2704 T( 0.123 , -4 ),
2705 ];
2706
foreach(elem;vals)2707 foreach (elem; vals)
2708 {
2709 assert(ilogb(elem[0]) == elem[1]);
2710 }
2711 }}
2712
2713 // min_normal and subnormals
2714 assert(ilogb(-float.min_normal) == -126);
2715 assert(ilogb(nextUp(-float.min_normal)) == -127);
2716 assert(ilogb(nextUp(-float(0.0))) == -149);
2717 assert(ilogb(-double.min_normal) == -1022);
2718 assert(ilogb(nextUp(-double.min_normal)) == -1023);
2719 assert(ilogb(nextUp(-double(0.0))) == -1074);
2720 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
2721 {
2722 assert(ilogb(-real.min_normal) == -16382);
2723 assert(ilogb(nextUp(-real.min_normal)) == -16383);
2724 assert(ilogb(nextUp(-real(0.0))) == -16445);
2725 }
2726 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2727 {
2728 assert(ilogb(-real.min_normal) == -1022);
2729 assert(ilogb(nextUp(-real.min_normal)) == -1023);
2730 assert(ilogb(nextUp(-real(0.0))) == -1074);
2731 }
2732
2733 // test integer types
2734 assert(ilogb(0) == FP_ILOGB0);
2735 assert(ilogb(int.max) == 30);
2736 assert(ilogb(int.min) == 31);
2737 assert(ilogb(uint.max) == 31);
2738 assert(ilogb(long.max) == 62);
2739 assert(ilogb(long.min) == 63);
2740 assert(ilogb(ulong.max) == 63);
2741 }
2742
2743 /*******************************************
2744 * Compute n * 2$(SUPERSCRIPT exp)
2745 * References: frexp
2746 */
2747
pragma(inline,true)2748 pragma(inline, true)
2749 real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2750 ///ditto
pragma(inline,true)2751 pragma(inline, true)
2752 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2753 ///ditto
pragma(inline,true)2754 pragma(inline, true)
2755 float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2756
2757 ///
2758 @nogc @safe pure nothrow unittest
2759 {
2760 import std.meta : AliasSeq;
2761 static foreach (T; AliasSeq!(float, double, real))
2762 {{
2763 T r;
2764
2765 r = ldexp(3.0L, 3);
2766 assert(r == 24);
2767
2768 r = ldexp(cast(T) 3.0, cast(int) 3);
2769 assert(r == 24);
2770
2771 T n = 3.0;
2772 int exp = 3;
2773 r = ldexp(n, exp);
2774 assert(r == 24);
2775 }}
2776 }
2777
2778 @safe pure nothrow @nogc unittest
2779 {
2780 import std.math : floatTraits, RealFormat;
2781
2782 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
2783 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
2784 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
2785 {
2786 assert(ldexp(1.0L, -16384) == 0x1p-16384L);
2787 assert(ldexp(1.0L, -16382) == 0x1p-16382L);
2788 int x;
2789 real n = frexp(0x1p-16384L, x);
2790 assert(n == 0.5L);
2791 assert(x==-16383);
2792 assert(ldexp(n, x)==0x1p-16384L);
2793 }
2794 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2795 {
2796 assert(ldexp(1.0L, -1024) == 0x1p-1024L);
2797 assert(ldexp(1.0L, -1022) == 0x1p-1022L);
2798 int x;
2799 real n = frexp(0x1p-1024L, x);
2800 assert(n == 0.5L);
2801 assert(x==-1023);
2802 assert(ldexp(n, x)==0x1p-1024L);
2803 }
2804 else static assert(false, "Floating point type real not supported");
2805 }
2806
2807 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718
2808 float parsing depends on platform strtold
2809 @safe pure nothrow @nogc unittest
2810 {
2811 assert(ldexp(1.0, -1024) == 0x1p-1024);
2812 assert(ldexp(1.0, -1022) == 0x1p-1022);
2813 int x;
2814 double n = frexp(0x1p-1024, x);
2815 assert(n == 0.5);
2816 assert(x==-1023);
2817 assert(ldexp(n, x)==0x1p-1024);
2818 }
2819
2820 @safe pure nothrow @nogc unittest
2821 {
2822 assert(ldexp(1.0f, -128) == 0x1p-128f);
2823 assert(ldexp(1.0f, -126) == 0x1p-126f);
2824 int x;
2825 float n = frexp(0x1p-128f, x);
2826 assert(n == 0.5f);
2827 assert(x==-127);
2828 assert(ldexp(n, x)==0x1p-128f);
2829 }
2830 */
2831
2832 @safe @nogc nothrow unittest
2833 {
2834 import std.math.operations : isClose;
2835
2836 static real[3][] vals = // value,exp,ldexp
2837 [
2838 [ 0, 0, 0],
2839 [ 1, 0, 1],
2840 [ -1, 0, -1],
2841 [ 1, 1, 2],
2842 [ 123, 10, 125952],
2843 [ real.max, int.max, real.infinity],
2844 [ real.max, -int.max, 0],
2845 [ real.min_normal, -int.max, 0],
2846 ];
2847 int i;
2848
2849 for (i = 0; i < vals.length; i++)
2850 {
2851 real x = vals[i][0];
2852 int exp = cast(int) vals[i][1];
2853 real z = vals[i][2];
2854 real l = ldexp(x, exp);
2855
2856 assert(isClose(z, l, 1e-6));
2857 }
2858
2859 real function(real, int) pldexp = &ldexp;
2860 assert(pldexp != null);
2861 }
2862
2863 private
2864 {
2865 import std.math : floatTraits, RealFormat;
2866
version(INLINE_YL2X)2867 version (INLINE_YL2X) {} else
2868 {
2869 static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
2870 {
2871 // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
2872 static immutable real[13] logCoeffsP = [
2873 1.313572404063446165910279910527789794488E4L,
2874 7.771154681358524243729929227226708890930E4L,
2875 2.014652742082537582487669938141683759923E5L,
2876 3.007007295140399532324943111654767187848E5L,
2877 2.854829159639697837788887080758954924001E5L,
2878 1.797628303815655343403735250238293741397E5L,
2879 7.594356839258970405033155585486712125861E4L,
2880 2.128857716871515081352991964243375186031E4L,
2881 3.824952356185897735160588078446136783779E3L,
2882 4.114517881637811823002128927449878962058E2L,
2883 2.321125933898420063925789532045674660756E1L,
2884 4.998469661968096229986658302195402690910E-1L,
2885 1.538612243596254322971797716843006400388E-6L
2886 ];
2887 static immutable real[13] logCoeffsQ = [
2888 3.940717212190338497730839731583397586124E4L,
2889 2.626900195321832660448791748036714883242E5L,
2890 7.777690340007566932935753241556479363645E5L,
2891 1.347518538384329112529391120390701166528E6L,
2892 1.514882452993549494932585972882995548426E6L,
2893 1.158019977462989115839826904108208787040E6L,
2894 6.132189329546557743179177159925690841200E5L,
2895 2.248234257620569139969141618556349415120E5L,
2896 5.605842085972455027590989944010492125825E4L,
2897 9.147150349299596453976674231612674085381E3L,
2898 9.104928120962988414618126155557301584078E2L,
2899 4.839208193348159620282142911143429644326E1L,
2900 1.0
2901 ];
2902
2903 // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
2904 // where z = 2(x-1)/(x+1)
2905 static immutable real[6] logCoeffsR = [
2906 1.418134209872192732479751274970992665513E5L,
2907 -8.977257995689735303686582344659576526998E4L,
2908 2.048819892795278657810231591630928516206E4L,
2909 -2.024301798136027039250415126250455056397E3L,
2910 8.057002716646055371965756206836056074715E1L,
2911 -8.828896441624934385266096344596648080902E-1L
2912 ];
2913 static immutable real[7] logCoeffsS = [
2914 1.701761051846631278975701529965589676574E6L,
2915 -1.332535117259762928288745111081235577029E6L,
2916 4.001557694070773974936904547424676279307E5L,
2917 -5.748542087379434595104154610899551484314E4L,
2918 3.998526750980007367835804959888064681098E3L,
2919 -1.186359407982897997337150403816839480438E2L,
2920 1.0
2921 ];
2922 }
2923 else
2924 {
2925 // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
2926 static immutable real[7] logCoeffsP = [
2927 2.0039553499201281259648E1L,
2928 5.7112963590585538103336E1L,
2929 6.0949667980987787057556E1L,
2930 2.9911919328553073277375E1L,
2931 6.5787325942061044846969E0L,
2932 4.9854102823193375972212E-1L,
2933 4.5270000862445199635215E-5L,
2934 ];
2935 static immutable real[7] logCoeffsQ = [
2936 6.0118660497603843919306E1L,
2937 2.1642788614495947685003E2L,
2938 3.0909872225312059774938E2L,
2939 2.2176239823732856465394E2L,
2940 8.3047565967967209469434E1L,
2941 1.5062909083469192043167E1L,
2942 1.0000000000000000000000E0L,
2943 ];
2944
2945 // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
2946 // where z = 2(x-1)/(x+1)
2947 static immutable real[4] logCoeffsR = [
2948 -3.5717684488096787370998E1L,
2949 1.0777257190312272158094E1L,
2950 -7.1990767473014147232598E-1L,
2951 1.9757429581415468984296E-3L,
2952 ];
2953 static immutable real[4] logCoeffsS = [
2954 -4.2861221385716144629696E2L,
2955 1.9361891836232102174846E2L,
2956 -2.6201045551331104417768E1L,
2957 1.0000000000000000000000E0L,
2958 ];
2959 }
2960 }
2961 }
2962
2963 /**************************************
2964 * Calculate the natural logarithm of x.
2965 *
2966 * $(TABLE_SV
2967 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?))
2968 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
2969 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
2970 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
2971 * )
2972 */
log(real x)2973 real log(real x) @safe pure nothrow @nogc
2974 {
2975 import std.math.constants : LN2, LOG2, SQRT1_2;
2976 import std.math.traits : isInfinity, isNaN, signbit;
2977 import std.math.algebraic : poly;
2978
2979 version (INLINE_YL2X)
2980 return core.math.yl2x(x, LN2);
2981 else
2982 {
2983 // C1 + C2 = LN2.
2984 enum real C1 = 6.93145751953125E-1L;
2985 enum real C2 = 1.428606820309417232121458176568075500134E-6L;
2986
2987 // Special cases.
2988 if (isNaN(x))
2989 return x;
2990 if (isInfinity(x) && !signbit(x))
2991 return x;
2992 if (x == 0.0)
2993 return -real.infinity;
2994 if (x < 0.0)
2995 return real.nan;
2996
2997 // Separate mantissa from exponent.
2998 // Note, frexp is used so that denormal numbers will be handled properly.
2999 real y, z;
3000 int exp;
3001
3002 x = frexp(x, exp);
3003
3004 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3005 // where z = 2(x - 1)/(x + 1)
3006 if ((exp > 2) || (exp < -2))
3007 {
3008 if (x < SQRT1_2)
3009 { // 2(2x - 1)/(2x + 1)
3010 exp -= 1;
3011 z = x - 0.5;
3012 y = 0.5 * z + 0.5;
3013 }
3014 else
3015 { // 2(x - 1)/(x + 1)
3016 z = x - 0.5;
3017 z -= 0.5;
3018 y = 0.5 * x + 0.5;
3019 }
3020 x = z / y;
3021 z = x * x;
3022 z = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
3023 z += exp * C2;
3024 z += x;
3025 z += exp * C1;
3026
3027 return z;
3028 }
3029
3030 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3031 if (x < SQRT1_2)
3032 {
3033 exp -= 1;
3034 x = 2.0 * x - 1.0;
3035 }
3036 else
3037 {
3038 x = x - 1.0;
3039 }
3040 z = x * x;
3041 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
3042 y += exp * C2;
3043 z = y - 0.5 * z;
3044
3045 // Note, the sum of above terms does not exceed x/4,
3046 // so it contributes at most about 1/4 lsb to the error.
3047 z += x;
3048 z += exp * C1;
3049
3050 return z;
3051 }
3052 }
3053
3054 ///
3055 @safe pure nothrow @nogc unittest
3056 {
3057 import std.math.operations : feqrel;
3058 import std.math.constants : E;
3059
3060 assert(feqrel(log(E), 1) >= real.mant_dig - 1);
3061 }
3062
3063 /**************************************
3064 * Calculate the base-10 logarithm of x.
3065 *
3066 * $(TABLE_SV
3067 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?))
3068 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3069 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
3070 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3071 * )
3072 */
log10(real x)3073 real log10(real x) @safe pure nothrow @nogc
3074 {
3075 import std.math.constants : LOG2, LN2, SQRT1_2;
3076 import std.math.algebraic : poly;
3077 import std.math.traits : isNaN, isInfinity, signbit;
3078
3079 version (INLINE_YL2X)
3080 return core.math.yl2x(x, LOG2);
3081 else
3082 {
3083 // log10(2) split into two parts.
3084 enum real L102A = 0.3125L;
3085 enum real L102B = -1.14700043360188047862611052755069732318101185E-2L;
3086
3087 // log10(e) split into two parts.
3088 enum real L10EA = 0.5L;
3089 enum real L10EB = -6.570551809674817234887108108339491770560299E-2L;
3090
3091 // Special cases are the same as for log.
3092 if (isNaN(x))
3093 return x;
3094 if (isInfinity(x) && !signbit(x))
3095 return x;
3096 if (x == 0.0)
3097 return -real.infinity;
3098 if (x < 0.0)
3099 return real.nan;
3100
3101 // Separate mantissa from exponent.
3102 // Note, frexp is used so that denormal numbers will be handled properly.
3103 real y, z;
3104 int exp;
3105
3106 x = frexp(x, exp);
3107
3108 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3109 // where z = 2(x - 1)/(x + 1)
3110 if ((exp > 2) || (exp < -2))
3111 {
3112 if (x < SQRT1_2)
3113 { // 2(2x - 1)/(2x + 1)
3114 exp -= 1;
3115 z = x - 0.5;
3116 y = 0.5 * z + 0.5;
3117 }
3118 else
3119 { // 2(x - 1)/(x + 1)
3120 z = x - 0.5;
3121 z -= 0.5;
3122 y = 0.5 * x + 0.5;
3123 }
3124 x = z / y;
3125 z = x * x;
3126 y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
3127 goto Ldone;
3128 }
3129
3130 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3131 if (x < SQRT1_2)
3132 {
3133 exp -= 1;
3134 x = 2.0 * x - 1.0;
3135 }
3136 else
3137 x = x - 1.0;
3138
3139 z = x * x;
3140 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
3141 y = y - 0.5 * z;
3142
3143 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3144 // This sequence of operations is critical and it may be horribly
3145 // defeated by some compiler optimizers.
3146 Ldone:
3147 z = y * L10EB;
3148 z += x * L10EB;
3149 z += exp * L102B;
3150 z += y * L10EA;
3151 z += x * L10EA;
3152 z += exp * L102A;
3153
3154 return z;
3155 }
3156 }
3157
3158 ///
3159 @safe pure nothrow @nogc unittest
3160 {
3161 import std.math.algebraic : fabs;
3162
3163 assert(fabs(log10(1000) - 3) < .000001);
3164 }
3165
3166 /**
3167 * Calculates the natural logarithm of 1 + x.
3168 *
3169 * For very small x, log1p(x) will be more accurate than
3170 * log(1 + x).
3171 *
3172 * $(TABLE_SV
3173 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?))
3174 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no))
3175 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3176 * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes))
3177 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3178 * )
3179 */
log1p(real x)3180 real log1p(real x) @safe pure nothrow @nogc
3181 {
3182 import std.math.traits : isNaN, isInfinity, signbit;
3183 import std.math.constants : LN2;
3184
3185 version (INLINE_YL2X)
3186 {
3187 // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
3188 // ie if -0.29 <= x <= 0.414
3189 return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
3190 }
3191 else
3192 {
3193 // Special cases.
3194 if (isNaN(x) || x == 0.0)
3195 return x;
3196 if (isInfinity(x) && !signbit(x))
3197 return x;
3198 if (x == -1.0)
3199 return -real.infinity;
3200 if (x < -1.0)
3201 return real.nan;
3202
3203 return log(x + 1.0);
3204 }
3205 }
3206
3207 ///
3208 @safe pure unittest
3209 {
3210 import std.math.traits : isIdentical, isNaN;
3211 import std.math.operations : feqrel;
3212
3213 assert(isIdentical(log1p(0.0), 0.0));
3214 assert(log1p(1.0).feqrel(0.69314) > 16);
3215
3216 assert(log1p(-1.0) == -real.infinity);
3217 assert(isNaN(log1p(-2.0)));
3218 assert(log1p(real.nan) is real.nan);
3219 assert(log1p(-real.nan) is -real.nan);
3220 assert(log1p(real.infinity) == real.infinity);
3221 }
3222
3223 /***************************************
3224 * Calculates the base-2 logarithm of x:
3225 * $(SUB log, 2)x
3226 *
3227 * $(TABLE_SV
3228 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?))
3229 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) )
3230 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) )
3231 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) )
3232 * )
3233 */
log2(real x)3234 real log2(real x) @safe pure nothrow @nogc
3235 {
3236 import std.math.traits : isNaN, isInfinity, signbit;
3237 import std.math.constants : SQRT1_2, LOG2E;
3238 import std.math.algebraic : poly;
3239
3240 version (INLINE_YL2X)
3241 return core.math.yl2x(x, 1.0L);
3242 else
3243 {
3244 // Special cases are the same as for log.
3245 if (isNaN(x))
3246 return x;
3247 if (isInfinity(x) && !signbit(x))
3248 return x;
3249 if (x == 0.0)
3250 return -real.infinity;
3251 if (x < 0.0)
3252 return real.nan;
3253
3254 // Separate mantissa from exponent.
3255 // Note, frexp is used so that denormal numbers will be handled properly.
3256 real y, z;
3257 int exp;
3258
3259 x = frexp(x, exp);
3260
3261 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3262 // where z = 2(x - 1)/(x + 1)
3263 if ((exp > 2) || (exp < -2))
3264 {
3265 if (x < SQRT1_2)
3266 { // 2(2x - 1)/(2x + 1)
3267 exp -= 1;
3268 z = x - 0.5;
3269 y = 0.5 * z + 0.5;
3270 }
3271 else
3272 { // 2(x - 1)/(x + 1)
3273 z = x - 0.5;
3274 z -= 0.5;
3275 y = 0.5 * x + 0.5;
3276 }
3277 x = z / y;
3278 z = x * x;
3279 y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
3280 goto Ldone;
3281 }
3282
3283 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3284 if (x < SQRT1_2)
3285 {
3286 exp -= 1;
3287 x = 2.0 * x - 1.0;
3288 }
3289 else
3290 x = x - 1.0;
3291
3292 z = x * x;
3293 y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
3294 y = y - 0.5 * z;
3295
3296 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3297 // This sequence of operations is critical and it may be horribly
3298 // defeated by some compiler optimizers.
3299 Ldone:
3300 z = y * (LOG2E - 1.0);
3301 z += x * (LOG2E - 1.0);
3302 z += y;
3303 z += x;
3304 z += exp;
3305
3306 return z;
3307 }
3308 }
3309
3310 ///
3311 @safe unittest
3312 {
3313 import std.math.operations : isClose;
3314
3315 assert(isClose(log2(1024.0L), 10));
3316 }
3317
3318 @safe @nogc nothrow unittest
3319 {
3320 import std.math.operations : isClose;
3321
3322 // check if values are equal to 19 decimal digits of precision
3323 assert(isClose(log2(1024.0L), 10, 1e-18));
3324 }
3325
3326 /*****************************************
3327 * Extracts the exponent of x as a signed integral value.
3328 *
3329 * If x is subnormal, it is treated as if it were normalized.
3330 * For a positive, finite x:
3331 *
3332 * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
3333 *
3334 * $(TABLE_SV
3335 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) )
3336 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
3337 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) )
3338 * )
3339 */
logb(real x)3340 real logb(real x) @trusted nothrow @nogc
3341 {
3342 version (InlineAsm_X87_MSVC)
3343 {
3344 version (X86_64)
3345 {
3346 asm pure nothrow @nogc
3347 {
3348 naked ;
3349 fld real ptr [RCX] ;
3350 fxtract ;
3351 fstp ST(0) ;
3352 ret ;
3353 }
3354 }
3355 else
3356 {
3357 asm pure nothrow @nogc
3358 {
3359 fld x ;
3360 fxtract ;
3361 fstp ST(0) ;
3362 }
3363 }
3364 }
3365 else
3366 return core.stdc.math.logbl(x);
3367 }
3368
3369 ///
3370 @safe @nogc nothrow unittest
3371 {
3372 assert(logb(1.0) == 0);
3373 assert(logb(100.0) == 6);
3374
3375 assert(logb(0.0) == -real.infinity);
3376 assert(logb(real.infinity) == real.infinity);
3377 assert(logb(-real.infinity) == real.infinity);
3378 }
3379
3380 /*************************************
3381 * Efficiently calculates x * 2$(SUPERSCRIPT n).
3382 *
3383 * scalbn handles underflow and overflow in
3384 * the same fashion as the basic arithmetic operators.
3385 *
3386 * $(TABLE_SV
3387 * $(TR $(TH x) $(TH scalb(x)))
3388 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) )
3389 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
3390 * )
3391 */
pragma(inline,true)3392 pragma(inline, true)
3393 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
3394
3395 /// ditto
pragma(inline,true)3396 pragma(inline, true)
3397 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
3398
3399 /// ditto
pragma(inline,true)3400 pragma(inline, true)
3401 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
3402
3403 ///
3404 @safe pure nothrow @nogc unittest
3405 {
3406 assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
3407 assert(scalbn(-real.infinity, 5) == -real.infinity);
3408 assert(scalbn(2.0,10) == 2048.0);
3409 assert(scalbn(2048.0f,-10) == 2.0f);
3410 }
3411
pragma(inline,true)3412 pragma(inline, true)
3413 private F _scalbn(F)(F x, int n)
3414 {
3415 import std.math.traits : isInfinity;
3416
3417 if (__ctfe)
3418 {
3419 // Handle special cases.
3420 if (x == F(0.0) || isInfinity(x))
3421 return x;
3422 }
3423 return core.math.ldexp(x, n);
3424 }
3425
3426 @safe pure nothrow @nogc unittest
3427 {
3428 // CTFE-able test
3429 static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
3430 static assert(scalbn(-real.infinity, 5) == -real.infinity);
3431 // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
3432 enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
3433 enum int n = resultExponent - initialExponent;
3434 enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent);
3435 enum staticResult = scalbn(x, n);
3436 static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent));
3437 assert(scalbn(x, n) == staticResult);
3438 }
3439
3440