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.traits : 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.traits : floatTraits, RealFormat, copysign, isInfinity, signbit;
679     import std.math.constants : PI_2, PI_4;
680     import std.math.algebraic : poly;
681 
682     // Coefficients for atan(x)
683     enum realFormat = floatTraits!T.realFormat;
684     static if (realFormat == RealFormat.ieeeQuadruple)
685     {
686         static immutable T[9] P = [
687             -6.880597774405940432145577545328795037141E2L,
688             -2.514829758941713674909996882101723647996E3L,
689             -3.696264445691821235400930243493001671932E3L,
690             -2.792272753241044941703278827346430350236E3L,
691             -1.148164399808514330375280133523543970854E3L,
692             -2.497759878476618348858065206895055957104E2L,
693             -2.548067867495502632615671450650071218995E1L,
694             -8.768423468036849091777415076702113400070E-1L,
695             -6.635810778635296712545011270011752799963E-4L
696         ];
697         static immutable T[9] Q = [
698             2.064179332321782129643673263598686441900E3L,
699             8.782996876218210302516194604424986107121E3L,
700             1.547394317752562611786521896296215170819E4L,
701             1.458510242529987155225086911411015961174E4L,
702             7.928572347062145288093560392463784743935E3L,
703             2.494680540950601626662048893678584497900E3L,
704             4.308348370818927353321556740027020068897E2L,
705             3.566239794444800849656497338030115886153E1L,
706             1.0
707         ];
708     }
709     else static if (realFormat == RealFormat.ieeeExtended)
710     {
711         static immutable T[5] P = [
712            -5.0894116899623603312185E1L,
713            -9.9988763777265819915721E1L,
714            -6.3976888655834347413154E1L,
715            -1.4683508633175792446076E1L,
716            -8.6863818178092187535440E-1L,
717         ];
718         static immutable T[6] Q = [
719             1.5268235069887081006606E2L,
720             3.9157570175111990631099E2L,
721             3.6144079386152023162701E2L,
722             1.4399096122250781605352E2L,
723             2.2981886733594175366172E1L,
724             1.0000000000000000000000E0L,
725         ];
726     }
727     else static if (realFormat == RealFormat.ieeeDouble)
728     {
729         static immutable T[5] P = [
730            -6.485021904942025371773E1L,
731            -1.228866684490136173410E2L,
732            -7.500855792314704667340E1L,
733            -1.615753718733365076637E1L,
734            -8.750608600031904122785E-1L,
735         ];
736         static immutable T[6] Q = [
737             1.945506571482613964425E2L,
738             4.853903996359136964868E2L,
739             4.328810604912902668951E2L,
740             1.650270098316988542046E2L,
741             2.485846490142306297962E1L,
742             1.000000000000000000000E0L,
743         ];
744 
745         enum T MOREBITS = 6.123233995736765886130E-17L;
746     }
747     else static if (realFormat == RealFormat.ieeeSingle)
748     {
749         static immutable T[4] P = [
750            -3.33329491539E-1,
751             1.99777106478E-1,
752            -1.38776856032E-1,
753             8.05374449538E-2,
754         ];
755     }
756     else
757         static assert(0, "no coefficients for atan()");
758 
759     // tan(PI/8)
760     enum T TAN_PI_8 = 0.414213562373095048801688724209698078569672L;
761     // tan(3 * PI/8)
762     enum T TAN3_PI_8 = 2.414213562373095048801688724209698078569672L;
763 
764     // Special cases.
765     if (x == cast(T) 0.0)
766         return x;
767     if (isInfinity(x))
768         return copysign(cast(T) PI_2, x);
769 
770     // Make argument positive but save the sign.
771     bool sign = false;
772     if (signbit(x))
773     {
774         sign = true;
775         x = -x;
776     }
777 
778     static if (realFormat == RealFormat.ieeeDouble) // special case for double precision
779     {
780         short flag = 0;
781         T y;
782         if (x > TAN3_PI_8)
783         {
784             y = PI_2;
785             flag = 1;
786             x = -(1.0 / x);
787         }
788         else if (x <= 0.66)
789         {
790             y = 0.0;
791         }
792         else
793         {
794             y = PI_4;
795             flag = 2;
796             x = (x - 1.0)/(x + 1.0);
797         }
798 
799         T z = x * x;
800         z = z * poly(z, P) / poly(z, Q);
801         z = x * z + x;
802         if (flag == 2)
803             z += 0.5 * MOREBITS;
804         else if (flag == 1)
805             z += MOREBITS;
806         y = y + z;
807     }
808     else
809     {
810         // Range reduction.
811         T y;
812         if (x > TAN3_PI_8)
813         {
814             y = PI_2;
815             x = -((cast(T) 1.0) / x);
816         }
817         else if (x > TAN_PI_8)
818         {
819             y = PI_4;
820             x = (x - cast(T) 1.0)/(x + cast(T) 1.0);
821         }
822         else
823             y = 0.0;
824 
825         // Rational form in x^^2.
826         const T z = x * x;
827         static if (realFormat == RealFormat.ieeeSingle)
828             y += poly(z, P) * z * x + x;
829         else
830             y = y + (poly(z, P) / poly(z, Q)) * z * x + x;
831     }
832 
833     return (sign) ? -y : y;
834 }
835 
836 @safe @nogc nothrow unittest
837 {
838     static void testAtan(T)()
839     {
840         import std.math.operations : CommonDefaultFor, isClose, NaN;
841         import std.math.traits : isIdentical;
842         import std.math.constants : PI_2, PI_4;
843 
844         // ±0
845         const T zero = 0.0;
846         assert(isIdentical(atan(zero), zero));
847         assert(isIdentical(atan(-zero), -zero));
848         // ±∞
849         const T inf = T.infinity;
850         assert(isClose(atan(inf), cast(T) PI_2));
851         assert(isClose(atan(-inf), cast(T) -PI_2));
852         // NaN
853         const T specialNaN = NaN(0x0123L);
854         assert(isIdentical(atan(specialNaN), specialNaN));
855 
856         static immutable T[2][] vals =
857         [
858             // x, atan(x)
859             [ 0.25, 0.244978663126864154172L ],
860             [ 0.5,  0.463647609000806116214L ],
861             [ 1,    PI_4                     ],
862             [ 1.5,  0.982793723247329067985L ],
863             [ 10,   1.471127674303734591852L ],
864         ];
865 
866         foreach (ref val; vals)
867         {
868             T x = val[0];
869             T r = val[1];
870             T a = atan(x);
871 
872             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
873             assert(isClose(r, a, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
874 
875             x = -x;
876             r = -r;
877             a = atan(x);
878             //printf("atan(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) a, cast(real) r);
879             assert(isClose(r, a, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T)));
880         }
881     }
882 
883     import std.meta : AliasSeq;
884     foreach (T; AliasSeq!(real, double, float))
885         testAtan!T();
886 
887     import std.math.operations : isClose;
888     import std.math.algebraic : sqrt;
889     import std.math.constants : PI;
890     assert(isClose(atan(sqrt(3.0L)), PI / 3, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
891 }
892 
893 /***************
894  * Calculates the arc tangent of y / x,
895  * returning a value ranging from -$(PI) to $(PI).
896  *
897  *      $(TABLE_SV
898  *      $(TR $(TH y)                 $(TH x)            $(TH atan(y, x)))
899  *      $(TR $(TD $(NAN))            $(TD anything)     $(TD $(NAN)) )
900  *      $(TR $(TD anything)          $(TD $(NAN))       $(TD $(NAN)) )
901  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(GT)0.0)     $(TD $(PLUSMN)0.0) )
902  *      $(TR $(TD $(PLUSMN)0.0)      $(TD +0.0)         $(TD $(PLUSMN)0.0) )
903  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(LT)0.0)     $(TD $(PLUSMN)$(PI)))
904  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -0.0)         $(TD $(PLUSMN)$(PI)))
905  *      $(TR $(TD $(GT)0.0)          $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
906  *      $(TR $(TD $(LT)0.0)          $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
907  *      $(TR $(TD $(GT)0.0)          $(TD $(INFIN))     $(TD $(PLUSMN)0.0) )
908  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything)     $(TD $(PLUSMN)$(PI)/2))
909  *      $(TR $(TD $(GT)0.0)          $(TD -$(INFIN))    $(TD $(PLUSMN)$(PI)) )
910  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN))     $(TD $(PLUSMN)$(PI)/4))
911  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN))    $(TD $(PLUSMN)3$(PI)/4))
912  *      )
913  */
914 pragma(inline, true)
915 real atan2(real y, real x) @trusted pure nothrow @nogc // TODO: @safe
916 {
917     version (InlineAsm_X87)
918     {
919         if (!__ctfe)
920             return atan2Asm(y, x);
921     }
922     return atan2Impl(y, x);
923 }
924 
925 /// ditto
926 pragma(inline, true)
927 double atan2(double y, double x) @safe pure nothrow @nogc
928 {
929     return __ctfe ? cast(double) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
930 }
931 
932 /// ditto
933 pragma(inline, true)
934 float atan2(float y, float x) @safe pure nothrow @nogc
935 {
936     return __ctfe ? cast(float) atan2(cast(real) y, cast(real) x) : atan2Impl(y, x);
937 }
938 
939 ///
940 @safe unittest
941 {
942     import std.math.operations : isClose;
943     import std.math.constants : PI;
944     import std.math.algebraic : sqrt;
945 
946     assert(atan2(1.0, sqrt(3.0)).isClose(PI / 6));
947 }
948 
949 version (InlineAsm_X87)
950 private real atan2Asm(real y, real x) @trusted pure nothrow @nogc
951 {
952     version (Win64)
953     {
954         asm pure nothrow @nogc {
955             naked;
956             fld real ptr [RDX]; // y
957             fld real ptr [RCX]; // x
958             fpatan;
959             ret;
960         }
961     }
962     else
963     {
964         asm pure nothrow @nogc {
965             fld y;
966             fld x;
967             fpatan;
968         }
969     }
970 }
971 
972 private T atan2Impl(T)(T y, T x) @safe pure nothrow @nogc
973 {
974     import std.math.traits : copysign, isInfinity, isNaN, signbit;
975     import std.math.constants : PI, PI_2, PI_4;
976 
977     // Special cases.
978     if (isNaN(x) || isNaN(y))
979         return T.nan;
980     if (y == cast(T) 0.0)
981     {
982         if (x >= 0 && !signbit(x))
983             return copysign(0, y);
984         else
985             return copysign(cast(T) PI, y);
986     }
987     if (x == cast(T) 0.0)
988         return copysign(cast(T) PI_2, y);
989     if (isInfinity(x))
990     {
991         if (signbit(x))
992         {
993             if (isInfinity(y))
994                 return copysign(3 * cast(T) PI_4, y);
995             else
996                 return copysign(cast(T) PI, y);
997         }
998         else
999         {
1000             if (isInfinity(y))
1001                 return copysign(cast(T) PI_4, y);
1002             else
1003                 return copysign(cast(T) 0.0, y);
1004         }
1005     }
1006     if (isInfinity(y))
1007         return copysign(cast(T) PI_2, y);
1008 
1009     // Call atan and determine the quadrant.
1010     T z = atan(y / x);
1011 
1012     if (signbit(x))
1013     {
1014         if (signbit(y))
1015             z = z - cast(T) PI;
1016         else
1017             z = z + cast(T) PI;
1018     }
1019 
1020     if (z == cast(T) 0.0)
1021         return copysign(z, y);
1022 
1023     return z;
1024 }
1025 
1026 @safe @nogc nothrow unittest
1027 {
1028     static void testAtan2(T)()
1029     {
1030         import std.math.operations : isClose;
1031         import std.math.traits : isIdentical, isNaN;
1032         import std.math.constants : PI, PI_2, PI_4;
1033 
1034         // NaN
1035         const T nan = T.nan;
1036         assert(isNaN(atan2(nan, cast(T) 1)));
1037         assert(isNaN(atan2(cast(T) 1, nan)));
1038 
1039         const T inf = T.infinity;
1040         static immutable T[3][] vals =
1041         [
1042             // y, x, atan2(y, x)
1043 
1044             // ±0
1045             [  0.0,  1.0,  0.0 ],
1046             [ -0.0,  1.0, -0.0 ],
1047             [  0.0,  0.0,  0.0 ],
1048             [ -0.0,  0.0, -0.0 ],
1049             [  0.0, -1.0,  PI ],
1050             [ -0.0, -1.0, -PI ],
1051             [  0.0, -0.0,  PI ],
1052             [ -0.0, -0.0, -PI ],
1053             [  1.0,  0.0,  PI_2 ],
1054             [  1.0, -0.0,  PI_2 ],
1055             [ -1.0,  0.0, -PI_2 ],
1056             [ -1.0, -0.0, -PI_2 ],
1057 
1058             // ±∞
1059             [  1.0,  inf,  0.0 ],
1060             [ -1.0,  inf, -0.0 ],
1061             [  1.0, -inf,  PI ],
1062             [ -1.0, -inf, -PI ],
1063             [  inf,  1.0,  PI_2 ],
1064             [  inf, -1.0,  PI_2 ],
1065             [ -inf,  1.0, -PI_2 ],
1066             [ -inf, -1.0, -PI_2 ],
1067             [  inf,  inf,  PI_4 ],
1068             [ -inf,  inf, -PI_4 ],
1069             [  inf, -inf,  3 * PI_4 ],
1070             [ -inf, -inf, -3 * PI_4 ],
1071 
1072             [  1.0,  1.0,  PI_4 ],
1073             [ -2.0,  2.0, -PI_4 ],
1074             [  3.0, -3.0,  3 * PI_4 ],
1075             [ -4.0, -4.0, -3 * PI_4 ],
1076 
1077             [  0.75,  0.25,   1.2490457723982544258299L ],
1078             [ -0.5,   0.375, -0.9272952180016122324285L ],
1079             [  0.5,  -0.125,  1.8157749899217607734034L ],
1080             [ -0.75, -0.5,   -2.1587989303424641704769L ],
1081         ];
1082 
1083         foreach (ref val; vals)
1084         {
1085             const T y = val[0];
1086             const T x = val[1];
1087             const T r = val[2];
1088             const T a = atan2(y, x);
1089 
1090             //printf("atan2(%Lg, %Lg) = %Lg, should be %Lg\n", cast(real) y, cast(real) x, cast(real) a, cast(real) r);
1091             if (r == 0)
1092                 assert(isIdentical(r, a)); // check sign
1093             else
1094                 assert(isClose(r, a));
1095         }
1096     }
1097 
1098     import std.meta : AliasSeq;
1099     foreach (T; AliasSeq!(real, double, float))
1100         testAtan2!T();
1101 
1102     import std.math.operations : isClose;
1103     import std.math.algebraic : sqrt;
1104     import std.math.constants : PI;
1105     assert(isClose(atan2(1.0L, sqrt(3.0L)), PI / 6, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1106 }
1107 
1108 /***********************************
1109  * Calculates the hyperbolic cosine of x.
1110  *
1111  *      $(TABLE_SV
1112  *      $(TR $(TH x)                 $(TH cosh(x))      $(TH invalid?))
1113  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
1114  *      )
1115  */
1116 real cosh(real x) @safe pure nothrow @nogc
1117 {
1118     import std.math.exponential : exp;
1119 
1120     //  cosh = (exp(x)+exp(-x))/2.
1121     // The naive implementation works correctly.
1122     const real y = exp(x);
1123     return (y + 1.0/y) * 0.5;
1124 }
1125 
1126 /// ditto
1127 double cosh(double x) @safe pure nothrow @nogc { return cosh(cast(real) x); }
1128 
1129 /// ditto
1130 float cosh(float x) @safe pure nothrow @nogc  { return cosh(cast(real) x); }
1131 
1132 ///
1133 @safe unittest
1134 {
1135     import std.math.constants : E;
1136     import std.math.operations : isClose;
1137 
1138     assert(cosh(0.0) == 1.0);
1139     assert(cosh(1.0).isClose((E + 1.0 / E) / 2));
1140 }
1141 
1142 @safe @nogc nothrow unittest
1143 {
1144     import std.math.constants : E;
1145     import std.math.operations : isClose;
1146 
1147     assert(isClose(cosh(1.0), (E + 1.0 / E) / 2, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1148 }
1149 
1150 /***********************************
1151  * Calculates the hyperbolic sine of x.
1152  *
1153  *      $(TABLE_SV
1154  *      $(TR $(TH x)                 $(TH sinh(x))           $(TH invalid?))
1155  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0)      $(TD no))
1156  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
1157  *      )
1158  */
1159 real sinh(real x) @safe pure nothrow @nogc { return _sinh(x); }
1160 
1161 /// ditto
1162 double sinh(double x) @safe pure nothrow @nogc { return _sinh(x); }
1163 
1164 /// ditto
1165 float sinh(float x) @safe pure nothrow @nogc { return _sinh(x); }
1166 
1167 ///
1168 @safe unittest
1169 {
1170     import std.math.constants : E;
1171     import std.math.operations : isClose;
1172     import std.math.traits : isIdentical;
1173 
1174     enum sinh1 = (E - 1.0 / E) / 2;
1175     import std.meta : AliasSeq;
1176     static foreach (F; AliasSeq!(float, double, real))
1177     {
1178         assert(isIdentical(sinh(F(0.0)), F(0.0)));
1179         assert(sinh(F(1.0)).isClose(F(sinh1)));
1180     }
1181 }
1182 
1183 private F _sinh(F)(F x)
1184 {
1185     import std.math.traits : copysign;
1186     import std.math.exponential : exp, expm1;
1187     import core.math : fabs;
1188     import std.math.constants : LN2;
1189 
1190     //  sinh(x) =  (exp(x)-exp(-x))/2;
1191     // Very large arguments could cause an overflow, but
1192     // the maximum value of x for which exp(x) + exp(-x)) != exp(x)
1193     // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
1194     if (fabs(x) > F.mant_dig * F(LN2))
1195     {
1196         return copysign(F(0.5) * exp(fabs(x)), x);
1197     }
1198 
1199     const y = expm1(x);
1200     return F(0.5) * y / (y+1) * (y+2);
1201 }
1202 
1203 @safe @nogc nothrow unittest
1204 {
1205     import std.math.constants : E;
1206     import std.math.operations : isClose;
1207 
1208     assert(isClose(sinh(1.0L), real((E - 1.0 / E) / 2), real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1209 }
1210 /***********************************
1211  * Calculates the hyperbolic tangent of x.
1212  *
1213  *      $(TABLE_SV
1214  *      $(TR $(TH x)                 $(TH tanh(x))      $(TH invalid?))
1215  *      $(TR $(TD $(PLUSMN)0.0)      $(TD $(PLUSMN)0.0) $(TD no) )
1216  *      $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
1217  *      )
1218  */
1219 real tanh(real x) @safe pure nothrow @nogc { return _tanh(x); }
1220 
1221 /// ditto
1222 double tanh(double x) @safe pure nothrow @nogc { return _tanh(x); }
1223 
1224 /// ditto
1225 float tanh(float x) @safe pure nothrow @nogc { return _tanh(x); }
1226 
1227 ///
1228 @safe unittest
1229 {
1230     import std.math.operations : isClose;
1231     import std.math.traits : isIdentical;
1232 
1233     assert(isIdentical(tanh(0.0), 0.0));
1234     assert(tanh(1.0).isClose(sinh(1.0) / cosh(1.0)));
1235 }
1236 
1237 private F _tanh(F)(F x)
1238 {
1239     import std.math.traits : copysign;
1240     import std.math.exponential : expm1;
1241     import core.math : fabs;
1242     import std.math.constants : LN2;
1243 
1244     //  tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
1245     if (fabs(x) > F.mant_dig * F(LN2))
1246     {
1247         return copysign(1, x);
1248     }
1249 
1250     const y = expm1(2*x);
1251     return y / (y + 2);
1252 }
1253 
1254 @safe @nogc nothrow unittest
1255 {
1256     import std.math.operations : isClose;
1257 
1258     assert(isClose(tanh(1.0L), sinh(1.0L) / cosh(1.0L), real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1259 }
1260 
1261 /***********************************
1262  * Calculates the inverse hyperbolic cosine of x.
1263  *
1264  *  Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
1265  *
1266  * $(TABLE_DOMRG
1267  *    $(DOMAIN 1..$(INFIN)),
1268  *    $(RANGE  0..$(INFIN))
1269  * )
1270  *
1271  *  $(TABLE_SV
1272  *    $(SVH  x,     acosh(x) )
1273  *    $(SV  $(NAN), $(NAN) )
1274  *    $(SV  $(LT)1,     $(NAN) )
1275  *    $(SV  1,      0       )
1276  *    $(SV  +$(INFIN),+$(INFIN))
1277  *  )
1278  */
1279 real acosh(real x) @safe pure nothrow @nogc { return _acosh(x); }
1280 
1281 /// ditto
1282 double acosh(double x) @safe pure nothrow @nogc { return _acosh(x); }
1283 
1284 /// ditto
1285 float acosh(float x) @safe pure nothrow @nogc { return _acosh(x); }
1286 
1287 ///
1288 @safe @nogc nothrow unittest
1289 {
1290     import std.math.traits : isIdentical, isNaN;
1291 
1292     assert(isNaN(acosh(0.9)));
1293     assert(isNaN(acosh(real.nan)));
1294     assert(isIdentical(acosh(1.0), 0.0));
1295     assert(acosh(real.infinity) == real.infinity);
1296     assert(isNaN(acosh(0.5)));
1297 }
1298 
1299 private F _acosh(F)(F x) @safe pure nothrow @nogc
1300 {
1301     import std.math.constants : LN2;
1302     import std.math.exponential : log;
1303     import core.math : sqrt;
1304 
1305     if (x > 1/F.epsilon)
1306         return F(LN2) + log(x);
1307     else
1308         return log(x + sqrt(x*x - 1));
1309 }
1310 
1311 @safe @nogc nothrow unittest
1312 {
1313     import std.math.operations : isClose;
1314 
1315     assert(isClose(acosh(cosh(3.0L)), 3.0L, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1316 }
1317 
1318 /***********************************
1319  * Calculates the inverse hyperbolic sine of x.
1320  *
1321  *  Mathematically,
1322  *  ---------------
1323  *  asinh(x) =  log( x + sqrt( x*x + 1 )) // if x >= +0
1324  *  asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
1325  *  -------------
1326  *
1327  *    $(TABLE_SV
1328  *    $(SVH x,                asinh(x)       )
1329  *    $(SV  $(NAN),           $(NAN)         )
1330  *    $(SV  $(PLUSMN)0,       $(PLUSMN)0      )
1331  *    $(SV  $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
1332  *    )
1333  */
1334 real asinh(real x) @safe pure nothrow @nogc { return _asinh(x); }
1335 
1336 /// ditto
1337 double asinh(double x) @safe pure nothrow @nogc { return _asinh(x); }
1338 
1339 /// ditto
1340 float asinh(float x) @safe pure nothrow @nogc { return _asinh(x); }
1341 
1342 ///
1343 @safe @nogc nothrow unittest
1344 {
1345     import std.math.traits : isIdentical, isNaN;
1346 
1347     assert(isIdentical(asinh(0.0), 0.0));
1348     assert(isIdentical(asinh(-0.0), -0.0));
1349     assert(asinh(real.infinity) == real.infinity);
1350     assert(asinh(-real.infinity) == -real.infinity);
1351     assert(isNaN(asinh(real.nan)));
1352 }
1353 
1354 private F _asinh(F)(F x)
1355 {
1356     import std.math.traits : copysign;
1357     import core.math : fabs, sqrt;
1358     import std.math.exponential : log, log1p;
1359     import std.math.constants : LN2;
1360 
1361     return (fabs(x) > 1 / F.epsilon)
1362         // beyond this point, x*x + 1 == x*x
1363         ? copysign(F(LN2) + log(fabs(x)), x)
1364         // sqrt(x*x + 1) ==  1 + x * x / ( 1 + sqrt(x*x + 1) )
1365         : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
1366 }
1367 
1368 @safe unittest
1369 {
1370     import std.math.operations : isClose;
1371 
1372     assert(isClose(asinh(sinh(3.0L)), 3.0L, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1373 }
1374 
1375 /***********************************
1376  * Calculates the inverse hyperbolic tangent of x,
1377  * returning a value from ranging from -1 to 1.
1378  *
1379  * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
1380  *
1381  * $(TABLE_DOMRG
1382  *    $(DOMAIN -$(INFIN)..$(INFIN)),
1383  *    $(RANGE  -1 .. 1)
1384  * )
1385  * $(BR)
1386  * $(TABLE_SV
1387  *    $(SVH  x,     acosh(x) )
1388  *    $(SV  $(NAN), $(NAN) )
1389  *    $(SV  $(PLUSMN)0, $(PLUSMN)0)
1390  *    $(SV  -$(INFIN), -0)
1391  * )
1392  */
1393 real atanh(real x) @safe pure nothrow @nogc
1394 {
1395     import std.math.exponential : log1p;
1396 
1397     // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
1398     return  0.5 * log1p( 2 * x / (1 - x) );
1399 }
1400 
1401 /// ditto
1402 double atanh(double x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
1403 
1404 /// ditto
1405 float atanh(float x) @safe pure nothrow @nogc { return atanh(cast(real) x); }
1406 
1407 ///
1408 @safe @nogc nothrow unittest
1409 {
1410     import std.math.traits : isIdentical, isNaN;
1411 
1412     assert(isIdentical(atanh(0.0), 0.0));
1413     assert(isIdentical(atanh(-0.0),-0.0));
1414     assert(isNaN(atanh(real.nan)));
1415     assert(isNaN(atanh(-real.infinity)));
1416     assert(atanh(0.0) == 0);
1417 }
1418 
1419 @safe unittest
1420 {
1421     import std.math.operations : isClose;
1422 
1423     assert(isClose(atanh(tanh(0.5L)), 0.5, real.sizeof > double.sizeof ? 1e-15 : 1e-14));
1424 }