?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
d: Merge upstream dmd, druntime 4ca4140e58, phobos 454dff14d.
[official-gcc.git] / libphobos / src / std / math / exponential.d
blobfd1ff244b626a1a2f95a43740f9ccaacd97e858a
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 = &plusmn;
27 INFIN = &infin;
28 PLUSMNINF = &plusmn;&infin;
29 LT = &lt;
30 GT = &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
40 version (DigitalMars)
41 {
42 version (OSX) { } // macOS 13 (M1) has issues emulating instruction
43 else version = INLINE_YL2X; // x87 has opcodes for these
44 }
45
46 version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
47 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
48
49 version (InlineAsm_X86_Any) version = InlineAsm_X87;
50 version (InlineAsm_X87)
51 {
52 static assert(real.mant_dig == 64);
53 version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
54 }
55
56 version (D_HardFloat)
57 {
58 // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport
59 version (IeeeFlagsSupport) version = FloatingPointControlSupport;
60 }
61
62 /**
63 * Compute the value of x $(SUPERSCRIPT n), where n is an integer
64 */
65 Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow
66 if (isFloatingPoint!(F) && isIntegral!(G))
67 {
68 import std.traits : Unsigned;
69
70 real p = 1.0, v = void;
71 Unsigned!(Unqual!G) m = n;
72
73 if (n < 0)
74 {
75 if (n == -1) return 1 / x;
76
77 m = cast(typeof(m))(0 - n);
78 v = p / x;
79 }
80 else
81 {
82 switch (n)
83 {
84 case 0:
85 return 1.0;
86 case 1:
87 return x;
88 case 2:
89 return x * x;
90 default:
91 }
92
93 v = x;
94 }
95
96 while (1)
97 {
98 if (m & 1)
99 p *= v;
100 m >>= 1;
101 if (!m)
102 break;
103 v *= v;
104 }
105 return p;
106 }
107
108 ///
109 @safe pure nothrow @nogc unittest
110 {
111 import std.math.operations : feqrel;
112
113 assert(pow(2.0, 5) == 32.0);
114 assert(pow(1.5, 9).feqrel(38.4433) > 16);
115 assert(pow(real.nan, 2) is real.nan);
116 assert(pow(real.infinity, 2) == real.infinity);
117 }
118
119 @safe pure nothrow @nogc unittest
120 {
121 import std.math.operations : isClose, feqrel;
122
123 // Make sure it instantiates and works properly on immutable values and
124 // with various integer and float types.
125 immutable real x = 46;
126 immutable float xf = x;
127 immutable double xd = x;
128 immutable uint one = 1;
129 immutable ushort two = 2;
130 immutable ubyte three = 3;
131 immutable ulong eight = 8;
132
133 immutable int neg1 = -1;
134 immutable short neg2 = -2;
135 immutable byte neg3 = -3;
136 immutable long neg8 = -8;
137
138
139 assert(pow(x,0) == 1.0);
140 assert(pow(xd,one) == x);
141 assert(pow(xf,two) == x * x);
142 assert(pow(x,three) == x * x * x);
143 assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
144
145 assert(pow(x, neg1) == 1 / x);
146
147 assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15));
148 assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15));
149
150 assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1);
151 }
152
153 @safe @nogc nothrow unittest
154 {
155 import std.math.operations : isClose;
156
157 assert(isClose(pow(2.0L, 10L), 1024, 1e-18));
158 }
159
160 // https://issues.dlang.org/show_bug.cgi?id=21601
161 @safe @nogc nothrow pure unittest
162 {
163 // When reals are large enough the results of pow(b, e) can be
164 // calculated correctly, if b is of type float or double and e is
165 // not too large.
166 static if (real.mant_dig >= 64)
167 {
168 // expected result: 3.790e-42
169 assert(pow(-513645318757045764096.0f, -2) > 0.0);
170
171 // expected result: 3.763915357831797e-309
172 assert(pow(-1.6299717435255677e+154, -2) > 0.0);
173 }
174 }
175
176 @safe @nogc nothrow unittest
177 {
178 import std.math.operations : isClose;
179 import std.math.traits : isInfinity;
180
181 static float f1 = 19100.0f;
182 static float f2 = 0.000012f;
183
184 assert(isClose(pow(f1,9), 3.3829868e+38f));
185 assert(isInfinity(pow(f1,10)));
186 assert(pow(f2,9) > 0.0f);
187 assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
188
189 static double d1 = 21800.0;
190 static double d2 = 0.000012;
191
192 assert(isClose(pow(d1,71), 1.0725339442974e+308));
193 assert(isInfinity(pow(d1,72)));
194 assert(pow(d2,65) > 0.0f);
195 assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
196
197 static if (real.mant_dig == 64) // x87
198 {
199 static real r1 = 21950.0L;
200 static real r2 = 0.000011L;
201
202 assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
203 assert(isInfinity(pow(r1,1137)));
204 assert(pow(r2,998) > 0.0L);
205 assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
206 }
207 }
208
209 @safe @nogc nothrow pure unittest
210 {
211 import std.math.operations : isClose;
212
213 enum f1 = 19100.0f;
214 enum f2 = 0.000012f;
215
216 static assert(isClose(pow(f1,9), 3.3829868e+38f));
217 static assert(pow(f1,10) > float.max);
218 static assert(pow(f2,9) > 0.0f);
219 static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal));
220
221 enum d1 = 21800.0;
222 enum d2 = 0.000012;
223
224 static assert(isClose(pow(d1,71), 1.0725339442974e+308));
225 static assert(pow(d1,72) > double.max);
226 static assert(pow(d2,65) > 0.0f);
227 static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal));
228
229 static if (real.mant_dig == 64) // x87
230 {
231 enum r1 = 21950.0L;
232 enum r2 = 0.000011L;
233
234 static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L));
235 static assert(pow(r1,1137) > real.max);
236 static assert(pow(r2,998) > 0.0L);
237 static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal));
238 }
239 }
240
241 /**
242 * Compute the power of two integral numbers.
243 *
244 * Params:
245 * x = base
246 * n = exponent
247 *
248 * Returns:
249 * x raised to the power of n. If n is negative the result is 1 / pow(x, -n),
250 * which is calculated as integer division with remainder. This may result in
251 * a division by zero error.
252 *
253 * If both x and n are 0, the result is 1.
254 *
255 * Throws:
256 * If x is 0 and n is negative, the result is the same as the result of a
257 * division by zero.
258 */
259 typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow
260 if (isIntegral!(F) && isIntegral!(G))
261 {
262 import std.traits : isSigned;
263
264 typeof(return) p, v = void;
265 Unqual!G m = n;
266
267 static if (isSigned!(F))
268 {
269 if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1);
270 }
271 static if (isSigned!(G))
272 {
273 if (x == 0 && m <= -1) return x / 0;
274 }
275 if (x == 1) return 1;
276 static if (isSigned!(G))
277 {
278 if (m < 0) return 0;
279 }
280
281 switch (m)
282 {
283 case 0:
284 p = 1;
285 break;
286
287 case 1:
288 p = x;
289 break;
290
291 case 2:
292 p = x * x;
293 break;
294
295 default:
296 v = x;
297 p = 1;
298 while (1)
299 {
300 if (m & 1)
301 p *= v;
302 m >>= 1;
303 if (!m)
304 break;
305 v *= v;
306 }
307 break;
308 }
309 return p;
310 }
311
312 ///
313 @safe pure nothrow @nogc unittest
314 {
315 assert(pow(2, 3) == 8);
316 assert(pow(3, 2) == 9);
317
318 assert(pow(2, 10) == 1_024);
319 assert(pow(2, 20) == 1_048_576);
320 assert(pow(2, 30) == 1_073_741_824);
321
322 assert(pow(0, 0) == 1);
323
324 assert(pow(1, -5) == 1);
325 assert(pow(1, -6) == 1);
326 assert(pow(-1, -5) == -1);
327 assert(pow(-1, -6) == 1);
328
329 assert(pow(-2, 5) == -32);
330 assert(pow(-2, -5) == 0);
331 assert(pow(cast(double) -2, -5) == -0.03125);
332 }
333
334 @safe pure nothrow @nogc unittest
335 {
336 immutable int one = 1;
337 immutable byte two = 2;
338 immutable ubyte three = 3;
339 immutable short four = 4;
340 immutable long ten = 10;
341
342 assert(pow(two, three) == 8);
343 assert(pow(two, ten) == 1024);
344 assert(pow(one, ten) == 1);
345 assert(pow(ten, four) == 10_000);
346 assert(pow(four, 10) == 1_048_576);
347 assert(pow(three, four) == 81);
348 }
349
350 // https://issues.dlang.org/show_bug.cgi?id=7006
351 @safe pure nothrow @nogc unittest
352 {
353 assert(pow(5, -1) == 0);
354 assert(pow(-5, -1) == 0);
355 assert(pow(5, -2) == 0);
356 assert(pow(-5, -2) == 0);
357 assert(pow(-1, int.min) == 1);
358 assert(pow(-2, int.min) == 0);
359
360 assert(pow(4294967290UL,2) == 18446744022169944100UL);
361 assert(pow(0,uint.max) == 0);
362 }
363
364 /**Computes integer to floating point powers.*/
365 real pow(I, F)(I x, F y) @nogc @trusted pure nothrow
366 if (isIntegral!I && isFloatingPoint!F)
367 {
368 return pow(cast(real) x, cast(Unqual!F) y);
369 }
370
371 ///
372 @safe pure nothrow @nogc unittest
373 {
374 assert(pow(2, 5.0) == 32.0);
375 assert(pow(7, 3.0) == 343.0);
376 assert(pow(2, real.nan) is real.nan);
377 assert(pow(2, real.infinity) == real.infinity);
378 }
379
380 /**
381 * Calculates x$(SUPERSCRIPT y).
382 *
383 * $(TABLE_SV
384 * $(TR $(TH x) $(TH y) $(TH pow(x, y))
385 * $(TH div 0) $(TH invalid?))
386 * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0)
387 * $(TD no) $(TD no) )
388 * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN))
389 * $(TD no) $(TD no) )
390 * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0)
391 * $(TD no) $(TD no) )
392 * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0)
393 * $(TD no) $(TD no) )
394 * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN))
395 * $(TD no) $(TD no) )
396 * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN))
397 * $(TD no) $(TD no) )
398 * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0)
399 * $(TD no) $(TD no) )
400 * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN))
401 * $(TD no) $(TD no) )
402 * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
403 * $(TD no) $(TD no))
404 * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0)
405 * $(TD no) $(TD no) )
406 * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
407 * $(TD no) $(TD no) )
408 * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN))
409 * $(TD no) $(TD yes) )
410 * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN))
411 * $(TD no) $(TD yes))
412 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF))
413 * $(TD yes) $(TD no) )
414 * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
415 * $(TD yes) $(TD no))
416 * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0)
417 * $(TD no) $(TD no) )
418 * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
419 * $(TD no) $(TD no) )
420 * )
421 */
422 Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow
423 if (isFloatingPoint!(F) && isFloatingPoint!(G))
424 {
425 import core.math : fabs, sqrt;
426 import std.math.traits : isInfinity, isNaN, signbit;
427
428 alias Float = typeof(return);
429
430 static real impl(real x, real y) @nogc pure nothrow
431 {
432 // Special cases.
433 if (isNaN(y))
434 return y;
435 if (isNaN(x) && y != 0.0)
436 return x;
437
438 // Even if x is NaN.
439 if (y == 0.0)
440 return 1.0;
441 if (y == 1.0)
442 return x;
443
444 if (isInfinity(y))
445 {
446 if (isInfinity(x))
447 {
448 if (!signbit(y) && !signbit(x))
449 return F.infinity;
450 else
451 return F.nan;
452 }
453 else if (fabs(x) > 1)
454 {
455 if (signbit(y))
456 return +0.0;
457 else
458 return F.infinity;
459 }
460 else if (fabs(x) == 1)
461 {
462 return F.nan;
463 }
464 else // < 1
465 {
466 if (signbit(y))
467 return F.infinity;
468 else
469 return +0.0;
470 }
471 }
472 if (isInfinity(x))
473 {
474 if (signbit(x))
475 {
476 long i = cast(long) y;
477 if (y > 0.0)
478 {
479 if (i == y && i & 1)
480 return -F.infinity;
481 else if (i == y)
482 return F.infinity;
483 else
484 return -F.nan;
485 }
486 else if (y < 0.0)
487 {
488 if (i == y && i & 1)
489 return -0.0;
490 else if (i == y)
491 return +0.0;
492 else
493 return F.nan;
494 }
495 }
496 else
497 {
498 if (y > 0.0)
499 return F.infinity;
500 else if (y < 0.0)
501 return +0.0;
502 }
503 }
504
505 if (x == 0.0)
506 {
507 if (signbit(x))
508 {
509 long i = cast(long) y;
510 if (y > 0.0)
511 {
512 if (i == y && i & 1)
513 return -0.0;
514 else
515 return +0.0;
516 }
517 else if (y < 0.0)
518 {
519 if (i == y && i & 1)
520 return -F.infinity;
521 else
522 return F.infinity;
523 }
524 }
525 else
526 {
527 if (y > 0.0)
528 return +0.0;
529 else if (y < 0.0)
530 return F.infinity;
531 }
532 }
533 if (x == 1.0)
534 return 1.0;
535
536 if (y >= F.max)
537 {
538 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
539 return 0.0;
540 if (x > 1.0 || x < -1.0)
541 return F.infinity;
542 }
543 if (y <= -F.max)
544 {
545 if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0))
546 return F.infinity;
547 if (x > 1.0 || x < -1.0)
548 return 0.0;
549 }
550
551 if (x >= F.max)
552 {
553 if (y > 0.0)
554 return F.infinity;
555 else
556 return 0.0;
557 }
558 if (x <= -F.max)
559 {
560 long i = cast(long) y;
561 if (y > 0.0)
562 {
563 if (i == y && i & 1)
564 return -F.infinity;
565 else
566 return F.infinity;
567 }
568 else if (y < 0.0)
569 {
570 if (i == y && i & 1)
571 return -0.0;
572 else
573 return +0.0;
574 }
575 }
576
577 // Integer power of x.
578 long iy = cast(long) y;
579 if (iy == y && fabs(y) < 32_768.0)
580 return pow(x, iy);
581
582 real sign = 1.0;
583 if (x < 0)
584 {
585 // Result is real only if y is an integer
586 // Check for a non-zero fractional part
587 enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
588 static if (maxOdd > ulong.max)
589 {
590 // Generic method, for any FP type
591 import std.math.rounding : floor;
592 if (floor(y) != y)
593 return sqrt(x); // Complex result -- create a NaN
594
595 const hy = 0.5 * y;
596 if (floor(hy) != hy)
597 sign = -1.0;
598 }
599 else
600 {
601 // Much faster, if ulong has enough precision
602 const absY = fabs(y);
603 if (absY <= maxOdd)
604 {
605 const uy = cast(ulong) absY;
606 if (uy != absY)
607 return sqrt(x); // Complex result -- create a NaN
608
609 if (uy & 1)
610 sign = -1.0;
611 }
612 }
613 x = -x;
614 }
615 version (INLINE_YL2X)
616 {
617 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
618 // TODO: This is not accurate in practice. A fast and accurate
619 // (though complicated) method is described in:
620 // "An efficient rounding boundary test for pow(x, y)
621 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
622 return sign * exp2( core.math.yl2x(x, y) );
623 }
624 else
625 {
626 // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
627 // TODO: This is not accurate in practice. A fast and accurate
628 // (though complicated) method is described in:
629 // "An efficient rounding boundary test for pow(x, y)
630 // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
631 Float w = exp2(y * log2(x));
632 return sign * w;
633 }
634 }
635 return impl(x, y);
636 }
637
638 ///
639 @safe pure nothrow @nogc unittest
640 {
641 import std.math.operations : isClose;
642
643 assert(isClose(pow(2.0, 3.0), 8.0));
644 assert(isClose(pow(1.5, 10.0), 57.6650390625));
645
646 // square root of 9
647 assert(isClose(pow(9.0, 0.5), 3.0));
648 // 10th root of 1024
649 assert(isClose(pow(1024.0, 0.1), 2.0));
650
651 assert(isClose(pow(-4.0, 3.0), -64.0));
652
653 // reciprocal of 4 ^^ 2
654 assert(isClose(pow(4.0, -2.0), 0.0625));
655 // reciprocal of (-2) ^^ 3
656 assert(isClose(pow(-2.0, -3.0), -0.125));
657
658 assert(isClose(pow(-2.5, 3.0), -15.625));
659 // reciprocal of 2.5 ^^ 3
660 assert(isClose(pow(2.5, -3.0), 0.064));
661 // reciprocal of (-2.5) ^^ 3
662 assert(isClose(pow(-2.5, -3.0), -0.064));
663
664 // reciprocal of square root of 4
665 assert(isClose(pow(4.0, -0.5), 0.5));
666
667 // per definition
668 assert(isClose(pow(0.0, 0.0), 1.0));
669 }
670
671 ///
672 @safe pure nothrow @nogc unittest
673 {
674 import std.math.operations : isClose;
675
676 // the result is a complex number
677 // which cannot be represented as floating point number
678 import std.math.traits : isNaN;
679 assert(isNaN(pow(-2.5, -1.5)));
680
681 // use the ^^-operator of std.complex instead
682 import std.complex : complex;
683 auto c1 = complex(-2.5, 0.0);
684 auto c2 = complex(-1.5, 0.0);
685 auto result = c1 ^^ c2;
686 // exact result apparently depends on `real` precision => increased tolerance
687 assert(isClose(result.re, -4.64705438e-17, 2e-4));
688 assert(isClose(result.im, 2.52982e-1, 2e-4));
689 }
690
691 @safe pure nothrow @nogc unittest
692 {
693 import std.math.traits : isNaN;
694
695 assert(pow(1.5, real.infinity) == real.infinity);
696 assert(pow(0.5, real.infinity) == 0.0);
697 assert(pow(1.5, -real.infinity) == 0.0);
698 assert(pow(0.5, -real.infinity) == real.infinity);
699 assert(pow(real.infinity, 1.0) == real.infinity);
700 assert(pow(real.infinity, -1.0) == 0.0);
701 assert(pow(real.infinity, real.infinity) == real.infinity);
702 assert(pow(-real.infinity, 1.0) == -real.infinity);
703 assert(pow(-real.infinity, 2.0) == real.infinity);
704 assert(pow(-real.infinity, -1.0) == -0.0);
705 assert(pow(-real.infinity, -2.0) == 0.0);
706 assert(isNaN(pow(1.0, real.infinity)));
707 assert(pow(0.0, -1.0) == real.infinity);
708 assert(pow(real.nan, 0.0) == 1.0);
709 assert(isNaN(pow(real.nan, 3.0)));
710 assert(isNaN(pow(3.0, real.nan)));
711 }
712
713 @safe @nogc nothrow unittest
714 {
715 import std.math.operations : isClose;
716
717 assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18));
718 }
719
720 @safe pure nothrow @nogc unittest
721 {
722 import std.math.operations : isClose;
723 import std.math.traits : isIdentical, isNaN;
724 import std.math.constants : PI;
725
726 // Test all the special values. These unittests can be run on Windows
727 // by temporarily changing the version (linux) to version (all).
728 immutable float zero = 0;
729 immutable real one = 1;
730 immutable double two = 2;
731 immutable float three = 3;
732 immutable float fnan = float.nan;
733 immutable double dnan = double.nan;
734 immutable real rnan = real.nan;
735 immutable dinf = double.infinity;
736 immutable rninf = -real.infinity;
737
738 assert(pow(fnan, zero) == 1);
739 assert(pow(dnan, zero) == 1);
740 assert(pow(rnan, zero) == 1);
741
742 assert(pow(two, dinf) == double.infinity);
743 assert(isIdentical(pow(0.2f, dinf), +0.0));
744 assert(pow(0.99999999L, rninf) == real.infinity);
745 assert(isIdentical(pow(1.000000001, rninf), +0.0));
746 assert(pow(dinf, 0.001) == dinf);
747 assert(isIdentical(pow(dinf, -0.001), +0.0));
748 assert(pow(rninf, 3.0L) == rninf);
749 assert(pow(rninf, 2.0L) == real.infinity);
750 assert(isIdentical(pow(rninf, -3.0), -0.0));
751 assert(isIdentical(pow(rninf, -2.0), +0.0));
752
753 // @@@BUG@@@ somewhere
754 version (OSX) {} else assert(isNaN(pow(one, dinf)));
755 version (OSX) {} else assert(isNaN(pow(-one, dinf)));
756 assert(isNaN(pow(-0.2, PI)));
757 // boundary cases. Note that epsilon == 2^^-n for some n,
758 // so 1/epsilon == 2^^n is always even.
759 assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
760 assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
761 assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
762 assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
763
764 assert(pow(0.0, -3.0) == double.infinity);
765 assert(pow(-0.0, -3.0) == -double.infinity);
766 assert(pow(0.0, -PI) == double.infinity);
767 assert(pow(-0.0, -PI) == double.infinity);
768 assert(isIdentical(pow(0.0, 5.0), 0.0));
769 assert(isIdentical(pow(-0.0, 5.0), -0.0));
770 assert(isIdentical(pow(0.0, 6.0), 0.0));
771 assert(isIdentical(pow(-0.0, 6.0), 0.0));
772
773 // https://issues.dlang.org/show_bug.cgi?id=14786 fixed
774 immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L;
775 assert(pow(-1.0L, maxOdd) == -1.0L);
776 assert(pow(-1.0L, -maxOdd) == -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 assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L);
781
782 // Now, actual numbers.
783 assert(isClose(pow(two, three), 8.0));
784 assert(isClose(pow(two, -2.5), 0.1767766953));
785
786 // Test integer to float power.
787 immutable uint twoI = 2;
788 assert(isClose(pow(twoI, three), 8.0));
789 }
790
791 // https://issues.dlang.org/show_bug.cgi?id=20508
792 @safe pure nothrow @nogc unittest
793 {
794 import std.math.traits : isNaN;
795
796 assert(isNaN(pow(-double.infinity, 0.5)));
797
798 assert(isNaN(pow(-real.infinity, real.infinity)));
799 assert(isNaN(pow(-real.infinity, -real.infinity)));
800 assert(isNaN(pow(-real.infinity, 1.234)));
801 assert(isNaN(pow(-real.infinity, -0.751)));
802 assert(pow(-real.infinity, 0.0) == 1.0);
803 }
804
805 /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`.
806 *
807 * Params:
808 * x = base
809 * n = exponent
810 * m = modulus
811 *
812 * Returns:
813 * `x` to the power `n`, modulo `m`.
814 * The return type is the largest of `x`'s and `m`'s type.
815 *
816 * The function requires that all values have unsigned types.
817 */
818 Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m)
819 if (isUnsigned!F && isUnsigned!G && isUnsigned!H)
820 {
821 import std.meta : AliasSeq;
822
823 alias T = Unqual!(Largest!(F, H));
824 static if (T.sizeof <= 4)
825 {
826 alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof];
827 }
828
829 static T mulmod(T a, T b, T c)
830 {
831 static if (T.sizeof == 8)
832 {
833 static T addmod(T a, T b, T c)
834 {
835 b = c - b;
836 if (a >= b)
837 return a - b;
838 else
839 return c - b + a;
840 }
841
842 T result = 0, tmp;
843
844 b %= c;
845 while (a > 0)
846 {
847 if (a & 1)
848 result = addmod(result, b, c);
849
850 a >>= 1;
851 b = addmod(b, b, c);
852 }
853
854 return result;
855 }
856 else
857 {
858 DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b);
859 return result % c;
860 }
861 }
862
863 T base = x, result = 1, modulus = m;
864 Unqual!G exponent = n;
865
866 while (exponent > 0)
867 {
868 if (exponent & 1)
869 result = mulmod(result, base, modulus);
870
871 base = mulmod(base, base, modulus);
872 exponent >>= 1;
873 }
874
875 return result;
876 }
877
878 ///
879 @safe pure nothrow @nogc unittest
880 {
881 assert(powmod(1U, 10U, 3U) == 1);
882 assert(powmod(3U, 2U, 6U) == 3);
883 assert(powmod(5U, 5U, 15U) == 5);
884 assert(powmod(2U, 3U, 5U) == 3);
885 assert(powmod(2U, 4U, 5U) == 1);
886 assert(powmod(2U, 5U, 5U) == 2);
887 }
888
889 @safe pure nothrow @nogc unittest
890 {
891 ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u;
892 assert(powmod(a, b, c) == 95367431640625u);
893 a = 100; b = 7919; c = 18446744073709551557u;
894 assert(powmod(a, b, c) == 18223853583554725198u);
895 a = 117; b = 7919; c = 18446744073709551557u;
896 assert(powmod(a, b, c) == 11493139548346411394u);
897 a = 134; b = 7919; c = 18446744073709551557u;
898 assert(powmod(a, b, c) == 10979163786734356774u);
899 a = 151; b = 7919; c = 18446744073709551557u;
900 assert(powmod(a, b, c) == 7023018419737782840u);
901 a = 168; b = 7919; c = 18446744073709551557u;
902 assert(powmod(a, b, c) == 58082701842386811u);
903 a = 185; b = 7919; c = 18446744073709551557u;
904 assert(powmod(a, b, c) == 17423478386299876798u);
905 a = 202; b = 7919; c = 18446744073709551557u;
906 assert(powmod(a, b, c) == 5522733478579799075u);
907 a = 219; b = 7919; c = 18446744073709551557u;
908 assert(powmod(a, b, c) == 15230218982491623487u);
909 a = 236; b = 7919; c = 18446744073709551557u;
910 assert(powmod(a, b, c) == 5198328724976436000u);
911
912 a = 0; b = 7919; c = 18446744073709551557u;
913 assert(powmod(a, b, c) == 0);
914 a = 123; b = 0; c = 18446744073709551557u;
915 assert(powmod(a, b, c) == 1);
916
917 immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u;
918 assert(powmod(a1, b1, c1) == 3883707345459248860u);
919
920 uint x = 100 ,y = 7919, z = 1844674407u;
921 assert(powmod(x, y, z) == 1613100340u);
922 x = 134; y = 7919; z = 1844674407u;
923 assert(powmod(x, y, z) == 734956622u);
924 x = 151; y = 7919; z = 1844674407u;
925 assert(powmod(x, y, z) == 1738696945u);
926 x = 168; y = 7919; z = 1844674407u;
927 assert(powmod(x, y, z) == 1247580927u);
928 x = 185; y = 7919; z = 1844674407u;
929 assert(powmod(x, y, z) == 1293855176u);
930 x = 202; y = 7919; z = 1844674407u;
931 assert(powmod(x, y, z) == 1566963682u);
932 x = 219; y = 7919; z = 1844674407u;
933 assert(powmod(x, y, z) == 181227807u);
934 x = 236; y = 7919; z = 1844674407u;
935 assert(powmod(x, y, z) == 217988321u);
936 x = 253; y = 7919; z = 1844674407u;
937 assert(powmod(x, y, z) == 1588843243u);
938
939 x = 0; y = 7919; z = 184467u;
940 assert(powmod(x, y, z) == 0);
941 x = 123; y = 0; z = 1844674u;
942 assert(powmod(x, y, z) == 1);
943
944 immutable ubyte x1 = 117;
945 immutable uint y1 = 7919;
946 immutable uint z1 = 1844674407u;
947 auto res = powmod(x1, y1, z1);
948 assert(is(typeof(res) == uint));
949 assert(res == 9479781u);
950
951 immutable ushort x2 = 123;
952 immutable uint y2 = 203;
953 immutable ubyte z2 = 113;
954 auto res2 = powmod(x2, y2, z2);
955 assert(is(typeof(res2) == ushort));
956 assert(res2 == 42u);
957 }
958
959 /**
960 * Calculates e$(SUPERSCRIPT x).
961 *
962 * $(TABLE_SV
963 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) )
964 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
965 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
966 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
967 * )
968 */
969 pragma(inline, true)
970 real exp(real x) @trusted pure nothrow @nogc // TODO: @safe
971 {
972 version (InlineAsm_X87)
973 {
974 import std.math.constants : LOG2E;
975
976 // e^^x = 2^^(LOG2E*x)
977 // (This is valid because the overflow & underflow limits for exp
978 // and exp2 are so similar).
979 if (!__ctfe)
980 return exp2Asm(LOG2E*x);
981 }
982 return expImpl(x);
983 }
984
985 /// ditto
986 pragma(inline, true)
987 double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); }
988
989 /// ditto
990 pragma(inline, true)
991 float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); }
992
993 ///
994 @safe unittest
995 {
996 import std.math.operations : feqrel;
997 import std.math.constants : E;
998
999 assert(exp(0.0) == 1.0);
1000 assert(exp(3.0).feqrel(E * E * E) > 16);
1001 }
1002
1003 private T expImpl(T)(T x) @safe pure nothrow @nogc
1004 {
1005 import std.math : floatTraits, RealFormat;
1006 import std.math.traits : isNaN;
1007 import std.math.rounding : floor;
1008 import std.math.algebraic : poly;
1009 import std.math.constants : LOG2E;
1010
1011 alias F = floatTraits!T;
1012 static if (F.realFormat == RealFormat.ieeeSingle)
1013 {
1014 static immutable T[6] P = [
1015 5.0000001201E-1,
1016 1.6666665459E-1,
1017 4.1665795894E-2,
1018 8.3334519073E-3,
1019 1.3981999507E-3,
1020 1.9875691500E-4,
1021 ];
1022
1023 enum T C1 = 0.693359375;
1024 enum T C2 = -2.12194440e-4;
1025
1026 // Overflow and Underflow limits.
1027 enum T OF = 88.72283905206835;
1028 enum T UF = -103.278929903431851103; // ln(2^-149)
1029 }
1030 else static if (F.realFormat == RealFormat.ieeeDouble)
1031 {
1032 // Coefficients for exp(x)
1033 static immutable T[3] P = [
1034 9.99999999999999999910E-1L,
1035 3.02994407707441961300E-2L,
1036 1.26177193074810590878E-4L,
1037 ];
1038 static immutable T[4] Q = [
1039 2.00000000000000000009E0L,
1040 2.27265548208155028766E-1L,
1041 2.52448340349684104192E-3L,
1042 3.00198505138664455042E-6L,
1043 ];
1044
1045 // C1 + C2 = LN2.
1046 enum T C1 = 6.93145751953125E-1;
1047 enum T C2 = 1.42860682030941723212E-6;
1048
1049 // Overflow and Underflow limits.
1050 enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024)
1051 enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
1052 }
1053 else static if (F.realFormat == RealFormat.ieeeExtended ||
1054 F.realFormat == RealFormat.ieeeExtended53)
1055 {
1056 // Coefficients for exp(x)
1057 static immutable T[3] P = [
1058 9.9999999999999999991025E-1L,
1059 3.0299440770744196129956E-2L,
1060 1.2617719307481059087798E-4L,
1061 ];
1062 static immutable T[4] Q = [
1063 2.0000000000000000000897E0L,
1064 2.2726554820815502876593E-1L,
1065 2.5244834034968410419224E-3L,
1066 3.0019850513866445504159E-6L,
1067 ];
1068
1069 // C1 + C2 = LN2.
1070 enum T C1 = 6.9314575195312500000000E-1L;
1071 enum T C2 = 1.4286068203094172321215E-6L;
1072
1073 // Overflow and Underflow limits.
1074 enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384)
1075 enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446)
1076 }
1077 else static if (F.realFormat == RealFormat.ieeeQuadruple)
1078 {
1079 // Coefficients for exp(x) - 1
1080 static immutable T[5] P = [
1081 9.999999999999999999999999999999999998502E-1L,
1082 3.508710990737834361215404761139478627390E-2L,
1083 2.708775201978218837374512615596512792224E-4L,
1084 6.141506007208645008909088812338454698548E-7L,
1085 3.279723985560247033712687707263393506266E-10L
1086 ];
1087 static immutable T[6] Q = [
1088 2.000000000000000000000000000000000000150E0,
1089 2.368408864814233538909747618894558968880E-1L,
1090 3.611828913847589925056132680618007270344E-3L,
1091 1.504792651814944826817779302637284053660E-5L,
1092 1.771372078166251484503904874657985291164E-8L,
1093 2.980756652081995192255342779918052538681E-12L
1094 ];
1095
1096 // C1 + C2 = LN2.
1097 enum T C1 = 6.93145751953125E-1L;
1098 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1099
1100 // Overflow and Underflow limits.
1101 enum T OF = 1.135583025911358400418251384584930671458833e4L;
1102 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1103 }
1104 else
1105 static assert(0, "Not implemented for this architecture");
1106
1107 // Special cases.
1108 if (isNaN(x))
1109 return x;
1110 if (x > OF)
1111 return real.infinity;
1112 if (x < UF)
1113 return 0.0;
1114
1115 // Express: e^^x = e^^g * 2^^n
1116 // = e^^g * e^^(n * LOG2E)
1117 // = e^^(g + n * LOG2E)
1118 T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5);
1119 const int n = cast(int) xx;
1120 x -= xx * C1;
1121 x -= xx * C2;
1122
1123 static if (F.realFormat == RealFormat.ieeeSingle)
1124 {
1125 xx = x * x;
1126 x = poly(x, P) * xx + x + 1.0f;
1127 }
1128 else
1129 {
1130 // Rational approximation for exponential of the fractional part:
1131 // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
1132 xx = x * x;
1133 const T px = x * poly(xx, P);
1134 x = px / (poly(xx, Q) - px);
1135 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
1136 }
1137
1138 // Scale by power of 2.
1139 x = core.math.ldexp(x, n);
1140
1141 return x;
1142 }
1143
1144 @safe @nogc nothrow unittest
1145 {
1146 import std.math : floatTraits, RealFormat;
1147 import std.math.operations : NaN, feqrel, isClose;
1148 import std.math.constants : E;
1149 import std.math.traits : isIdentical;
1150 import std.math.algebraic : abs;
1151
1152 version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags;
1153 version (FloatingPointControlSupport)
1154 {
1155 import std.math.hardware : FloatingPointControl;
1156
1157 FloatingPointControl ctrl;
1158 if (FloatingPointControl.hasExceptionTraps)
1159 ctrl.disableExceptions(FloatingPointControl.allExceptions);
1160 ctrl.rounding = FloatingPointControl.roundToNearest;
1161 }
1162
1163 static void testExp(T)()
1164 {
1165 enum realFormat = floatTraits!T.realFormat;
1166 static if (realFormat == RealFormat.ieeeQuadruple)
1167 {
1168 static immutable T[2][] exptestpoints =
1169 [ // x exp(x)
1170 [ 1.0L, E ],
1171 [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ],
1172 [ 3.0L, E*E*E ],
1173 [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow
1174 [ 0x1.7p+13L, T.infinity ], // close overflow
1175 [ 0x1p+80L, T.infinity ], // far overflow
1176 [ T.infinity, T.infinity ],
1177 [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow
1178 [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto
1179 [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal
1180 [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto
1181 [-0x1.655p+13L, 0 ], // close underflow
1182 [-0x1p+30L, 0 ], // far underflow
1183 ];
1184 }
1185 else static if (realFormat == RealFormat.ieeeExtended ||
1186 realFormat == RealFormat.ieeeExtended53)
1187 {
1188 static immutable T[2][] exptestpoints =
1189 [ // x exp(x)
1190 [ 1.0L, E ],
1191 [ 0.5L, 0x1.a61298e1e069bc97p+0L ],
1192 [ 3.0L, E*E*E ],
1193 [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow
1194 [ 0x1.7p+13L, T.infinity ], // close overflow
1195 [ 0x1p+80L, T.infinity ], // far overflow
1196 [ T.infinity, T.infinity ],
1197 [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow
1198 [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto
1199 [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal
1200 [-0x1.643p+13L, 0x1p-16444L ], // ditto
1201 [-0x1.645p+13L, 0 ], // close underflow
1202 [-0x1p+30L, 0 ], // far underflow
1203 ];
1204 }
1205 else static if (realFormat == RealFormat.ieeeDouble)
1206 {
1207 static immutable T[2][] exptestpoints =
1208 [ // x, exp(x)
1209 [ 1.0L, E ],
1210 [ 0.5L, 0x1.a61298e1e069cp+0L ],
1211 [ 3.0L, E*E*E ],
1212 [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow
1213 [ 0x1.7p+9L, T.infinity ], // close overflow
1214 [ 0x1p+80L, T.infinity ], // far overflow
1215 [ T.infinity, T.infinity ],
1216 [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow
1217 [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal
1218 [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto
1219 [-0x1.8p+9L, 0 ], // close underflow
1220 [-0x1p+30L, 0 ], // far underflow
1221 ];
1222 }
1223 else static if (realFormat == RealFormat.ieeeSingle)
1224 {
1225 static immutable T[2][] exptestpoints =
1226 [ // x, exp(x)
1227 [ 1.0L, E ],
1228 [ 0.5L, 0x1.a61299p+0L ],
1229 [ 3.0L, E*E*E ],
1230 [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow
1231 [ 0x1.7p+6L, T.infinity ], // close overflow
1232 [ 0x1p+80L, T.infinity ], // far overflow
1233 [ T.infinity, T.infinity ],
1234 [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow
1235 [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal
1236 [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto
1237 [-0x1.ap+6L, 0 ], // close underflow
1238 [-0x1p+30L, 0 ], // far underflow
1239 ];
1240 }
1241 else
1242 static assert(0, "No exp() tests for real type!");
1243
1244 const minEqualMantissaBits = T.mant_dig - 13;
1245 T x;
1246 version (IeeeFlagsSupport) IeeeFlags f;
1247 foreach (ref pair; exptestpoints)
1248 {
1249 version (IeeeFlagsSupport) resetIeeeFlags();
1250 x = exp(pair[0]);
1251 //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]);
1252 assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
1253 }
1254
1255 // Ideally, exp(0) would not set the inexact flag.
1256 // Unfortunately, fldl2e sets it!
1257 // So it's not realistic to avoid setting it.
1258 assert(exp(cast(T) 0.0) == 1.0);
1259
1260 // NaN propagation. Doesn't set flags, bcos was already NaN.
1261 version (IeeeFlagsSupport)
1262 {
1263 resetIeeeFlags();
1264 x = exp(T.nan);
1265 f = ieeeFlags;
1266 assert(isIdentical(abs(x), T.nan));
1267 assert(f.flags == 0);
1268
1269 resetIeeeFlags();
1270 x = exp(-T.nan);
1271 f = ieeeFlags;
1272 assert(isIdentical(abs(x), T.nan));
1273 assert(f.flags == 0);
1274 }
1275 else
1276 {
1277 x = exp(T.nan);
1278 assert(isIdentical(abs(x), T.nan));
1279
1280 x = exp(-T.nan);
1281 assert(isIdentical(abs(x), T.nan));
1282 }
1283
1284 x = exp(NaN(0x123));
1285 assert(isIdentical(x, NaN(0x123)));
1286 }
1287
1288 import std.meta : AliasSeq;
1289 foreach (T; AliasSeq!(real, double, float))
1290 testExp!T();
1291
1292 // High resolution test (verified against GNU MPFR/Mathematica).
1293 assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L);
1294
1295 assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1296 }
1297
1298 /**
1299 * Calculates the value of the natural logarithm base (e)
1300 * raised to the power of x, minus 1.
1301 *
1302 * For very small x, expm1(x) is more accurate
1303 * than exp(x)-1.
1304 *
1305 * $(TABLE_SV
1306 * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) )
1307 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
1308 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1309 * $(TR $(TD -$(INFIN)) $(TD -1.0) )
1310 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1311 * )
1312 */
1313 pragma(inline, true)
1314 real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe
1315 {
1316 version (InlineAsm_X87)
1317 {
1318 if (!__ctfe)
1319 return expm1Asm(x);
1320 }
1321 return expm1Impl(x);
1322 }
1323
1324 /// ditto
1325 pragma(inline, true)
1326 double expm1(double x) @safe pure nothrow @nogc
1327 {
1328 return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x);
1329 }
1330
1331 /// ditto
1332 pragma(inline, true)
1333 float expm1(float x) @safe pure nothrow @nogc
1334 {
1335 // no single-precision version in Cephes => use double precision
1336 return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x);
1337 }
1338
1339 ///
1340 @safe unittest
1341 {
1342 import std.math.traits : isIdentical;
1343 import std.math.operations : feqrel;
1344
1345 assert(isIdentical(expm1(0.0), 0.0));
1346 assert(expm1(1.0).feqrel(1.71828) > 16);
1347 assert(expm1(2.0).feqrel(6.3890) > 16);
1348 }
1349
1350 version (InlineAsm_X87)
1351 private real expm1Asm(real x) @trusted pure nothrow @nogc
1352 {
1353 version (X86)
1354 {
1355 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1356 asm pure nothrow @nogc
1357 {
1358 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1359 * Author: Don Clugston.
1360 *
1361 * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
1362 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
1363 * and 2ym1 = (2^^(y-rndint(y))-1).
1364 * If 2rndy < 0.5*real.epsilon, result is -1.
1365 * Implementation is otherwise the same as for exp2()
1366 */
1367 naked;
1368 fld real ptr [ESP+4] ; // x
1369 mov AX, [ESP+4+8]; // AX = exponent and sign
1370 sub ESP, 12+8; // Create scratch space on the stack
1371 // [ESP,ESP+2] = scratchint
1372 // [ESP+4..+6, +8..+10, +10] = scratchreal
1373 // set scratchreal mantissa = 1.0
1374 mov dword ptr [ESP+8], 0;
1375 mov dword ptr [ESP+8+4], 0x80000000;
1376 and AX, 0x7FFF; // drop sign bit
1377 cmp AX, 0x401D; // avoid InvalidException in fist
1378 jae L_extreme;
1379 fldl2e;
1380 fmulp ST(1), ST; // y = x*log2(e)
1381 fist dword ptr [ESP]; // scratchint = rndint(y)
1382 fisub dword ptr [ESP]; // y - rndint(y)
1383 // and now set scratchreal exponent
1384 mov EAX, [ESP];
1385 add EAX, 0x3fff;
1386 jle short L_largenegative;
1387 cmp EAX,0x8000;
1388 jge short L_largepositive;
1389 mov [ESP+8+8],AX;
1390 f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1391 fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1392 fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1
1393 fld1;
1394 fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1395 faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1396 add ESP,12+8;
1397 ret PARAMSIZE;
1398
1399 L_extreme: // Extreme exponent. X is very large positive, very
1400 // large negative, infinity, or NaN.
1401 fxam;
1402 fstsw AX;
1403 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1404 jz L_was_nan; // if x is NaN, returns x
1405 test AX, 0x0200;
1406 jnz L_largenegative;
1407 L_largepositive:
1408 // Set scratchreal = real.max.
1409 // squaring it will create infinity, and set overflow flag.
1410 mov word ptr [ESP+8+8], 0x7FFE;
1411 fstp ST(0);
1412 fld real ptr [ESP+8]; // load scratchreal
1413 fmul ST(0), ST; // square it, to create havoc!
1414 L_was_nan:
1415 add ESP,12+8;
1416 ret PARAMSIZE;
1417 L_largenegative:
1418 fstp ST(0);
1419 fld1;
1420 fchs; // return -1. Underflow flag is not set.
1421 add ESP,12+8;
1422 ret PARAMSIZE;
1423 }
1424 }
1425 else version (X86_64)
1426 {
1427 asm pure nothrow @nogc
1428 {
1429 naked;
1430 }
1431 version (Win64)
1432 {
1433 asm pure nothrow @nogc
1434 {
1435 fld real ptr [RCX]; // x
1436 mov AX,[RCX+8]; // AX = exponent and sign
1437 }
1438 }
1439 else
1440 {
1441 asm pure nothrow @nogc
1442 {
1443 fld real ptr [RSP+8]; // x
1444 mov AX,[RSP+8+8]; // AX = exponent and sign
1445 }
1446 }
1447 asm pure nothrow @nogc
1448 {
1449 /* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1450 * Author: Don Clugston.
1451 *
1452 * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1453 * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1454 * and 2ym1 = (2^(y-rndint(y))-1).
1455 * If 2rndy < 0.5*real.epsilon, result is -1.
1456 * Implementation is otherwise the same as for exp2()
1457 */
1458 sub RSP, 24; // Create scratch space on the stack
1459 // [RSP,RSP+2] = scratchint
1460 // [RSP+4..+6, +8..+10, +10] = scratchreal
1461 // set scratchreal mantissa = 1.0
1462 mov dword ptr [RSP+8], 0;
1463 mov dword ptr [RSP+8+4], 0x80000000;
1464 and AX, 0x7FFF; // drop sign bit
1465 cmp AX, 0x401D; // avoid InvalidException in fist
1466 jae L_extreme;
1467 fldl2e;
1468 fmul ; // y = x*log2(e)
1469 fist dword ptr [RSP]; // scratchint = rndint(y)
1470 fisub dword ptr [RSP]; // y - rndint(y)
1471 // and now set scratchreal exponent
1472 mov EAX, [RSP];
1473 add EAX, 0x3fff;
1474 jle short L_largenegative;
1475 cmp EAX,0x8000;
1476 jge short L_largepositive;
1477 mov [RSP+8+8],AX;
1478 f2xm1; // 2^(y-rndint(y)) -1
1479 fld real ptr [RSP+8] ; // 2^rndint(y)
1480 fmul ST(1), ST;
1481 fld1;
1482 fsubp ST(1), ST;
1483 fadd;
1484 add RSP,24;
1485 ret;
1486
1487 L_extreme: // Extreme exponent. X is very large positive, very
1488 // large negative, infinity, or NaN.
1489 fxam;
1490 fstsw AX;
1491 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1492 jz L_was_nan; // if x is NaN, returns x
1493 test AX, 0x0200;
1494 jnz L_largenegative;
1495 L_largepositive:
1496 // Set scratchreal = real.max.
1497 // squaring it will create infinity, and set overflow flag.
1498 mov word ptr [RSP+8+8], 0x7FFE;
1499 fstp ST(0);
1500 fld real ptr [RSP+8]; // load scratchreal
1501 fmul ST(0), ST; // square it, to create havoc!
1502 L_was_nan:
1503 add RSP,24;
1504 ret;
1505
1506 L_largenegative:
1507 fstp ST(0);
1508 fld1;
1509 fchs; // return -1. Underflow flag is not set.
1510 add RSP,24;
1511 ret;
1512 }
1513 }
1514 else
1515 static assert(0);
1516 }
1517
1518 private T expm1Impl(T)(T x) @safe pure nothrow @nogc
1519 {
1520 import std.math : floatTraits, RealFormat;
1521 import std.math.rounding : floor;
1522 import std.math.algebraic : poly;
1523 import std.math.constants : LN2;
1524
1525 // Coefficients for exp(x) - 1 and overflow/underflow limits.
1526 enum realFormat = floatTraits!T.realFormat;
1527 static if (realFormat == RealFormat.ieeeQuadruple)
1528 {
1529 static immutable T[8] P = [
1530 2.943520915569954073888921213330863757240E8L,
1531 -5.722847283900608941516165725053359168840E7L,
1532 8.944630806357575461578107295909719817253E6L,
1533 -7.212432713558031519943281748462837065308E5L,
1534 4.578962475841642634225390068461943438441E4L,
1535 -1.716772506388927649032068540558788106762E3L,
1536 4.401308817383362136048032038528753151144E1L,
1537 -4.888737542888633647784737721812546636240E-1L
1538 ];
1539
1540 static immutable T[9] Q = [
1541 1.766112549341972444333352727998584753865E9L,
1542 -7.848989743695296475743081255027098295771E8L,
1543 1.615869009634292424463780387327037251069E8L,
1544 -2.019684072836541751428967854947019415698E7L,
1545 1.682912729190313538934190635536631941751E6L,
1546 -9.615511549171441430850103489315371768998E4L,
1547 3.697714952261803935521187272204485251835E3L,
1548 -8.802340681794263968892934703309274564037E1L,
1549 1.0
1550 ];
1551
1552 enum T OF = 1.1356523406294143949491931077970764891253E4L;
1553 enum T UF = -1.143276959615573793352782661133116431383730e4L;
1554 }
1555 else static if (realFormat == RealFormat.ieeeExtended)
1556 {
1557 static immutable T[5] P = [
1558 -1.586135578666346600772998894928250240826E4L,
1559 2.642771505685952966904660652518429479531E3L,
1560 -3.423199068835684263987132888286791620673E2L,
1561 1.800826371455042224581246202420972737840E1L,
1562 -5.238523121205561042771939008061958820811E-1L,
1563 ];
1564 static immutable T[6] Q = [
1565 -9.516813471998079611319047060563358064497E4L,
1566 3.964866271411091674556850458227710004570E4L,
1567 -7.207678383830091850230366618190187434796E3L,
1568 7.206038318724600171970199625081491823079E2L,
1569 -4.002027679107076077238836622982900945173E1L,
1570 1.0
1571 ];
1572
1573 enum T OF = 1.1356523406294143949492E4L;
1574 enum T UF = -4.5054566736396445112120088E1L;
1575 }
1576 else static if (realFormat == RealFormat.ieeeDouble)
1577 {
1578 static immutable T[3] P = [
1579 9.9999999999999999991025E-1,
1580 3.0299440770744196129956E-2,
1581 1.2617719307481059087798E-4,
1582 ];
1583 static immutable T[4] Q = [
1584 2.0000000000000000000897E0,
1585 2.2726554820815502876593E-1,
1586 2.5244834034968410419224E-3,
1587 3.0019850513866445504159E-6,
1588 ];
1589 }
1590 else
1591 static assert(0, "no coefficients for expm1()");
1592
1593 static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
1594 {
1595 if (x < -0.5 || x > 0.5)
1596 return exp(x) - 1.0;
1597 if (x == 0.0)
1598 return x;
1599
1600 const T xx = x * x;
1601 x = x * poly(xx, P);
1602 x = x / (poly(xx, Q) - x);
1603 return x + x;
1604 }
1605 else
1606 {
1607 // C1 + C2 = LN2.
1608 enum T C1 = 6.9314575195312500000000E-1L;
1609 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
1610
1611 // Special cases.
1612 if (x > OF)
1613 return real.infinity;
1614 if (x == cast(T) 0.0)
1615 return x;
1616 if (x < UF)
1617 return -1.0;
1618
1619 // Express x = LN2 (n + remainder), remainder not exceeding 1/2.
1620 int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2);
1621 x -= n * C1;
1622 x -= n * C2;
1623
1624 // Rational approximation:
1625 // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x)
1626 T px = x * poly(x, P);
1627 T qx = poly(x, Q);
1628 const T xx = x * x;
1629 qx = x + ((cast(T) 0.5) * xx + xx * px / qx);
1630
1631 // We have qx = exp(remainder LN2) - 1, so:
1632 // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1.
1633 px = core.math.ldexp(cast(T) 1.0, n);
1634 x = px * qx + (px - cast(T) 1.0);
1635
1636 return x;
1637 }
1638 }
1639
1640 @safe @nogc nothrow unittest
1641 {
1642 import std.math.traits : isNaN;
1643 import std.math.operations : isClose, CommonDefaultFor;
1644
1645 static void testExpm1(T)()
1646 {
1647 // NaN
1648 assert(isNaN(expm1(cast(T) T.nan)));
1649
1650 static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ];
1651 foreach (x; xs)
1652 {
1653 const T e = expm1(x);
1654 const T r = exp(x) - 1;
1655
1656 //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
1657 assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
1658 }
1659 }
1660
1661 import std.meta : AliasSeq;
1662 foreach (T; AliasSeq!(real, double))
1663 testExpm1!T();
1664 }
1665
1666 /**
1667 * Calculates 2$(SUPERSCRIPT x).
1668 *
1669 * $(TABLE_SV
1670 * $(TR $(TH x) $(TH exp2(x)) )
1671 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1672 * $(TR $(TD -$(INFIN)) $(TD +0.0) )
1673 * $(TR $(TD $(NAN)) $(TD $(NAN)) )
1674 * )
1675 */
1676 pragma(inline, true)
1677 real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe
1678 {
1679 version (InlineAsm_X87)
1680 {
1681 if (!__ctfe)
1682 return exp2Asm(x);
1683 }
1684 return exp2Impl(x);
1685 }
1686
1687 /// ditto
1688 pragma(inline, true)
1689 double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); }
1690
1691 /// ditto
1692 pragma(inline, true)
1693 float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); }
1694
1695 ///
1696 @safe unittest
1697 {
1698 import std.math.traits : isIdentical;
1699 import std.math.operations : feqrel;
1700
1701 assert(isIdentical(exp2(0.0), 1.0));
1702 assert(exp2(2.0).feqrel(4.0) > 16);
1703 assert(exp2(8.0).feqrel(256.0) > 16);
1704 }
1705
1706 @safe unittest
1707 {
1708 version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented
1709 {
1710 assert( core.stdc.math.exp2f(0.0f) == 1 );
1711 assert( core.stdc.math.exp2 (0.0) == 1 );
1712 assert( core.stdc.math.exp2l(0.0L) == 1 );
1713 }
1714 }
1715
1716 version (InlineAsm_X87)
1717 private real exp2Asm(real x) @nogc @trusted pure nothrow
1718 {
1719 version (X86)
1720 {
1721 enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4
1722
1723 asm pure nothrow @nogc
1724 {
1725 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1726 * Author: Don Clugston.
1727 *
1728 * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1729 * The trick for high performance is to avoid the fscale(28cycles on core2),
1730 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1731 *
1732 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1733 * because it will set the Invalid Operation flag if overflow or NaN occurs.
1734 * Fortunately, whenever this happens the result would be zero or infinity.
1735 *
1736 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1737 * work for the (very rare) cases where the result is subnormal. So we fall back
1738 * to the slow method in that case.
1739 */
1740 naked;
1741 fld real ptr [ESP+4] ; // x
1742 mov AX, [ESP+4+8]; // AX = exponent and sign
1743 sub ESP, 12+8; // Create scratch space on the stack
1744 // [ESP,ESP+2] = scratchint
1745 // [ESP+4..+6, +8..+10, +10] = scratchreal
1746 // set scratchreal mantissa = 1.0
1747 mov dword ptr [ESP+8], 0;
1748 mov dword ptr [ESP+8+4], 0x80000000;
1749 and AX, 0x7FFF; // drop sign bit
1750 cmp AX, 0x401D; // avoid InvalidException in fist
1751 jae L_extreme;
1752 fist dword ptr [ESP]; // scratchint = rndint(x)
1753 fisub dword ptr [ESP]; // x - rndint(x)
1754 // and now set scratchreal exponent
1755 mov EAX, [ESP];
1756 add EAX, 0x3fff;
1757 jle short L_subnormal;
1758 cmp EAX,0x8000;
1759 jge short L_overflow;
1760 mov [ESP+8+8],AX;
1761 L_normal:
1762 f2xm1;
1763 fld1;
1764 faddp ST(1), ST; // 2^^(x-rndint(x))
1765 fld real ptr [ESP+8] ; // 2^^rndint(x)
1766 add ESP,12+8;
1767 fmulp ST(1), ST;
1768 ret PARAMSIZE;
1769
1770 L_subnormal:
1771 // Result will be subnormal.
1772 // In this rare case, the simple poking method doesn't work.
1773 // The speed doesn't matter, so use the slow fscale method.
1774 fild dword ptr [ESP]; // scratchint
1775 fld1;
1776 fscale;
1777 fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1778 fstp ST(0); // drop scratchint
1779 jmp L_normal;
1780
1781 L_extreme: // Extreme exponent. X is very large positive, very
1782 // large negative, infinity, or NaN.
1783 fxam;
1784 fstsw AX;
1785 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1786 jz L_was_nan; // if x is NaN, returns x
1787 // set scratchreal = real.min_normal
1788 // squaring it will return 0, setting underflow flag
1789 mov word ptr [ESP+8+8], 1;
1790 test AX, 0x0200;
1791 jnz L_waslargenegative;
1792 L_overflow:
1793 // Set scratchreal = real.max.
1794 // squaring it will create infinity, and set overflow flag.
1795 mov word ptr [ESP+8+8], 0x7FFE;
1796 L_waslargenegative:
1797 fstp ST(0);
1798 fld real ptr [ESP+8]; // load scratchreal
1799 fmul ST(0), ST; // square it, to create havoc!
1800 L_was_nan:
1801 add ESP,12+8;
1802 ret PARAMSIZE;
1803 }
1804 }
1805 else version (X86_64)
1806 {
1807 asm pure nothrow @nogc
1808 {
1809 naked;
1810 }
1811 version (Win64)
1812 {
1813 asm pure nothrow @nogc
1814 {
1815 fld real ptr [RCX]; // x
1816 mov AX,[RCX+8]; // AX = exponent and sign
1817 }
1818 }
1819 else
1820 {
1821 asm pure nothrow @nogc
1822 {
1823 fld real ptr [RSP+8]; // x
1824 mov AX,[RSP+8+8]; // AX = exponent and sign
1825 }
1826 }
1827 asm pure nothrow @nogc
1828 {
1829 /* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1830 * Author: Don Clugston.
1831 *
1832 * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1833 * The trick for high performance is to avoid the fscale(28cycles on core2),
1834 * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1835 *
1836 * We can do frndint by using fist. BUT we can't use it for huge numbers,
1837 * because it will set the Invalid Operation flag is overflow or NaN occurs.
1838 * Fortunately, whenever this happens the result would be zero or infinity.
1839 *
1840 * We can perform fscale by directly poking into the exponent. BUT this doesn't
1841 * work for the (very rare) cases where the result is subnormal. So we fall back
1842 * to the slow method in that case.
1843 */
1844 sub RSP, 24; // Create scratch space on the stack
1845 // [RSP,RSP+2] = scratchint
1846 // [RSP+4..+6, +8..+10, +10] = scratchreal
1847 // set scratchreal mantissa = 1.0
1848 mov dword ptr [RSP+8], 0;
1849 mov dword ptr [RSP+8+4], 0x80000000;
1850 and AX, 0x7FFF; // drop sign bit
1851 cmp AX, 0x401D; // avoid InvalidException in fist
1852 jae L_extreme;
1853 fist dword ptr [RSP]; // scratchint = rndint(x)
1854 fisub dword ptr [RSP]; // x - rndint(x)
1855 // and now set scratchreal exponent
1856 mov EAX, [RSP];
1857 add EAX, 0x3fff;
1858 jle short L_subnormal;
1859 cmp EAX,0x8000;
1860 jge short L_overflow;
1861 mov [RSP+8+8],AX;
1862 L_normal:
1863 f2xm1;
1864 fld1;
1865 fadd; // 2^(x-rndint(x))
1866 fld real ptr [RSP+8] ; // 2^rndint(x)
1867 add RSP,24;
1868 fmulp ST(1), ST;
1869 ret;
1870
1871 L_subnormal:
1872 // Result will be subnormal.
1873 // In this rare case, the simple poking method doesn't work.
1874 // The speed doesn't matter, so use the slow fscale method.
1875 fild dword ptr [RSP]; // scratchint
1876 fld1;
1877 fscale;
1878 fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1879 fstp ST(0); // drop scratchint
1880 jmp L_normal;
1881
1882 L_extreme: // Extreme exponent. X is very large positive, very
1883 // large negative, infinity, or NaN.
1884 fxam;
1885 fstsw AX;
1886 test AX, 0x0400; // NaN_or_zero, but we already know x != 0
1887 jz L_was_nan; // if x is NaN, returns x
1888 // set scratchreal = real.min
1889 // squaring it will return 0, setting underflow flag
1890 mov word ptr [RSP+8+8], 1;
1891 test AX, 0x0200;
1892 jnz L_waslargenegative;
1893 L_overflow:
1894 // Set scratchreal = real.max.
1895 // squaring it will create infinity, and set overflow flag.
1896 mov word ptr [RSP+8+8], 0x7FFE;
1897 L_waslargenegative:
1898 fstp ST(0);
1899 fld real ptr [RSP+8]; // load scratchreal
1900 fmul ST(0), ST; // square it, to create havoc!
1901 L_was_nan:
1902 add RSP,24;
1903 ret;
1904 }
1905 }
1906 else
1907 static assert(0);
1908 }
1909
1910 private T exp2Impl(T)(T x) @nogc @safe pure nothrow
1911 {
1912 import std.math : floatTraits, RealFormat;
1913 import std.math.traits : isNaN;
1914 import std.math.rounding : floor;
1915 import std.math.algebraic : poly;
1916
1917 // Coefficients for exp2(x)
1918 enum realFormat = floatTraits!T.realFormat;
1919 static if (realFormat == RealFormat.ieeeQuadruple)
1920 {
1921 static immutable T[5] P = [
1922 9.079594442980146270952372234833529694788E12L,
1923 1.530625323728429161131811299626419117557E11L,
1924 5.677513871931844661829755443994214173883E8L,
1925 6.185032670011643762127954396427045467506E5L,
1926 1.587171580015525194694938306936721666031E2L
1927 ];
1928
1929 static immutable T[6] Q = [
1930 2.619817175234089411411070339065679229869E13L,
1931 1.490560994263653042761789432690793026977E12L,
1932 1.092141473886177435056423606755843616331E10L,
1933 2.186249607051644894762167991800811827835E7L,
1934 1.236602014442099053716561665053645270207E4L,
1935 1.0
1936 ];
1937 }
1938 else static if (realFormat == RealFormat.ieeeExtended)
1939 {
1940 static immutable T[3] P = [
1941 2.0803843631901852422887E6L,
1942 3.0286971917562792508623E4L,
1943 6.0614853552242266094567E1L,
1944 ];
1945 static immutable T[4] Q = [
1946 6.0027204078348487957118E6L,
1947 3.2772515434906797273099E5L,
1948 1.7492876999891839021063E3L,
1949 1.0000000000000000000000E0L,
1950 ];
1951 }
1952 else static if (realFormat == RealFormat.ieeeDouble)
1953 {
1954 static immutable T[3] P = [
1955 1.51390680115615096133E3L,
1956 2.02020656693165307700E1L,
1957 2.30933477057345225087E-2L,
1958 ];
1959 static immutable T[3] Q = [
1960 4.36821166879210612817E3L,
1961 2.33184211722314911771E2L,
1962 1.00000000000000000000E0L,
1963 ];
1964 }
1965 else static if (realFormat == RealFormat.ieeeSingle)
1966 {
1967 static immutable T[6] P = [
1968 6.931472028550421E-001L,
1969 2.402264791363012E-001L,
1970 5.550332471162809E-002L,
1971 9.618437357674640E-003L,
1972 1.339887440266574E-003L,
1973 1.535336188319500E-004L,
1974 ];
1975 }
1976 else
1977 static assert(0, "no coefficients for exp2()");
1978
1979 // Overflow and Underflow limits.
1980 enum T OF = T.max_exp;
1981 enum T UF = T.min_exp - 1;
1982
1983 // Special cases.
1984 if (isNaN(x))
1985 return x;
1986 if (x > OF)
1987 return real.infinity;
1988 if (x < UF)
1989 return 0.0;
1990
1991 static if (realFormat == RealFormat.ieeeSingle) // special case for single precision
1992 {
1993 // The following is necessary because range reduction blows up.
1994 if (x == 0.0f)
1995 return 1.0f;
1996
1997 // Separate into integer and fractional parts.
1998 const T i = floor(x);
1999 int n = cast(int) i;
2000 x -= i;
2001 if (x > 0.5f)
2002 {
2003 n += 1;
2004 x -= 1.0f;
2005 }
2006
2007 // Rational approximation:
2008 // exp2(x) = 1.0 + x P(x)
2009 x = 1.0f + x * poly(x, P);
2010 }
2011 else
2012 {
2013 // Separate into integer and fractional parts.
2014 const T i = floor(x + cast(T) 0.5);
2015 int n = cast(int) i;
2016 x -= i;
2017
2018 // Rational approximation:
2019 // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2))
2020 const T xx = x * x;
2021 const T px = x * poly(xx, P);
2022 x = px / (poly(xx, Q) - px);
2023 x = (cast(T) 1.0) + (cast(T) 2.0) * x;
2024 }
2025
2026 // Scale by power of 2.
2027 x = core.math.ldexp(x, n);
2028
2029 return x;
2030 }
2031
2032 @safe @nogc nothrow unittest
2033 {
2034 import std.math.operations : feqrel, NaN, isClose;
2035 import std.math.traits : isIdentical;
2036 import std.math.constants : SQRT2;
2037
2038 assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1);
2039 assert(exp2(8.0L) == 256.0);
2040 assert(exp2(-9.0L)== 1.0L/512.0);
2041
2042 static void testExp2(T)()
2043 {
2044 // NaN
2045 const T specialNaN = NaN(0x0123L);
2046 assert(isIdentical(exp2(specialNaN), specialNaN));
2047
2048 // over-/underflow
2049 enum T OF = T.max_exp;
2050 enum T UF = T.min_exp - T.mant_dig;
2051 assert(isIdentical(exp2(OF + 1), cast(T) T.infinity));
2052 assert(isIdentical(exp2(UF - 1), cast(T) 0.0));
2053
2054 static immutable T[2][] vals =
2055 [
2056 // x, exp2(x)
2057 [ 0.0, 1.0 ],
2058 [ -0.0, 1.0 ],
2059 [ 0.5, SQRT2 ],
2060 [ 8.0, 256.0 ],
2061 [ -9.0, 1.0 / 512 ],
2062 ];
2063
2064 foreach (ref val; vals)
2065 {
2066 const T x = val[0];
2067 const T r = val[1];
2068 const T e = exp2(x);
2069
2070 //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r);
2071 assert(isClose(r, e));
2072 }
2073 }
2074
2075 import std.meta : AliasSeq;
2076 foreach (T; AliasSeq!(real, double, float))
2077 testExp2!T();
2078 }
2079
2080 /*********************************************************************
2081 * Separate floating point value into significand and exponent.
2082 *
2083 * Returns:
2084 * Calculate and return $(I x) and $(I exp) such that
2085 * value =$(I x)*2$(SUPERSCRIPT exp) and
2086 * .5 $(LT)= |$(I x)| $(LT) 1.0
2087 *
2088 * $(I x) has same sign as value.
2089 *
2090 * $(TABLE_SV
2091 * $(TR $(TH value) $(TH returns) $(TH exp))
2092 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0))
2093 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max))
2094 * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min))
2095 * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
2096 * )
2097 */
2098 T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc
2099 if (isFloatingPoint!T)
2100 {
2101 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2102 import std.math.traits : isSubnormal;
2103
2104 if (__ctfe)
2105 {
2106 // Handle special cases.
2107 if (value == 0) { exp = 0; return value; }
2108 if (value == T.infinity) { exp = int.max; return value; }
2109 if (value == -T.infinity || value != value) { exp = int.min; return value; }
2110 // Handle ordinary cases.
2111 // In CTFE there is no performance advantage for having separate
2112 // paths for different floating point types.
2113 T absValue = value < 0 ? -value : value;
2114 int expCount;
2115 static if (T.mant_dig > double.mant_dig)
2116 {
2117 for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L)
2118 expCount += 1024;
2119 for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L)
2120 expCount -= 1021;
2121 }
2122 const double dval = cast(double) absValue;
2123 int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2;
2124 dexp++;
2125 expCount += dexp;
2126 absValue *= 2.0 ^^ -dexp;
2127 // If the original value was subnormal or if it was a real
2128 // then absValue can still be outside the [0.5, 1.0) range.
2129 if (absValue < 0.5)
2130 {
2131 assert(T.mant_dig > double.mant_dig || isSubnormal(value));
2132 do
2133 {
2134 absValue += absValue;
2135 expCount--;
2136 } while (absValue < 0.5);
2137 }
2138 else
2139 {
2140 assert(absValue < 1 || T.mant_dig > double.mant_dig);
2141 for (; absValue >= 1; absValue *= T(0.5))
2142 expCount++;
2143 }
2144 exp = expCount;
2145 return value < 0 ? -absValue : absValue;
2146 }
2147
2148 Unqual!T vf = value;
2149 ushort* vu = cast(ushort*)&vf;
2150 static if (is(immutable T == immutable float))
2151 int* vi = cast(int*)&vf;
2152 else
2153 long* vl = cast(long*)&vf;
2154 int ex;
2155 alias F = floatTraits!T;
2156
2157 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2158 static if (F.realFormat == RealFormat.ieeeExtended ||
2159 F.realFormat == RealFormat.ieeeExtended53)
2160 {
2161 if (ex)
2162 { // If exponent is non-zero
2163 if (ex == F.EXPMASK) // infinity or NaN
2164 {
2165 if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2166 {
2167 *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ
2168 exp = int.min;
2169 }
2170 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2171 exp = int.min;
2172 else // positive infinity
2173 exp = int.max;
2174
2175 }
2176 else
2177 {
2178 exp = ex - F.EXPBIAS;
2179 vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2180 }
2181 }
2182 else if (!*vl)
2183 {
2184 // vf is +-0.0
2185 exp = 0;
2186 }
2187 else
2188 {
2189 // subnormal
2190
2191 vf *= F.RECIP_EPSILON;
2192 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2193 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2194 vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE;
2195 }
2196 return vf;
2197 }
2198 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2199 {
2200 if (ex) // If exponent is non-zero
2201 {
2202 if (ex == F.EXPMASK)
2203 {
2204 // infinity or NaN
2205 if (vl[MANTISSA_LSB] |
2206 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2207 {
2208 // convert NaNS to NaNQ
2209 vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000;
2210 exp = int.min;
2211 }
2212 else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity
2213 exp = int.min;
2214 else // positive infinity
2215 exp = int.max;
2216 }
2217 else
2218 {
2219 exp = ex - F.EXPBIAS;
2220 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2221 }
2222 }
2223 else if ((vl[MANTISSA_LSB] |
2224 (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2225 {
2226 // vf is +-0.0
2227 exp = 0;
2228 }
2229 else
2230 {
2231 // subnormal
2232 vf *= F.RECIP_EPSILON;
2233 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2234 exp = ex - F.EXPBIAS - T.mant_dig + 1;
2235 vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]);
2236 }
2237 return vf;
2238 }
2239 else static if (F.realFormat == RealFormat.ieeeDouble)
2240 {
2241 if (ex) // If exponent is non-zero
2242 {
2243 if (ex == F.EXPMASK) // infinity or NaN
2244 {
2245 if (*vl == 0x7FF0_0000_0000_0000) // positive infinity
2246 {
2247 exp = int.max;
2248 }
2249 else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity
2250 exp = int.min;
2251 else
2252 { // NaN
2253 *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ
2254 exp = int.min;
2255 }
2256 }
2257 else
2258 {
2259 exp = (ex - F.EXPBIAS) >> 4;
2260 vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2261 }
2262 }
2263 else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF))
2264 {
2265 // vf is +-0.0
2266 exp = 0;
2267 }
2268 else
2269 {
2270 // subnormal
2271 vf *= F.RECIP_EPSILON;
2272 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2273 exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1;
2274 vu[F.EXPPOS_SHORT] =
2275 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0);
2276 }
2277 return vf;
2278 }
2279 else static if (F.realFormat == RealFormat.ieeeSingle)
2280 {
2281 if (ex) // If exponent is non-zero
2282 {
2283 if (ex == F.EXPMASK) // infinity or NaN
2284 {
2285 if (*vi == 0x7F80_0000) // positive infinity
2286 {
2287 exp = int.max;
2288 }
2289 else if (*vi == 0xFF80_0000) // negative infinity
2290 exp = int.min;
2291 else
2292 { // NaN
2293 *vi |= 0x0040_0000; // convert NaNS to NaNQ
2294 exp = int.min;
2295 }
2296 }
2297 else
2298 {
2299 exp = (ex - F.EXPBIAS) >> 7;
2300 vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00);
2301 }
2302 }
2303 else if (!(*vi & 0x7FFF_FFFF))
2304 {
2305 // vf is +-0.0
2306 exp = 0;
2307 }
2308 else
2309 {
2310 // subnormal
2311 vf *= F.RECIP_EPSILON;
2312 ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
2313 exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1;
2314 vu[F.EXPPOS_SHORT] =
2315 cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00);
2316 }
2317 return vf;
2318 }
2319 else // static if (F.realFormat == RealFormat.ibmExtended)
2320 {
2321 assert(0, "frexp not implemented");
2322 }
2323 }
2324
2325 ///
2326 @safe unittest
2327 {
2328 import std.math.operations : isClose;
2329
2330 int exp;
2331 real mantissa = frexp(123.456L, exp);
2332
2333 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L));
2334
2335 assert(frexp(-real.nan, exp) && exp == int.min);
2336 assert(frexp(real.nan, exp) && exp == int.min);
2337 assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min);
2338 assert(frexp(real.infinity, exp) == real.infinity && exp == int.max);
2339 assert(frexp(-0.0, exp) == -0.0 && exp == 0);
2340 assert(frexp(0.0, exp) == 0.0 && exp == 0);
2341 }
2342
2343 @safe @nogc nothrow unittest
2344 {
2345 import std.math.operations : isClose;
2346
2347 int exp;
2348 real mantissa = frexp(123.456L, exp);
2349
2350 // check if values are equal to 19 decimal digits of precision
2351 assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18));
2352 }
2353
2354 @safe unittest
2355 {
2356 import std.math : floatTraits, RealFormat;
2357 import std.math.traits : isIdentical;
2358 import std.meta : AliasSeq;
2359 import std.typecons : tuple, Tuple;
2360
2361 static foreach (T; AliasSeq!(real, double, float))
2362 {{
2363 Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2364 tuple(T(0.0), T( 0.0 ), 0),
2365 tuple(T(-0.0), T( -0.0), 0),
2366 tuple(T(1.0), T( .5 ), 1),
2367 tuple(T(-1.0), T( -.5 ), 1),
2368 tuple(T(2.0), T( .5 ), 2),
2369 tuple(T(float.min_normal/2.0f), T(.5), -126),
2370 tuple(T.infinity, T.infinity, int.max),
2371 tuple(-T.infinity, -T.infinity, int.min),
2372 tuple(T.nan, T.nan, int.min),
2373 tuple(-T.nan, -T.nan, int.min),
2374
2375 // https://issues.dlang.org/show_bug.cgi?id=16026:
2376 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2377 ];
2378
2379 foreach (elem; vals)
2380 {
2381 T x = elem[0];
2382 T e = elem[1];
2383 int exp = elem[2];
2384 int eptr;
2385 T v = frexp(x, eptr);
2386 assert(isIdentical(e, v));
2387 assert(exp == eptr);
2388 }
2389
2390 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2391 {
2392 static T[3][] extendedvals = [ // x,frexp,exp
2393 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2394 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2395 [T.min_normal, .5, -16381],
2396 [T.min_normal/2.0L, .5, -16382] // subnormal
2397 ];
2398 foreach (elem; extendedvals)
2399 {
2400 T x = elem[0];
2401 T e = elem[1];
2402 int exp = cast(int) elem[2];
2403 int eptr;
2404 T v = frexp(x, eptr);
2405 assert(isIdentical(e, v));
2406 assert(exp == eptr);
2407 }
2408 }
2409 }}
2410
2411 // CTFE
2412 alias CtfeFrexpResult= Tuple!(real, int);
2413 static CtfeFrexpResult ctfeFrexp(T)(const T value)
2414 {
2415 int exp;
2416 auto significand = frexp(value, exp);
2417 return CtfeFrexpResult(significand, exp);
2418 }
2419 static foreach (T; AliasSeq!(real, double, float))
2420 {{
2421 enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp
2422 tuple(T(0.0), T( 0.0 ), 0),
2423 tuple(T(-0.0), T( -0.0), 0),
2424 tuple(T(1.0), T( .5 ), 1),
2425 tuple(T(-1.0), T( -.5 ), 1),
2426 tuple(T(2.0), T( .5 ), 2),
2427 tuple(T(float.min_normal/2.0f), T(.5), -126),
2428 tuple(T.infinity, T.infinity, int.max),
2429 tuple(-T.infinity, -T.infinity, int.min),
2430 tuple(T.nan, T.nan, int.min),
2431 tuple(-T.nan, -T.nan, int.min),
2432
2433 // https://issues.dlang.org/show_bug.cgi?id=16026:
2434 tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2)
2435 ];
2436
2437 static foreach (elem; vals)
2438 {
2439 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2]));
2440 }
2441
2442 static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended)
2443 {
2444 enum T[3][] extendedvals = [ // x,frexp,exp
2445 [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
2446 [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
2447 [T.min_normal, .5, -16381],
2448 [T.min_normal/2.0L, .5, -16382] // subnormal
2449 ];
2450 static foreach (elem; extendedvals)
2451 {
2452 static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2]));
2453 }
2454 }
2455 }}
2456 }
2457
2458 @safe unittest
2459 {
2460 import std.meta : AliasSeq;
2461 void foo() {
2462 static foreach (T; AliasSeq!(real, double, float))
2463 {{
2464 int exp;
2465 const T a = 1;
2466 immutable T b = 2;
2467 auto c = frexp(a, exp);
2468 auto d = frexp(b, exp);
2469 }}
2470 }
2471 }
2472
2473 /******************************************
2474 * Extracts the exponent of x as a signed integral value.
2475 *
2476 * If x is not a special value, the result is the same as
2477 * $(D cast(int) logb(x)).
2478 *
2479 * $(TABLE_SV
2480 * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?))
2481 * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes))
2482 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no))
2483 * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no))
2484 * )
2485 */
2486 int ilogb(T)(const T x) @trusted pure nothrow @nogc
2487 if (isFloatingPoint!T)
2488 {
2489 import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB;
2490
2491 import core.bitop : bsr;
2492 alias F = floatTraits!T;
2493
2494 union floatBits
2495 {
2496 T rv;
2497 ushort[T.sizeof/2] vu;
2498 uint[T.sizeof/4] vui;
2499 static if (T.sizeof >= 8)
2500 ulong[T.sizeof/8] vul;
2501 }
2502 floatBits y = void;
2503 y.rv = x;
2504
2505 int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
2506 static if (F.realFormat == RealFormat.ieeeExtended ||
2507 F.realFormat == RealFormat.ieeeExtended53)
2508 {
2509 if (ex)
2510 {
2511 // If exponent is non-zero
2512 if (ex == F.EXPMASK) // infinity or NaN
2513 {
2514 if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN
2515 return FP_ILOGBNAN;
2516 else // +-infinity
2517 return int.max;
2518 }
2519 else
2520 {
2521 return ex - F.EXPBIAS - 1;
2522 }
2523 }
2524 else if (!y.vul[0])
2525 {
2526 // vf is +-0.0
2527 return FP_ILOGB0;
2528 }
2529 else
2530 {
2531 // subnormal
2532 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]);
2533 }
2534 }
2535 else static if (F.realFormat == RealFormat.ieeeQuadruple)
2536 {
2537 if (ex) // If exponent is non-zero
2538 {
2539 if (ex == F.EXPMASK)
2540 {
2541 // infinity or NaN
2542 if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN
2543 return FP_ILOGBNAN;
2544 else // +- infinity
2545 return int.max;
2546 }
2547 else
2548 {
2549 return ex - F.EXPBIAS - 1;
2550 }
2551 }
2552 else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0)
2553 {
2554 // vf is +-0.0
2555 return FP_ILOGB0;
2556 }
2557 else
2558 {
2559 // subnormal
2560 const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF;
2561 const ulong lsb = y.vul[MANTISSA_LSB];
2562 if (msb)
2563 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64;
2564 else
2565 return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb);
2566 }
2567 }
2568 else static if (F.realFormat == RealFormat.ieeeDouble)
2569 {
2570 if (ex) // If exponent is non-zero
2571 {
2572 if (ex == F.EXPMASK) // infinity or NaN
2573 {
2574 if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity
2575 return int.max;
2576 else // NaN
2577 return FP_ILOGBNAN;
2578 }
2579 else
2580 {
2581 return ((ex - F.EXPBIAS) >> 4) - 1;
2582 }
2583 }
2584 else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF))
2585 {
2586 // vf is +-0.0
2587 return FP_ILOGB0;
2588 }
2589 else
2590 {
2591 // subnormal
2592 enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF;
2593 return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64);
2594 }
2595 }
2596 else static if (F.realFormat == RealFormat.ieeeSingle)
2597 {
2598 if (ex) // If exponent is non-zero
2599 {
2600 if (ex == F.EXPMASK) // infinity or NaN
2601 {
2602 if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity
2603 return int.max;
2604 else // NaN
2605 return FP_ILOGBNAN;
2606 }
2607 else
2608 {
2609 return ((ex - F.EXPBIAS) >> 7) - 1;
2610 }
2611 }
2612 else if (!(y.vui[0] & 0x7FFF_FFFF))
2613 {
2614 // vf is +-0.0
2615 return FP_ILOGB0;
2616 }
2617 else
2618 {
2619 // subnormal
2620 const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT;
2621 return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa);
2622 }
2623 }
2624 else // static if (F.realFormat == RealFormat.ibmExtended)
2625 {
2626 assert(0, "ilogb not implemented");
2627 }
2628 }
2629 /// ditto
2630 int ilogb(T)(const T x) @safe pure nothrow @nogc
2631 if (isIntegral!T && isUnsigned!T)
2632 {
2633 import core.bitop : bsr;
2634 if (x == 0)
2635 return FP_ILOGB0;
2636 else
2637 {
2638 static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation");
2639 return bsr(x);
2640 }
2641 }
2642 /// ditto
2643 int ilogb(T)(const T x) @safe pure nothrow @nogc
2644 if (isIntegral!T && isSigned!T)
2645 {
2646 import std.traits : Unsigned;
2647 // Note: abs(x) can not be used because the return type is not Unsigned and
2648 // the return value would be wrong for x == int.min
2649 Unsigned!T absx = x >= 0 ? x : -x;
2650 return ilogb(absx);
2651 }
2652
2653 ///
2654 @safe pure unittest
2655 {
2656 assert(ilogb(1) == 0);
2657 assert(ilogb(3) == 1);
2658 assert(ilogb(3.0) == 1);
2659 assert(ilogb(100_000_000) == 26);
2660
2661 assert(ilogb(0) == FP_ILOGB0);
2662 assert(ilogb(0.0) == FP_ILOGB0);
2663 assert(ilogb(double.nan) == FP_ILOGBNAN);
2664 assert(ilogb(double.infinity) == int.max);
2665 }
2666
2667 /**
2668 Special return values of $(LREF ilogb).
2669 */
2670 alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0;
2671 /// ditto
2672 alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN;
2673
2674 ///
2675 @safe pure unittest
2676 {
2677 assert(ilogb(0) == FP_ILOGB0);
2678 assert(ilogb(0.0) == FP_ILOGB0);
2679 assert(ilogb(double.nan) == FP_ILOGBNAN);
2680 }
2681
2682 @safe nothrow @nogc unittest
2683 {
2684 import std.math : floatTraits, RealFormat;
2685 import std.math.operations : nextUp;
2686 import std.meta : AliasSeq;
2687 import std.typecons : Tuple;
2688 static foreach (F; AliasSeq!(float, double, real))
2689 {{
2690 alias T = Tuple!(F, int);
2691 T[13] vals = // x, ilogb(x)
2692 [
2693 T( F.nan , FP_ILOGBNAN ),
2694 T( -F.nan , FP_ILOGBNAN ),
2695 T( F.infinity, int.max ),
2696 T( -F.infinity, int.max ),
2697 T( 0.0 , FP_ILOGB0 ),
2698 T( -0.0 , FP_ILOGB0 ),
2699 T( 2.0 , 1 ),
2700 T( 2.0001 , 1 ),
2701 T( 1.9999 , 0 ),
2702 T( 0.5 , -1 ),
2703 T( 123.123 , 6 ),
2704 T( -123.123 , 6 ),
2705 T( 0.123 , -4 ),
2706 ];
2707
2708 foreach (elem; vals)
2709 {
2710 assert(ilogb(elem[0]) == elem[1]);
2711 }
2712 }}
2713
2714 // min_normal and subnormals
2715 assert(ilogb(-float.min_normal) == -126);
2716 assert(ilogb(nextUp(-float.min_normal)) == -127);
2717 assert(ilogb(nextUp(-float(0.0))) == -149);
2718 assert(ilogb(-double.min_normal) == -1022);
2719 assert(ilogb(nextUp(-double.min_normal)) == -1023);
2720 assert(ilogb(nextUp(-double(0.0))) == -1074);
2721 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
2722 {
2723 assert(ilogb(-real.min_normal) == -16382);
2724 assert(ilogb(nextUp(-real.min_normal)) == -16383);
2725 assert(ilogb(nextUp(-real(0.0))) == -16445);
2726 }
2727 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2728 {
2729 assert(ilogb(-real.min_normal) == -1022);
2730 assert(ilogb(nextUp(-real.min_normal)) == -1023);
2731 assert(ilogb(nextUp(-real(0.0))) == -1074);
2732 }
2733
2734 // test integer types
2735 assert(ilogb(0) == FP_ILOGB0);
2736 assert(ilogb(int.max) == 30);
2737 assert(ilogb(int.min) == 31);
2738 assert(ilogb(uint.max) == 31);
2739 assert(ilogb(long.max) == 62);
2740 assert(ilogb(long.min) == 63);
2741 assert(ilogb(ulong.max) == 63);
2742 }
2743
2744 /*******************************************
2745 * Compute n * 2$(SUPERSCRIPT exp)
2746 * References: frexp
2747 */
2748
2749 pragma(inline, true)
2750 real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2751 ///ditto
2752 pragma(inline, true)
2753 double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2754 ///ditto
2755 pragma(inline, true)
2756 float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); }
2757
2758 ///
2759 @nogc @safe pure nothrow unittest
2760 {
2761 import std.meta : AliasSeq;
2762 static foreach (T; AliasSeq!(float, double, real))
2763 {{
2764 T r;
2765
2766 r = ldexp(3.0L, 3);
2767 assert(r == 24);
2768
2769 r = ldexp(cast(T) 3.0, cast(int) 3);
2770 assert(r == 24);
2771
2772 T n = 3.0;
2773 int exp = 3;
2774 r = ldexp(n, exp);
2775 assert(r == 24);
2776 }}
2777 }
2778
2779 @safe pure nothrow @nogc unittest
2780 {
2781 import std.math : floatTraits, RealFormat;
2782
2783 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
2784 floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
2785 floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
2786 {
2787 assert(ldexp(1.0L, -16384) == 0x1p-16384L);
2788 assert(ldexp(1.0L, -16382) == 0x1p-16382L);
2789 int x;
2790 real n = frexp(0x1p-16384L, x);
2791 assert(n == 0.5L);
2792 assert(x==-16383);
2793 assert(ldexp(n, x)==0x1p-16384L);
2794 }
2795 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
2796 {
2797 assert(ldexp(1.0L, -1024) == 0x1p-1024L);
2798 assert(ldexp(1.0L, -1022) == 0x1p-1022L);
2799 int x;
2800 real n = frexp(0x1p-1024L, x);
2801 assert(n == 0.5L);
2802 assert(x==-1023);
2803 assert(ldexp(n, x)==0x1p-1024L);
2804 }
2805 else static assert(false, "Floating point type real not supported");
2806 }
2807
2808 /* workaround https://issues.dlang.org/show_bug.cgi?id=14718
2809 float parsing depends on platform strtold
2810 @safe pure nothrow @nogc unittest
2811 {
2812 assert(ldexp(1.0, -1024) == 0x1p-1024);
2813 assert(ldexp(1.0, -1022) == 0x1p-1022);
2814 int x;
2815 double n = frexp(0x1p-1024, x);
2816 assert(n == 0.5);
2817 assert(x==-1023);
2818 assert(ldexp(n, x)==0x1p-1024);
2819 }
2820
2821 @safe pure nothrow @nogc unittest
2822 {
2823 assert(ldexp(1.0f, -128) == 0x1p-128f);
2824 assert(ldexp(1.0f, -126) == 0x1p-126f);
2825 int x;
2826 float n = frexp(0x1p-128f, x);
2827 assert(n == 0.5f);
2828 assert(x==-127);
2829 assert(ldexp(n, x)==0x1p-128f);
2830 }
2831 */
2832
2833 @safe @nogc nothrow unittest
2834 {
2835 import std.math.operations : isClose;
2836
2837 static real[3][] vals = // value,exp,ldexp
2838 [
2839 [ 0, 0, 0],
2840 [ 1, 0, 1],
2841 [ -1, 0, -1],
2842 [ 1, 1, 2],
2843 [ 123, 10, 125952],
2844 [ real.max, int.max, real.infinity],
2845 [ real.max, -int.max, 0],
2846 [ real.min_normal, -int.max, 0],
2847 ];
2848 int i;
2849
2850 for (i = 0; i < vals.length; i++)
2851 {
2852 real x = vals[i][0];
2853 int exp = cast(int) vals[i][1];
2854 real z = vals[i][2];
2855 real l = ldexp(x, exp);
2856
2857 assert(isClose(z, l, 1e-6));
2858 }
2859
2860 real function(real, int) pldexp = &ldexp;
2861 assert(pldexp != null);
2862 }
2863
2864 private
2865 {
2866 // Coefficients shared across log(), log2(), log10(), log1p().
2867 template LogCoeffs(T)
2868 {
2869 import std.math : floatTraits, RealFormat;
2870
2871 static if (floatTraits!T.realFormat == RealFormat.ieeeQuadruple)
2872 {
2873 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2874 // Theoretical peak relative error = 5.3e-37
2875 static immutable real[13] logP = [
2876 1.313572404063446165910279910527789794488E4L,
2877 7.771154681358524243729929227226708890930E4L,
2878 2.014652742082537582487669938141683759923E5L,
2879 3.007007295140399532324943111654767187848E5L,
2880 2.854829159639697837788887080758954924001E5L,
2881 1.797628303815655343403735250238293741397E5L,
2882 7.594356839258970405033155585486712125861E4L,
2883 2.128857716871515081352991964243375186031E4L,
2884 3.824952356185897735160588078446136783779E3L,
2885 4.114517881637811823002128927449878962058E2L,
2886 2.321125933898420063925789532045674660756E1L,
2887 4.998469661968096229986658302195402690910E-1L,
2888 1.538612243596254322971797716843006400388E-6L
2889 ];
2890 static immutable real[13] logQ = [
2891 3.940717212190338497730839731583397586124E4L,
2892 2.626900195321832660448791748036714883242E5L,
2893 7.777690340007566932935753241556479363645E5L,
2894 1.347518538384329112529391120390701166528E6L,
2895 1.514882452993549494932585972882995548426E6L,
2896 1.158019977462989115839826904108208787040E6L,
2897 6.132189329546557743179177159925690841200E5L,
2898 2.248234257620569139969141618556349415120E5L,
2899 5.605842085972455027590989944010492125825E4L,
2900 9.147150349299596453976674231612674085381E3L,
2901 9.104928120962988414618126155557301584078E2L,
2902 4.839208193348159620282142911143429644326E1L,
2903 1.0
2904 ];
2905
2906 // log2 uses the same coefficients as log.
2907 alias log2P = logP;
2908 alias log2Q = logQ;
2909
2910 // log10 uses the same coefficients as log.
2911 alias log10P = logP;
2912 alias log10Q = logQ;
2913
2914 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
2915 // where z = 2(x-1)/(x+1)
2916 // Theoretical peak relative error = 1.1e-35
2917 static immutable real[6] logR = [
2918 1.418134209872192732479751274970992665513E5L,
2919 -8.977257995689735303686582344659576526998E4L,
2920 2.048819892795278657810231591630928516206E4L,
2921 -2.024301798136027039250415126250455056397E3L,
2922 8.057002716646055371965756206836056074715E1L,
2923 -8.828896441624934385266096344596648080902E-1L
2924 ];
2925 static immutable real[7] logS = [
2926 1.701761051846631278975701529965589676574E6L,
2927 -1.332535117259762928288745111081235577029E6L,
2928 4.001557694070773974936904547424676279307E5L,
2929 -5.748542087379434595104154610899551484314E4L,
2930 3.998526750980007367835804959888064681098E3L,
2931 -1.186359407982897997337150403816839480438E2L,
2932 1.0
2933 ];
2934 }
2935 else static if (floatTraits!T.realFormat == RealFormat.ieeeExtended ||
2936 floatTraits!T.realFormat == RealFormat.ieeeExtended53)
2937 {
2938 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2939 // Theoretical peak relative error = 2.32e-20
2940 static immutable real[7] logP = [
2941 2.0039553499201281259648E1L,
2942 5.7112963590585538103336E1L,
2943 6.0949667980987787057556E1L,
2944 2.9911919328553073277375E1L,
2945 6.5787325942061044846969E0L,
2946 4.9854102823193375972212E-1L,
2947 4.5270000862445199635215E-5L,
2948 ];
2949 static immutable real[7] logQ = [
2950 6.0118660497603843919306E1L,
2951 2.1642788614495947685003E2L,
2952 3.0909872225312059774938E2L,
2953 2.2176239823732856465394E2L,
2954 8.3047565967967209469434E1L,
2955 1.5062909083469192043167E1L,
2956 1.0000000000000000000000E0L,
2957 ];
2958
2959 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
2960 // Theoretical peak relative error = 6.2e-22
2961 static immutable real[7] log2P = [
2962 1.0747524399916215149070E2L,
2963 3.4258224542413922935104E2L,
2964 4.2401812743503691187826E2L,
2965 2.5620629828144409632571E2L,
2966 7.7671073698359539859595E1L,
2967 1.0767376367209449010438E1L,
2968 4.9962495940332550844739E-1L,
2969 ];
2970 static immutable real[8] log2Q = [
2971 3.2242573199748645407652E2L,
2972 1.2695660352705325274404E3L,
2973 2.0307734695595183428202E3L,
2974 1.6911722418503949084863E3L,
2975 7.7952888181207260646090E2L,
2976 1.9444210022760132894510E2L,
2977 2.3479774160285863271658E1L,
2978 1.0000000000000000000000E0,
2979 ];
2980
2981 // log10 uses the same coefficients as log2.
2982 alias log10P = log2P;
2983 alias log10Q = log2Q;
2984
2985 // Coefficients for log(x) = z + z^^3 R(z^^2)/S(z^^2)
2986 // where z = 2(x-1)/(x+1)
2987 // Theoretical peak relative error = 6.16e-22
2988 static immutable real[4] logR = [
2989 -3.5717684488096787370998E1L,
2990 1.0777257190312272158094E1L,
2991 -7.1990767473014147232598E-1L,
2992 1.9757429581415468984296E-3L,
2993 ];
2994 static immutable real[4] logS = [
2995 -4.2861221385716144629696E2L,
2996 1.9361891836232102174846E2L,
2997 -2.6201045551331104417768E1L,
2998 1.0000000000000000000000E0L,
2999 ];
3000 }
3001 else static if (floatTraits!T.realFormat == RealFormat.ieeeDouble)
3002 {
3003 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3004 static immutable double[6] logP = [
3005 7.70838733755885391666E0,
3006 1.79368678507819816313E1,
3007 1.44989225341610930846E1,
3008 4.70579119878881725854E0,
3009 4.97494994976747001425E-1,
3010 1.01875663804580931796E-4,
3011 ];
3012 static immutable double[6] logQ = [
3013 2.31251620126765340583E1,
3014 7.11544750618563894466E1,
3015 8.29875266912776603211E1,
3016 4.52279145837532221105E1,
3017 1.12873587189167450590E1,
3018 1.00000000000000000000E0,
3019 ];
3020
3021 // log2 uses the same coefficients as log.
3022 alias log2P = logP;
3023 alias log2Q = logQ;
3024
3025 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3026 static immutable double[7] logp1P = [
3027 2.0039553499201281259648E1,
3028 5.7112963590585538103336E1,
3029 6.0949667980987787057556E1,
3030 2.9911919328553073277375E1,
3031 6.5787325942061044846969E0,
3032 4.9854102823193375972212E-1,
3033 4.5270000862445199635215E-5,
3034 ];
3035 static immutable double[7] logp1Q = [
3036 1.0000000000000000000000E0,
3037 6.0118660497603843919306E1,
3038 2.1642788614495947685003E2,
3039 3.0909872225312059774938E2,
3040 2.2176239823732856465394E2,
3041 8.3047565967967209469434E1,
3042 1.5062909083469192043167E1,
3043 ];
3044
3045 static immutable double[7] log10P = [
3046 1.98892446572874072159E1,
3047 5.67349287391754285487E1,
3048 6.06127134467767258030E1,
3049 2.97877425097986925891E1,
3050 6.56312093769992875930E0,
3051 4.98531067254050724270E-1,
3052 4.58482948458143443514E-5,
3053 ];
3054 static immutable double[7] log10Q = [
3055 5.96677339718622216300E1,
3056 2.14955586696422947765E2,
3057 3.07254189979530058263E2,
3058 2.20664384982121929218E2,
3059 8.27410449222435217021E1,
3060 1.50314182634250003249E1,
3061 1.00000000000000000000E0,
3062 ];
3063
3064 // Coefficients for log(x) = z + z^^3 R(z)/S(z)
3065 // where z = 2(x-1)/(x+1)
3066 static immutable double[3] logR = [
3067 -6.41409952958715622951E1,
3068 1.63866645699558079767E1,
3069 -7.89580278884799154124E-1,
3070 ];
3071 static immutable double[4] logS = [
3072 -7.69691943550460008604E2,
3073 3.12093766372244180303E2,
3074 -3.56722798256324312549E1,
3075 1.00000000000000000000E0,
3076 ];
3077 }
3078 else static if (floatTraits!T.realFormat == RealFormat.ieeeSingle)
3079 {
3080 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)
3081 static immutable float[9] logP = [
3082 3.3333331174E-1,
3083 -2.4999993993E-1,
3084 2.0000714765E-1,
3085 -1.6668057665E-1,
3086 1.4249322787E-1,
3087 -1.2420140846E-1,
3088 1.1676998740E-1,
3089 -1.1514610310E-1,
3090 7.0376836292E-2,
3091 ];
3092
3093 // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
3094 static immutable float[7] logp1P = [
3095 2.0039553499E1,
3096 5.7112963590E1,
3097 6.0949667980E1,
3098 2.9911919328E1,
3099 6.5787325942E0,
3100 4.9854102823E-1,
3101 4.5270000862E-5,
3102 ];
3103 static immutable float[7] logp1Q = [
3104 1.00000000000E0,
3105 6.01186604976E1,
3106 2.16427886144E2,
3107 3.09098722253E2,
3108 2.21762398237E2,
3109 8.30475659679E1,
3110 1.50629090834E1,
3111 ];
3112
3113 // log2 and log10 uses the same coefficients as log.
3114 alias log2P = logP;
3115 alias log10P = logP;
3116 }
3117 else
3118 static assert(0, "no coefficients for log()");
3119 }
3120 }
3121
3122 /**************************************
3123 * Calculate the natural logarithm of x.
3124 *
3125 * $(TABLE_SV
3126 * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?))
3127 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3128 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
3129 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3130 * )
3131 */
3132 pragma(inline, true)
3133 real log(real x) @safe pure nothrow @nogc
3134 {
3135 version (INLINE_YL2X)
3136 {
3137 import std.math.constants : LN2;
3138 return core.math.yl2x(x, LN2);
3139 }
3140 else
3141 return logImpl(x);
3142 }
3143
3144 /// ditto
3145 pragma(inline, true)
3146 double log(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log(cast(real) x) : logImpl(x); }
3147
3148 /// ditto
3149 pragma(inline, true)
3150 float log(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log(cast(real) x) : logImpl(x); }
3151
3152 // @@@DEPRECATED_[2.112.0]@@@
3153 deprecated("`std.math.exponential.log` called with argument types `(int)` matches both "
3154 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3155 real log(int x) @safe pure nothrow @nogc { return log(cast(real) x); }
3156 // @@@DEPRECATED_[2.112.0]@@@
3157 deprecated("`std.math.exponential.log` called with argument types `(uint)` matches both "
3158 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3159 real log(uint x) @safe pure nothrow @nogc { return log(cast(real) x); }
3160 // @@@DEPRECATED_[2.112.0]@@@
3161 deprecated("`std.math.exponential.log` called with argument types `(long)` matches both "
3162 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3163 real log(long x) @safe pure nothrow @nogc { return log(cast(real) x); }
3164 // @@@DEPRECATED_[2.112.0]@@@
3165 deprecated("`std.math.exponential.log` called with argument types `(ulong)` matches both "
3166 ~ "`log(real)`, `log(double)`, and `log(float)`. Cast argument to floating point type instead.")
3167 real log(ulong x) @safe pure nothrow @nogc { return log(cast(real) x); }
3168
3169 ///
3170 @safe pure nothrow @nogc unittest
3171 {
3172 import std.math.operations : feqrel;
3173 import std.math.constants : E;
3174
3175 assert(feqrel(log(E), 1) >= real.mant_dig - 1);
3176 }
3177
3178 private T logImpl(T, bool LOG1P = false)(T x) @safe pure nothrow @nogc
3179 {
3180 import std.math.constants : SQRT1_2;
3181 import std.math.algebraic : poly;
3182 import std.math.traits : isInfinity, isNaN, signbit;
3183 import std.math : floatTraits, RealFormat;
3184
3185 alias coeffs = LogCoeffs!T;
3186 alias F = floatTraits!T;
3187
3188 static if (LOG1P)
3189 {
3190 const T xm1 = x;
3191 x = x + 1.0;
3192 }
3193
3194 static if (F.realFormat == RealFormat.ieeeExtended ||
3195 F.realFormat == RealFormat.ieeeExtended53 ||
3196 F.realFormat == RealFormat.ieeeQuadruple)
3197 {
3198 // C1 + C2 = LN2.
3199 enum T C1 = 6.93145751953125E-1L;
3200 enum T C2 = 1.428606820309417232121458176568075500134E-6L;
3201 }
3202 else static if (F.realFormat == RealFormat.ieeeDouble)
3203 {
3204 enum T C1 = 0.693359375;
3205 enum T C2 = -2.121944400546905827679e-4;
3206 }
3207 else static if (F.realFormat == RealFormat.ieeeSingle)
3208 {
3209 enum T C1 = 0.693359375;
3210 enum T C2 = -2.12194440e-4;
3211 }
3212 else
3213 static assert(0, "Not implemented for this architecture");
3214
3215 // Special cases.
3216 if (isNaN(x))
3217 return x;
3218 if (isInfinity(x) && !signbit(x))
3219 return x;
3220 if (x == 0.0)
3221 return -T.infinity;
3222 if (x < 0.0)
3223 return T.nan;
3224
3225 // Separate mantissa from exponent.
3226 // Note, frexp is used so that denormal numbers will be handled properly.
3227 T y, z;
3228 int exp;
3229
3230 x = frexp(x, exp);
3231
3232 static if (F.realFormat == RealFormat.ieeeDouble ||
3233 F.realFormat == RealFormat.ieeeExtended ||
3234 F.realFormat == RealFormat.ieeeExtended53 ||
3235 F.realFormat == RealFormat.ieeeQuadruple)
3236 {
3237 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3238 // where z = 2(x - 1)/(x + 1)
3239 if ((exp > 2) || (exp < -2))
3240 {
3241 if (x < SQRT1_2)
3242 { // 2(2x - 1)/(2x + 1)
3243 exp -= 1;
3244 z = x - 0.5;
3245 y = 0.5 * z + 0.5;
3246 }
3247 else
3248 { // 2(x - 1)/(x + 1)
3249 z = x - 0.5;
3250 z -= 0.5;
3251 y = 0.5 * x + 0.5;
3252 }
3253 x = z / y;
3254 z = x * x;
3255 z = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3256 z += exp * C2;
3257 z += x;
3258 z += exp * C1;
3259
3260 return z;
3261 }
3262 }
3263
3264 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3265 if (x < SQRT1_2)
3266 {
3267 exp -= 1;
3268 static if (LOG1P)
3269 {
3270 if (exp != 0)
3271 x = 2.0 * x - 1.0;
3272 else
3273 x = xm1;
3274 }
3275 else
3276 x = 2.0 * x - 1.0;
3277
3278 }
3279 else
3280 {
3281 static if (LOG1P)
3282 {
3283 if (exp != 0)
3284 x = x - 1.0;
3285 else
3286 x = xm1;
3287 }
3288 else
3289 x = x - 1.0;
3290 }
3291 z = x * x;
3292 static if (F.realFormat == RealFormat.ieeeSingle)
3293 y = x * (z * poly(x, coeffs.logP));
3294 else
3295 y = x * (z * poly(x, coeffs.logP) / poly(x, coeffs.logQ));
3296 y += exp * C2;
3297 z = y - 0.5 * z;
3298
3299 // Note, the sum of above terms does not exceed x/4,
3300 // so it contributes at most about 1/4 lsb to the error.
3301 z += x;
3302 z += exp * C1;
3303
3304 return z;
3305 }
3306
3307 @safe @nogc nothrow unittest
3308 {
3309 import std.math : floatTraits, RealFormat;
3310 import std.meta : AliasSeq;
3311
3312 static void testLog(T)(T[2][] vals)
3313 {
3314 import std.math.operations : isClose;
3315 import std.math.traits : isNaN;
3316 foreach (ref pair; vals)
3317 {
3318 if (isNaN(pair[1]))
3319 assert(isNaN(log(pair[0])));
3320 else
3321 assert(isClose(log(pair[0]), pair[1]));
3322 }
3323 }
3324 static foreach (F; AliasSeq!(float, double, real))
3325 {{
3326 F[2][24] vals = [
3327 [F(1), F(0x0p+0)], [F(2), F(0x1.62e42fefa39ef358p-1)],
3328 [F(4), F(0x1.62e42fefa39ef358p+0)], [F(8), F(0x1.0a2b23f3bab73682p+1)],
3329 [F(16), F(0x1.62e42fefa39ef358p+1)], [F(32), F(0x1.bb9d3beb8c86b02ep+1)],
3330 [F(64), F(0x1.0a2b23f3bab73682p+2)], [F(128), F(0x1.3687a9f1af2b14ecp+2)],
3331 [F(256), F(0x1.62e42fefa39ef358p+2)], [F(512), F(0x1.8f40b5ed9812d1c2p+2)],
3332 [F(1024), F(0x1.bb9d3beb8c86b02ep+2)], [F(2048), F(0x1.e7f9c1e980fa8e98p+2)],
3333 [F(3), F(0x1.193ea7aad030a976p+0)], [F(5), F(0x1.9c041f7ed8d336bp+0)],
3334 [F(7), F(0x1.f2272ae325a57546p+0)], [F(15), F(0x1.5aa16394d481f014p+1)],
3335 [F(17), F(0x1.6aa6bc1fa7f79cfp+1)], [F(31), F(0x1.b78ce48912b59f12p+1)],
3336 [F(33), F(0x1.bf8d8f4d5b8d1038p+1)], [F(63), F(0x1.09291e8e3181b20ep+2)],
3337 [F(65), F(0x1.0b292939429755ap+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
3338 [F(10000), F(0x1.26bb1bbb5551582ep+3)],
3339 ];
3340 testLog(vals);
3341 }}
3342 {
3343 float[2][16] vals = [
3344 [float.nan, float.nan],[-float.nan, float.nan],
3345 [float.infinity, float.infinity], [-float.infinity, float.nan],
3346 [float.min_normal, -0x1.5d58ap+6f], [-float.min_normal, float.nan],
3347 [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3348 [float.min_normal / 2, -0x1.601e68p+6f], [-float.min_normal / 2, float.nan],
3349 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3350 [float.min_normal / 3, -0x1.61bd9ap+6f], [-float.min_normal / 3, float.nan],
3351 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3352 ];
3353 testLog(vals);
3354 }
3355 {
3356 double[2][16] vals = [
3357 [double.nan, double.nan],[-double.nan, double.nan],
3358 [double.infinity, double.infinity], [-double.infinity, double.nan],
3359 [double.min_normal, -0x1.6232bdd7abcd2p+9], [-double.min_normal, double.nan],
3360 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3361 [double.min_normal / 2, -0x1.628b76e3a7b61p+9], [-double.min_normal / 2, double.nan],
3362 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3363 [double.min_normal / 3, -0x1.62bf5d2b81354p+9], [-double.min_normal / 3, double.nan],
3364 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3365 ];
3366 testLog(vals);
3367 }
3368 alias F = floatTraits!real;
3369 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3370 {{
3371 real[2][16] vals = [
3372 [real.nan, real.nan],[-real.nan, real.nan],
3373 [real.infinity, real.infinity], [-real.infinity, real.nan],
3374 [real.min_normal, -0x1.62d918ce2421d66p+13L], [-real.min_normal, real.nan],
3375 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3376 [real.min_normal / 2, -0x1.62dea45ee3e064dcp+13L], [-real.min_normal / 2, real.nan],
3377 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3378 [real.min_normal / 3, -0x1.62e1e2c3617857e6p+13L], [-real.min_normal / 3, real.nan],
3379 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3380 ];
3381 testLog(vals);
3382 }}
3383 }
3384
3385 /**************************************
3386 * Calculate the base-10 logarithm of x.
3387 *
3388 * $(TABLE_SV
3389 * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?))
3390 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3391 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
3392 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3393 * )
3394 */
3395 pragma(inline, true)
3396 real log10(real x) @safe pure nothrow @nogc
3397 {
3398 version (INLINE_YL2X)
3399 {
3400 import std.math.constants : LOG2;
3401 return core.math.yl2x(x, LOG2);
3402 }
3403 else
3404 return log10Impl(x);
3405 }
3406
3407 /// ditto
3408 pragma(inline, true)
3409 double log10(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log10(cast(real) x) : log10Impl(x); }
3410
3411 /// ditto
3412 pragma(inline, true)
3413 float log10(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log10(cast(real) x) : log10Impl(x); }
3414
3415 // @@@DEPRECATED_[2.112.0]@@@
3416 deprecated("`std.math.exponential.log10` called with argument types `(int)` matches both "
3417 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3418 real log10(int x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3419 // @@@DEPRECATED_[2.112.0]@@@
3420 deprecated("`std.math.exponential.log10` called with argument types `(uint)` matches both "
3421 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3422 real log10(uint x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3423 // @@@DEPRECATED_[2.112.0]@@@
3424 deprecated("`std.math.exponential.log10` called with argument types `(long)` matches both "
3425 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3426 real log10(long x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3427 // @@@DEPRECATED_[2.112.0]@@@
3428 deprecated("`std.math.exponential.log10` called with argument types `(ulong)` matches both "
3429 ~ "`log10(real)`, `log10(double)`, and `log10(float)`. Cast argument to floating point type instead.")
3430 real log10(ulong x) @safe pure nothrow @nogc { return log10(cast(real) x); }
3431
3432 ///
3433 @safe pure nothrow @nogc unittest
3434 {
3435 import std.math.algebraic : fabs;
3436
3437 assert(fabs(log10(1000.0L) - 3) < .000001);
3438 }
3439
3440 @safe pure nothrow @nogc unittest
3441 {
3442 import std.math.algebraic : fabs;
3443
3444 assert(fabs(log10(1000.0) - 3) < .000001);
3445 assert(fabs(log10(1000.0f) - 3) < .000001);
3446 }
3447
3448 private T log10Impl(T)(T x) @safe pure nothrow @nogc
3449 {
3450 import std.math.constants : SQRT1_2;
3451 import std.math.algebraic : poly;
3452 import std.math.traits : isNaN, isInfinity, signbit;
3453 import std.math : floatTraits, RealFormat;
3454
3455 alias coeffs = LogCoeffs!T;
3456 alias F = floatTraits!T;
3457
3458 static if (F.realFormat == RealFormat.ieeeExtended ||
3459 F.realFormat == RealFormat.ieeeExtended53 ||
3460 F.realFormat == RealFormat.ieeeQuadruple)
3461 {
3462 // log10(2) split into two parts.
3463 enum T L102A = 0.3125L;
3464 enum T L102B = -1.14700043360188047862611052755069732318101185E-2L;
3465
3466 // log10(e) split into two parts.
3467 enum T L10EA = 0.5L;
3468 enum T L10EB = -6.570551809674817234887108108339491770560299E-2L;
3469 }
3470 else static if (F.realFormat == RealFormat.ieeeDouble ||
3471 F.realFormat == RealFormat.ieeeSingle)
3472 {
3473 enum T L102A = 3.0078125E-1;
3474 enum T L102B = 2.48745663981195213739E-4;
3475
3476 enum T L10EA = 4.3359375E-1;
3477 enum T L10EB = 7.00731903251827651129E-4;
3478 }
3479 else
3480 static assert(0, "Not implemented for this architecture");
3481
3482 // Special cases are the same as for log.
3483 if (isNaN(x))
3484 return x;
3485 if (isInfinity(x) && !signbit(x))
3486 return x;
3487 if (x == 0.0)
3488 return -T.infinity;
3489 if (x < 0.0)
3490 return T.nan;
3491
3492 // Separate mantissa from exponent.
3493 // Note, frexp is used so that denormal numbers will be handled properly.
3494 T y, z;
3495 int exp;
3496
3497 x = frexp(x, exp);
3498
3499 static if (F.realFormat == RealFormat.ieeeExtended ||
3500 F.realFormat == RealFormat.ieeeExtended53 ||
3501 F.realFormat == RealFormat.ieeeQuadruple)
3502 {
3503 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3504 // where z = 2(x - 1)/(x + 1)
3505 if ((exp > 2) || (exp < -2))
3506 {
3507 if (x < SQRT1_2)
3508 { // 2(2x - 1)/(2x + 1)
3509 exp -= 1;
3510 z = x - 0.5;
3511 y = 0.5 * z + 0.5;
3512 }
3513 else
3514 { // 2(x - 1)/(x + 1)
3515 z = x - 0.5;
3516 z -= 0.5;
3517 y = 0.5 * x + 0.5;
3518 }
3519 x = z / y;
3520 z = x * x;
3521 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3522 goto Ldone;
3523 }
3524 }
3525
3526 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3527 if (x < SQRT1_2)
3528 {
3529 exp -= 1;
3530 x = 2.0 * x - 1.0;
3531 }
3532 else
3533 x = x - 1.0;
3534
3535 z = x * x;
3536 static if (F.realFormat == RealFormat.ieeeSingle)
3537 y = x * (z * poly(x, coeffs.log10P));
3538 else
3539 y = x * (z * poly(x, coeffs.log10P) / poly(x, coeffs.log10Q));
3540 y = y - 0.5 * z;
3541
3542 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3543 // This sequence of operations is critical and it may be horribly
3544 // defeated by some compiler optimizers.
3545 Ldone:
3546 z = y * L10EB;
3547 z += x * L10EB;
3548 z += exp * L102B;
3549 z += y * L10EA;
3550 z += x * L10EA;
3551 z += exp * L102A;
3552
3553 return z;
3554 }
3555
3556 @safe @nogc nothrow unittest
3557 {
3558 import std.math : floatTraits, RealFormat;
3559 import std.meta : AliasSeq;
3560
3561 static void testLog10(T)(T[2][] vals)
3562 {
3563 import std.math.operations : isClose;
3564 import std.math.traits : isNaN;
3565 foreach (ref pair; vals)
3566 {
3567 if (isNaN(pair[1]))
3568 assert(isNaN(log10(pair[0])));
3569 else
3570 assert(isClose(log10(pair[0]), pair[1]));
3571 }
3572 }
3573 static foreach (F; AliasSeq!(float, double, real))
3574 {{
3575 F[2][24] vals = [
3576 [F(1), F(0x0p+0)], [F(2), F(0x1.34413509f79fef32p-2)],
3577 [F(4), F(0x1.34413509f79fef32p-1)], [F(8), F(0x1.ce61cf8ef36fe6cap-1)],
3578 [F(16), F(0x1.34413509f79fef32p+0)], [F(32), F(0x1.8151824c7587eafep+0)],
3579 [F(64), F(0x1.ce61cf8ef36fe6cap+0)], [F(128), F(0x1.0db90e68b8abf14cp+1)],
3580 [F(256), F(0x1.34413509f79fef32p+1)], [F(512), F(0x1.5ac95bab3693ed18p+1)],
3581 [F(1024), F(0x1.8151824c7587eafep+1)], [F(2048), F(0x1.a7d9a8edb47be8e4p+1)],
3582 [F(3), F(0x1.e8927964fd5fd08cp-2)], [F(5), F(0x1.65df657b04300868p-1)],
3583 [F(7), F(0x1.b0b0b0b78cc3f296p-1)], [F(15), F(0x1.2d145116c16ff856p+0)],
3584 [F(17), F(0x1.3afeb354b7d9731ap+0)], [F(31), F(0x1.7dc9e145867e62eap+0)],
3585 [F(33), F(0x1.84bd545e4baeddp+0)], [F(63), F(0x1.cca1950e4511e192p+0)],
3586 [F(65), F(0x1.d01b16f9433cf7b8p+0)], [F(-0), -F.infinity], [F(0), -F.infinity],
3587 [F(10000), F(0x1p+2)],
3588 ];
3589 testLog10(vals);
3590 }}
3591 {
3592 float[2][16] vals = [
3593 [float.nan, float.nan],[-float.nan, float.nan],
3594 [float.infinity, float.infinity], [-float.infinity, float.nan],
3595 [float.min_normal, -0x1.2f703p+5f], [-float.min_normal, float.nan],
3596 [float.max, 0x1.344136p+5f], [-float.max, float.nan],
3597 [float.min_normal / 2, -0x1.31d8b2p+5f], [-float.min_normal / 2, float.nan],
3598 [float.max / 2, 0x1.31d8b2p+5f], [-float.max / 2, float.nan],
3599 [float.min_normal / 3, -0x1.334156p+5f], [-float.min_normal / 3, float.nan],
3600 [float.max / 3, 0x1.30701p+5f], [-float.max / 3, float.nan],
3601 ];
3602 testLog10(vals);
3603 }
3604 {
3605 double[2][16] vals = [
3606 [double.nan, double.nan],[-double.nan, double.nan],
3607 [double.infinity, double.infinity], [-double.infinity, double.nan],
3608 [double.min_normal, -0x1.33a7146f72a42p+8], [-double.min_normal, double.nan],
3609 [double.max, 0x1.34413509f79ffp+8], [-double.max, double.nan],
3610 [double.min_normal / 2, -0x1.33f424bcb522p+8], [-double.min_normal / 2, double.nan],
3611 [double.max / 2, 0x1.33f424bcb522p+8], [-double.max / 2, double.nan],
3612 [double.min_normal / 3, -0x1.3421390dcbe37p+8], [-double.min_normal / 3, double.nan],
3613 [double.max / 3, 0x1.33c7106b9e609p+8], [-double.max / 3, double.nan],
3614 ];
3615 testLog10(vals);
3616 }
3617 alias F = floatTraits!real;
3618 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3619 {{
3620 real[2][16] vals = [
3621 [real.nan, real.nan],[-real.nan, real.nan],
3622 [real.infinity, real.infinity], [-real.infinity, real.nan],
3623 [real.min_normal, -0x1.343793004f503232p+12L], [-real.min_normal, real.nan],
3624 [real.max, 0x1.34413509f79fef32p+12L], [-real.max, real.nan],
3625 [real.min_normal / 2, -0x1.343c6405237810b2p+12L], [-real.min_normal / 2, real.nan],
3626 [real.max / 2, 0x1.343c6405237810b2p+12L], [-real.max / 2, real.nan],
3627 [real.min_normal / 3, -0x1.343f354a34e427bp+12L], [-real.min_normal / 3, real.nan],
3628 [real.max / 3, 0x1.343992c0120bf9b2p+12L], [-real.max / 3, real.nan],
3629 ];
3630 testLog10(vals);
3631 }}
3632 }
3633
3634 /**
3635 * Calculates the natural logarithm of 1 + x.
3636 *
3637 * For very small x, log1p(x) will be more accurate than
3638 * log(1 + x).
3639 *
3640 * $(TABLE_SV
3641 * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?))
3642 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no))
3643 * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
3644 * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes))
3645 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
3646 * )
3647 */
3648 pragma(inline, true)
3649 real log1p(real x) @safe pure nothrow @nogc
3650 {
3651 version (INLINE_YL2X)
3652 {
3653 // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
3654 // ie if -0.29 <= x <= 0.414
3655 import std.math.constants : LN2;
3656 return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
3657 }
3658 else
3659 return log1pImpl(x);
3660 }
3661
3662 /// ditto
3663 pragma(inline, true)
3664 double log1p(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log1p(cast(real) x) : log1pImpl(x); }
3665
3666 /// ditto
3667 pragma(inline, true)
3668 float log1p(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log1p(cast(real) x) : log1pImpl(x); }
3669
3670 // @@@DEPRECATED_[2.112.0]@@@
3671 deprecated("`std.math.exponential.log1p` called with argument types `(int)` matches both "
3672 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3673 real log1p(int x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3674 // @@@DEPRECATED_[2.112.0]@@@
3675 deprecated("`std.math.exponential.log1p` called with argument types `(uint)` matches both "
3676 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3677 real log1p(uint x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3678 // @@@DEPRECATED_[2.112.0]@@@
3679 deprecated("`std.math.exponential.log1p` called with argument types `(long)` matches both "
3680 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3681 real log1p(long x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3682 // @@@DEPRECATED_[2.112.0]@@@
3683 deprecated("`std.math.exponential.log1p` called with argument types `(ulong)` matches both "
3684 ~ "`log1p(real)`, `log1p(double)`, and `log1p(float)`. Cast argument to floating point type instead.")
3685 real log1p(ulong x) @safe pure nothrow @nogc { return log1p(cast(real) x); }
3686
3687 ///
3688 @safe pure unittest
3689 {
3690 import std.math.traits : isIdentical, isNaN;
3691 import std.math.operations : feqrel;
3692
3693 assert(isIdentical(log1p(0.0), 0.0));
3694 assert(log1p(1.0).feqrel(0.69314) > 16);
3695
3696 assert(log1p(-1.0) == -real.infinity);
3697 assert(isNaN(log1p(-2.0)));
3698 assert(log1p(real.nan) is real.nan);
3699 assert(log1p(-real.nan) is -real.nan);
3700 assert(log1p(real.infinity) == real.infinity);
3701 }
3702
3703 private T log1pImpl(T)(T x) @safe pure nothrow @nogc
3704 {
3705 import std.math.traits : isNaN, isInfinity, signbit;
3706 import std.math.algebraic : poly;
3707 import std.math.constants : SQRT1_2, SQRT2;
3708 import std.math : floatTraits, RealFormat;
3709
3710 // Special cases.
3711 if (isNaN(x) || x == 0.0)
3712 return x;
3713 if (isInfinity(x) && !signbit(x))
3714 return x;
3715 if (x == -1.0)
3716 return -T.infinity;
3717 if (x < -1.0)
3718 return T.nan;
3719
3720 alias F = floatTraits!T;
3721 static if (F.realFormat == RealFormat.ieeeSingle ||
3722 F.realFormat == RealFormat.ieeeDouble)
3723 {
3724 // When the input is within the range 1/sqrt(2) <= x+1 <= sqrt(2), compute
3725 // log1p inline. Forwarding to log() would otherwise result in inaccuracies.
3726 const T xp1 = x + 1.0;
3727 if (xp1 >= SQRT1_2 && xp1 <= SQRT2)
3728 {
3729 alias coeffs = LogCoeffs!T;
3730
3731 T px = poly(x, coeffs.logp1P);
3732 T qx = poly(x, coeffs.logp1Q);
3733 const T xx = x * x;
3734 qx = x + ((cast(T) -0.5) * xx + x * (xx * px / qx));
3735 return qx;
3736 }
3737 }
3738
3739 return logImpl!(T, true)(x);
3740 }
3741
3742 @safe @nogc nothrow unittest
3743 {
3744 import std.math : floatTraits, RealFormat;
3745 import std.meta : AliasSeq;
3746
3747 static void testLog1p(T)(T[2][] vals)
3748 {
3749 import std.math.operations : isClose;
3750 import std.math.traits : isNaN;
3751 foreach (ref pair; vals)
3752 {
3753 if (isNaN(pair[1]))
3754 assert(isNaN(log1p(pair[0])));
3755 else
3756 assert(isClose(log1p(pair[0]), pair[1]));
3757 }
3758 }
3759 static foreach (F; AliasSeq!(float, double, real))
3760 {{
3761 F[2][24] vals = [
3762 [F(1), F(0x1.62e42fefa39ef358p-1)], [F(2), F(0x1.193ea7aad030a976p+0)],
3763 [F(4), F(0x1.9c041f7ed8d336bp+0)], [F(8), F(0x1.193ea7aad030a976p+1)],
3764 [F(16), F(0x1.6aa6bc1fa7f79cfp+1)], [F(32), F(0x1.bf8d8f4d5b8d1038p+1)],
3765 [F(64), F(0x1.0b292939429755ap+2)], [F(128), F(0x1.37072a9b5b6cb31p+2)],
3766 [F(256), F(0x1.63241004e9010ad8p+2)], [F(512), F(0x1.8f60adf041bde2a8p+2)],
3767 [F(1024), F(0x1.bbad39ebe1cc08b6p+2)], [F(2048), F(0x1.e801c1698ba4395cp+2)],
3768 [F(3), F(0x1.62e42fefa39ef358p+0)], [F(5), F(0x1.cab0bfa2a2002322p+0)],
3769 [F(7), F(0x1.0a2b23f3bab73682p+1)], [F(15), F(0x1.62e42fefa39ef358p+1)],
3770 [F(17), F(0x1.71f7b3a6b918664cp+1)], [F(31), F(0x1.bb9d3beb8c86b02ep+1)],
3771 [F(33), F(0x1.c35fc81b90df59c6p+1)], [F(63), F(0x1.0a2b23f3bab73682p+2)],
3772 [F(65), F(0x1.0c234da4a23a6686p+2)], [F(-0), F(-0x0p+0)], [F(0), F(0x0p+0)],
3773 [F(10000), F(0x1.26bbed6fbd84182ep+3)],
3774 ];
3775 testLog1p(vals);
3776 }}
3777 {
3778 float[2][16] vals = [
3779 [float.nan, float.nan],[-float.nan, float.nan],
3780 [float.infinity, float.infinity], [-float.infinity, float.nan],
3781 [float.min_normal, 0x1p-126f], [-float.min_normal, -0x1p-126f],
3782 [float.max, 0x1.62e43p+6f], [-float.max, float.nan],
3783 [float.min_normal / 2, 0x0.8p-126f], [-float.min_normal / 2, -0x0.8p-126f],
3784 [float.max / 2, 0x1.601e68p+6f], [-float.max / 2, float.nan],
3785 [float.min_normal / 3, 0x0.555556p-126f], [-float.min_normal / 3, -0x0.555556p-126f],
3786 [float.max / 3, 0x1.5e7f36p+6f], [-float.max / 3, float.nan],
3787 ];
3788 testLog1p(vals);
3789 }
3790 {
3791 double[2][16] vals = [
3792 [double.nan, double.nan],[-double.nan, double.nan],
3793 [double.infinity, double.infinity], [-double.infinity, double.nan],
3794 [double.min_normal, 0x1p-1022], [-double.min_normal, -0x1p-1022],
3795 [double.max, 0x1.62e42fefa39efp+9], [-double.max, double.nan],
3796 [double.min_normal / 2, 0x0.8p-1022], [-double.min_normal / 2, -0x0.8p-1022],
3797 [double.max / 2, 0x1.628b76e3a7b61p+9], [-double.max / 2, double.nan],
3798 [double.min_normal / 3, 0x0.5555555555555p-1022], [-double.min_normal / 3, -0x0.5555555555555p-1022],
3799 [double.max / 3, 0x1.6257909bce36ep+9], [-double.max / 3, double.nan],
3800 ];
3801 testLog1p(vals);
3802 }
3803 alias F = floatTraits!real;
3804 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
3805 {{
3806 real[2][16] vals = [
3807 [real.nan, real.nan],[-real.nan, real.nan],
3808 [real.infinity, real.infinity], [-real.infinity, real.nan],
3809 [real.min_normal, 0x1p-16382L], [-real.min_normal, -0x1p-16382L],
3810 [real.max, 0x1.62e42fefa39ef358p+13L], [-real.max, real.nan],
3811 [real.min_normal / 2, 0x0.8p-16382L], [-real.min_normal / 2, -0x0.8p-16382L],
3812 [real.max / 2, 0x1.62dea45ee3e064dcp+13L], [-real.max / 2, real.nan],
3813 [real.min_normal / 3, 0x0.5555555555555556p-16382L], [-real.min_normal / 3, -0x0.5555555555555556p-16382L],
3814 [real.max / 3, 0x1.62db65fa664871d2p+13L], [-real.max / 3, real.nan],
3815 ];
3816 testLog1p(vals);
3817 }}
3818 }
3819
3820 /***************************************
3821 * Calculates the base-2 logarithm of x:
3822 * $(SUB log, 2)x
3823 *
3824 * $(TABLE_SV
3825 * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?))
3826 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) )
3827 * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) )
3828 * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) )
3829 * )
3830 */
3831 pragma(inline, true)
3832 real log2(real x) @safe pure nothrow @nogc
3833 {
3834 version (INLINE_YL2X)
3835 return core.math.yl2x(x, 1.0L);
3836 else
3837 return log2Impl(x);
3838 }
3839
3840 /// ditto
3841 pragma(inline, true)
3842 double log2(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) log2(cast(real) x) : log2Impl(x); }
3843
3844 /// ditto
3845 pragma(inline, true)
3846 float log2(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) log2(cast(real) x) : log2Impl(x); }
3847
3848 // @@@DEPRECATED_[2.112.0]@@@
3849 deprecated("`std.math.exponential.log2` called with argument types `(int)` matches both "
3850 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3851 real log2(int x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3852 // @@@DEPRECATED_[2.112.0]@@@
3853 deprecated("`std.math.exponential.log2` called with argument types `(uint)` matches both "
3854 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3855 real log2(uint x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3856 // @@@DEPRECATED_[2.112.0]@@@
3857 deprecated("`std.math.exponential.log2` called with argument types `(long)` matches both "
3858 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3859 real log2(long x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3860 // @@@DEPRECATED_[2.112.0]@@@
3861 deprecated("`std.math.exponential.log2` called with argument types `(ulong)` matches both "
3862 ~ "`log2(real)`, `log2(double)`, and `log2(float)`. Cast argument to floating point type instead.")
3863 real log2(ulong x) @safe pure nothrow @nogc { return log2(cast(real) x); }
3864
3865 ///
3866 @safe unittest
3867 {
3868 import std.math.operations : isClose;
3869
3870 assert(isClose(log2(1024.0L), 10));
3871 }
3872
3873 @safe @nogc nothrow unittest
3874 {
3875 import std.math.operations : isClose;
3876
3877 // check if values are equal to 19 decimal digits of precision
3878 assert(isClose(log2(1024.0L), 10, 1e-18));
3879 }
3880
3881 private T log2Impl(T)(T x) @safe pure nothrow @nogc
3882 {
3883 import std.math.traits : isNaN, isInfinity, signbit;
3884 import std.math.constants : SQRT1_2, LOG2E;
3885 import std.math.algebraic : poly;
3886 import std.math : floatTraits, RealFormat;
3887
3888 alias coeffs = LogCoeffs!T;
3889 alias F = floatTraits!T;
3890
3891 // Special cases are the same as for log.
3892 if (isNaN(x))
3893 return x;
3894 if (isInfinity(x) && !signbit(x))
3895 return x;
3896 if (x == 0.0)
3897 return -T.infinity;
3898 if (x < 0.0)
3899 return T.nan;
3900
3901 // Separate mantissa from exponent.
3902 // Note, frexp is used so that denormal numbers will be handled properly.
3903 T y, z;
3904 int exp;
3905
3906 x = frexp(x, exp);
3907
3908 static if (F.realFormat == RealFormat.ieeeDouble ||
3909 F.realFormat == RealFormat.ieeeExtended ||
3910 F.realFormat == RealFormat.ieeeExtended53 ||
3911 F.realFormat == RealFormat.ieeeQuadruple)
3912 {
3913 // Logarithm using log(x) = z + z^^3 R(z) / S(z),
3914 // where z = 2(x - 1)/(x + 1)
3915 if ((exp > 2) || (exp < -2))
3916 {
3917 if (x < SQRT1_2)
3918 { // 2(2x - 1)/(2x + 1)
3919 exp -= 1;
3920 z = x - 0.5;
3921 y = 0.5 * z + 0.5;
3922 }
3923 else
3924 { // 2(x - 1)/(x + 1)
3925 z = x - 0.5;
3926 z -= 0.5;
3927 y = 0.5 * x + 0.5;
3928 }
3929 x = z / y;
3930 z = x * x;
3931 y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
3932 goto Ldone;
3933 }
3934 }
3935
3936 // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
3937 if (x < SQRT1_2)
3938 {
3939 exp -= 1;
3940 x = 2.0 * x - 1.0;
3941 }
3942 else
3943 x = x - 1.0;
3944
3945 z = x * x;
3946 static if (F.realFormat == RealFormat.ieeeSingle)
3947 y = x * (z * poly(x, coeffs.log2P));
3948 else
3949 y = x * (z * poly(x, coeffs.log2P) / poly(x, coeffs.log2Q));
3950 y = y - 0.5 * z;
3951
3952 // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
3953 // This sequence of operations is critical and it may be horribly
3954 // defeated by some compiler optimizers.
3955 Ldone:
3956 z = y * (LOG2E - 1.0);
3957 z += x * (LOG2E - 1.0);
3958 z += y;
3959 z += x;
3960 z += exp;
3961
3962 return z;
3963 }
3964
3965 @safe @nogc nothrow unittest
3966 {
3967 import std.math : floatTraits, RealFormat;
3968 import std.meta : AliasSeq;
3969
3970 static void testLog2(T)(T[2][] vals)
3971 {
3972 import std.math.operations : isClose;
3973 import std.math.traits : isNaN;
3974 foreach (ref pair; vals)
3975 {
3976 if (isNaN(pair[1]))
3977 assert(isNaN(log2(pair[0])));
3978 else
3979 assert(isClose(log2(pair[0]), pair[1]));
3980 }
3981 }
3982 static foreach (F; AliasSeq!(float, double, real))
3983 {{
3984 F[2][24] vals = [
3985 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
3986 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
3987 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
3988 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
3989 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
3990 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
3991 [F(3), F(0x1.95c01a39fbd687ap+0)], [F(5), F(0x1.2934f0979a3715fcp+1)],
3992 [F(7), F(0x1.675767f54042cd9ap+1)], [F(15), F(0x1.f414fdb4982259ccp+1)],
3993 [F(17), F(0x1.0598fdbeb244c5ap+2)], [F(31), F(0x1.3d118d66c4d4e554p+2)],
3994 [F(33), F(0x1.42d75a6eb1dfb0e6p+2)], [F(63), F(0x1.7e8bc1179e0caa9cp+2)],
3995 [F(65), F(0x1.816e79685c2d2298p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
3996 [F(10000), F(0x1.a934f0979a3715fcp+3)],
3997 ];
3998 testLog2(vals);
3999 }}
4000 {
4001 float[2][16] vals = [
4002 [float.nan, float.nan],[-float.nan, float.nan],
4003 [float.infinity, float.infinity], [-float.infinity, float.nan],
4004 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, float.nan],
4005 [float.max, 0x1p+7f], [-float.max, float.nan],
4006 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, float.nan],
4007 [float.max / 2, 0x1.fcp+6f], [-float.max / 2, float.nan],
4008 [float.min_normal / 3, -0x1.fe57p+6f], [-float.min_normal / 3, float.nan],
4009 [float.max / 3, 0x1.f9a9p+6f], [-float.max / 3, float.nan],
4010 ];
4011 testLog2(vals);
4012 }
4013 {
4014 double[2][16] vals = [
4015 [double.nan, double.nan],[-double.nan, double.nan],
4016 [double.infinity, double.infinity], [-double.infinity, double.nan],
4017 [double.min_normal, -0x1.ffp+9], [-double.min_normal, double.nan],
4018 [double.max, 0x1p+10], [-double.max, double.nan],
4019 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, double.nan],
4020 [double.max / 2, 0x1.ff8p+9], [-double.max / 2, double.nan],
4021 [double.min_normal / 3, -0x1.ffcae00d1cfdfp+9], [-double.min_normal / 3, double.nan],
4022 [double.max / 3, 0x1.ff351ff2e3021p+9], [-double.max / 3, double.nan],
4023 ];
4024 testLog2(vals);
4025 }
4026 alias F = floatTraits!real;
4027 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4028 {{
4029 real[2][16] vals = [
4030 [real.nan, real.nan],[-real.nan, real.nan],
4031 [real.infinity, real.infinity], [-real.infinity, real.nan],
4032 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, real.nan],
4033 [real.max, 0x1p+14L], [-real.max, real.nan],
4034 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, real.nan],
4035 [real.max / 2, 0x1.fff8p+13L], [-real.max / 2, real.nan],
4036 [real.min_normal / 3, -0x1.fffcae00d1cfdeb4p+13L], [-real.min_normal / 3, real.nan],
4037 [real.max / 3, 0x1.fff351ff2e30214cp+13L], [-real.max / 3, real.nan],
4038 ];
4039 testLog2(vals);
4040 }}
4041 }
4042
4043 /*****************************************
4044 * Extracts the exponent of x as a signed integral value.
4045 *
4046 * If x is subnormal, it is treated as if it were normalized.
4047 * For a positive, finite x:
4048 *
4049 * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX
4050 *
4051 * $(TABLE_SV
4052 * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) )
4053 * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
4054 * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) )
4055 * )
4056 */
4057 pragma(inline, true)
4058 real logb(real x) @trusted pure nothrow @nogc
4059 {
4060 version (InlineAsm_X87_MSVC)
4061 return logbAsm(x);
4062 else
4063 return logbImpl(x);
4064 }
4065
4066 /// ditto
4067 pragma(inline, true)
4068 double logb(double x) @trusted pure nothrow @nogc { return logbImpl(x); }
4069
4070 /// ditto
4071 pragma(inline, true)
4072 float logb(float x) @trusted pure nothrow @nogc { return logbImpl(x); }
4073
4074 ///
4075 @safe @nogc nothrow unittest
4076 {
4077 assert(logb(1.0) == 0);
4078 assert(logb(100.0) == 6);
4079
4080 assert(logb(0.0) == -real.infinity);
4081 assert(logb(real.infinity) == real.infinity);
4082 assert(logb(-real.infinity) == real.infinity);
4083 }
4084
4085 @safe @nogc nothrow unittest
4086 {
4087 import std.meta : AliasSeq;
4088 import std.typecons : Tuple;
4089 import std.math.traits : isNaN;
4090 static foreach (F; AliasSeq!(float, double, real))
4091 {{
4092 alias T = Tuple!(F, F);
4093 T[17] vals = // x, logb(x)
4094 [
4095 T(1.0 , 0 ),
4096 T(100.0 , 6 ),
4097 T(0.0 , -F.infinity),
4098 T(-0.0 , -F.infinity),
4099 T(1024 , 10 ),
4100 T(-2000 , 10 ),
4101 T(0x0.1p-127 , -131 ),
4102 T(0x0.01p-127 , -135 ),
4103 T(0x0.011p-127 , -135 ),
4104 T(F.nan , F.nan ),
4105 T(-F.nan , F.nan ),
4106 T(F.infinity , F.infinity ),
4107 T(-F.infinity , F.infinity ),
4108 T(F.min_normal , F.min_exp-1),
4109 T(-F.min_normal, F.min_exp-1),
4110 T(F.max , F.max_exp-1),
4111 T(-F.max , F.max_exp-1),
4112 ];
4113
4114 foreach (elem; vals)
4115 {
4116 if (isNaN(elem[1]))
4117 assert(isNaN(logb(elem[1])));
4118 else
4119 assert(logb(elem[0]) == elem[1]);
4120 }
4121 }}
4122 }
4123
4124 version (InlineAsm_X87_MSVC)
4125 private T logbAsm(T)(T x) @trusted pure nothrow @nogc
4126 {
4127 version (X86_64)
4128 {
4129 asm pure nothrow @nogc
4130 {
4131 naked ;
4132 fld real ptr [RCX] ;
4133 fxtract ;
4134 fstp ST(0) ;
4135 ret ;
4136 }
4137 }
4138 else
4139 {
4140 asm pure nothrow @nogc
4141 {
4142 fld x ;
4143 fxtract ;
4144 fstp ST(0) ;
4145 }
4146 }
4147 }
4148
4149 private T logbImpl(T)(T x) @trusted pure nothrow @nogc
4150 {
4151 import std.math.traits : isFinite;
4152
4153 // Handle special cases.
4154 if (!isFinite(x))
4155 return x * x;
4156 if (x == 0)
4157 return -1 / (x * x);
4158
4159 return ilogb(x);
4160 }
4161
4162 @safe @nogc nothrow unittest
4163 {
4164 import std.math : floatTraits, RealFormat;
4165 import std.meta : AliasSeq;
4166
4167 static void testLogb(T)(T[2][] vals)
4168 {
4169 import std.math.operations : isClose;
4170 import std.math.traits : isNaN;
4171 foreach (ref pair; vals)
4172 {
4173 if (isNaN(pair[1]))
4174 assert(isNaN(logb(pair[0])));
4175 else
4176 assert(isClose(logb(pair[0]), pair[1]));
4177 }
4178 }
4179 static foreach (F; AliasSeq!(float, double, real))
4180 {{
4181 F[2][24] vals = [
4182 [F(1), F(0x0p+0)], [F(2), F(0x1p+0)],
4183 [F(4), F(0x1p+1)], [F(8), F(0x1.8p+1)],
4184 [F(16), F(0x1p+2)], [F(32), F(0x1.4p+2)],
4185 [F(64), F(0x1.8p+2)], [F(128), F(0x1.cp+2)],
4186 [F(256), F(0x1p+3)], [F(512), F(0x1.2p+3)],
4187 [F(1024), F(0x1.4p+3)], [F(2048), F(0x1.6p+3)],
4188 [F(3), F(0x1p+0)], [F(5), F(0x1p+1)],
4189 [F(7), F(0x1p+1)], [F(15), F(0x1.8p+1)],
4190 [F(17), F(0x1p+2)], [F(31), F(0x1p+2)],
4191 [F(33), F(0x1.4p+2)], [F(63), F(0x1.4p+2)],
4192 [F(65), F(0x1.8p+2)], [F(-0), -F.infinity], [F(0), -F.infinity],
4193 [F(10000), F(0x1.ap+3)],
4194 ];
4195 testLogb(vals);
4196 }}
4197 {
4198 float[2][16] vals = [
4199 [float.nan, float.nan],[-float.nan, float.nan],
4200 [float.infinity, float.infinity], [-float.infinity, float.infinity],
4201 [float.min_normal, -0x1.f8p+6f], [-float.min_normal, -0x1.f8p+6f],
4202 [float.max, 0x1.fcp+6f], [-float.max, 0x1.fcp+6f],
4203 [float.min_normal / 2, -0x1.fcp+6f], [-float.min_normal / 2, -0x1.fcp+6f],
4204 [float.max / 2, 0x1.f8p+6f], [-float.max / 2, 0x1.f8p+6f],
4205 [float.min_normal / 3, -0x1p+7f], [-float.min_normal / 3, -0x1p+7f],
4206 [float.max / 3, 0x1.f8p+6f], [-float.max / 3, 0x1.f8p+6f],
4207 ];
4208 testLogb(vals);
4209 }
4210 {
4211 double[2][16] vals = [
4212 [double.nan, double.nan],[-double.nan, double.nan],
4213 [double.infinity, double.infinity], [-double.infinity, double.infinity],
4214 [double.min_normal, -0x1.ffp+9], [-double.min_normal, -0x1.ffp+9],
4215 [double.max, 0x1.ff8p+9], [-double.max, 0x1.ff8p+9],
4216 [double.min_normal / 2, -0x1.ff8p+9], [-double.min_normal / 2, -0x1.ff8p+9],
4217 [double.max / 2, 0x1.ffp+9], [-double.max / 2, 0x1.ffp+9],
4218 [double.min_normal / 3, -0x1p+10], [-double.min_normal / 3, -0x1p+10],
4219 [double.max / 3, 0x1.ffp+9], [-double.max / 3, 0x1.ffp+9],
4220 ];
4221 testLogb(vals);
4222 }
4223 alias F = floatTraits!real;
4224 static if (F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeQuadruple)
4225 {{
4226 real[2][16] vals = [
4227 [real.nan, real.nan],[-real.nan, real.nan],
4228 [real.infinity, real.infinity], [-real.infinity, real.infinity],
4229 [real.min_normal, -0x1.fffp+13L], [-real.min_normal, -0x1.fffp+13L],
4230 [real.max, 0x1.fff8p+13L], [-real.max, 0x1.fff8p+13L],
4231 [real.min_normal / 2, -0x1.fff8p+13L], [-real.min_normal / 2, -0x1.fff8p+13L],
4232 [real.max / 2, 0x1.fffp+13L], [-real.max / 2, 0x1.fffp+13L],
4233 [real.min_normal / 3, -0x1p+14L], [-real.min_normal / 3, -0x1p+14L],
4234 [real.max / 3, 0x1.fffp+13L], [-real.max / 3, 0x1.fffp+13L],
4235 ];
4236 testLogb(vals);
4237 }}
4238 }
4239
4240 /*************************************
4241 * Efficiently calculates x * 2$(SUPERSCRIPT n).
4242 *
4243 * scalbn handles underflow and overflow in
4244 * the same fashion as the basic arithmetic operators.
4245 *
4246 * $(TABLE_SV
4247 * $(TR $(TH x) $(TH scalb(x)))
4248 * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) )
4249 * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
4250 * )
4251 */
4252 pragma(inline, true)
4253 real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4254
4255 /// ditto
4256 pragma(inline, true)
4257 double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4258
4259 /// ditto
4260 pragma(inline, true)
4261 float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); }
4262
4263 ///
4264 @safe pure nothrow @nogc unittest
4265 {
4266 assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4267 assert(scalbn(-real.infinity, 5) == -real.infinity);
4268 assert(scalbn(2.0,10) == 2048.0);
4269 assert(scalbn(2048.0f,-10) == 2.0f);
4270 }
4271
4272 pragma(inline, true)
4273 private F _scalbn(F)(F x, int n)
4274 {
4275 import std.math.traits : isInfinity;
4276
4277 if (__ctfe)
4278 {
4279 // Handle special cases.
4280 if (x == F(0.0) || isInfinity(x))
4281 return x;
4282 }
4283 return core.math.ldexp(x, n);
4284 }
4285
4286 @safe pure nothrow @nogc unittest
4287 {
4288 // CTFE-able test
4289 static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L);
4290 static assert(scalbn(-real.infinity, 5) == -real.infinity);
4291 // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not.
4292 enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2;
4293 enum int n = resultExponent - initialExponent;
4294 enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent);
4295 enum staticResult = scalbn(x, n);
4296 static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent));
4297 assert(scalbn(x, n) == staticResult);
4298 }
4299
Mirror of the offical gcc git repository hosted at gcc.gnu.org