1 /** This module contains the $(LREF Dual) type, which is used to represent
2     dual numbers, along with related mathematical operations and functions.
3 
4     The intention of the this module is to provide a dual number type, which
5     can be used for automatic differentation. This means you can calculate
6     with it the same way you would with floating point numbers. The real
7     part will be identically to the result of a calculation with floating
8     point numbers. Comparison operations (opCmp) ignore the value of the dual part.
9 
10     To calculate the derivate with respect to a variable, set the dual part
11     of this variable to 1. All other dual parts should be 0. After the
12     calculation the dual part contains the derivate.
13 
14     See_Also:
15         https://en.wikipedia.org/wiki/Dual_number
16         https://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers
17 
18     Current_limitations:
19         "^^" operator works only for integral numbers.
20 
21     Authors:    René Heldmaier and others (see individual files)
22     Copyright:  Copyright (c) 2019, René Heldmaier and others (see individual files)
23     License:    $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0)
24     Source:     https://github.com/the5avage/dualnumbers
25 */
26 
27 /* The functions
28 
29    auto dual(R)(const R re),
30    auto dual(R, D)(const R re, const D du),
31    string toString() const,
32    void toString(Writer, Char)(scope Writer w, scope const ref FormatSpec!Char formatSpec) const
33 
34    and the corresponding unittests are derived from phobos std.complex
35    (github.com/dlang/phobos/blob/master/std/complex.d).
36    See comments for more details
37 */
38 
39 module dualnumbers;
40 
41 import std.traits;
42 import std.math;
43 
44 /* The following two functions and the corresponding unittest are derived from
45    phobos std.complex (github.com/dlang/phobos/blob/master/std/complex.d)
46    Authors:     Lars Tandle Kyllingstad, Don Clugston
47    Copyright:   Copyright (c) 2010, Lars T. Kyllingstad.
48    License:     $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0)
49    Modified by: René Heldmaier
50 */
51 /** Helper function that returns a dual number with the specified
52     real and dual parts.
53     Params:
54         R = (template parameter) type of real part of dual number
55         D = (template parameter) type of dual part of dual number
56         re = real part of complex number to be constructed
57         du = (optional) dual part of complex number, 0 if omitted.
58     Returns:
59         `Dual` instance with real and dual parts set
60         to the values provided as input.  If neither `re` nor
61         `du` are floating-point numbers, the return type will
62         be `Dual!double`.  Otherwise, the return type is
63         deduced using $(D std.traits.CommonType!(R, D)).
64 */
65 auto dual(R)(const R re)  @safe pure nothrow @nogc
66 if (is(R : double))
67 {
68     static if (isFloatingPoint!R)
69         return Dual!R(re, 0);
70     else
71         return Dual!double(re, 0);
72 }
73 
74 /// ditto
75 auto dual(R, D)(const R re, const D du)  @safe pure nothrow @nogc
76 if (is(R : double) && is(D : double))
77 {
78     static if (isFloatingPoint!R || isFloatingPoint!D)
79         return Dual!(CommonType!(R, D))(re, du);
80     else
81         return Dual!double(re, du);
82 }
83 
84 ///
85 @safe pure nothrow unittest
86 {
87     auto a = dual(1.0);
88     static assert(is(typeof(a) == Dual!double));
89     assert(a.re == 1.0);
90     assert(a.du == 0.0);
91 
92     auto b = dual(2.0L);
93     static assert(is(typeof(b) == Dual!real));
94     assert(b.re == 2.0L);
95     assert(b.du == 0.0L);
96 
97     auto c = dual(1.0, 2.0);
98     static assert(is(typeof(c) == Dual!double));
99     assert(c.re == 1.0);
100     assert(c.du == 2.0);
101 
102     auto d = dual(3.0, 4.0L);
103     static assert(is(typeof(d) == Dual!real));
104     assert(d.re == 3.0);
105     assert(d.du == 4.0L);
106 
107     auto e = dual(1);
108     static assert(is(typeof(e) == Dual!double));
109     assert(e.re == 1);
110     assert(e.du == 0);
111 
112     auto f = dual(1L, 2);
113     static assert(is(typeof(f) == Dual!double));
114     assert(f.re == 1L);
115     assert(f.du == 2);
116 
117     auto g = dual(3, 4.0L);
118     static assert(is(typeof(g) == Dual!real));
119     assert(g.re == 3);
120     assert(g.du == 4.0L);
121 }
122 
123 /// Dual number for automatic differentation
124 struct Dual(T) if (isFloatingPoint!T)
125 {
126     import std.format : FormatSpec;
127     import std.range.primitives : isOutputRange;
128 
129     /// The real part
130     T re;
131     /// The dual part
132     T du;
133 
134     /* The following two functions and the corresponding unittest are derived from
135        phobos std.complex (github.com/dlang/phobos/blob/master/std/complex.d)
136        Authors:     Lars Tandle Kyllingstad, Don Clugston
137        Copyright:   Copyright (c) 2010, Lars T. Kyllingstad.
138        License:     $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0)
139        Modified by: René Heldmaier
140     */
141     /** Converts the dual number to a string representation.
142         The second form of this function is usually not called directly;
143         instead, it is used via $(REF format, std,string), as shown in the examples
144         below.  Supported format characters are 'e', 'f', 'g', 'a', and 's'.
145         See the $(MREF std, format) and $(REF format, std,string)
146         documentation for more information.
147     */
148     string toString() const @safe /* TODO: pure nothrow */
149     {
150             import std.exception : assumeUnique;
151             char[] buf;
152             buf.reserve(100);
153             auto fmt = FormatSpec!char("%s");
154             toString((const(char)[] s) { buf ~= s; }, fmt);
155             static trustedAssumeUnique(T)(T t) @trusted { return assumeUnique(t); }
156             return trustedAssumeUnique(buf);
157     }
158 
159     ///
160     @safe unittest
161     {
162             auto c = dual(1.2, 3.4);
163 
164             // Vanilla toString formatting:
165             assert(c.toString() == "1.2+3.4ε");
166 
167             // Formatting with std.string.format specs: the precision and width
168             // specifiers apply to both the real and imaginary parts of the
169             // complex number.
170             import std.format : format;
171             assert(format("%.2f", c)  == "1.20+3.40ε");
172             assert(format("%4.1f", c) == " 1.2+ 3.4ε");
173     }
174 
175     /// ditto
176     void toString(Writer, Char)(scope Writer w, scope const ref FormatSpec!Char formatSpec) const
177         if (isOutputRange!(Writer, const(Char)[]))
178         {
179             import std.format : formatValue;
180             import std.math : signbit;
181             import std.range.primitives : put;
182             formatValue(w, re, formatSpec);
183             if (signbit(du) == 0)
184                 put(w, "+");
185             formatValue(w, du, formatSpec);
186             put(w, "ε");
187     }
188 
189     @safe pure nothrow @nogc:
190 
191     /**
192     Construct a dual number with the specified real and
193     dual parts. In the case where a single argument is passed
194     that is not a dual number, the imaginary part of the result will be
195     zero.
196     */
197     this(R : T)(const Dual!R d)
198     {
199         re = d.re;
200         du = d.du;
201     }
202 
203     /// ditto
204     this(R : T)(const R re)
205     {
206         this.re = re;
207         du = 0.0;
208     }
209 
210     /// ditto
211     this(Rr : T, Rd : T)(const Rr re, const Rd du)
212     {
213         this.re = re;
214         this.du = du;
215     }
216 
217     /* assignment operators */
218 
219     // this = numeric
220     ref Dual opAssign(R : T)(const R rhs)
221     {
222         this.re = rhs;
223         this.du = 0.0;
224         return this;
225     }
226 
227     // this = dual
228     ref Dual opAssign(R : T)(const Dual!R rhs)
229     {
230         re = rhs.re;
231         du = rhs.du;
232         return this;
233     }
234 
235     /* unary operators */
236 
237     // unary + and -
238     Dual opUnary(string op)() const
239         if (op == "+" || op == "-")
240     {
241         return mixin("Dual(" ~op~ "re," ~op~ "du)");
242     }
243 
244     /* binary operators */
245 
246     // dual + dual, dual - dual
247     Dual!(CommonType!(T, R)) opBinary(string op, R)(const Dual!R rhs) const
248         if (op == "+" || op == "-")
249     {
250         alias D = typeof(return);
251         return mixin("D(re" ~ op ~ "rhs.re, du" ~ op ~ "rhs.du)");
252     }
253 
254     // dual * dual
255     Dual!(CommonType!(T, R)) opBinary(string op, R)(const Dual!R rhs) const
256         if (op == "*")
257     {
258         alias D = typeof(return);
259         return D(re * rhs.re, du * rhs.re + re * rhs.du);
260     }
261 
262     // dual / dual
263     Dual!(CommonType!(T, R)) opBinary(string op, R)(const Dual!R rhs) const
264         if (op == "/")
265     {
266         alias D = typeof(return);
267         return D(re/rhs.re, (rhs.re*du - rhs.du*re)/rhs.re/rhs.re);
268     }
269 
270     // dual + numeric, dual - numeric
271     Dual!(CommonType!(T, R)) opBinary(string op, R)(const R rhs) const
272         if ((op == "+" || op == "-") && isNumeric!R)
273     {
274         alias D = typeof(return);
275         return mixin("D(re" ~op~ "rhs, du)");
276     }
277 
278     // dual * numeric, dual / numeric
279     Dual!(CommonType!(T, R)) opBinary(string op, R)(const R rhs) const
280         if ((op == "*" || op == "/") && isNumeric!R)
281     {
282         alias D = typeof(return);
283         return mixin("D(re" ~op~ "rhs, du" ~op~ "rhs)");
284     }
285 
286     // numeric * dual, numeric + dual
287     Dual!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const
288         if ((op == "+" || op == "*") && isNumeric!R)
289     {
290         return opBinary!op(lhs);
291     }
292 
293     // numeric - dual
294     Dual!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const
295         if ((op == "-") && isNumeric!R)
296     {
297         alias D = typeof(return);
298         return D(lhs - re, -du);
299     }
300 
301     // numeric / dual
302     Dual!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const
303         if ((op == "/") && isNumeric!R)
304     {
305         alias D = typeof(return);
306         return D(lhs/re, -du*lhs/re/re);
307     }
308 
309     // dual ^^ integer
310     Dual!T opBinary(string op, R)(R rhs) const
311         if (op == "^^" && isIntegral!R)
312     {
313         if (rhs > 0) {
314             Dual!T result = this;
315             while (--rhs) {
316                 result *= this;
317             }
318             return result;
319         } else if (rhs < 0) {
320             auto result = Dual!T(1.0);
321             while (rhs++) {
322                 result /= this;
323             }
324             return result;
325         }
326         return Dual!T(1.0);
327     }
328 
329     /* op-assign operators */
330 
331     // dual += dual, dual -= dual, dual *= dual, dual /= dual
332     ref Dual opOpAssign(string op, D)(const D rhs)
333         if ((op == "+" || op == "-" || op == "*" || op == "/") && is (D R == Dual!R))
334     {
335         mixin("this = this " ~op~ "rhs;");
336         return this;
337     }
338 
339     // dual ^^= numeric
340     ref Dual opOpAssign(string op, R : T)(R rhs)
341         if (op == "^^")
342     {
343             this = this ^^ rhs;
344             return this;
345     }
346 
347     // dual += numeric, dual -= numeric, dual *= numeric, dual /= numeric
348     ref Dual opOpAssign(string op, R : T)(const R rhs)
349         if (op == "+" || op == "-" || op == "*" || op =="/")
350     {
351         mixin("this = this " ~op~ "rhs;");
352         return this;
353     }
354 
355     /* comparison operators
356     Ignore dual part because dual numbers shall behave identically
357     to normal numbers (they are made for automatic differentation).
358     */
359 
360     // this == dual
361     bool opEquals(R : T)(const Dual!R z) const
362     {
363         return re == z.re && du == z.du;
364     }
365 
366     // this == numeric
367     bool opEquals(R : T)(const R r) const
368     {
369         return re == r && du == 0;
370     }
371 
372     // this == dual
373     int opCmp(R : T)(const Dual!R z) const
374     {
375         if (re < z.re)
376             return -1;
377         else if (re > z.re)
378             return 1;
379         else
380             return 0;
381     }
382 
383     // this == dual
384     int opCmp(R : T)(const R r) const
385     {
386             if (re < r)
387                 return -1;
388             else if (re > r)
389                 return 1;
390             else
391                 return 0;
392     }
393 }
394 
395 /** To calculate the derivate of a function, set the dual part of the variable
396     by which you want to derive to 1. Other dual numbers should have dual part of 0.
397     Then just calculate the function like you would with normal numbers.
398     The real part is always identically to the result with normal numbers.
399     After the calculation the dual part holds the derivate of the function.
400 */
401 unittest
402 {
403     import std.math;
404     // f(x) = x⁵ f'(x) = 5x⁴ for x = 2
405     const x = Dual!double(2.0, 1.0);
406     auto f2 = x^^5;
407     assert(approxEqual(f2.re, 2.0^^5));
408     assert(approxEqual(f2.du, 5.0 * 2.0^^4));
409 
410     // f(x) = 3x² f'(x) = 6x for x = 2
411     f2 = 3.0 * x * x;
412     assert(approxEqual(f2.re, 3.0 * 2.0 * 2.0));
413     assert(approxEqual(f2.du, 6.0 * 2.0));
414 
415     // f(x) = 3/(1-x) f'(x) = 3/(1-x)^2 for x = 2
416     f2 = 3.0/(1.0 - x);
417     assert(approxEqual(f2.re, -3.0));
418     assert(approxEqual(f2.du, 3.0));
419 
420     // f(x) = 3exp(2x), f'(x) = 6exp(2x) for x = 2
421     f2 = 3 * exp(2 * x);
422     assert(approxEqual(f2.re, 3 * std.math.exp(4.0)));
423     assert(approxEqual(f2.du, 6 * std.math.exp(4.0)));
424 }
425 
426 /** Exponential function on dual numbers.
427 
428     Params: x = A dual number.
429     Returns: exp(x)
430 */
431 Dual!T exp(T)(Dual!T x)  @safe pure nothrow @nogc
432 {
433     import std.math: exp;
434     Dual!T result = void;
435     result.re = exp(x.re);
436     result.du = result.re * x.du;
437     return result;
438 }
439 
440 ///
441 unittest
442 {
443     import std.math;
444     // f(x) = exp(x), f'(x) = exp(x)
445     auto x = dual(5.0, 1.0);
446     auto res = x.exp();
447     assert(approxEqual(res.re, std.math.exp(5.0)));
448     assert(approxEqual(res.du, std.math.exp(5.0)));
449 
450     // f(x) = exp(3x), f'(x) = 3*exp(3x)
451     res = exp(3 * x);
452     assert(approxEqual(res.re, std.math.exp(15.0)));
453     assert(approxEqual(res.du, 3 * std.math.exp(15.0)));
454 }
455 
456 /** Abs function on dual numbers.
457 
458     Params: x = A dual number.
459     Returns: |x|
460 */
461 Dual!T abs(T)(Dual!T x)  @safe pure nothrow @nogc
462 {
463     import std.math: signbit;
464     if (signbit(x.re))
465         return -x;
466     else
467         return x;
468 }
469 
470 ///
471 unittest
472 {
473     import std.math: approxEqual;
474     // f(x) = |x|, f'(x) = signum(x)
475     auto x = dual(2.0, 1.0);
476     auto result = abs(x);
477     assert(approxEqual(result.re, 2.0)); // f(2) = 2
478     assert(approxEqual(result.du, 1.0)); // f'(2) = 1
479 
480     x = dual(-2.0, 1.0);
481     result = abs(x);
482     assert(approxEqual(result.re, 2.0)); // f(-2) = 2
483     assert(approxEqual(result.du, -1.0)); // f'(-2) = -1
484 
485     // because floating point numbers have -0 and +0 f'(x) is defined for x = 0
486     x = dual(0.0, 1.0);
487     result = abs(x);
488     assert(approxEqual(result.du, 1.0));
489     x = dual(-0.0, 1.0);
490     result = abs(x);
491     assert(approxEqual(result.du, -1.0));
492 }
493 
494 /* Trigonometric functions */
495 
496 /** Sinus function on dual numbers.
497 
498     Params: x = A dual number.
499     Returns: sin(x)
500 */
501 Dual!T sin(T)(Dual!T x)  @safe pure nothrow @nogc
502 {
503     import std.math: sin, cos;
504     Dual!T result = void;
505     result.re = sin(x.re);
506     result.du = x.du * cos(x.re);
507     return result;
508 }
509 
510 ///
511 unittest
512 {
513     import std.math;
514     // f(x) = sin(x), f'(x) = cos(x)
515     auto x = dual(0.0, 1.0);
516     auto result = sin(x);
517     assert(approxEqual(result.re, 0.0)); // sin(0) = 0
518     assert(approxEqual(result.du, 1.0)); // cos(0) = 1
519 
520     x = dual(PI / 2.0, 1.0);
521     result = sin(x);
522     assert(approxEqual(result.re, 1.0)); // sin(π/2) = 1
523     assert(approxEqual(result.du, 0.0)); // cos(π/2) = 0
524 
525     // f(x) = sin(3x), f'(x) = 3cos(3x)
526     x = dual(5.0, 1.0);
527     result = sin(3 * x);
528     assert(approxEqual(result.re, std.math.sin(15.0)));
529     assert(approxEqual(result.du, 3.0 * std.math.cos(15.0)));
530 }
531 
532 /** Cosinus function on dual numbers.
533 
534     Params: x = A dual number.
535     Returns: cos(x)
536 */
537 Dual!T cos(T)(Dual!T x)  @safe pure nothrow @nogc
538 {
539     import std.math: sin, cos;
540     Dual!T result = void;
541     result.re = cos(x.re);
542     result.du = -x.du * sin(x.re);
543     return result;
544 }
545 
546 ///
547 unittest
548 {
549     import std.math;
550     // f(x) = cos(x), f'(x) = -sin(x)
551     auto x = dual(0.0, 1.0);
552     auto result = cos(x);
553     assert(approxEqual(result.re, 1.0)); // cos(0) = 1
554     assert(approxEqual(result.du, 0.0)); // -sin(0) = 0
555 
556     x = dual(PI / 2.0, 1.0);
557     result = cos(x);
558     assert(approxEqual(result.re, 0.0)); // cos(π/2) = 0
559     assert(approxEqual(result.du, -1.0)); // -sin(π/2) = -1
560 
561     // f(x) = cos(3x), f'(x) = -3sin(3x)
562     x = dual(5.0, 1.0);
563     result = cos(3 * x);
564     assert(approxEqual(result.re, std.math.cos(15.0)));
565     assert(approxEqual(result.du, -3.0 * std.math.sin(15.0)));
566 }
567 
568 /* Test constructors */
569 
570 // Works with floating point types
571 unittest
572 {
573     Dual!float df;
574     Dual!double dd;
575     Dual!real dr;
576     assert(df.sizeof == 2*float.sizeof);
577     assert(dd.sizeof == 2*double.sizeof);
578     assert(dr.sizeof == 2*real.sizeof);
579 }
580 
581 // Constructor with one parameter
582 unittest
583 {
584     import std.math: approxEqual;
585     Dual!double d = Dual!double(33.3);
586     assert(d.re == 33.3 && d.du == 0.0);
587     Dual!real dr = Dual!real(44.4);
588     assert(dr.re.approxEqual(44.4) && dr.du.approxEqual(0.0));
589 }
590 
591 // Constructor with two parameters
592 unittest
593 {
594     import std.math: approxEqual;
595     Dual!double d = Dual!double(1.1, 2.2);
596     assert(d.re.approxEqual(1.1) && d.du.approxEqual(2.2));
597 }
598 
599 // Construct from other dual
600 unittest
601 {
602     import std.math: approxEqual;
603     const d = Dual!double(1.0, 2.0);
604     Dual!float f = d;
605     assert(d.re.approxEqual(f.re) && d.du.approxEqual(f.du));
606 }
607 
608 /* Test operators */
609 
610 // dual = numeric
611 unittest
612 {
613     import std.math: approxEqual;
614     immutable f = 4.0;
615     auto d = Dual!double(0.9, 0.9);
616     auto d2 = d = f;
617     assert(d2.re.approxEqual(4.0) && d2.du.approxEqual(0.0));
618 }
619 
620 // dual = dual
621 unittest
622 {
623     auto d = Dual!double(1.1, 1.1);
624     const  d2 = Dual!double(2.2, 3.3);
625     d = d2;
626     assert(d.re == 2.2 && d.du == 3.3);
627 }
628 
629 // Chain assignment
630 unittest
631 {
632     import std.math: approxEqual;
633     Dual!double dd = Dual!double(1.1, 2.2);
634     Dual!float df;
635     df = dd = Dual!real(1.1, 2.2);
636     assert(df.re.approxEqual(1.1) && df.du.approxEqual(2.2));
637 }
638 
639 // unary + and -
640 unittest
641 {
642     import std.math: approxEqual;
643     const d = Dual!double(1.1, 2.2);
644     Dual!double dp = +d;
645     assert(approxEqual(dp.re, 1.1) && approxEqual(dp.du, 2.2));
646     Dual!double dm = -d;
647     assert(approxEqual(dm.re, -1.1) && approxEqual(dm.du, -2.2));
648 }
649 
650 // dual + dual
651 unittest
652 {
653     import std.math: approxEqual;
654     auto d = Dual!double(1.1, 2.2);
655     const d2 = Dual!float(3.3, 4.4);
656     d = d2 + d;
657     assert(approxEqual(d.re, 4.4) && approxEqual(d.du, 6.6));
658 }
659 
660 // dual - dual
661 unittest
662 {
663     import std.math: approxEqual;
664     auto d = Dual!double(1.1, 2.2);
665     const d2 = Dual!float(3.3, 5.5);
666     d = d2 - d;
667     assert(approxEqual(d.re, 2.2) && approxEqual(d.du, 3.3));
668 }
669 
670 // dual + numeric, numeric + dual
671 unittest
672 {
673     import std.math: approxEqual;
674     const d = Dual!double(1.1, 2.2);
675     auto res = d + 5.0;
676     assert(approxEqual(res.re, 6.1) && approxEqual(res.du, 2.2));
677     res = 1.0 + d;
678     assert(approxEqual(res.re, 2.1) && approxEqual(res.du, 2.2));
679 }
680 
681 // dual - numeric, numeric - dual
682 unittest
683 {
684     import std.math: approxEqual;
685     const d = Dual!double(1.1, 2.2);
686     auto res = d - 5.0;
687     assert(approxEqual(res.re, -3.9) && approxEqual(res.du, 2.2));
688     res = 1.0 - d;
689     assert(approxEqual(res.re, -0.1) && approxEqual(res.du, -2.2));
690 }
691 
692 // dual * dual
693 unittest
694 {
695     import std.math: approxEqual;
696     const a = Dual!double(1.1, 2.2);
697     const b = Dual!float(3.3, 4.4);
698     Dual!double d = a * b;
699     assert(approxEqual(d.re, 1.1 * 3.3) && approxEqual(d.du, 1.1*4.4 + 2.2*3.3));
700 }
701 
702 // dual / dual
703 unittest
704 {
705     import std.math: approxEqual;
706     const a = Dual!double(1.1, 2.2);
707     const b = Dual!float(3.3, 4.4);
708     Dual!double d = a / b;
709     assert(approxEqual(d.re, 1.1 / 3.3) && approxEqual(d.du, (3.3*2.2 - 1.1*4.4)/3.3/3.3));
710 }
711 
712 // dual * numeric
713 unittest
714 {
715     import std.math: approxEqual;
716     Dual!double d = Dual!double(1.1, 2.2);
717     d = d * 3.0;
718     assert(approxEqual(d.re, 3.3) && approxEqual(d.du, 6.6));
719 }
720 
721 // numeric * dual
722 unittest
723 {
724     import std.math: approxEqual;
725     const d = Dual!double(1.1, 2.2);
726     const res = 3.0 * d;
727     assert(approxEqual(res.re, 3.3) && approxEqual(res.du, 6.6));
728 }
729 
730 // dual / numeric
731 unittest
732 {
733     import std.math: approxEqual;
734     const d = Dual!double(1.1, 2.2);
735     const res = d / 3.0;
736     assert(approxEqual(res.re, 1.1/3.0) && approxEqual(res.du, 2.2/3.0));
737 }
738 
739 // numeric / dual
740 unittest
741 {
742     import std.math: approxEqual;
743     const d = Dual!double(1.1, 2.2);
744     Dual!double expect = Dual!double(3.0, 0.0);
745     Dual!double res = 3.0 / d;
746     expect = expect / d;
747     assert(approxEqual(res.re, expect.re) && approxEqual(res.du, expect.du));
748 }
749 
750 // dual += dual
751 unittest
752 {
753     import std.math: approxEqual;
754     const d = Dual!double(1.1, 2.2);
755     Dual!double res = d;
756     res += d;
757     assert(res.re.approxEqual(2.2) && res.du.approxEqual(4.4));
758 }
759 
760 // dual -= dual
761 unittest
762 {
763     import std.math: approxEqual;
764     const d = Dual!double(1.1, 2.2);
765     auto res = Dual!double(0.5, 0.7);
766     res -= d;
767     assert(res.re.approxEqual(-0.6) && res.du.approxEqual(-1.5));
768 }
769 
770 // dual *= dual
771 unittest
772 {
773     import std.math: approxEqual;
774     const d = Dual!double(1.1, 2.2);
775     auto res = Dual!double(0.5, 0.7);
776     auto expect = res * d;
777     res *= d;
778     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
779 }
780 
781 // dual /= dual
782 unittest
783 {
784     import std.math: approxEqual;
785     immutable d = Dual!double(1.1, 2.2);
786     auto res = Dual!double(0.5, 0.7);
787     auto expect = res / d;
788     res /= d;
789     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
790 }
791 
792 // dual += numeric
793 unittest
794 {
795     import std.math: approxEqual;
796     auto res = Dual!double(0.5, 0.7);
797     res += 5.0;
798     assert(res.re.approxEqual(5.5) && res.du.approxEqual(0.7));
799 }
800 
801 // dual -= numeric
802 unittest
803 {
804     import std.math: approxEqual;
805     auto res = Dual!double(0.5, 0.7);
806     res -= 5.0;
807     assert(res.re.approxEqual(-4.5) && res.du.approxEqual(0.7));
808 }
809 
810 // dual *= numeric
811 unittest
812 {
813     import std.math: approxEqual;
814     auto res = Dual!double(0.5, 0.7);
815     auto expect = res * 5.0;
816     res *= 5.0;
817     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
818 }
819 
820 // dual /= numeric
821 unittest
822 {
823     import std.math: approxEqual;
824     auto res = Dual!double(0.5, 0.7);
825     auto expect = res / 5.0;
826     res /= 5.0;
827     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
828 }
829 
830 // dual ^^= numeric
831 unittest
832 {
833     import std.math: approxEqual;
834     auto res = Dual!double(0.5, 0.7);
835     auto expect = res ^^ 5;
836     res ^^= 5;
837     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
838 }
839 
840 // dual ^^ positive integer
841 unittest
842 {
843     import std.math: approxEqual;
844     long exponent = 4;
845     const d = Dual!double(2.2, 3.3);
846     const expect = d * d * d * d;
847     const res = d ^^ exponent;
848     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
849 }
850 
851 // dual ^^ 0
852 unittest
853 {
854     import std.math: approxEqual;
855     short exponent = 0;
856     const d = Dual!double(3.3, 2.2);
857     const res = d ^^ exponent;
858     assert(res.re.approxEqual(1.0) && res.du.approxEqual(0.0));
859 }
860 
861 // dual ^^ negative integer
862 unittest
863 {
864     import std.math: approxEqual;
865     int exponent = -4;
866     const d = Dual!double(2.2, 3.3);
867     const expect = 1.0 / d / d / d / d;
868     const res = d ^^ exponent;
869     assert(res.re.approxEqual(expect.re) && res.du.approxEqual(expect.du));
870 }
871 
872 // opEquals
873 unittest
874 {
875     const a = Dual!double(1.1, 2.2);
876     const b = Dual!double(1.1, 3.3);
877     const c = a;
878     const d = Dual!float(1.1f, 0.0f);
879     assert(a != b && a == c);
880     assert(d == 1.1f && a != 1.1);
881 }
882 
883 // opCmp
884 unittest
885 {
886     const a = Dual!double(1.1, 100.0);
887     const b = Dual!double(2.0, -30.0);
888     assert(a < b && a <= b);
889     assert(b > a && b >= a);
890     assert(a <= a && a >= a);
891     assert(a < 2.0 && a <= 2.0);
892     assert(2.0 > a && 2.0 >= a);
893     assert(a <= 1.1 && a >= 1.1);
894 }