1 /**
2  * Several simple math functions.
3  *
4  * They are rather dedicated to the X86_64 architecture but implemented for
5  * X86 too, even if parameters loading is less adequate and slower under this
6  * architecture.
7  */
8 module iz.math;
9 
10 import
11     std.traits, std.complex;
12 
13 /**
14  * When the conditional version identifer "ensure_sseRoundingMode" is
15  * set then the current rounding mode is reset to the default value.
16  */
17 version(ensure_sseRoundingMode) static this()
18 {
19     roundToNearest;
20 }
21 
22 /**
23  * Saves the rounding mode used in iz.math
24  * (always the 13th and 14th bit of SSE MXCSR).
25  *
26  * Returns:
27  *      The current rounding mode, as a read-only integer value.
28  */
29 const(int) saveIzRoundingMode() @trusted pure nothrow
30 {
31     version(X86) asm pure nothrow
32     {
33         naked;
34         sub     ESP, 4;
35         stmxcsr dword ptr [ESP-4];
36         mov     EAX, dword ptr [ESP-4];
37         add     ESP, 4;
38         ret;
39     }
40     else version(X86_64) asm pure nothrow
41     {
42         naked;
43         sub     RSP, 4;
44         stmxcsr dword ptr [RSP-4];
45         mov     EAX, dword ptr [RSP-4];
46         add     RSP, 4;
47         ret;
48     }
49     else static assert(0, "unsupported architecture");
50 }
51 
52 /**
53  * Restores the rounding mode used in iz.math (MXCSR).
54  *
55  * iz.math provides the functions to set all the possible modes and this function
56  * should only be used to restore a value saved by the roundToXXX functions.
57  * When used freely, this function will only work when the 13th and the 14th bit
58  * modify the current MXCSR register.
59  */
60 void setIzRoundingMode(int value) @trusted pure nothrow
61 {
62     int curr;
63     asm pure nothrow
64     {
65         stmxcsr curr;
66     }
67     import iz.sugar: maskBit;
68     curr = maskBit!13(curr);
69     curr = maskBit!14(curr);
70     curr |= value;
71     asm pure nothrow
72     {
73         ldmxcsr curr;
74     }
75 }
76 
77 /**
78  * Resets the default rounding mode.
79  *
80  * After the call, floor() and ceil() are conform to the specifications.
81  *
82  * Returns:
83  *   The previous rounding mode, which can be restored with setIzRoundingMode().
84  */
85 const(int) roundToNearest() @trusted pure nothrow
86 {
87     const int result = saveIzRoundingMode;
88     import iz.sugar: maskBit;
89     int newMode = result;
90     newMode = maskBit!13(newMode);
91     newMode = maskBit!14(newMode);
92     setIzRoundingMode(newMode);
93     return result;
94 }
95 
96 /**
97  * Sets round() to behave like floor().
98  *
99  * After the call, floor() and ceil() are not anymore reliable.
100  *
101  * Returns:
102  *   The previous rounding mode, which can be restored with setIzRoundingMode().
103  */
104 const(int) roundToPositive() @trusted pure nothrow
105 {
106     const int result = saveIzRoundingMode;
107     import iz.sugar: maskBit;
108     int newMode = result;
109     newMode = maskBit!13(newMode);
110     newMode |= 1 << 14;
111     setIzRoundingMode(newMode);
112     return result;
113 }
114 
115 /**
116  * Sets round() to behave like ceil().
117  *
118  * After the call, floor() and ceil() are not anymore reliable.
119  *
120  * Returns:
121  *   The previous rounding mode, which can be restored with setIzRoundingMode().
122  */
123 const(int) roundToNegative() @trusted pure nothrow
124 {
125     const int result = saveIzRoundingMode;
126     import iz.sugar: maskBit;
127     int newMode = result;
128     newMode = maskBit!14(newMode);
129     newMode |= 1 << 13;
130     setIzRoundingMode(newMode);
131     return result;
132 }
133 
134 /**
135  * Sets round() to behave like trunc().
136  *
137  * After the call, floor() and ceil() are not anymore reliable.
138  *
139  * Returns:
140  *      The previous rounding mode, which can be restored with setIzRoundingMode().
141  */
142 const(int) roundToZero() @trusted pure nothrow
143 {
144     const int result = saveIzRoundingMode;
145     int newMode = result;
146     newMode |= 1 << 13;
147     newMode |= 1 << 14;
148     setIzRoundingMode(newMode);
149     return result;
150 }
151 ///
152 @safe pure nothrow unittest
153 {
154     auto sav = roundToNearest;
155     assert(round(0.4) == 0);
156     auto dnc0 = roundToZero;
157     assert(round(0.8) == 0);
158     assert(round(-0.8) == 0);
159     auto dnc1 = roundToNegative;
160     assert(round(1.8) == 1);
161     assert(round(-0.8) == -1);
162     auto dnc2 = roundToPositive;
163     assert(round(1.8) == 2);
164     assert(round(-0.8) == 0);
165     setIzRoundingMode(sav);
166 }
167 
168 /**
169  * Converts a floating point value to an integer.
170  *
171  * This function requires the SSE2 instruction set and it can't be inlined.
172  *
173  * Params:
174  *      value = Either a float or a double.
175  *
176  * Returns:
177  *      An integer equal to the nearest integral value.
178  */
179 extern(C) int round(T)(T value) @trusted pure nothrow
180 {
181     version(X86_64)
182     {
183         static if (is(T==float)) asm pure nothrow
184         {
185             naked;
186             cvtss2si EAX, XMM0;
187             ret;
188         }
189         else static if (is(T==double)) asm pure nothrow
190         {
191             naked;
192             cvtsd2si EAX, XMM0;
193             ret;
194         }
195         else static assert(0, "unsupported FP type");
196     }
197     else version(X86)
198     {
199         static if (is(T==float)) asm pure nothrow
200         {
201             naked;
202             push EBP;
203             mov EBP, ESP;
204             movss XMM0, value;
205             cvtss2si EAX, XMM0;
206             mov ESP, EBP;
207             pop EBP;
208             ret;
209         }
210         else static if (is(T==double)) asm pure nothrow
211         {
212             naked;
213             push EBP;
214             mov EBP, ESP;
215             movsd XMM0, value;
216             cvtsd2si EAX, XMM0;
217             mov ESP, EBP;
218             pop EBP;
219             ret;
220         }
221         else static assert(0, "unsupported FP type");          
222     }
223     else static assert(0, "unsupported architecture");
224 }
225 ///
226 @safe pure nothrow unittest
227 {
228     assert(round(0.2f) == 0);
229     assert(round(0.8f) == 1);
230     assert(round(-0.2f) == 0);
231     assert(round(-0.8f) == -1);
232 }
233 
234 /**
235  * Converts a floating point value to an integer.
236  *
237  * This function requires the SSE2 instruction set and it can't be inlined.
238  *
239  * Params:
240  *      value = Either a float or a double.
241  *
242  * Returns:
243  *      The largest integral value that is not greater than $(D_PARAM value).
244  */
245 int floor(T)(T value) @trusted pure nothrow
246 if (isFloatingPoint!T)
247 {
248     const T offs = -0.5;
249     return round(value + offs);
250 }
251 ///
252 @safe pure nothrow unittest
253 {
254     assert(floor(0.2f) == 0);
255     assert(floor(0.8f) == 0);
256     assert(floor(-0.2f) == -1);
257     assert(floor(-0.8f) == -1);
258 }
259 
260 /**
261  * Converts a floating point value to an integer.
262  *
263  * This function requires the SSE2 instruction set and it can't be inlined.
264  *
265  * Params:
266  *      value = Either a float or a double.
267  *
268  * Returns:
269  *      The smallest integral value that is not less than $(D_PARAM value).
270  */
271 int ceil(T)(T value) @trusted pure nothrow
272 if (isFloatingPoint!T)
273 {
274     const T offs = 0.5;
275     return round(value + offs);
276 }
277 ///
278 @safe pure nothrow unittest
279 {
280     assert(ceil(0.2f) == 1);
281     assert(ceil(0.8f) == 1);
282     assert(ceil(-0.2f) == 0);
283     assert(ceil(-0.8f) == 0);
284 }
285 
286 /**
287  * Converts a floating point value to an integer.
288  *
289  * This function requires the SSE2 instruction set and it can't be inlined.
290  *
291  * Params:
292  *      value = Either a float or a double.
293  *
294  * Returns:
295  *       An integral value equal to the $(D_PARAM value) nearest integral toward 0.
296  */
297 extern(C) int trunc(T)(T value) @trusted pure nothrow
298 {
299     version(X86_64)
300     {
301         static if (is(T==float)) asm pure nothrow
302         {
303             naked;
304             cvttss2si EAX, XMM0;
305             ret;
306         }
307         else static if (is(T==double)) asm pure nothrow
308         {
309             naked;
310             cvttsd2si EAX, XMM0;
311             ret;
312         }
313         else static assert(0, "unsupported FP type");
314     }
315     else version(X86)
316     {
317         static if (is(T==float)) asm pure nothrow
318         {
319             naked;
320             push EBP;
321             mov EBP, ESP;
322             movss XMM0, value;
323             cvttss2si EAX, XMM0;
324             mov ESP, EBP;
325             pop EBP;
326             ret;
327         }
328         else static if (is(T==double)) asm pure nothrow
329         {
330             naked;
331             push EBP;
332             mov EBP, ESP;
333             movsd XMM0, value;
334             cvttsd2si EAX, XMM0;
335             mov ESP, EBP;
336             pop EBP;
337             ret;
338         }
339         else static assert(0, "unsupported FP type");          
340     }
341     else static assert(0, "unsupported architecture");
342 }
343 ///
344 @safe pure nothrow unittest
345 {
346     assert(trunc(0.2f) == 0);
347     assert(trunc(0.8f) == 0);
348     assert(trunc(-0.2f) == 0);
349     assert(trunc(-8.8f) == -8);
350 }
351 
352 /**
353  * Converts a floating point value to an integer.
354  *
355  * This function relies on the D behavior when casting a floating point value
356  * to an integer. It uses SSE2 on X86_64 and the FPU on X86 but it can always be 
357  * inlined.
358  *
359  * Params:
360  *      value = Either a float or a double.
361  *
362  * Returns:
363  *      An integral value equal to the $(D_PARAM value) nearest integral toward 0.
364  */
365 int dtrunc(T)(T value) @trusted pure nothrow
366 {
367     static if (is(T==float) || is(T==double))
368         return cast(int) value;
369     else
370         static assert(0, "unsupported FP type");
371 }
372 ///
373 @safe pure nothrow unittest
374 {
375     assert(dtrunc(0.2f) == 0);
376     assert(dtrunc(0.8f) == 0);
377     assert(dtrunc(-0.2f) == 0);
378     assert(dtrunc(-8.8f) == -8);
379 }
380 
381 /**
382  * Computes the hypothenus of two FP numbers.
383  *
384  * This function requires the SSE2 instruction set and it can't be inlined.
385  *
386  * Params:
387  *      x = Either a float or a double.
388  *      y = Either a float or a double.
389  *
390  * Returns:
391  *      A floating point value of type T.
392  */
393 extern(C) T hypot(T)(T x, T y) pure @trusted nothrow
394 {
395     version(X86_64)
396     {
397         static if (is(T==float)) asm pure nothrow
398         {
399             naked;
400             mulps   XMM0, XMM0;
401             mulps   XMM1, XMM1;
402             addps   XMM0, XMM1;
403             sqrtps  XMM0, XMM0;
404             ret;
405         }
406         else static if (is(T==double)) asm pure nothrow
407         {
408             naked;
409             mulpd   XMM0, XMM0;
410             mulpd   XMM1, XMM1;
411             addpd   XMM0, XMM1;
412             sqrtpd  XMM0, XMM0;
413             ret;
414         }
415         else static assert(0, "unsupported FP type");
416     }
417     else version(X86)
418     {
419         static if (is(T==float)) asm pure nothrow
420         {
421             naked;
422             push    EBP;
423             mov     EBP, ESP;
424             sub     ESP, 4;
425             movss   XMM0, x;
426             movss   XMM1, y;
427             mulps   XMM0, XMM0;
428             mulps   XMM1, XMM1;
429             addps   XMM0, XMM1;
430             sqrtps  XMM0, XMM0;
431             movss   dword ptr [ESP-4], XMM0;
432             fld     dword ptr [ESP-4];
433             add     ESP, 4;
434             mov     ESP, EBP;
435             pop     EBP;
436             ret;
437         }
438         else static if (is(T==double)) asm pure nothrow
439         {
440             naked;
441             push    EBP;
442             mov     EBP, ESP;
443             sub     ESP, 8;
444             movsd   XMM0, x;
445             movsd   XMM1, y;
446             mulpd   XMM0, XMM0;
447             mulpd   XMM1, XMM1;
448             addpd   XMM0, XMM1;
449             sqrtpd  XMM0, XMM0;
450             movsd   qword ptr [ESP-8], XMM0;
451             fld     qword ptr [ESP-8];
452             add     ESP, 8;
453             mov     ESP, EBP;
454             pop     EBP;
455             ret;
456         }
457         else static assert(0, "unsupported FP type");
458     }
459     else static assert(0, "unsupported architecture");
460 }
461 ///
462 pure @safe nothrow unittest
463 {
464     assert(hypot(3.0,4.0) == 5.0);
465     assert(hypot(3.0f,4.0f) == 5.0f);
466 }
467 
468 /**
469  * Returns the magnitude of a complex number.
470  *
471  * Convenience function that calls hypot() either with a clfloat, a cdouble or
472  * a std.complex.Complex.
473  */
474 auto magn(T)(T t)
475 if (is(T==cfloat) || is(T==cdouble) || is(T==Complex!float) || is(T==Complex!double))
476 {
477     return hypot(t.re, t.im);
478 }
479 ///
480 nothrow pure @safe unittest
481 {
482     assert(magn(cdouble(3.0+4.0i)) == 5.0);
483     assert(magn(cfloat(3.0f+4.0fi)) == 5.0f);
484     assert(magn(complex(3.0f,4.0f)) == 5.0f);
485     assert(magn(complex(3.0,4.0)) == 5.0);
486 }
487 
488 /**
489  * Wraps a numeric value between 0 and a max value.
490  *
491  * Params:
492  *      bound = a string that indicates if the max is excluded ($(D ")")), the default
493  *          or if the max is included ($(D "]")).
494  *      value = the value to wrap.
495  *      max the maximal value.
496  */
497 T wrap(string bound = ")", T)(T value, T max)
498 if (isNumeric!T && (bound == ")" || bound == "]"))
499 {
500     static if (bound == ")")
501     {
502         if (value > max)
503             return value - max;
504     }
505     static if (bound == "]")
506     {
507         if (value >= max)
508             return value - max;
509     }
510     if (value < 0)
511         return max + value;
512     else
513         return value;
514 }
515 ///
516 nothrow @nogc @safe pure unittest
517 {
518     import std.math: approxEqual;
519     import std.math: modf;
520     // wrap from max
521     assert(wrap(1,1) == 1);
522     assert(wrap(3,2) == 1);
523     assert(wrap(-1,3) == 2);
524     assert(wrap(1.5,1).approxEqual(0.5));
525     assert(wrap(1.01,1).approxEqual(0.01));
526     assert(wrap(-0.5,2).approxEqual(1.5));
527     // wrap past max
528     assert(wrap!"]"(1,1) == 0);
529     assert(wrap!"]"(3,2) == 1);
530     assert(wrap!"]"(3,3) == 0);
531     assert(wrap!"]"(1.0,1.0) == 0.0);
532     assert(wrap!"]"(-0.5,1.0) == 0.5);
533     // two phases in sync
534     double incr = 0.0125;
535     double sync = 0.25;
536     double phase1 = 0, phase2 = 0;
537     foreach(i; 0 .. 100000)
538     {
539         phase1 = wrap(phase1 + incr, 1.0);
540         phase2 = wrap(phase1 - sync, 1.0);
541         assert(phase1 < 1.0 && phase2 < 1.0);
542     }
543 }
544 
545 /**
546  * Allows to represent fractions of PI without using the usual suffixes such as
547  * "two", "half", etc.
548  */
549 template Pi(int a, int b = 1)
550 {
551     import std.math: PI;
552     enum Pi = PI * a / b;
553 }
554 ///
555 unittest
556 {
557     static assert(Pi!(2,1) == Pi!(4,2));
558     static assert(Pi!(1,2) == Pi!(1,8) * 4);
559     static assert(Pi!(4,2) == Pi!2);
560 }
561 
562 /// Log(0.5) as a double;
563 enum logHalf = -0.69314718055994531;
564 
565 /**
566  * Retrieves the exponent part on the result of pow().
567  *
568  * Params:
569  *      x = The first pow() argument.
570  *      y = The pow() result.
571  * Returns:
572  *      The exponent N that verifies y = pow(x, N).
573  */
574 double logN(double x, double y)
575 {
576     import std.math;
577     return log(y) / log(x);
578 }
579 ///
580 unittest
581 {
582     import std.math;
583     double[] exps = [3, 1.4, 2.2, 0.0001, 0.9999];
584     foreach(e; exps)
585     {
586         double y = pow(0.5, e);
587         double n = logN(0.5, y);
588         assert(n.approxEqual(e));
589     }
590 }
591 
592 /**
593  * A simple and fast parametric easing function.
594  *
595  * Its name comes from the fact that the plot of f(x) in [0..1] and with a
596  * control of 3 forms a regular parabolic curve.
597  */
598 struct VariableParabol
599 {
600     @disable this(this);
601 
602     /**
603      * Applies the standard transformation.
604      *
605      * Params:
606      *      x = The X coordinate, between 0.0 and 1.0
607      *      c = The control of the speed, between 0 (similar to InExpo)
608      *      and 3.0 (similar to OutExpo)
609      *
610      * Returns:
611      *      The Y coordinate, a value between 0.0 and 1.0.
612      */
613     static double fx(int NC = 1)(double x, double c) pure @safe @nogc
614     {
615         import iz.logicver;
616         static if (verX86_64 & verDigitalMars)
617             enum PureD = false;
618         else // X86 with FPU, better codegen with LDC, etc
619             enum PureD = true;
620         static if (PureD)
621         {
622             assert(0 <= x && x <= 1.0);
623             assert(0 <= c && c <= 3.0);
624             return x*x*x - x*x*c + x*c;
625         }
626         else asm pure nothrow @nogc @safe
627         {
628             naked;
629             movapd  XMM2, XMM1; // saves x
630             mulsd   XMM2, XMM2; // x²
631             movapd  XMM3, XMM2; // saves x²
632             mulsd   XMM3, XMM1; // x³
633             mulsd   XMM2, XMM0; // cx²
634             subsd   XMM3, XMM2; // x³ - cx²
635             mulsd   XMM0, XMM1; // cx
636             addsd   XMM0, XMM3; // x³ - cx² + cx
637             ret;
638         }
639     }
640 
641     /**
642      * Retrieves the control point coefficient.
643      *
644      * Params:
645      *      y = The value as obtained by fx(0.5, control).
646      *
647      * Returns:
648      *      a value between 0.0 and 3.0.
649      */
650     static double control(int N = 0)(double y) pure @safe @nogc
651     in
652     {
653         assert(0 <= y && y <= 1.0);
654     }
655     out (c)
656     {
657         assert(0 <= c && c <= 3.0);
658     }
659     body
660     {
661         return controlClip(4 * (y - 0.125));
662     }
663 
664     /**
665      * Computes the Y coordinate of the control point.
666      *
667      * Returns:
668      *      The same as fx(0.5, c) but faster.
669      */
670     static double controlFx(int N = 0)(double c) pure @safe @nogc
671     in
672     {
673         assert(0 <= c && c <= 3.0);
674     }
675     out (y)
676     {
677         assert(0 <= y && y <= 1.0);
678     }
679     body
680     {
681         return 0.125 + 0.25 * c;
682     }
683 
684     /**
685      * Clips the control point coefficient.
686      *
687      * Params:
688      *      c = The control point coefficient.
689      * Returns:
690      *      The control point, validated for fx().
691      */
692     static double controlClip(double c) pure @safe @nogc
693     {
694         import std.algorithm.comparison: clamp;
695         return c.clamp(0.0, 3.0);
696     }
697 
698     /**
699      * Indicates that 1 control point is used.
700      */
701     enum numControls = 1;
702 
703     /// Uniform Easing API.
704     static shared VariableParabol ease;
705 }
706 ///
707 @nogc pure @safe unittest
708 {
709     double x = 0;
710     double x0 = 120;
711     double x1 = 480;
712     double h = 200;
713     const int numPoints = 1000;
714     const double increment = 1.0 / numPoints;
715 
716     // draw a curve...
717     foreach(i; 0..numPoints)
718     {
719         double y = variableParabol.fx(x, 0.1);
720         // lineTo(x * (480 - 120) + 120, y * h);
721         x += increment;
722     }
723     import std.math: approxEqual;
724     assert(x.approxEqual(1.0));
725 
726     enum inc = 0.001;
727     double control = 0;
728     foreach(i; 0..1000)
729     {
730         auto cy = variableParabol.fx(0.5, control);
731         // Operation to get the control value as function of the cursor
732         assert(variableParabol.control(cy) == control);
733         control += 3 / 1000;
734     }
735 
736     assert(variableParabol.fx(0.0,3.0) == 0.0);
737     assert(variableParabol.fx(1.0,3.0) == 1.0);
738     assert(variableParabol.fx(0.0,2.0) == 0.0);
739     assert(variableParabol.fx(1.0,2.0) == 1.0);
740 }
741 
742 ///
743 alias variableParabol = VariableParabol;
744 
745 
746 /**
747  * Parametric easing function using pow().
748  */
749 struct VariablePow
750 {
751     @disable this(this);
752 
753     /**
754      * Applies the standard transformation.
755      *
756      * Params:
757      *      x = The X coordinate, between 0.0 and 1.0
758      *      c = The control of the speed, between 0 (InExpo)
759      *      and 9.0 (OutExpo)
760      *
761      * Returns:
762      *      The Y coordinate, a value between 0.0 and 1.0.
763      */
764     static double fx(int NC = 1)(double x, double c) pure @safe @nogc
765     in
766     {
767         assert(0 <= x && x <= 1.0);
768         assert(0.05 <= c && c <= 9.0);
769     }
770     body
771     {
772         import std.math: pow;
773         return pow(x, c);
774     }
775 
776     /**
777      * Retrieves the control point coefficient.
778      *
779      * Params:
780      *      y = The value as obtained by fx(0.5, control).
781      *
782      * Returns:
783      *      a value between 0.0 and 9.0.
784      */
785     static double control(int N = 0)(double y) pure @safe @nogc
786     in
787     {
788         assert(0 <= y && y <= 1.0);
789     }
790     out (c)
791     {
792         assert(0.05 <= c && c <= 9.0);
793     }
794     body
795     {
796         import std.math: log;
797         return controlClip(log(y) / logHalf);
798     }
799 
800     /**
801      * Computes the Y coordinate of the control point.
802      *
803      * Returns:
804      *      The same as fx(0.5, c).
805      */
806     static double controlFx(int N = 0)(double c) pure @safe @nogc
807     in
808     {
809         assert(0.05 <= c && c <= 9.0);
810     }
811     out (y)
812     {
813         assert(0 <= y && y <= 1.0);
814     }
815     body
816     {
817         return fx(0.5, c);
818     }
819 
820     /**
821      * Clips the control point coefficient.
822      *
823      * Params:
824      *      c = The control point coefficient.
825      * Returns:
826      *      The control point, validated for fx().
827      */
828     static double controlClip(double c) pure @safe @nogc
829     {
830         import std.algorithm.comparison: clamp;
831         return c.clamp(0.05, 9.0);
832     }
833 
834     /**
835      * Indicates that 1 control point is used.
836      */
837     enum numControls = 1;
838 
839     /// Uniform Easing API.
840     static shared VariablePow ease;
841 }
842 
843 ///
844 alias variablePow = VariablePow;
845 
846 /**
847  * An parametric easing function based on the Supper Ellipse.
848  */
849 struct VariableEllipse
850 {
851     @disable this(this);
852 
853     /**
854      * Applies the standard transformation.
855      *
856      * Params:
857      *      x = The X coordinate, between 0.0 and 1.0
858      *      c = The control of the speed, between 0 (similar to InExpo)
859      *      and 3.0 (similar to OutExpo)
860      *
861      * Returns:
862      *      The Y coordinate, a value between 0.0 and 1.0.
863      */
864     static double fx(int NC = 1)(double x, double c) pure @safe @nogc
865     in
866     {
867         assert(0 <= x && x <= 1.0);
868         assert(0.1 <= c && x <= 8.0);
869     }
870     body
871     {
872         import std.math: pow;
873         return 1.0 - pow(1.0 - pow(x, 2.0/c), c * 0.5);
874     }
875 
876     /**
877      * Retrieves the control point coefficient.
878      *
879      * Params:
880      *      y = The value as obtained by fx(0.5, control).
881      *
882      * Returns:
883      *      A value between 0.0 and 8.0.
884      */
885     static double control(int N = 0)(double y) pure @safe @nogc
886     in
887     {
888         assert(0 <= y && y <= 1.0);
889     }
890     out (c)
891     {
892         assert(0.1 <= c && c <= 8.0);
893     }
894     body
895     {
896         return controlClip(variablePow.control(1 - y));
897     }
898 
899     /**
900      * Computes the Y coordinate of the control point.
901      *
902      * Returns:
903      *      The same as fx(0.5, c).
904      */
905     static double controlFx(int N = 0)(double c) pure @safe @nogc
906     in
907     {
908         assert(0.1 <= c && c <= 8.0);
909     }
910     out (y)
911     {
912         assert(0 <= y && y <= 1.0);
913     }
914     body
915     {
916         return fx(0.5, c);
917     }
918 
919     /**
920      * Clips the control point coefficient.
921      *
922      * Params:
923      *      c = The control point coefficient.
924      * Returns:
925      *      The control point, validated for fx().
926      */
927     static double controlClip(double c) pure @safe @nogc
928     {
929         import std.algorithm.comparison: clamp;
930         return c.clamp(0.1, 8.0);
931     }
932 
933     /**
934      * Indicates that 1 control point is used.
935      */
936     enum numControls = 1;
937 
938     /// Uniform Easing API.
939     static shared VariableEllipse ease;
940 }
941 
942 ///
943 alias variableEllipse = VariableEllipse;
944