std.math.trigonometry source code

1 // Written in the D programming language.
2 
3 /**
4 This is a submodule of $(MREF std, math).
5 
6 It contains several trigonometric functions.
7 
8 Copyright: Copyright The D Language Foundation 2000 - 2011.
9            D implementations of tan, atan, and atan2 functions are based on the
10            CEPHES math library, which is Copyright (C) 2001 Stephen L. Moshier
11            $(LT)steve@moshier.net$(GT) and are incorporated herein by permission
12            of the author. The author reserves the right to distribute this
13            material elsewhere under different copying permissions.
14            These modifications are distributed here under the following terms:
15 License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
16 Authors:   $(HTTP digitalmars.com, Walter Bright), Don Clugston,
17            Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
18 Source: $(PHOBOSSRC std/math/trigonometry.d)
19 
20 Macros:
21     TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
22                <caption>Special Values</caption>
23                $0</table>
24     SVH = $(TR $(TH $1) $(TH $2))
25     SV  = $(TR $(TD $1) $(TD $2))
26     TH3 = $(TR $(TH $1) $(TH $2) $(TH $3))
27     TD3 = $(TR $(TD $1) $(TD $2) $(TD $3))
28     TABLE_DOMRG = <table border="1" cellpadding="4" cellspacing="0">
29                   $(SVH Domain X, Range Y)
30                   $(SV $1, $2)
31                   </table>
32     DOMAIN=$1
33     RANGE=$1
34     POWER = $1<sup>$2</sup>
35     NAN = $(RED NAN)
36     PLUSMN = &plusmn;
37     INFIN = &infin;
38     PLUSMNINF = &plusmn;&infin;
39  */
40 
41 module std.math.trigonometry;
42 
43 static import core.math;
44 
45 version (D_InlineAsm_X86)    version = InlineAsm_X86_Any;
46 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
47 
48 version (InlineAsm_X86_Any) version = InlineAsm_X87;
49 version (InlineAsm_X87)
50 {
51     static assert(real.mant_dig == 64);
52     version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC;
53 }
54 
55 /***********************************
56  * Returns cosine of x. x is in radians.
57  *
58  *      $(TABLE_SV
59  *      $(TR $(TH x)                 $(TH cos(x)) $(TH invalid?))
60  *      $(TR $(TD $(NAN))            $(TD $(NAN)) $(TD yes)     )
61  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)     )
62  *      )
63  * Bugs:
64  *      Results are undefined if |x| >= $(POWER 2,64).
65  */
66 pragma(inline, true)
67 real cos(real x) @safe pure nothrow @nogc { return core.math.cos(x); }
68 ///ditto
69 pragma(inline, true)
70 double cos(double x) @safe pure nothrow @nogc { return core.math.cos(x); }
71 ///ditto
72 pragma(inline, true)
73 float cos(float x) @safe pure nothrow @nogc { return core.math.cos(x); }
74 
75 ///
76 @safe unittest
77 {
78     import std.math.operations : isClose;
79 
80     assert(cos(0.0) == 1.0);
81     assert(cos(1.0).isClose(0.5403023059));
82     assert(cos(3.0).isClose(-0.9899924966));
83 }
84 
85 @safe unittest
86 {
87     real function(real) pcos = &cos;
88     assert(pcos != null);
89 }
90 
91 @safe pure nothrow @nogc unittest
92 {
93     import std.math.algebraic : fabs;
94 
95     float f = cos(-2.0f);
96     assert(fabs(f - -0.416147f) < .00001);
97 
98     double d = cos(-2.0);
99     assert(fabs(d - -0.416147f) < .00001);
100 
101     real r = cos(-2.0L);
102     assert(fabs(r - -0.416147f) < .00001);
103 }
104 
105 /***********************************
106  * Returns $(HTTP en.wikipedia.org/wiki/Sine, sine) of x. x is in $(HTTP en.wikipedia.org/wiki/Radian, radians).
107  *
108  *      $(TABLE_SV
109  *      $(TH3 x           ,  sin(x)      ,  invalid?)
110  *      $(TD3 $(NAN)      ,  $(NAN)      ,  yes     )
111  *      $(TD3 $(PLUSMN)0.0,  $(PLUSMN)0.0,  no      )
112  *      $(TD3 $(PLUSMNINF),  $(NAN)      ,  yes     )
113  *      )
114  *
115  * Params:
116  *      x = angle in radians (not degrees)
117  * Returns:
118  *      sine of x
119  * See_Also:
120  *      $(MYREF cos), $(MYREF tan), $(MYREF asin)
121  * Bugs:
122  *      Results are undefined if |x| >= $(POWER 2,64).
123  */
124 pragma(inline, true)
125 real sin(real x) @safe pure nothrow @nogc { return core.math.sin(x); }
126 ///ditto
127 pragma(inline, true)
128 double sin(double x) @safe pure nothrow @nogc { return core.math.sin(x); }
129 ///ditto
130 pragma(inline, true)
131 float sin(float x) @safe pure nothrow @nogc { return core.math.sin(x); }
132 
133 ///
134 @safe unittest
135 {
136     import std.math.constants : PI;
137     import std.stdio : writefln;
138 
139     void someFunc()
140     {
141       real x = 30.0;
142       auto result = sin(x * (PI / 180)); // convert degrees to radians
143       writefln("The sine of %s degrees is %s", x, result);
144     }
145 }
146 
147 @safe unittest
148 {
149     real function(real) psin = &sin;
150     assert(psin != null);
151 }
152 
153 @safe pure nothrow @nogc unittest
154 {
155     import std.math.algebraic : fabs;
156 
157     float f = sin(-2.0f);
158     assert(fabs(f - -0.909297f) < .00001);
159 
160     double d = sin(-2.0);
161     assert(fabs(d - -0.909297f) < .00001);
162 
163     real r = sin(-2.0L);
164     assert(fabs(r - -0.909297f) < .00001);
165 }
166 
167 /****************************************************************************
168  * Returns tangent of x. x is in radians.
169  *
170  *      $(TABLE_SV
171  *      $(TR $(TH x)             $(TH tan(x))       $(TH invalid?))
172  *      $(TR $(TD $(NAN))        $(TD $(NAN))       $(TD yes))
173  *      $(TR $(TD $(PLUSMN)0.0)  $(TD $(PLUSMN)0.0) $(TD no))
174  *      $(TR $(TD $(PLUSMNINF))  $(TD $(NAN))       $(TD yes))
175  *      )
176  */
177 pragma(inline, true)
178 real tan(real x) @safe pure nothrow @nogc
179 {
180     version (InlineAsm_X87)
181     {
182         if (!__ctfe)
183             return tanAsm(x);
184     }
185     return tanImpl(x);
186 }
187 
188 /// ditto
189 pragma(inline, true)
190 double tan(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) tan(cast(real) x) : tanImpl(x); }
191 
192 /// ditto
193 pragma(inline, true)
194 float tan(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) tan(cast(real) x) : tanImpl(x); }
195 
196 ///
197 @safe unittest
198 {
199     import std.math.operations : isClose;
200     import std.math.traits : isIdentical;
201     import std.math.constants : PI;
202     import std.math.algebraic : sqrt;
203 
204     assert(isIdentical(tan(0.0), 0.0));
205     assert(tan(PI).isClose(0, 0.0, 1e-10));
206     assert(tan(PI / 3).isClose(sqrt(3.0)));
207 }
208 
209 version (InlineAsm_X87)
210 private real tanAsm(real x) @trusted pure nothrow @nogc
211 {
212     // Separating `return real.nan` from the asm block on LDC produces unintended
213     // behaviour as additional instructions are generated, invalidating the asm
214     // logic inside the previous block. To circumvent this, we can push rnan
215     // manually by creating an immutable variable in the stack.
216     immutable rnan = real.nan;
217 
218     version (X86)
219     {
220     asm pure nothrow @nogc
221     {
222         fld     x[EBP]                  ; // load theta
223         fxam                            ; // test for oddball values
224         fstsw   AX                      ;
225         sahf                            ;
226         jc      trigerr                 ; // x is NAN, infinity, or empty
227                                           // 387's can handle subnormals
228 SC18:   fptan                           ;
229         fstsw   AX                      ;
230         sahf                            ;
231         jnp     Clear1                  ; // C2 = 1 (x is out of range)
232 
233         // Do argument reduction to bring x into range
234         fldpi                           ;
235         fxch                            ;
236 SC17:   fprem1                          ;
237         fstsw   AX                      ;
238         sahf                            ;
239         jp      SC17                    ;
240         fstp    ST(1)                   ; // remove pi from stack
241         jmp     SC18                    ;
242 
243 trigerr:
244         jnp     Lret                    ; // if theta is NAN, return theta
245         fstp    ST(0)                   ; // dump theta
246         fld     rnan                    ; // return rnan
247         jmp     Lret                    ;
248 Clear1:
249         fstp    ST(0)                   ; // dump X, which is always 1
250 Lret:
251         ;
252     }
253     }
254     else version (X86_64)
255     {
256         version (Win64)
257         {
258             asm pure nothrow @nogc
259             {
260                 fld     real ptr [RCX]  ; // load theta
261             }
262         }
263         else
264         {
265             asm pure nothrow @nogc
266             {
267                 fld     x[RBP]          ; // load theta
268             }
269         }
270     asm pure nothrow @nogc
271     {
272         fxam                            ; // test for oddball values
273         fstsw   AX                      ;
274         test    AH,1                    ;
275         jnz     trigerr                 ; // x is NAN, infinity, or empty
276                                           // 387's can handle subnormals
277 SC18:   fptan                           ;
278         fstsw   AX                      ;
279         test    AH,4                    ;
280         jz      Clear1                  ; // C2 = 1 (x is out of range)
281 
282         // Do argument reduction to bring x into range
283         fldpi                           ;
284         fxch                            ;
285 SC17:   fprem1                          ;
286         fstsw   AX                      ;
287         test    AH,4                    ;
288         jnz     SC17                    ;
289         fstp    ST(1)                   ; // remove pi from stack
290         jmp     SC18                    ;
291 
292 trigerr:
293         test    AH,4                    ;
294         jz      Lret                    ; // if theta is NAN, return theta
295         fstp    ST(0)                   ; // dump theta
296         fld     rnan                    ; // return rnan
297         jmp     Lret                    ;
298 Clear1:
299         fstp    ST(0)                   ; // dump X, which is always 1
300 Lret:
301         ;
302     }
303     }
304     else
305         static assert(0);
306 }
307 
308 private T tanImpl(T)(T x) @safe pure nothrow @nogc
309 {
310     import std.math : floatTraits, RealFormat;
311     import std.math.constants : PI, PI_4;
312     import std.math.rounding : floor;
313     import std.math.algebraic : poly;
314     import std.math.traits : isInfinity, isNaN, signbit;
315 
316     // Coefficients for tan(x) and PI/4 split into three parts.
317     enum realFormat = floatTraits!T.realFormat;
318     static if (realFormat == RealFormat.ieeeQuadruple)
319     {
320         static immutable T[6] P = [
321             2.883414728874239697964612246732416606301E10L,
322             -2.307030822693734879744223131873392503321E9L,
323             5.160188250214037865511600561074819366815E7L,
324             -4.249691853501233575668486667664718192660E5L,
325             1.272297782199996882828849455156962260810E3L,
326             -9.889929415807650724957118893791829849557E-1L
327         ];
328         static immutable T[7] Q = [
329             8.650244186622719093893836740197250197602E10L,
330             -4.152206921457208101480801635640958361612E10L,
331             2.758476078803232151774723646710890525496E9L,
332             -5.733709132766856723608447733926138506824E7L,
333             4.529422062441341616231663543669583527923E5L,
334             -1.317243702830553658702531997959756728291E3L,
335             1.0
336         ];
337 
338         enum T P1 =
339             7.853981633974483067550664827649598009884357452392578125E-1L;
340         enum T P2 =
341             2.8605943630549158983813312792950660807511260829685741796657E-18L;
342         enum T P3 =
343             2.1679525325309452561992610065108379921905808E-35L;
344     }
345     else static if (realFormat == RealFormat.ieeeExtended ||
346                     realFormat == RealFormat.ieeeDouble)
347     {
348         static immutable T[3] P = [
349            -1.7956525197648487798769E7L,
350             1.1535166483858741613983E6L,
351            -1.3093693918138377764608E4L,
352         ];
353         static immutable T[5] Q = [
354            -5.3869575592945462988123E7L,
355             2.5008380182335791583922E7L,
356            -1.3208923444021096744731E6L,
357             1.3681296347069295467845E4L,
358             1.0000000000000000000000E0L,
359         ];
360 
361         enum T P1 = 7.853981554508209228515625E-1L;
362         enum T P2 = 7.946627356147928367136046290398E-9L;
363         enum T P3 = 3.061616997868382943065164830688E-17L;
364     }
365     else static if (realFormat == RealFormat.ieeeSingle)
366     {
367         static immutable T[6] P = [
368             3.33331568548E-1,
369             1.33387994085E-1,
370             5.34112807005E-2,
371             2.44301354525E-2,
372             3.11992232697E-3,
373             9.38540185543E-3,
374         ];
375 
376         enum T P1 = 0.78515625;
377         enum T P2 = 2.4187564849853515625E-4;
378         enum T P3 = 3.77489497744594108E-8;
379     }
380     else
381         static assert(0, "no coefficients for tan()");
382 
383     // Special cases.
384     if (x == cast(T) 0.0 || isNaN(x))
385         return x;
386     if (isInfinity(x))
387         return T.nan;
388 
389     // Make argument positive but save the sign.
390     bool sign = false;
391     if (signbit(x))
392     {
393         sign = true;
394         x = -x;
395     }
396 
397     // Compute x mod PI/4.
398     static if (realFormat == RealFormat.ieeeSingle)
399     {
400         enum T FOPI = 4 / PI;
401         int j = cast(int) (FOPI * x);
402         T y = j;
403         T z;
404     }
405     else
406     {
407         T y = floor(x / cast(T) PI_4);
408         // Strip high bits of integer part.
409         enum T highBitsFactor = (realFormat == RealFormat.ieeeDouble ? 0x1p3 : 0x1p4);
410         enum T highBitsInv = 1.0 / highBitsFactor;
411         T z = y * highBitsInv;
412         // Compute y - 2^numHighBits * (y / 2^numHighBits).
413         z = y - highBitsFactor * floor(z);
414 
415         // Integer and fraction part modulo one octant.
416         int j = cast(int)(z);
417     }
418 
419     // Map zeros and singularities to origin.
420     if (j & 1)
421     {
422         j += 1;
423         y += cast(T) 1.0;
424     }
425 
426     z = ((x - y * P1) - y * P2) - y * P3;
427     const T zz = z * z;
428 
429     enum T zzThreshold = (realFormat == RealFormat.ieeeSingle ? 1.0e-4L :
430                           realFormat == RealFormat.ieeeDouble ? 1.0e-14L : 1.0e-20L);
431     if (zz > zzThreshold)
432     {
433         static if (realFormat == RealFormat.ieeeSingle)
434             y = z + z * (zz * poly(zz, P));
435         else
436             y = z + z * (zz * poly(zz, P) / poly(zz, Q));
437     }
438     else
439         y = z;
440 
441     if (j & 2)
442         y = (cast(T) -1.0) / y;
443 
444     return (sign) ? -y : y;
445 }
446 
447 @safe @nogc nothrow unittest
448 {
449     static void testTan(T)()
450     {
451         import std.math.operations : CommonDefaultFor, isClose, NaN;
452         import std.math.traits : isIdentical, isNaN;
453         import std.math.constants : PI, PI_4;
454 
455         // ±0
456         const T zero = 0.0;
457         assert(isIdentical(tan(zero), zero));
458         assert(isIdentical(tan(-zero), -zero));
459         // ±∞
460         const T inf = T.infinity;
461         assert(isNaN(tan(inf)));
462         assert(isNaN(tan(-inf)));
463         // NaN
464         const T specialNaN = NaN(0x0123L);
465         assert(isIdentical(tan(specialNaN), specialNaN));
466 
467         static immutable T[2][] vals =
468         [
469             // angle, tan
470             [   .5,  .546302489843790513255L],
471             [   1,   1.55740772465490223050L],
472             [   1.5, 14.1014199471717193876L],
473             [   2,  -2.18503986326151899164L],
474             [   2.5,-.747022297238660279355L],
475             [   3,  -.142546543074277805295L],
476             [   3.5, .374585640158594666330L],
477             [   4,   1.15782128234957758313L],
478             [   4.5, 4.63733205455118446831L],
479             [   5,  -3.38051500624658563698L],
480             [   5.5,-.995584052213885017701L],
481             [   6,  -.291006191384749157053L],
482             [   6.5, .220277200345896811825L],
483             [   10,  .648360827459086671259L],
484 
485             // special angles
486             [   PI_4,   1],
487             //[   PI_2,   T.infinity], // PI_2 is not _exactly_ pi/2.
488             [   3*PI_4, -1],
489             [   PI,     0],
490             [   5*PI_4, 1],
491             //[   3*PI_2, -T.infinity],
492             [   7*PI_4, -1],
493             [   2*PI,   0],
494          ];
495 
496         foreach (ref val; vals)
497         {
498             T x = val[0];
499             T r = val[1];
500             T t = tan(x);
501 
502             //printf("tan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) t, cast(real) r);
503             assert(isClose(r, t, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
504 
505             x = -x;
506             r = -r;
507             t = tan(x);
508             //printf("tan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) t, cast(real) r);
509             assert(isClose(r, t, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
510         }
511     }
512 
513     import std.meta : AliasSeq;
514     foreach (T; AliasSeq!(real, double, float))
515         testTan!T();
516 
517     import std.math.operations : isClose;
518     import std.math.constants : PI;
519     import std.math.algebraic : sqrt;
520     assert(isClose(tan(PI / 3), sqrt(3.0L), real.sizeof > double.sizeof ? 1e-15 : 1e-14));
521 }
522 
523 @safe pure nothrow @nogc unittest
524 {
525     import std.math.algebraic : fabs;
526     import std.math.traits : isNaN;
527 
528     float f = tan(-2.0f);
529     assert(fabs(f - 2.18504f) < .00001);
530 
531     double d = tan(-2.0);
532     assert(fabs(d - 2.18504f) < .00001);
533 
534     real r = tan(-2.0L);
535     assert(fabs(r - 2.18504f) < .00001);
536 
537     // Verify correct behavior for large inputs
538     assert(!isNaN(tan(0x1p63)));
539     assert(!isNaN(tan(-0x1p63)));
540     static if (real.mant_dig >= 64)
541     {
542         assert(!isNaN(tan(0x1p300L)));
543         assert(!isNaN(tan(-0x1p300L)));
544     }
545 }
546 
547 /***************
548  * Calculates the arc cosine of x,
549  * returning a value ranging from 0 to $(PI).
550  *
551  *      $(TABLE_SV
552  *      $(TR $(TH x)         $(TH acos(x)) $(TH invalid?))
553  *      $(TR $(TD $(GT)1.0)  $(TD $(NAN))  $(TD yes))
554  *      $(TR $(TD $(LT)-1.0) $(TD $(NAN))  $(TD yes))
555  *      $(TR $(TD $(NAN))    $(TD $(NAN))  $(TD yes))
556  *  )
557  */
558 real acos(real x) @safe pure nothrow @nogc
559 {
560     import core.math : sqrt;
561 
562     return atan2(sqrt(1-x*x), x);
563 }
564 
565 /// ditto
566 double acos(double x) @safe pure nothrow @nogc { return acos(cast(real) x); }
567 
568 /// ditto
569 float acos(float x) @safe pure nothrow @nogc  { return acos(cast(real) x); }
570 
571 ///
572 @safe unittest
573 {
574     import std.math.operations : isClose;
575     import std.math.traits : isNaN;
576     import std.math.constants : PI;
577 
578     assert(acos(0.0).isClose(1.570796327));
579     assert(acos(0.5).isClose(PI / 3));
580     assert(acos(PI).isNaN);
581 }
582 
583 @safe @nogc nothrow unittest
584 {
585     import std.math.operations : isClose;
586     import std.math.constants : PI;
587 
588     assert(isClose(acos(0.5), PI / 3, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
589 }
590 
591 /***************
592  * Calculates the arc sine of x,
593  * returning a value ranging from -$(PI)/2 to $(PI)/2.
594  *
595  *      $(TABLE_SV
596  *      $(TR $(TH x)            $(TH asin(x))      $(TH invalid?))
597  *      $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
598  *      $(TR $(TD $(GT)1.0)     $(TD $(NAN))       $(TD yes))
599  *      $(TR $(TD $(LT)-1.0)    $(TD $(NAN))       $(TD yes))
600  *  )
601  */
602 real asin(real x) @safe pure nothrow @nogc
603 {
604     import core.math : sqrt;
605 
606     return atan2(x, sqrt(1-x*x));
607 }
608 
609 /// ditto
610 double asin(double x) @safe pure nothrow @nogc { return asin(cast(real) x); }
611 
612 /// ditto
613 float asin(float x) @safe pure nothrow @nogc  { return asin(cast(real) x); }
614 
615 ///
616 @safe unittest
617 {
618     import std.math.operations : isClose;
619     import std.math.traits : isIdentical, isNaN;
620     import std.math.constants : PI;
621 
622     assert(isIdentical(asin(0.0), 0.0));
623     assert(asin(0.5).isClose(PI / 6));
624     assert(asin(PI).isNaN);
625 }
626 
627 @safe @nogc nothrow unittest
628 {
629     import std.math.operations : isClose;
630     import std.math.constants : PI;
631 
632     assert(isClose(asin(0.5), PI / 6, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
633 }
634 
635 /***************
636  * Calculates the arc tangent of x,
637  * returning a value ranging from -$(PI)/2 to $(PI)/2.
638  *
639  *      $(TABLE_SV
640  *      $(TR $(TH x)                 $(TH atan(x))      $(TH invalid?))
641  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no))
642  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))       $(TD yes))
643  *  )
644  */
645 pragma(inline, true)
646 real atan(real x) @safe pure nothrow @nogc
647 {
648     version (InlineAsm_X87)
649     {
650         if (!__ctfe)
651             return atan2Asm(x, 1.0L);
652     }
653     return atanImpl(x);
654 }
655 
656 /// ditto
657 pragma(inline, true)
658 double atan(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) atan(cast(real) x) : atanImpl(x); }
659 
660 /// ditto
661 pragma(inline, true)
662 float atan(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) atan(cast(real) x) : atanImpl(x); }
663 
664 ///
665 @safe unittest
666 {
667     import std.math.operations : isClose;
668     import std.math.traits : isIdentical;
669     import std.math.constants : PI;
670     import std.math.algebraic : sqrt;
671 
672     assert(isIdentical(atan(0.0), 0.0));
673     assert(atan(sqrt(3.0)).isClose(PI / 3));
674 }
675 
676 private T atanImpl(T)(T x) @safe pure nothrow @nogc
677 {
678     import std.math : floatTraits, RealFormat;
679     import std.math.traits : copysign, isInfinity, signbit;
680     import std.math.constants : PI_2, PI_4;
681     import std.math.algebraic : poly;
682 
683     // Coefficients for atan(x)
684     enum realFormat = floatTraits!T.realFormat;
685     static if (realFormat == RealFormat.ieeeQuadruple)
686     {
687         static immutable T[9] P = [
688             -6.880597774405940432145577545328795037141E2L,
689             -2.514829758941713674909996882101723647996E3L,
690             -3.696264445691821235400930243493001671932E3L,
691             -2.792272753241044941703278827346430350236E3L,
692             -1.148164399808514330375280133523543970854E3L,
693             -2.497759878476618348858065206895055957104E2L,
694             -2.548067867495502632615671450650071218995E1L,
695             -8.768423468036849091777415076702113400070E-1L,
696             -6.635810778635296712545011270011752799963E-4L
697         ];
698         static immutable T[9] Q = [
699             2.064179332321782129643673263598686441900E3L,
700             8.782996876218210302516194604424986107121E3L,
701             1.547394317752562611786521896296215170819E4L,
702             1.458510242529987155225086911411015961174E4L,
703             7.928572347062145288093560392463784743935E3L,
704             2.494680540950601626662048893678584497900E3L,
705             4.308348370818927353321556740027020068897E2L,
706             3.566239794444800849656497338030115886153E1L,
707             1.0
708         ];
709     }
710     else static if (realFormat == RealFormat.ieeeExtended)
711     {
712         static immutable T[5] P = [
713            -5.0894116899623603312185E1L,
714            -9.9988763777265819915721E1L,
715            -6.3976888655834347413154E1L,
716            -1.4683508633175792446076E1L,
717            -8.6863818178092187535440E-1L,
718         ];
719         static immutable T[6] Q = [
720             1.5268235069887081006606E2L,
721             3.9157570175111990631099E2L,
722             3.6144079386152023162701E2L,
723             1.4399096122250781605352E2L,
724             2.2981886733594175366172E1L,
725             1.0000000000000000000000E0L,
726         ];
727     }
728     else static if (realFormat == RealFormat.ieeeDouble)
729     {
730         static immutable T[5] P = [
731            -6.485021904942025371773E1L,
732            -1.228866684490136173410E2L,
733            -7.500855792314704667340E1L,
734            -1.615753718733365076637E1L,
735            -8.750608600031904122785E-1L,
736         ];
737         static immutable T[6] Q = [
738             1.945506571482613964425E2L,
739             4.853903996359136964868E2L,
740             4.328810604912902668951E2L,
741             1.650270098316988542046E2L,
742             2.485846490142306297962E1L,
743             1.000000000000000000000E0L,
744         ];
745 
746         enum T MOREBITS = 6.123233995736765886130E-17L;
747     }
748     else static if (realFormat == RealFormat.ieeeSingle)
749     {
750         static immutable T[4] P = [
751            -3.33329491539E-1,
752             1.99777106478E-1,
753            -1.38776856032E-1,
754             8.05374449538E-2,
755         ];
756     }
757     else
758         static assert(0, "no coefficients for atan()");
759 
760     // tan(PI/8)
761     enum T TAN_PI_8 = 0.414213562373095048801688724209698078569672L;
762     // tan(3 * PI/8)
763     enum T TAN3_PI_8 = 2.414213562373095048801688724209698078569672L;
764 
765     // Special cases.
766     if (x == cast(T) 0.0)
767         return x;
768     if (isInfinity(x))
769         return copysign(cast(T) PI_2, x);
770 
771     // Make argument positive but save the sign.
772     bool sign = false;
773     if (signbit(x))
774     {
775         sign = true;
776         x = -x;
777     }
778 
779     static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
780     {
781         short flag = 0;
782         T y;
783         if (x > TAN3_PI_8)
784         {
785             y = PI_2;
786             flag = 1;
787             x = -(1.0 / x);
788         }
789         else if (x <= 0.66)
790         {
791             y = 0.0;
792         }
793         else
794         {
795             y = PI_4;
796             flag = 2;
797             x = (x - 1.0)/(x + 1.0);
798         }
799 
800         T z = x * x;
801         z = z * poly(z, P) / poly(z, Q);
802         z = x * z + x;
803         if (flag == 2)
804             z += 0.5 * MOREBITS;
805         else if (flag == 1)
806             z += MOREBITS;
807         y = y + z;
808     }
809     else
810     {
811         // Range reduction.
812         T y;
813         if (x > TAN3_PI_8)
814         {
815             y = PI_2;
816             x = -((cast(T) 1.0) / x);
817         }
818         else if (x > TAN_PI_8)
819         {
820             y = PI_4;
821             x = (x - cast(T) 1.0)/(x + cast(T) 1.0);
822         }
823         else
824             y = 0.0;
825 
826         // Rational form in x^^2.
827         const T z = x * x;
828         static if (realFormat == RealFormat.ieeeSingle)
829             y += poly(z, P) * z * x + x;
830         else
831             y = y + (poly(z, P) / poly(z, Q)) * z * x + x;
832     }
833 
834     return (sign) ? -y : y;
835 }
836 
837 @safe @nogc nothrow unittest
838 {
839     static void testAtan(T)()
840     {
841         import std.math.operations : CommonDefaultFor, isClose, NaN;
842         import std.math.traits : isIdentical;
843         import std.math.constants : PI_2, PI_4;
844 
845         // ±0
846         const T zero = 0.0;
847         assert(isIdentical(atan(zero), zero));
848         assert(isIdentical(atan(-zero), -zero));
849         // ±∞
850         const T inf = T.infinity;
851         assert(isClose(atan(inf), cast(T) PI_2));
852         assert(isClose(atan(-inf), cast(T) -PI_2));
853         // NaN
854         const T specialNaN = NaN(0x0123L);
855         assert(isIdentical(atan(specialNaN), specialNaN));
856 
857         static immutable T[2][] vals =
858         [
859             // x, atan(x)
860             [ 0.25, 0.244978663126864154172L ],
861             [ 0.5,  0.463647609000806116214L ],
862             [ 1,    PI_4                     ],
863             [ 1.5,  0.982793723247329067985L ],
864             [ 10,   1.471127674303734591852L ],
865         ];
866 
867         foreach (ref val; vals)
868         {
869             T x = val[0];
870             T r = val[1];
871             T a = atan(x);
872 
873             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
874             assert(isClose(r, a, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
875 
876             x = -x;
877             r = -r;
878             a = atan(x);
879             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
880             assert(isClose(r, a, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
881         }
882     }
883 
884     import std.meta : AliasSeq;
885     foreach (T; AliasSeq!(real, double, float))
886         testAtan!T();
887 
888     import std.math.operations : isClose;
889     import std.math.algebraic : sqrt;
890     import std.math.constants : PI;
891     assert(isClose(atan(sqrt(3.0L)), PI / 3, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
892 }
893 
894 /***************
895  * Calculates the arc tangent of y / x,
896  * returning a value ranging from -$(PI) to $(PI).
897  *
898  *      $(TABLE_SV
899  *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
900  *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
901  *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
902  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
903  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
904  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
905  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
906  *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
907  *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
908  *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
909  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
910  *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
911  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
912  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
913  *      )
914  */
915 pragma(inline, true)
916 real atan2(real y, real x) @trusted pure nothrow @nogc // TODO: @safe
917 {
918     version (InlineAsm_X87)
919     {
920         if (!__ctfe)
921             return atan2Asm(y, x);
922     }
923     return atan2Impl(y, x);
924 }
925 
926 /// ditto
927 pragma(inline, true)
928 double atan2(double y, double x) @safe pure nothrow @nogc
929 {
930     return __ctfe ? cast(double) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
931 }
932 
933 /// ditto
934 pragma(inline, true)
935 float atan2(float y, float x) @safe pure nothrow @nogc
936 {
937     return __ctfe ? cast(float) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
938 }
939 
940 ///
941 @safe unittest
942 {
943     import std.math.operations : isClose;
944     import std.math.constants : PI;
945     import std.math.algebraic : sqrt;
946 
947     assert(atan2(1.0, sqrt(3.0)).isClose(PI / 6));
948 }
949 
950 version (InlineAsm_X87)
951 private real atan2Asm(real y, real x) @trusted pure nothrow @nogc
952 {
953     version (Win64)
954     {
955         asm pure nothrow @nogc {
956             naked;
957             fld real ptr [RDX]; // y
958             fld real ptr [RCX]; // x
959             fpatan;
960             ret;
961         }
962     }
963     else
964     {
965         asm pure nothrow @nogc {
966             fld y;
967             fld x;
968             fpatan;
969         }
970     }
971 }
972 
973 private T atan2Impl(T)(T y, T x) @safe pure nothrow @nogc
974 {
975     import std.math.traits : copysign, isInfinity, isNaN, signbit;
976     import std.math.constants : PI, PI_2, PI_4;
977 
978     // Special cases.
979     if (isNaN(x) || isNaN(y))
980         return T.nan;
981     if (y == cast(T) 0.0)
982     {
983         if (x >= 0 && !signbit(x))
984             return copysign(0, y);
985         else
986             return copysign(cast(T) PI, y);
987     }
988     if (x == cast(T) 0.0)
989         return copysign(cast(T) PI_2, y);
990     if (isInfinity(x))
991     {
992         if (signbit(x))
993         {
994             if (isInfinity(y))
995                 return copysign(3 * cast(T) PI_4, y);
996             else
997                 return copysign(cast(T) PI, y);
998         }
999         else
1000         {
1001             if (isInfinity(y))
1002                 return copysign(cast(T) PI_4, y);
1003             else
1004                 return copysign(cast(T) 0.0, y);
1005         }
1006     }
1007     if (isInfinity(y))
1008         return copysign(cast(T) PI_2, y);
1009 
1010     // Call atan and determine the quadrant.
1011     T z = atan(y / x);
1012 
1013     if (signbit(x))
1014     {
1015         if (signbit(y))
1016             z = z - cast(T) PI;
1017         else
1018             z = z + cast(T) PI;
1019     }
1020 
1021     if (z == cast(T) 0.0)
1022         return copysign(z, y);
1023 
1024     return z;
1025 }
1026 
1027 @safe @nogc nothrow unittest
1028 {
1029     static void testAtan2(T)()
1030     {
1031         import std.math.operations : isClose;
1032         import std.math.traits : isIdentical, isNaN;
1033         import std.math.constants : PI, PI_2, PI_4;
1034 
1035         // NaN
1036         const T nan = T.nan;
1037         assert(isNaN(atan2(nan, cast(T) 1)));
1038         assert(isNaN(atan2(cast(T) 1, nan)));
1039 
1040         const T inf = T.infinity;
1041         static immutable T[3][] vals =
1042         [
1043             // y, x, atan2(y, x)
1044 
1045             // ±0
1046             [  0.0,  1.0,  0.0 ],
1047             [ -0.0,  1.0, -0.0 ],
1048             [  0.0,  0.0,  0.0 ],
1049             [ -0.0,  0.0, -0.0 ],
1050             [  0.0, -1.0,  PI ],
1051             [ -0.0, -1.0, -PI ],
1052             [  0.0, -0.0,  PI ],
1053             [ -0.0, -0.0, -PI ],
1054             [  1.0,  0.0,  PI_2 ],
1055             [  1.0, -0.0,  PI_2 ],
1056             [ -1.0,  0.0, -PI_2 ],
1057             [ -1.0, -0.0, -PI_2 ],
1058 
1059             // ±∞
1060             [  1.0,  inf,  0.0 ],
1061             [ -1.0,  inf, -0.0 ],
1062             [  1.0, -inf,  PI ],
1063             [ -1.0, -inf, -PI ],
1064             [  inf,  1.0,  PI_2 ],
1065             [  inf, -1.0,  PI_2 ],
1066             [ -inf,  1.0, -PI_2 ],
1067             [ -inf, -1.0, -PI_2 ],
1068             [  inf,  inf,  PI_4 ],
1069             [ -inf,  inf, -PI_4 ],
1070             [  inf, -inf,  3 * PI_4 ],
1071             [ -inf, -inf, -3 * PI_4 ],
1072 
1073             [  1.0,  1.0,  PI_4 ],
1074             [ -2.0,  2.0, -PI_4 ],
1075             [  3.0, -3.0,  3 * PI_4 ],
1076             [ -4.0, -4.0, -3 * PI_4 ],
1077 
1078             [  0.75,  0.25,   1.2490457723982544258299L ],
1079             [ -0.5,   0.375, -0.9272952180016122324285L ],
1080             [  0.5,  -0.125,  1.8157749899217607734034L ],
1081             [ -0.75, -0.5,   -2.1587989303424641704769L ],
1082         ];
1083 
1084         foreach (ref val; vals)
1085         {
1086             const T y = val[0];
1087             const T x = val[1];
1088             const T r = val[2];
1089             const T a = atan2(y, x);
1090 
1091             //printf("atan2(%Lg, %Lg) = %Lg, should be %Lg\n", cast(real) y, cast(real) x, cast(real) a, cast(real) r);
1092             if (r == 0)
1093                 assert(isIdentical(r, a)); // check sign
1094             else
1095                 assert(isClose(r, a));
1096         }
1097     }
1098 
1099     import std.meta : AliasSeq;
1100     foreach (T; AliasSeq!(real, double, float))
1101         testAtan2!T();
1102 
1103     import std.math.operations : isClose;
1104     import std.math.algebraic : sqrt;
1105     import std.math.constants : PI;
1106     assert(isClose(atan2(1.0L, sqrt(3.0L)), PI / 6, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1107 }
1108 
1109 /***********************************
1110  * Calculates the hyperbolic cosine of x.
1111  *
1112  *      $(TABLE_SV
1113  *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
1114  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
1115  *      )
1116  */
1117 real cosh(real x) @safe pure nothrow @nogc
1118 {
1119     import std.math.exponential : exp;
1120 
1121     //  cosh = (exp(x)+exp(-x))/2.
1122     // The naive implementation works correctly.
1123     const real y = exp(x);
1124     return (y + 1.0/y) * 0.5;
1125 }
1126 
1127 /// ditto
1128 double cosh(double x) @safe pure nothrow @nogc { return cosh(cast(real) x); }
1129 
1130 /// ditto
1131 float cosh(float x) @safe pure nothrow @nogc  { return cosh(cast(real) x); }
1132 
1133 ///
1134 @safe unittest
1135 {
1136     import std.math.constants : E;
1137     import std.math.operations : isClose;
1138 
1139     assert(cosh(0.0) == 1.0);
1140     assert(cosh(1.0).isClose((E + 1.0 / E) / 2));
1141 }
1142 
1143 @safe @nogc nothrow unittest
1144 {
1145     import std.math.constants : E;
1146     import std.math.operations : isClose;
1147 
1148     assert(isClose(cosh(1.0), (E + 1.0 / E) / 2, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1149 }
1150 
1151 /***********************************
1152  * Calculates the hyperbolic sine of x.
1153  *
1154  *      $(TABLE_SV
1155  *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
1156  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
1157  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
1158  *      )
1159  */
1160 real sinh(real x) @safe pure nothrow @nogc { return _sinh(x); }
1161 
1162 /// ditto
1163 double sinh(double x) @safe pure nothrow @nogc { return _sinh(x); }
1164 
1165 /// ditto
1166 float sinh(float x) @safe pure nothrow @nogc { return _sinh(x); }
1167 
1168 ///
1169 @safe unittest
1170 {
1171     import std.math.constants : E;
1172     import std.math.operations : isClose;
1173     import std.math.traits : isIdentical;
1174 
1175     enum sinh1 = (E - 1.0 / E) / 2;
1176     import std.meta : AliasSeq;
1177     static foreach (F; AliasSeq!(float, double, real))
1178     {
1179         assert(isIdentical(sinh(F(0.0)), F(0.0)));
1180         assert(sinh(F(1.0)).isClose(F(sinh1)));
1181     }
1182 }
1183 
1184 private F _sinh(F)(F x)
1185 {
1186     import std.math.traits : copysign;
1187     import std.math.exponential : exp, expm1;
1188     import core.math : fabs;
1189     import std.math.constants : LN2;
1190 
1191     //  sinh(x) =  (exp(x)-exp(-x))/2;
1192     // Very large arguments could cause an overflow, but
1193     // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
1194     // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
1195     if (fabs(x) > F.mant_dig * F(LN2))
1196     {
1197         return copysign(F(0.5) * exp(fabs(x)), x);
1198     }
1199 
1200     const y = expm1(x);
1201     return F(0.5) * y / (y+1) * (y+2);
1202 }
1203 
1204 @safe @nogc nothrow unittest
1205 {
1206     import std.math.constants : E;
1207     import std.math.operations : isClose;
1208 
1209     assert(isClose(sinh(1.0L), real((E - 1.0 / E) / 2), real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1210 }
1211 /***********************************
1212  * Calculates the hyperbolic tangent of x.
1213  *
1214  *      $(TABLE_SV
1215  *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
1216  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
1217  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
1218  *      )
1219  */
1220 real tanh(real x) @safe pure nothrow @nogc { return _tanh(x); }
1221 
1222 /// ditto
1223 double tanh(double x) @safe pure nothrow @nogc { return _tanh(x); }
1224 
1225 /// ditto
1226 float tanh(float x) @safe pure nothrow @nogc { return _tanh(x); }
1227 
1228 ///
1229 @safe unittest
1230 {
1231     import std.math.operations : isClose;
1232     import std.math.traits : isIdentical;
1233 
1234     assert(isIdentical(tanh(0.0), 0.0));
1235     assert(tanh(1.0).isClose(sinh(1.0) / cosh(1.0)));
1236 }
1237 
1238 private F _tanh(F)(F x)
1239 {
1240     import std.math.traits : copysign;
1241     import std.math.exponential : expm1;
1242     import core.math : fabs;
1243     import std.math.constants : LN2;
1244 
1245     //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
1246     if (fabs(x) > F.mant_dig * F(LN2))
1247     {
1248         return copysign(1, x);
1249     }
1250 
1251     const y = expm1(2*x);
1252     return y / (y + 2);
1253 }
1254 
1255 @safe @nogc nothrow unittest
1256 {
1257     import std.math.operations : isClose;
1258 
1259     assert(isClose(tanh(1.0L), sinh(1.0L) / cosh(1.0L), real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1260 }
1261 
1262 /***********************************
1263  * Calculates the inverse hyperbolic cosine of x.
1264  *
1265  *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
1266  *
1267  * $(TABLE_DOMRG
1268  *    $(DOMAIN 1..$(INFIN)),
1269  *    $(RANGE  0..$(INFIN))
1270  * )
1271  *
1272  *  $(TABLE_SV
1273  *    $(SVH  x,     acosh(x) )
1274  *    $(SV  $(NAN), $(NAN) )
1275  *    $(SV  $(LT)1,     $(NAN) )
1276  *    $(SV  1,      0       )
1277  *    $(SV  +$(INFIN),+$(INFIN))
1278  *  )
1279  */
1280 real acosh(real x) @safe pure nothrow @nogc { return _acosh(x); }
1281 
1282 /// ditto
1283 double acosh(double x) @safe pure nothrow @nogc { return _acosh(x); }
1284 
1285 /// ditto
1286 float acosh(float x) @safe pure nothrow @nogc { return _acosh(x); }
1287 
1288 ///
1289 @safe @nogc nothrow unittest
1290 {
1291     import std.math.traits : isIdentical, isNaN;
1292 
1293     assert(isNaN(acosh(0.9)));
1294     assert(isNaN(acosh(real.nan)));
1295     assert(isIdentical(acosh(1.0), 0.0));
1296     assert(acosh(real.infinity) == real.infinity);
1297     assert(isNaN(acosh(0.5)));
1298 }
1299 
1300 private F _acosh(F)(F x) @safe pure nothrow @nogc
1301 {
1302     import std.math.constants : LN2;
1303     import std.math.exponential : log;
1304     import core.math : sqrt;
1305 
1306     if (x > 1/F.epsilon)
1307         return F(LN2) + log(x);
1308     else
1309         return log(x + sqrt(x*x - 1));
1310 }
1311 
1312 @safe @nogc nothrow unittest
1313 {
1314     import std.math.operations : isClose;
1315 
1316     assert(isClose(acosh(cosh(3.0L)), 3.0L, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1317 }
1318 
1319 /***********************************
1320  * Calculates the inverse hyperbolic sine of x.
1321  *
1322  *  Mathematically,
1323  *  ---------------
1324  *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
1325  *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
1326  *  -------------
1327  *
1328  *    $(TABLE_SV
1329  *    $(SVH x,                asinh(x)       )
1330  *    $(SV  $(NAN),           $(NAN)         )
1331  *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
1332  *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
1333  *    )
1334  */
1335 real asinh(real x) @safe pure nothrow @nogc { return _asinh(x); }
1336 
1337 /// ditto
1338 double asinh(double x) @safe pure nothrow @nogc { return _asinh(x); }
1339 
1340 /// ditto
1341 float asinh(float x) @safe pure nothrow @nogc { return _asinh(x); }
1342 
1343 ///
1344 @safe @nogc nothrow unittest
1345 {
1346     import std.math.traits : isIdentical, isNaN;
1347 
1348     assert(isIdentical(asinh(0.0), 0.0));
1349     assert(isIdentical(asinh(-0.0), -0.0));
1350     assert(asinh(real.infinity) == real.infinity);
1351     assert(asinh(-real.infinity) == -real.infinity);
1352     assert(isNaN(asinh(real.nan)));
1353 }
1354 
1355 private F _asinh(F)(F x)
1356 {
1357     import std.math.traits : copysign;
1358     import core.math : fabs, sqrt;
1359     import std.math.exponential : log, log1p;
1360     import std.math.constants : LN2;
1361 
1362     return (fabs(x) > 1 / F.epsilon)
1363         // beyond this point, x*x + 1 == x*x
1364         ? copysign(F(LN2) + log(fabs(x)), x)
1365         // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
1366         : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
1367 }
1368 
1369 @safe unittest
1370 {
1371     import std.math.operations : isClose;
1372 
1373     assert(isClose(asinh(sinh(3.0L)), 3.0L, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1374 }
1375 
1376 /***********************************
1377  * Calculates the inverse hyperbolic tangent of x,
1378  * returning a value from ranging from -1 to 1.
1379  *
1380  * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
1381  *
1382  * $(TABLE_DOMRG
1383  *    $(DOMAIN -$(INFIN)..$(INFIN)),
1384  *    $(RANGE  -1 .. 1)
1385  * )
1386  * $(BR)
1387  * $(TABLE_SV
1388  *    $(SVH  x,     acosh(x) )
1389  *    $(SV  $(NAN), $(NAN) )
1390  *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
1391  *    $(SV  -$(INFIN), -0)
1392  * )
1393  */
1394 real atanh(real x) @safe pure nothrow @nogc
1395 {
1396     import std.math.exponential : log1p;
1397 
1398     // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
1399     return  0.5 * log1p( 2 * x / (1 - x) );
1400 }
1401 
1402 /// ditto
1403 double atanh(double x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
1404 
1405 /// ditto
1406 float atanh(float x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
1407 
1408 ///
1409 @safe @nogc nothrow unittest
1410 {
1411     import std.math.traits : isIdentical, isNaN;
1412 
1413     assert(isIdentical(atanh(0.0), 0.0));
1414     assert(isIdentical(atanh(-0.0),-0.0));
1415     assert(isNaN(atanh(real.nan)));
1416     assert(isNaN(atanh(-real.infinity)));
1417     assert(atanh(0.0) == 0);
1418 }
1419 
1420 @safe unittest
1421 {
1422     import std.math.operations : isClose;
1423 
1424     assert(isClose(atanh(tanh(0.5L)), 0.5, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1425 }