1 /*
2     SuperCollider real time audio synthesis system
3     Copyright (c) 2002 James McCartney. All rights reserved.
4     http://www.audiosynth.com
5 
6     This program is free software; you can redistribute it and/or modify
7     it under the terms of the GNU General Public License as published by
8     the Free Software Foundation; either version 2 of the License, or
9     (at your option) any later version.
10 
11     This program is distributed in the hope that it will be useful,
12     but WITHOUT ANY WARRANTY; without even the implied warranty of
13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14     GNU General Public License for more details.
15 
16     You should have received a copy of the GNU General Public License
17     along with this program; if not, write to the Free Software
18     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
19 */
20 
21 #include "PyrKernel.h"
22 #include "PyrPrimitive.h"
23 #include "PyrMathPrim.h"
24 #include "MiscInlineMath.h"
25 #include "SC_InlineUnaryOp.h"
26 #include "SC_InlineBinaryOp.h"
27 #include "PyrSignal.h"
28 #include "PyrParseNode.h"
29 #include "PyrMessage.h"
30 #include "clz.h"
31 #include <stdlib.h>
32 #include <string.h>
33 #include <math.h>
34 #include "SC_Endian.h"
35 #include "SCBase.h"
36 
37 #define BOOST_MATH_DISABLE_FLOAT128 1
38 #include "boost/math/special_functions.hpp"
39 #include "boost/math/special_functions/spherical_harmonic.hpp"
40 #include "boost/math/special_functions/legendre.hpp"
41 #include "PyrSymbol.h"
42 #include <vector>
43 
44 const int INT_MAX_BY_PyrSlot = INT_MAX / sizeof(PyrSlot);
45 
IsSignal(PyrSlot * slot)46 inline bool IsSignal(PyrSlot* slot) { return (IsObj(slot) && slotRawObject(slot)->classptr == class_signal); }
NotSignal(PyrSlot * slot)47 inline bool NotSignal(PyrSlot* slot) { return (NotObj(slot) || slotRawObject(slot)->classptr != class_signal); }
48 
49 
50 /* functors for dispatching template code */
51 
52 struct addNum {
runaddNum53     static inline double run(double lhs, double rhs) { return lhs + rhs; }
runaddNum54     static inline int run(int lhs, int rhs) { return lhs + rhs; }
signal_xfaddNum55     static inline PyrObject* signal_xf(VMGlobals* g, PyrObject* ina, float inb) { return signal_add_xf(g, ina, inb); }
signal_fxaddNum56     static inline PyrObject* signal_fx(VMGlobals* g, float ina, PyrObject* inb) { return signal_xf(g, inb, ina); }
signal_xxaddNum57     static inline PyrObject* signal_xx(VMGlobals* g, PyrObject* ina, PyrObject* inb) {
58         return signal_add_xx(g, ina, inb);
59     }
60 };
61 
62 struct mulNum {
runmulNum63     static inline double run(double lhs, double rhs) { return lhs * rhs; }
runmulNum64     static inline int run(int lhs, int rhs) { return lhs * rhs; }
signal_xfmulNum65     static inline PyrObject* signal_xf(VMGlobals* g, PyrObject* ina, float inb) { return signal_mul_xf(g, ina, inb); }
signal_fxmulNum66     static inline PyrObject* signal_fx(VMGlobals* g, float ina, PyrObject* inb) { return signal_xf(g, inb, ina); }
signal_xxmulNum67     static inline PyrObject* signal_xx(VMGlobals* g, PyrObject* ina, PyrObject* inb) {
68         return signal_mul_xx(g, ina, inb);
69     }
70 };
71 
72 struct subNum {
runsubNum73     static inline double run(double lhs, double rhs) { return lhs - rhs; }
runsubNum74     static inline int run(int lhs, int rhs) { return lhs - rhs; }
signal_xfsubNum75     static inline PyrObject* signal_xf(VMGlobals* g, PyrObject* ina, float inb) { return signal_sub_xf(g, ina, inb); }
signal_fxsubNum76     static inline PyrObject* signal_fx(VMGlobals* g, float ina, PyrObject* inb) { return signal_sub_fx(g, ina, inb); }
signal_xxsubNum77     static inline PyrObject* signal_xx(VMGlobals* g, PyrObject* ina, PyrObject* inb) {
78         return signal_sub_xx(g, ina, inb);
79     }
80 };
81 
prOpNum(VMGlobals * g,int numArgsPushed)82 template <typename Functor> inline int prOpNum(VMGlobals* g, int numArgsPushed) {
83     PyrSlot *a, *b;
84     PyrSymbol* msg;
85 
86     a = g->sp - 1;
87     b = g->sp;
88 
89     switch (GetTag(a)) {
90     case tagInt:
91         switch (GetTag(b)) {
92         case tagInt:
93             SetRaw(a, Functor::run(slotRawInt(a), slotRawInt(b)));
94             break;
95         case tagChar:
96         case tagPtr:
97         case tagNil:
98         case tagFalse:
99         case tagTrue:
100             goto send_normal_2;
101         case tagSym:
102             SetSymbol(a, slotRawSymbol(b));
103             break;
104         case tagObj:
105             if (isKindOf(slotRawObject(b), class_signal))
106                 SetObject(a, Functor::signal_fx(g, slotRawInt(a), slotRawObject(b)));
107             else
108                 goto send_normal_2;
109             break;
110         default:
111             SetFloat(a, slotRawInt(a) + slotRawFloat(b));
112             break;
113         }
114         break;
115     case tagChar:
116     case tagPtr:
117     case tagNil:
118     case tagFalse:
119     case tagTrue:
120         goto send_normal_2;
121     case tagSym:
122         // leave self in 'a'
123         break;
124     case tagObj:
125         if (isKindOf(slotRawObject(a), class_signal)) {
126             switch (GetTag(b)) {
127             case tagInt:
128                 SetRaw(a, Functor::signal_xf(g, slotRawObject(a), slotRawInt(b)));
129                 break;
130             case tagChar:
131             case tagPtr:
132             case tagNil:
133             case tagFalse:
134             case tagTrue:
135                 goto send_normal_2;
136             case tagSym:
137                 SetSymbol(a, slotRawSymbol(b));
138                 break;
139             case tagObj:
140                 if (isKindOf(slotRawObject(b), class_signal)) {
141                     SetRaw(a, Functor::signal_xx(g, slotRawObject(a), slotRawObject(b)));
142                 } else
143                     goto send_normal_2;
144                 break;
145             default: // double
146                 SetRaw(a, Functor::signal_xf(g, slotRawObject(a), slotRawFloat(b)));
147                 break;
148             }
149         } else
150             goto send_normal_2;
151         break;
152     default: // double
153         switch (GetTag(b)) {
154         case tagInt:
155             SetRaw(a, Functor::run(slotRawFloat(a), (double)slotRawInt(b)));
156             break;
157         case tagChar:
158         case tagPtr:
159         case tagNil:
160         case tagFalse:
161         case tagTrue:
162             goto send_normal_2;
163         case tagSym:
164             SetSymbol(a, slotRawSymbol(b));
165             break;
166         case tagObj:
167             if (isKindOf(slotRawObject(b), class_signal))
168                 SetObject(a, Functor::signal_fx(g, slotRawFloat(a), slotRawObject(b)));
169             else
170                 goto send_normal_2;
171             break;
172         default: // double
173             SetRaw(a, Functor::run(slotRawFloat(a), slotRawFloat(b)));
174             break;
175         }
176         break;
177     }
178     g->sp--; // drop
179     g->numpop = 0;
180 #if TAILCALLOPTIMIZE
181     g->tailCall = 0;
182 #endif
183     return errNone;
184 
185 send_normal_2:
186     if (numArgsPushed != -1) // special case flag meaning it is a primitive
187         return errFailed; // arguments remain on the stack
188 
189     msg = gSpecialBinarySelectors[g->primitiveIndex];
190     sendMessage(g, msg, 2);
191     return errNone;
192 }
193 
prOpInt(VMGlobals * g,int numArgsPushed)194 template <typename Functor> inline int prOpInt(VMGlobals* g, int numArgsPushed) {
195     PyrSlot *a, *b;
196     PyrSymbol* msg;
197 
198     a = g->sp - 1;
199     b = g->sp;
200 
201     switch (GetTag(b)) {
202     case tagInt:
203         SetRaw(a, Functor::run(slotRawInt(a), slotRawInt(b)));
204         break;
205     case tagChar:
206     case tagPtr:
207     case tagNil:
208     case tagFalse:
209     case tagTrue:
210         goto send_normal_2;
211     case tagSym:
212         SetSymbol(a, slotRawSymbol(b));
213         break;
214     case tagObj:
215         if (isKindOf(slotRawObject(b), class_signal))
216             SetObject(a, Functor::signal_fx(g, slotRawInt(a), slotRawObject(b)));
217         else
218             goto send_normal_2;
219         break;
220     default:
221         SetFloat(a, Functor::run((double)slotRawInt(a), slotRawFloat(b)));
222         break;
223     }
224     g->sp--; // drop
225     g->numpop = 0;
226 #if TAILCALLOPTIMIZE
227     g->tailCall = 0;
228 #endif
229     return errNone;
230 
231 send_normal_2:
232     if (numArgsPushed != -1) // special case flag meaning it is a primitive
233         return errFailed; // arguments remain on the stack
234 
235     msg = gSpecialBinarySelectors[g->primitiveIndex];
236     sendMessage(g, msg, 2);
237     return errNone;
238 }
239 
prOpFloat(VMGlobals * g,int numArgsPushed)240 template <typename Functor> inline int prOpFloat(VMGlobals* g, int numArgsPushed) {
241     PyrSlot *a, *b;
242     PyrSymbol* msg;
243 
244     a = g->sp - 1;
245     b = g->sp;
246 
247     switch (GetTag(b)) {
248     case tagInt:
249         SetRaw(a, Functor::run(slotRawFloat(a), (double)slotRawInt(b)));
250         break;
251     case tagChar:
252     case tagPtr:
253     case tagNil:
254     case tagFalse:
255     case tagTrue:
256         goto send_normal_2;
257     case tagSym:
258         SetSymbol(a, slotRawSymbol(b));
259         break;
260     case tagObj:
261         if (isKindOf(slotRawObject(b), class_signal))
262             SetObject(a, Functor::signal_fx(g, slotRawFloat(a), slotRawObject(b)));
263         else
264             goto send_normal_2;
265         break;
266     default:
267         SetRaw(a, Functor::run(slotRawFloat(a), slotRawFloat(b)));
268         break;
269     }
270     g->sp--; // drop
271     g->numpop = 0;
272 #if TAILCALLOPTIMIZE
273     g->tailCall = 0;
274 #endif
275     return errNone;
276 
277 send_normal_2:
278     if (numArgsPushed != -1) // special case flag meaning it is a primitive
279         return errFailed; // arguments remain on the stack
280 
281     msg = gSpecialBinarySelectors[g->primitiveIndex];
282     sendMessage(g, msg, 2);
283     return errNone;
284 }
285 
prAddNum(VMGlobals * g,int numArgsPushed)286 int prAddNum(VMGlobals* g, int numArgsPushed) { return prOpNum<addNum>(g, numArgsPushed); }
287 
prSubNum(VMGlobals * g,int numArgsPushed)288 int prSubNum(VMGlobals* g, int numArgsPushed) { return prOpNum<subNum>(g, numArgsPushed); }
289 
prMulNum(VMGlobals * g,int numArgsPushed)290 int prMulNum(VMGlobals* g, int numArgsPushed) { return prOpNum<mulNum>(g, numArgsPushed); }
291 
292 
prAddFloat(VMGlobals * g,int numArgsPushed)293 int prAddFloat(VMGlobals* g, int numArgsPushed) { return prOpFloat<addNum>(g, numArgsPushed); }
294 
prSubFloat(VMGlobals * g,int numArgsPushed)295 int prSubFloat(VMGlobals* g, int numArgsPushed) { return prOpFloat<subNum>(g, numArgsPushed); }
296 
prMulFloat(VMGlobals * g,int numArgsPushed)297 int prMulFloat(VMGlobals* g, int numArgsPushed) { return prOpFloat<mulNum>(g, numArgsPushed); }
298 
prAddInt(VMGlobals * g,int numArgsPushed)299 int prAddInt(VMGlobals* g, int numArgsPushed) { return prOpInt<addNum>(g, numArgsPushed); }
300 
prSubInt(VMGlobals * g,int numArgsPushed)301 int prSubInt(VMGlobals* g, int numArgsPushed) { return prOpInt<subNum>(g, numArgsPushed); }
302 
prMulInt(VMGlobals * g,int numArgsPushed)303 int prMulInt(VMGlobals* g, int numArgsPushed) { return prOpInt<mulNum>(g, numArgsPushed); }
304 
305 
306 int prNthPrime(VMGlobals* g, int numArgsPushed);
prNthPrime(VMGlobals * g,int numArgsPushed)307 int prNthPrime(VMGlobals* g, int numArgsPushed) {
308     PyrSlot* a;
309     int n, p;
310 
311     a = g->sp;
312     n = slotRawInt(a);
313     p = nthPrime(n);
314     if (p == 0) {
315         SetNil(a);
316     } else {
317         SetInt(a, p);
318     }
319     return errNone;
320 }
321 
322 int prPrevPrime(VMGlobals* g, int numArgsPushed);
prPrevPrime(VMGlobals * g,int numArgsPushed)323 int prPrevPrime(VMGlobals* g, int numArgsPushed) {
324     PyrSlot* a;
325     int n, p, i;
326 
327     a = g->sp;
328     n = slotRawInt(a);
329     i = prevPrime(n);
330     p = nthPrime(i);
331     if (p == 0) {
332         SetNil(a);
333     } else {
334         SetInt(a, p);
335     }
336     return errNone;
337 }
338 
339 int prNextPrime(VMGlobals* g, int numArgsPushed);
prNextPrime(VMGlobals * g,int numArgsPushed)340 int prNextPrime(VMGlobals* g, int numArgsPushed) {
341     PyrSlot* a;
342     int n, p, i;
343 
344     a = g->sp;
345     n = slotRawInt(a);
346     i = nextPrime(n);
347     p = nthPrime(i);
348     if (p == 0) {
349         SetNil(a);
350     } else {
351         SetInt(a, p);
352     }
353     return errNone;
354 }
355 
356 
357 int prIsPrime(VMGlobals* g, int numArgsPushed);
prIsPrime(VMGlobals * g,int numArgsPushed)358 int prIsPrime(VMGlobals* g, int numArgsPushed) {
359     PyrSlot* a;
360     int n, p, sqrtn, i;
361 
362     a = g->sp;
363     n = slotRawInt(a);
364     SetNil(a);
365     if (n <= 2) {
366         if (n == 2) {
367             SetTrue(a);
368         } else {
369             SetFalse(a);
370         }
371     } else if (n <= nthPrime(NUMPRIMES - 1)) {
372         // do a search of the primes table
373         i = findPrime(n);
374         if (i >= 0) {
375             SetTrue(a);
376         } else {
377             SetFalse(a);
378         }
379     } else {
380 #ifdef _WIN32
381         sqrtn = (int)sqrt(static_cast<double>(n));
382 #else
383         sqrtn = (int)sqrt(n);
384 #endif
385         for (i = 0; i < NUMPRIMES; ++i) {
386             p = nthPrime(i);
387             if (n % p == 0) {
388                 SetFalse(a);
389                 break;
390             }
391             if (p >= sqrtn) {
392                 SetTrue(a);
393                 break;
394             }
395         }
396     }
397     return errNone;
398 }
399 
400 int prIndexOfPrime(VMGlobals* g, int numArgsPushed);
prIndexOfPrime(VMGlobals * g,int numArgsPushed)401 int prIndexOfPrime(VMGlobals* g, int numArgsPushed) {
402     PyrSlot* a;
403     int n, p;
404 
405     a = g->sp;
406     n = slotRawInt(a);
407     if (n <= 2) {
408         if (n == 2) {
409             SetInt(a, 0);
410         } else {
411             SetNil(a);
412         }
413     } else if (n <= nthPrime(NUMPRIMES - 1)) {
414         p = findPrime(n);
415         if (p < 0) {
416             SetNil(a);
417         } else {
418             SetInt(a, p);
419         }
420     } else {
421         SetNil(a);
422     }
423     return errNone;
424 }
425 
426 
427 int prAs32Bits(VMGlobals* g, int numArgsPushed);
prAs32Bits(VMGlobals * g,int numArgsPushed)428 int prAs32Bits(VMGlobals* g, int numArgsPushed) {
429     PyrSlot* a = g->sp;
430     // return an integer that is a bit pattern for the 32 bit float representation
431     union {
432         float f;
433         int32 i;
434     } u;
435     u.f = slotRawFloat(a);
436     SetInt(a, u.i);
437     return errNone;
438 }
439 
440 int prHigh32Bits(VMGlobals* g, int numArgsPushed);
prHigh32Bits(VMGlobals * g,int numArgsPushed)441 int prHigh32Bits(VMGlobals* g, int numArgsPushed) {
442     PyrSlot* a = g->sp;
443 
444 #if BYTE_ORDER == BIG_ENDIAN
445     union {
446         struct {
447             uint32 hi, lo;
448         } i;
449         double f;
450     } du;
451 #else
452     union {
453         struct {
454             uint32 lo, hi;
455         } i;
456         double f;
457     } du;
458 #endif
459 
460     du.f = slotRawFloat(a);
461     SetInt(a, du.i.hi);
462     return errNone;
463 }
464 
465 int prLow32Bits(VMGlobals* g, int numArgsPushed);
prLow32Bits(VMGlobals * g,int numArgsPushed)466 int prLow32Bits(VMGlobals* g, int numArgsPushed) {
467     PyrSlot* a = g->sp;
468 
469 #if BYTE_ORDER == BIG_ENDIAN
470     union {
471         struct {
472             uint32 hi, lo;
473         } i;
474         double f;
475     } du;
476 #else
477     union {
478         struct {
479             uint32 lo, hi;
480         } i;
481         double f;
482     } du;
483 #endif
484 
485     du.f = slotRawFloat(a);
486     SetInt(a, du.i.lo);
487     return errNone;
488 }
489 
490 int prFrom32Bits(VMGlobals* g, int numArgsPushed);
prFrom32Bits(VMGlobals * g,int numArgsPushed)491 int prFrom32Bits(VMGlobals* g, int numArgsPushed) {
492     PyrSlot* a = g->sp - 1;
493     PyrSlot* b = g->sp;
494     int err, word;
495 
496     err = slotIntVal(b, &word);
497     if (err)
498         return err;
499 
500     union {
501         float f;
502         int32 i;
503     } u;
504     u.i = word;
505     SetFloat(a, u.f);
506 
507     return errNone;
508 }
509 
510 int prFrom64Bits(VMGlobals* g, int numArgsPushed);
prFrom64Bits(VMGlobals * g,int numArgsPushed)511 int prFrom64Bits(VMGlobals* g, int numArgsPushed) {
512     PyrSlot* a = g->sp - 2;
513     PyrSlot* b = g->sp - 1;
514     PyrSlot* c = g->sp;
515     int err, hi, lo;
516 
517     err = slotIntVal(b, &hi);
518     if (err)
519         return err;
520 
521     err = slotIntVal(c, &lo);
522     if (err)
523         return err;
524 
525 #if BYTE_ORDER == BIG_ENDIAN
526     union {
527         struct {
528             uint32 hi, lo;
529         } i;
530         double f;
531     } du;
532 #else
533     union {
534         struct {
535             uint32 lo, hi;
536         } i;
537         double f;
538     } du;
539 #endif
540     du.i.hi = hi;
541     du.i.lo = lo;
542     SetFloat(a, du.f);
543 
544     return errNone;
545 }
546 
mathClipInt(struct VMGlobals * g,int numArgsPushed)547 int mathClipInt(struct VMGlobals* g, int numArgsPushed) {
548     PyrSlot *a, *b, *c;
549     double lo, hi;
550     int err;
551 
552     a = g->sp - 2;
553     b = g->sp - 1;
554     c = g->sp;
555 
556     if (IsSym(b)) {
557         *a = *b;
558     } else if (IsSym(c)) {
559         *a = *c;
560     } else if (IsInt(b) && IsInt(c)) {
561         SetRaw(a, sc_clip(slotRawInt(a), slotRawInt(b), slotRawInt(c)));
562     } else {
563         err = slotDoubleVal(b, &lo);
564         if (err)
565             return err;
566         err = slotDoubleVal(c, &hi);
567         if (err)
568             return err;
569         SetFloat(a, sc_clip((double)slotRawInt(a), lo, hi));
570     }
571     return errNone;
572 }
573 
mathClipFloat(struct VMGlobals * g,int numArgsPushed)574 int mathClipFloat(struct VMGlobals* g, int numArgsPushed) {
575     PyrSlot *a, *b, *c;
576     double lo, hi;
577     int err;
578 
579     a = g->sp - 2;
580     b = g->sp - 1;
581     c = g->sp;
582 
583     if (IsSym(b)) {
584         *a = *b;
585     } else if (IsSym(c)) {
586         *a = *c;
587     } else {
588         err = slotDoubleVal(b, &lo);
589         if (err)
590             return err;
591         err = slotDoubleVal(c, &hi);
592         if (err)
593             return err;
594         SetRaw(a, sc_clip(slotRawFloat(a), lo, hi));
595     }
596     return errNone;
597 }
598 
mathClipSignal(struct VMGlobals * g,int numArgsPushed)599 int mathClipSignal(struct VMGlobals* g, int numArgsPushed) {
600     PyrSlot *a, *b, *c;
601     float lo, hi;
602     int err;
603     PyrObject* sig;
604 
605     a = g->sp - 2;
606     b = g->sp - 1;
607     c = g->sp;
608 
609     if (IsSym(b)) {
610         *a = *b;
611     } else if (IsSym(c)) {
612         *a = *c;
613     } else if (IsSignal(b) && IsSignal(c)) {
614         sig = signal_clip_x(g, slotRawObject(a), slotRawObject(b), slotRawObject(c));
615         SetObject(a, sig);
616     } else {
617         err = slotFloatVal(b, &lo);
618         if (err)
619             return err;
620         err = slotFloatVal(c, &hi);
621         if (err)
622             return err;
623         sig = signal_clip_f(g, slotRawObject(a), lo, hi);
624         SetObject(a, sig);
625     }
626     return errNone;
627 }
628 
mathWrapInt(struct VMGlobals * g,int numArgsPushed)629 int mathWrapInt(struct VMGlobals* g, int numArgsPushed) {
630     PyrSlot *a, *b, *c;
631     int err;
632 
633     a = g->sp - 2;
634     b = g->sp - 1;
635     c = g->sp;
636 
637     if (IsSym(b)) {
638         *a = *b;
639     } else if (IsSym(c)) {
640         *a = *c;
641     } else if (IsInt(b) && IsInt(c)) {
642         SetRaw(a,
643                sc_mod((int)(slotRawInt(a) - slotRawInt(b)), (int)(slotRawInt(c) - slotRawInt(b) + 1)) + slotRawInt(b));
644     } else {
645         double x, lo, hi;
646         x = slotRawInt(a);
647         err = slotDoubleVal(b, &lo);
648         if (err)
649             return err;
650         err = slotDoubleVal(c, &hi);
651         if (err)
652             return err;
653         SetFloat(a, sc_mod(x - lo, hi - lo) + lo);
654     }
655     return errNone;
656 }
657 
mathWrapFloat(struct VMGlobals * g,int numArgsPushed)658 int mathWrapFloat(struct VMGlobals* g, int numArgsPushed) {
659     PyrSlot *a, *b, *c;
660     double lo, hi;
661     int err;
662 
663     a = g->sp - 2;
664     b = g->sp - 1;
665     c = g->sp;
666 
667     if (IsSym(b)) {
668         *a = *b;
669     } else if (IsSym(c)) {
670         *a = *c;
671     } else {
672         err = slotDoubleVal(b, &lo);
673         if (err)
674             return err;
675         err = slotDoubleVal(c, &hi);
676         if (err)
677             return err;
678         SetRaw(a, sc_mod(slotRawFloat(a) - lo, hi - lo) + lo);
679     }
680     return errNone;
681 }
682 
mathWrapSignal(struct VMGlobals * g,int numArgsPushed)683 int mathWrapSignal(struct VMGlobals* g, int numArgsPushed) {
684     PyrSlot *a, *b, *c;
685     float lo, hi;
686     int err;
687     PyrObject* sig;
688 
689     a = g->sp - 2;
690     b = g->sp - 1;
691     c = g->sp;
692 
693     if (IsSym(b)) {
694         *a = *b;
695     } else if (IsSym(c)) {
696         *a = *c;
697     } else if (IsSignal(b) && IsSignal(c)) {
698         sig = signal_wrap_x(g, slotRawObject(a), slotRawObject(b), slotRawObject(c));
699         SetObject(a, sig);
700     } else {
701         err = slotFloatVal(b, &lo);
702         if (err)
703             return err;
704         err = slotFloatVal(c, &hi);
705         if (err)
706             return err;
707         sig = signal_wrap_f(g, slotRawObject(a), lo, hi);
708         SetObject(a, sig);
709     }
710     return errNone;
711 }
712 
mathFoldInt(struct VMGlobals * g,int numArgsPushed)713 int mathFoldInt(struct VMGlobals* g, int numArgsPushed) {
714     PyrSlot *a, *b, *c;
715     int err;
716 
717     a = g->sp - 2;
718     b = g->sp - 1;
719     c = g->sp;
720 
721     if (IsSym(b)) {
722         *a = *b;
723     } else if (IsSym(c)) {
724         *a = *c;
725     } else if (IsInt(b) && IsInt(c)) {
726         SetRaw(a, sc_fold(slotRawInt(a), slotRawInt(b), slotRawInt(c)));
727     } else {
728         double x, lo, hi;
729         x = slotRawInt(a);
730         err = slotDoubleVal(b, &lo);
731         if (err)
732             return err;
733         err = slotDoubleVal(c, &hi);
734         if (err)
735             return err;
736         SetFloat(a, sc_fold(x, lo, hi));
737     }
738     return errNone;
739 }
740 
mathFoldFloat(struct VMGlobals * g,int numArgsPushed)741 int mathFoldFloat(struct VMGlobals* g, int numArgsPushed) {
742     PyrSlot *a, *b, *c;
743     double lo, hi;
744     int err;
745 
746     a = g->sp - 2;
747     b = g->sp - 1;
748     c = g->sp;
749 
750     if (IsSym(b)) {
751         *a = *b;
752     } else if (IsSym(c)) {
753         *a = *c;
754     } else {
755         err = slotDoubleVal(b, &lo);
756         if (err)
757             return err;
758         err = slotDoubleVal(c, &hi);
759         if (err)
760             return err;
761         SetRaw(a, sc_fold(slotRawFloat(a), lo, hi));
762     }
763     return errNone;
764 }
765 
mathFoldSignal(struct VMGlobals * g,int numArgsPushed)766 int mathFoldSignal(struct VMGlobals* g, int numArgsPushed) {
767     PyrSlot *a, *b, *c;
768     float lo, hi;
769     int err;
770     PyrObject* sig;
771 
772     a = g->sp - 2;
773     b = g->sp - 1;
774     c = g->sp;
775 
776     if (IsSym(b)) {
777         *a = *b;
778     } else if (IsSym(c)) {
779         *a = *c;
780     } else if (IsSignal(b) && IsSignal(c)) {
781         sig = signal_fold_x(g, slotRawObject(a), slotRawObject(b), slotRawObject(c));
782         SetObject(a, sig);
783     } else {
784         err = slotFloatVal(b, &lo);
785         if (err)
786             return err;
787         err = slotFloatVal(c, &hi);
788         if (err)
789             return err;
790         sig = signal_fold_f(g, slotRawObject(a), lo, hi);
791         SetObject(a, sig);
792     }
793     return errNone;
794 }
795 
prSimpleNumberSeries(struct VMGlobals * g,int numArgsPushed)796 int prSimpleNumberSeries(struct VMGlobals* g, int numArgsPushed) {
797     PyrSlot* a = g->sp - 2;
798     PyrSlot* b = g->sp - 1;
799     PyrSlot* c = g->sp;
800 
801     int err, size;
802 
803     if (IsInt(a) && (IsInt(b) || IsNil(b)) && IsInt(c)) {
804         int first, second, last, step;
805         first = slotRawInt(a);
806         last = slotRawInt(c);
807         second = IsInt(b) ? slotRawInt(b) : (first < last ? first + 1 : first - 1);
808         step = second - first;
809 
810         if (step == 0)
811             size = 1;
812         else
813             size = ((last - first) / step) + 1;
814 
815         if ((size < 1) || ((step >= 0) && (last < first)) || ((step <= 0) && (last > first))) {
816             post("prSimpleNumberSeries: arguments do not form an arithmetic progression:\n\tfirst: %i, step: %i, last "
817                  "%i\n",
818                  first, step, last);
819             return errFailed;
820         }
821         if (size > INT_MAX_BY_PyrSlot) {
822             post("prSimpleNumberSeries: array size %i exceeds limit (%i)\n", size, INT_MAX_BY_PyrSlot);
823             return errFailed;
824         }
825 
826         PyrObject* obj = newPyrArray(g->gc, size, 0, true);
827         obj->size = size;
828         PyrSlot* slots = obj->slots;
829         if (step == 1) {
830             // Faster iteration for common case
831             if (first == 0) {
832                 for (int i = 0; i < size; ++i) {
833                     SetInt(slots + i, i);
834                 }
835             } else {
836                 for (int i = 0; i < size; ++i) {
837                     SetInt(slots + i, first++);
838                 }
839             }
840         } else {
841             int val = first;
842             for (int i = 0; i < size; ++i) {
843                 SetInt(slots + i, val);
844                 val += step;
845             }
846         }
847         SetObject(a, obj);
848     } else {
849         double first, second, last, step;
850         err = slotDoubleVal(a, &first);
851         if (err)
852             return err;
853         err = slotDoubleVal(c, &last);
854         if (err)
855             return err;
856         err = slotDoubleVal(b, &second);
857         if (err) {
858             if (first < last)
859                 second = first + 1.;
860             else
861                 second = first - 1.;
862         }
863 
864         step = second - first;
865         if (step == 0.f) {
866             size = 1;
867         } else {
868             size = (int)((last - first) / step) + 1;
869         }
870 
871         if ((size < 1) || ((step >= 0) && (last < first)) || ((step <= 0) && (last > first))) {
872             post("prSimpleNumberSeries: arguments do not form an arithmetic progression:\n\tfirst: %f, step: %f, last "
873                  "%f\n",
874                  first, step, last);
875             return errFailed;
876         }
877         if (size > INT_MAX_BY_PyrSlot) {
878             post("prSimpleNumberSeries: array size %i exceeds limit (%i)\n", size, INT_MAX_BY_PyrSlot);
879             return errFailed;
880         }
881         PyrObject* obj = newPyrArray(g->gc, size, 0, true);
882         obj->size = size;
883         PyrSlot* slots = obj->slots;
884         if (first == 0. && step == 1.) {
885             // Faster iteration for common case
886             for (long i = 0; i < size; ++i) {
887                 SetFloat(slots + i, i);
888             }
889         } else {
890             double val = first;
891             for (long i = 0; i < size; ++i) {
892                 val = first + step * i;
893                 SetFloat(slots + i, val);
894             }
895         }
896         SetObject(a, obj);
897     }
898     return errNone;
899 }
900 
901 /*
902 
903     asFraction {|maxDenominator=100|
904         var mediant, lower, upper, temp;
905         var n,d, k, k1;
906         if (this < 0) {
907             #n, d = this.neg.asFraction(maxDenominator);
908             ^[n.neg, d]
909         };
910         if (this < 1.0) {
911             upper = [1.0, this.reciprocal.floor];
912             lower = [1.0, upper[1]+1.0];
913         }{
914             lower = [this.floor, 1.0];
915             upper = [lower[0]+1.0, 1.0];
916         };
917         mediant = [lower[0] + upper[0], lower[1] + upper[1]];
918         loop {
919             mediant = [lower[0] + upper[0], lower[1] + upper[1]];
920             case
921             { (this * mediant[1]) > mediant[0] }
922             {
923                 if (maxDenominator < mediant[1]) {^upper};
924                 d = upper[0] - (this * upper[1]);
925                 if (d == 0) {^upper};
926                 lower = mediant;
927                 k = floor(((this * lower[1]) - lower[0]) / d);
928                 k1 = k + 1;
929                 temp = [lower[0] + (k1 * upper[0]), lower[1] + (k1 * upper[1])];
930                 lower = [lower[0] + (k * upper[0]), lower[1] + (k * upper[1])];
931                 upper = temp;
932             }
933             { (this * mediant[1]) == mediant[0] }
934             {
935                 if (maxDenominator >= mediant[1]) {^mediant};
936                 if (lower[1] < upper[1]) {^lower};
937                 ^upper
938             }
939             {
940                 if (maxDenominator < mediant[1]) {^lower};
941                 d = lower[0] - (this * lower[1]);
942                 if (d == 0) {^lower};
943                 upper = mediant;
944                 k = floor(((this * upper[1]) - upper[0]) / d);
945                 k1 = k + 1;
946                 temp = [(k1 * lower[0]) + upper[0], (k1 * lower[1]) + upper[1]];
947                 upper = [(k * lower[0]) + upper[0], (k * lower[1]) + upper[1]];
948                 lower = temp;
949             };
950         }
951     }
952 
953 */
954 
prAsFraction(struct VMGlobals * g,int numArgsPushed)955 int prAsFraction(struct VMGlobals* g, int numArgsPushed) {
956     PyrSlot* a = g->sp - 2;
957     PyrSlot* b = g->sp - 1;
958     PyrSlot* c = g->sp;
959 
960     double mediant_num, lower_num, upper_num, temp_num;
961     double mediant_den, lower_den, upper_den, temp_den;
962     double x, d;
963     int k, k1;
964     int maxDenominator;
965     int err;
966     bool neg = false;
967 
968     err = slotDoubleVal(a, &x);
969     if (err)
970         return err;
971 
972     err = slotIntVal(b, &maxDenominator);
973     if (err)
974         return err;
975 
976     bool faster = IsTrue(c);
977 
978     PyrObject* obj = newPyrArray(g->gc, 2, 0, true);
979     obj->size = 2;
980     PyrSlot* slots = obj->slots;
981     SetObject(a, obj);
982 
983     if (x < 0.0) {
984         x = -x;
985         neg = true;
986     }
987 
988     if (x == 0.0) {
989         SetInt(slots + 0, 0);
990         SetInt(slots + 1, 1);
991         return errNone;
992     }
993 
994     if (x < 1.0) {
995         upper_num = 1.0;
996         upper_den = floor(1. / x);
997         lower_num = 1.0;
998         lower_den = upper_den + 1.;
999     } else {
1000         lower_num = floor(x);
1001         lower_den = 1.0;
1002         upper_num = lower_num + 1.;
1003         upper_den = 1.0;
1004     }
1005 
1006     while (true) {
1007         mediant_num = lower_num + upper_num;
1008         mediant_den = lower_den + upper_den;
1009         // post("  md %g %g    %g %g    %g %g\n", mediant_num, mediant_den, lower_num, lower_den, upper_num, upper_den);
1010 
1011         if (x * mediant_den > mediant_num) {
1012             d = upper_num - (x * upper_den);
1013             if (maxDenominator < mediant_den || fabs(d) < 1e-5) {
1014                 if (neg)
1015                     upper_num = -upper_num;
1016                 SetInt(slots + 0, (int)upper_num);
1017                 SetInt(slots + 1, (int)upper_den);
1018                 return errNone;
1019             }
1020             lower_num = mediant_num;
1021             lower_den = mediant_den;
1022             if (faster) {
1023                 k = (int)floor(((x * lower_den) - lower_num) / d);
1024                 if (k < 10000) {
1025                     k1 = k + 1;
1026                     temp_num = lower_num + (k1 * upper_num);
1027                     temp_den = lower_den + (k1 * upper_den);
1028                     lower_num = lower_num + (k * upper_num);
1029                     lower_den = lower_den + (k * upper_den);
1030                     upper_num = temp_num;
1031                     upper_den = temp_den;
1032                 }
1033             }
1034         } else if (x * mediant_den == mediant_num) {
1035             if (maxDenominator >= mediant_den) {
1036                 if (neg)
1037                     mediant_num = -mediant_num;
1038                 SetInt(slots + 0, (int)mediant_num);
1039                 SetInt(slots + 1, (int)mediant_den);
1040                 return errNone;
1041             } else if (lower_den < upper_den) {
1042                 if (neg)
1043                     lower_num = -lower_num;
1044                 SetInt(slots + 0, (int)lower_num);
1045                 SetInt(slots + 1, (int)lower_den);
1046                 return errNone;
1047             } else {
1048                 if (neg)
1049                     upper_num = -upper_num;
1050                 SetInt(slots + 0, (int)upper_num);
1051                 SetInt(slots + 1, (int)upper_den);
1052                 return errNone;
1053             }
1054         } else {
1055             d = lower_num - (x * lower_den);
1056             if (maxDenominator < mediant_den || fabs(d) < 1e-5) {
1057                 if (neg)
1058                     lower_num = -lower_num;
1059                 SetInt(slots + 0, (int)lower_num);
1060                 SetInt(slots + 1, (int)lower_den);
1061                 return errNone;
1062             }
1063             upper_num = mediant_num;
1064             upper_den = mediant_den;
1065             if (faster) {
1066                 k = (int)floor(((x * upper_den) - upper_num) / d);
1067                 if (k < 10000) {
1068                     k1 = k + 1;
1069                     temp_num = (k1 * lower_num) + upper_num;
1070                     temp_den = (k1 * lower_den) + upper_den;
1071                     upper_num = (k * lower_num) + upper_num;
1072                     upper_den = (k * lower_den) + upper_den;
1073                     lower_num = temp_num;
1074                     lower_den = temp_den;
1075                 }
1076             }
1077         }
1078     }
1079 }
1080 
prSphericalHarmonic(struct VMGlobals * g,int numArgsPushed)1081 int prSphericalHarmonic(struct VMGlobals* g, int numArgsPushed) {
1082     PyrSlot* a = g->sp - 3; // n (spherical order)
1083     PyrSlot* b = g->sp - 2; // m (index within the order)
1084     PyrSlot* c = g->sp - 1; // theta
1085     PyrSlot* d = g->sp; // phi
1086 
1087     int err;
1088 
1089     int n = 0;
1090     int m = 0;
1091     double theta = 0.0;
1092     double phi = 0.0;
1093 
1094     err = slotDoubleVal(d, &phi);
1095     if (err)
1096         return err;
1097     err = slotDoubleVal(c, &theta);
1098     if (err)
1099         return err;
1100     err = slotIntVal(b, &m);
1101     if (err)
1102         return err;
1103     err = slotIntVal(a, &n);
1104     if (err)
1105         return err;
1106 
1107     auto res_cmplx = boost::math::spherical_harmonic(n, m, theta, phi);
1108 
1109     PyrObject* obj = instantiateObject(gMainVMGlobals->gc, getsym("Complex")->u.classobj, 0, true, true);
1110     SetObject(a, obj);
1111 
1112     PyrSlot* slots = obj->slots;
1113     SetFloat(slots + 0, res_cmplx.real());
1114     SetFloat(slots + 1, res_cmplx.imag());
1115 
1116     return errNone;
1117 }
1118 
1119 template <typename Arg1T, typename Arg2T, std::complex<double> BoostFunctionT(Arg1T, Arg2T)>
prBoostTwoArgRetComplex(struct VMGlobals * g,int numArgsPushed)1120 int prBoostTwoArgRetComplex(struct VMGlobals* g, int numArgsPushed) {
1121     PyrSlot* a = g->sp - 1;
1122     PyrSlot* b = g->sp;
1123 
1124     Arg1T arg1 = {};
1125     Arg2T arg2 = {};
1126 
1127     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1128     if (err)
1129         return err;
1130     err = slotVal(b, &arg2);
1131     if (err)
1132         return err;
1133 
1134     auto res_cmplx = BoostFunctionT(arg1, arg2);
1135 
1136     PyrObject* obj = instantiateObject(gMainVMGlobals->gc, getsym("Complex")->u.classobj, 0, true, true);
1137     SetObject(a, obj);
1138 
1139     PyrSlot* slots = obj->slots;
1140     SetFloat(slots + 0, res_cmplx.real());
1141     SetFloat(slots + 1, res_cmplx.imag());
1142 
1143     return errNone;
1144 }
1145 
1146 template <typename Arg1T, std::vector<double> BoostFunctionT(Arg1T)>
prBoostOneArgRetArray(struct VMGlobals * g,int numArgsPushed)1147 int prBoostOneArgRetArray(struct VMGlobals* g, int numArgsPushed) {
1148     PyrSlot* a = g->sp;
1149 
1150     Arg1T arg1;
1151 
1152     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1153     if (err)
1154         return err;
1155 
1156     auto res = BoostFunctionT(arg1);
1157 
1158     PyrObject* array = newPyrArray(g->gc, res.size(), 0, true);
1159     SetObject(a, array);
1160 
1161     PyrSlot* s = array->slots;
1162     for (auto item : res) {
1163         SetFloat(s, item);
1164         ++array->size;
1165         ++s;
1166     }
1167 
1168     return errNone;
1169 }
1170 
1171 template <typename RetT, RetT BoostFunctionT(unsigned, double const&)>
prBoostCheby(struct VMGlobals * g,int numArgsPushed)1172 int prBoostCheby(struct VMGlobals* g, int numArgsPushed) {
1173     PyrSlot* a = g->sp - 1;
1174     PyrSlot* b = g->sp;
1175 
1176     int arg1 = 0;
1177     double arg2 = 0.0;
1178 
1179     int err = slotIntVal(a, &arg1);
1180     if (err)
1181         return err;
1182     err = slotDoubleVal(b, &arg2);
1183     if (err)
1184         return err;
1185 
1186     RetT res = BoostFunctionT(arg1, arg2);
1187     setSlotVal(a, res);
1188 
1189     return errNone;
1190 }
1191 
1192 template <typename RetT, RetT BoostFunctionT(double const&)>
prBoostsqrt1pm1(struct VMGlobals * g,int numArgsPushed)1193 int prBoostsqrt1pm1(struct VMGlobals* g, int numArgsPushed) {
1194     PyrSlot* a = g->sp;
1195 
1196     double arg1 = 0;
1197 
1198     int err = slotDoubleVal(a, &arg1);
1199     if (err)
1200         return err;
1201 
1202     RetT res = BoostFunctionT(arg1);
1203     setSlotVal(a, res);
1204 
1205     return errNone;
1206 }
1207 
1208 template <typename RetT, typename Arg1T, RetT BoostFunctionT(Arg1T)>
prBoostOneArg(struct VMGlobals * g,int numArgsPushed)1209 int prBoostOneArg(struct VMGlobals* g, int numArgsPushed) {
1210     PyrSlot* a = g->sp;
1211 
1212     Arg1T arg1;
1213 
1214     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1215     if (err)
1216         return err;
1217 
1218     RetT res = BoostFunctionT(arg1);
1219     setSlotVal(a, res);
1220 
1221     return errNone;
1222 }
1223 
1224 template <typename RetT, typename Arg1T, typename Arg2T, RetT BoostFunctionT(Arg1T, Arg2T)>
prBoostTwoArg(struct VMGlobals * g,int numArgsPushed)1225 int prBoostTwoArg(struct VMGlobals* g, int numArgsPushed) {
1226     PyrSlot* a = g->sp - 1;
1227     PyrSlot* b = g->sp;
1228 
1229     Arg1T arg1 = {};
1230     Arg2T arg2 = {};
1231 
1232     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1233     if (err)
1234         return err;
1235     err = slotVal(b, &arg2);
1236     if (err)
1237         return err;
1238 
1239     RetT res = BoostFunctionT(arg1, arg2);
1240     setSlotVal(a, res);
1241 
1242     return errNone;
1243 }
1244 
1245 template <typename RetT, typename Arg1T, typename Arg2T, typename Arg3T, RetT BoostFunctionT(Arg1T, Arg2T, Arg3T)>
prBoostThreeArg(struct VMGlobals * g,int numArgsPushed)1246 int prBoostThreeArg(struct VMGlobals* g, int numArgsPushed) {
1247     PyrSlot* a = g->sp - 2;
1248     PyrSlot* b = g->sp - 1;
1249     PyrSlot* c = g->sp;
1250 
1251     Arg1T arg1 = {};
1252     Arg2T arg2 = {};
1253     Arg3T arg3 = {};
1254 
1255     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1256     if (err)
1257         return err;
1258     err = slotVal(b, &arg2);
1259     if (err)
1260         return err;
1261     err = slotVal(c, &arg3);
1262     if (err)
1263         return err;
1264 
1265     RetT res = BoostFunctionT(arg1, arg2, arg3);
1266     setSlotVal(a, res);
1267 
1268     return errNone;
1269 }
1270 
1271 template <typename RetT, typename Arg1T, typename Arg2T, typename Arg3T, typename Arg4T,
1272           RetT BoostFunctionT(Arg1T, Arg2T, Arg3T, Arg4T)>
prBoostFourArg(struct VMGlobals * g,int numArgsPushed)1273 int prBoostFourArg(struct VMGlobals* g, int numArgsPushed) {
1274     PyrSlot* a = g->sp - 3;
1275     PyrSlot* b = g->sp - 2;
1276     PyrSlot* c = g->sp - 1;
1277     PyrSlot* d = g->sp;
1278 
1279     Arg1T arg1 = {};
1280     Arg2T arg2 = {};
1281     Arg3T arg3 = {};
1282     Arg4T arg4 = {};
1283 
1284     int err = slotVal(a, &arg1); // slotVal only unwraps/extracts numeric values
1285     if (err)
1286         return err;
1287     err = slotVal(b, &arg2);
1288     if (err)
1289         return err;
1290     err = slotVal(c, &arg3);
1291     if (err)
1292         return err;
1293     err = slotVal(d, &arg4);
1294     if (err)
1295         return err;
1296 
1297     RetT res = BoostFunctionT(arg1, arg2, arg3, arg4);
1298     setSlotVal(a, res);
1299 
1300     return errNone;
1301 }
1302 
1303 
initMathPrimitives()1304 void initMathPrimitives() {
1305     int base, index;
1306 
1307     base = nextPrimitiveIndex();
1308     index = 0;
1309     definePrimitive(base, index++, "_AddInt", prAddInt, 2, 0);
1310     definePrimitive(base, index++, "_SubInt", prSubInt, 2, 0);
1311     definePrimitive(base, index++, "_MulInt", prMulInt, 2, 0);
1312     definePrimitive(base, index++, "_AddFloat", prAddFloat, 2, 0);
1313     definePrimitive(base, index++, "_SubFloat", prSubFloat, 2, 0);
1314     definePrimitive(base, index++, "_MulFloat", prMulFloat, 2, 0);
1315     definePrimitive(base, index++, "_NthPrime", prNthPrime, 1, 0);
1316     definePrimitive(base, index++, "_PrevPrime", prPrevPrime, 1, 0);
1317     definePrimitive(base, index++, "_NextPrime", prNextPrime, 1, 0);
1318     definePrimitive(base, index++, "_IsPrime", prIsPrime, 1, 0);
1319     definePrimitive(base, index++, "_IndexOfPrime", prIndexOfPrime, 1, 0);
1320     definePrimitive(base, index++, "_As32Bits", prAs32Bits, 1, 0);
1321     definePrimitive(base, index++, "_High32Bits", prHigh32Bits, 1, 0);
1322     definePrimitive(base, index++, "_Low32Bits", prLow32Bits, 1, 0);
1323     definePrimitive(base, index++, "_From32Bits", prFrom32Bits, 2, 0);
1324     definePrimitive(base, index++, "_From64Bits", prFrom64Bits, 3, 0);
1325 
1326     definePrimitive(base, index++, "_ClipInt", mathClipInt, 3, 0);
1327     definePrimitive(base, index++, "_ClipFloat", mathClipFloat, 3, 0);
1328     definePrimitive(base, index++, "_ClipSignal", mathClipSignal, 3, 0);
1329     definePrimitive(base, index++, "_WrapInt", mathWrapInt, 3, 0);
1330     definePrimitive(base, index++, "_WrapFloat", mathWrapFloat, 3, 0);
1331     definePrimitive(base, index++, "_WrapSignal", mathWrapSignal, 3, 0);
1332     definePrimitive(base, index++, "_FoldInt", mathFoldInt, 3, 0);
1333     definePrimitive(base, index++, "_FoldFloat", mathFoldFloat, 3, 0);
1334     definePrimitive(base, index++, "_FoldSignal", mathFoldSignal, 3, 0);
1335 
1336     definePrimitive(base, index++, "_SimpleNumberSeries", prSimpleNumberSeries, 3, 0);
1337     definePrimitive(base, index++, "_AsFraction", prAsFraction, 3, 0);
1338 
1339     //  Number Series:
1340     definePrimitive(base, index++, "_BernouliB2n", prBoostOneArg<double, int, boost::math::bernoulli_b2n<double>>, 1,
1341                     0);
1342     definePrimitive(base, index++, "_TangentT2n", prBoostOneArg<double, int, boost::math::tangent_t2n<double>>, 1, 0);
1343 
1344     //  Gamma:
1345     definePrimitive(base, index++, "_TGamma", prBoostOneArg<double, double, boost::math::tgamma<double>>, 1, 0);
1346     definePrimitive(base, index++, "_TGamma1pm1", prBoostOneArg<double, double, boost::math::tgamma1pm1<double>>, 1, 0);
1347     definePrimitive(base, index++, "_LGamma", prBoostOneArg<double, double, boost::math::lgamma<double>>, 1, 0);
1348     definePrimitive(base, index++, "_Digamma", prBoostOneArg<double, double, boost::math::digamma<double>>, 1, 0);
1349     definePrimitive(base, index++, "_Trigamma", prBoostOneArg<double, double, boost::math::trigamma<double>>, 1, 0);
1350     definePrimitive(base, index++, "_Polygamma", prBoostTwoArg<double, int, double, boost::math::polygamma<double>>, 2,
1351                     0);
1352     definePrimitive(base, index++, "_TGammaRatio",
1353                     prBoostTwoArg<double, double, double, boost::math::tgamma_ratio<double>>, 2, 0);
1354     definePrimitive(base, index++, "_TGammaDeltaRatio",
1355                     prBoostTwoArg<double, double, double, boost::math::tgamma_delta_ratio<double>>, 2, 0);
1356     // Incomplete Gamma Functions
1357     definePrimitive(base, index++, "_GammaP", prBoostTwoArg<double, double, double, boost::math::gamma_p<double>>, 2,
1358                     0);
1359     definePrimitive(base, index++, "_GammaQ", prBoostTwoArg<double, double, double, boost::math::gamma_q<double>>, 2,
1360                     0);
1361     definePrimitive(base, index++, "_TGammaLower",
1362                     prBoostTwoArg<double, double, double, boost::math::tgamma_lower<double>>, 2, 0);
1363     definePrimitive(base, index++, "_TGammaI", prBoostTwoArg<double, double, double, boost::math::tgamma<double>>, 2,
1364                     0);
1365     // Incomplete Gamma Function Inverses
1366     definePrimitive(base, index++, "_GammaPInv",
1367                     prBoostTwoArg<double, double, double, boost::math::gamma_p_inv<double>>, 2, 0);
1368     definePrimitive(base, index++, "_GammaQInv",
1369                     prBoostTwoArg<double, double, double, boost::math::gamma_q_inv<double>>, 2, 0);
1370     definePrimitive(base, index++, "_GammaPInvA",
1371                     prBoostTwoArg<double, double, double, boost::math::gamma_p_inva<double>>, 2, 0);
1372     definePrimitive(base, index++, "_GammaQInvA",
1373                     prBoostTwoArg<double, double, double, boost::math::gamma_q_inva<double>>, 2, 0);
1374     // Incomplete Gamma Function Derivative
1375     definePrimitive(base, index++, "_GammaPDerivative",
1376                     prBoostTwoArg<double, double, double, boost::math::gamma_p_derivative<double>>, 2, 0);
1377 
1378     //	Factorials and Binomial Coefficients:
1379     definePrimitive(base, index++, "_Factorial", prBoostOneArg<double, unsigned, boost::math::factorial<double>>, 1, 0);
1380     definePrimitive(base, index++, "_DoubleFactorial",
1381                     prBoostOneArg<double, unsigned int, boost::math::double_factorial<double>>, 1, 0);
1382     definePrimitive(base, index++, "_RisingFactorial",
1383                     prBoostTwoArg<double, double, int, boost::math::rising_factorial<double>>, 2, 0);
1384     definePrimitive(base, index++, "_FallingFactorial",
1385                     prBoostTwoArg<double, double, unsigned, boost::math::falling_factorial<double>>, 2, 0);
1386     definePrimitive(base, index++, "_BinomialCoefficient",
1387                     prBoostTwoArg<double, unsigned, unsigned, boost::math::binomial_coefficient<double>>, 2, 0);
1388 
1389     //	Beta Functions:
1390     definePrimitive(base, index++, "_Beta", prBoostTwoArg<double, double, double, boost::math::beta<double>>, 2, 0);
1391     // Incomplete Betas, normalized and non-normalized
1392     definePrimitive(base, index++, "_IBeta",
1393                     prBoostThreeArg<double, double, double, double, boost::math::ibeta<double>>, 3, 0);
1394     definePrimitive(base, index++, "_IBetaC",
1395                     prBoostThreeArg<double, double, double, double, boost::math::ibetac<double>>, 3, 0);
1396     definePrimitive(base, index++, "_BetaFull",
1397                     prBoostThreeArg<double, double, double, double, boost::math::beta<double>>, 3, 0);
1398     definePrimitive(base, index++, "_BetaCFull",
1399                     prBoostThreeArg<double, double, double, double, boost::math::betac<double>>, 3, 0);
1400     // Incomplete Betas, inverse
1401     definePrimitive(base, index++, "_IBetaInv",
1402                     prBoostThreeArg<double, double, double, double, boost::math::ibeta_inv<double>>, 3, 0);
1403     definePrimitive(base, index++, "_IBetaCInv",
1404                     prBoostThreeArg<double, double, double, double, boost::math::ibetac_inv<double>>, 3, 0);
1405     definePrimitive(base, index++, "_IBetaInvA",
1406                     prBoostThreeArg<double, double, double, double, boost::math::ibeta_inva<double>>, 3, 0);
1407     definePrimitive(base, index++, "_IBetaCInvA",
1408                     prBoostThreeArg<double, double, double, double, boost::math::ibetac_inva<double>>, 3, 0);
1409     definePrimitive(base, index++, "_IBetaInvB",
1410                     prBoostThreeArg<double, double, double, double, boost::math::ibeta_invb<double>>, 3, 0);
1411     definePrimitive(base, index++, "_IBetaCInvB",
1412                     prBoostThreeArg<double, double, double, double, boost::math::ibetac_invb<double>>, 3, 0);
1413     // 	Incomplete Betas, derivative
1414     definePrimitive(base, index++, "_IBetaDerivative",
1415                     prBoostThreeArg<double, double, double, double, boost::math::ibeta_derivative<double>>, 3, 0);
1416 
1417     //  Error functions:
1418     definePrimitive(base, index++, "_Erf", prBoostOneArg<double, double, boost::math::erf<double>>, 1, 0);
1419     definePrimitive(base, index++, "_ErfC", prBoostOneArg<double, double, boost::math::erfc<double>>, 1, 0);
1420     definePrimitive(base, index++, "_ErfInv", prBoostOneArg<double, double, boost::math::erf_inv<double>>, 1, 0);
1421     definePrimitive(base, index++, "_ErfCInv", prBoostOneArg<double, double, boost::math::erfc_inv<double>>, 1, 0);
1422 
1423     //  Polynomials:
1424     // Legendre (and Associated)
1425     definePrimitive(base, index++, "_LegendreP", prBoostTwoArg<double, int, double, boost::math::legendre_p<double>>, 2,
1426                     0);
1427     definePrimitive(base, index++, "_LegendrePPrime",
1428                     prBoostTwoArg<double, int, double, boost::math::legendre_p_prime<double>>, 2, 0);
1429     definePrimitive(base, index++, "_LegendrePZeros", prBoostOneArgRetArray<int, boost::math::legendre_p_zeros<double>>,
1430                     1, 0); // TODO: generalize return array form?
1431     definePrimitive(base, index++, "_LegendrePAssoc",
1432                     prBoostThreeArg<double, int, int, double, boost::math::legendre_p<double>>, 3, 0);
1433     definePrimitive(base, index++, "_LegendreQ",
1434                     prBoostTwoArg<double, unsigned, double, boost::math::legendre_q<double>>, 2, 0);
1435     // Laguerre (and Associated)
1436     definePrimitive(base, index++, "_Laguerre", prBoostTwoArg<double, unsigned, double, boost::math::laguerre<double>>,
1437                     2, 0);
1438     definePrimitive(base, index++, "_LaguerreAssoc",
1439                     prBoostThreeArg<double, unsigned, double, double, boost::math::laguerre<double>>, 3, 0);
1440     // Hermite
1441     definePrimitive(base, index++, "_Hermite", prBoostTwoArg<double, unsigned, double, boost::math::hermite<double>>, 2,
1442                     0);
1443     // Chebyshev
1444 
1445     definePrimitive(base, index++, "_ChebyshevT", prBoostCheby<double, boost::math::chebyshev_t<double>>, 2, 0);
1446     definePrimitive(base, index++, "_ChebyshevU", prBoostCheby<double, boost::math::chebyshev_u<double>>, 2, 0);
1447     definePrimitive(base, index++, "_ChebyshevTPrime", prBoostCheby<double, boost::math::chebyshev_t_prime<double>>, 2,
1448                     0);
1449     // Spherical Harmonics
1450     definePrimitive(base, index++, "_SphericalHarmonic", prSphericalHarmonic, 4, 0);
1451     definePrimitive(base, index++, "_SphericalHarmonicR",
1452                     prBoostFourArg<double, unsigned, int, double, double, boost::math::spherical_harmonic_r<double>>, 4,
1453                     0);
1454     definePrimitive(base, index++, "_SphericalHarmonicI",
1455                     prBoostFourArg<double, unsigned, int, double, double, boost::math::spherical_harmonic_i<double>>, 4,
1456                     0);
1457 
1458     //	Bessel Functions:
1459     // First and Second Kinds
1460     definePrimitive(base, index++, "_CylBesselJ",
1461                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_j<double>>, 2, 0);
1462     definePrimitive(base, index++, "_CylNeumann",
1463                     prBoostTwoArg<double, double, double, boost::math::cyl_neumann<double>>, 2, 0);
1464     // Zero finder
1465     definePrimitive(base, index++, "_CylBesselJZero",
1466                     prBoostTwoArg<double, double, int, boost::math::cyl_bessel_j_zero<double>>, 2, 0);
1467     definePrimitive(base, index++, "_CylNeumannZero",
1468                     prBoostTwoArg<double, double, int, boost::math::cyl_neumann_zero<double>>, 2, 0);
1469     // Modified, First and Second Kinds
1470     definePrimitive(base, index++, "_CylBesselI",
1471                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_i<double>>, 2, 0);
1472     definePrimitive(base, index++, "_CylBesselK",
1473                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_k<double>>, 2, 0);
1474     // Spherical, First and Second Kinds
1475     definePrimitive(base, index++, "_SphBessel",
1476                     prBoostTwoArg<double, unsigned, double, boost::math::sph_bessel<double>>, 2, 0);
1477     definePrimitive(base, index++, "_SphNeumann",
1478                     prBoostTwoArg<double, unsigned, double, boost::math::sph_neumann<double>>, 2, 0);
1479     // Derivatives
1480     definePrimitive(base, index++, "_CylBesselJPrime",
1481                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_j_prime<double>>, 2, 0);
1482     definePrimitive(base, index++, "_CylNeumannPrime",
1483                     prBoostTwoArg<double, double, double, boost::math::cyl_neumann_prime<double>>, 2, 0);
1484     definePrimitive(base, index++, "_CylBesselIPrime",
1485                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_i_prime<double>>, 2, 0);
1486     definePrimitive(base, index++, "_CylBesselKPrime",
1487                     prBoostTwoArg<double, double, double, boost::math::cyl_bessel_k_prime<double>>, 2, 0);
1488     definePrimitive(base, index++, "_SphBesselPrime",
1489                     prBoostTwoArg<double, unsigned, double, boost::math::sph_bessel_prime<double>>, 2, 0);
1490     definePrimitive(base, index++, "_SphNeumannPrime",
1491                     prBoostTwoArg<double, unsigned, double, boost::math::sph_neumann_prime<double>>, 2, 0);
1492 
1493     //  Hankel Functions:
1494     // Cyclic
1495     definePrimitive(base, index++, "_CylHankel1",
1496                     prBoostTwoArgRetComplex<double, double, boost::math::cyl_hankel_1<double>>, 2, 0);
1497     definePrimitive(base, index++, "_CylHankel2",
1498                     prBoostTwoArgRetComplex<double, double, boost::math::cyl_hankel_2<double>>, 2, 0);
1499     // Spherical
1500     definePrimitive(base, index++, "_SphHankel1",
1501                     prBoostTwoArgRetComplex<double, double, boost::math::sph_hankel_1<double>>, 2, 0);
1502     definePrimitive(base, index++, "_SphHankel2",
1503                     prBoostTwoArgRetComplex<double, double, boost::math::sph_hankel_2<double>>, 2, 0);
1504 
1505     //  Airy Functions:
1506     definePrimitive(base, index++, "_AiryAi", prBoostOneArg<double, double, boost::math::airy_ai<double>>, 1, 0);
1507     definePrimitive(base, index++, "_AiryBi", prBoostOneArg<double, double, boost::math::airy_bi<double>>, 1, 0);
1508     definePrimitive(base, index++, "_AiryAiPrime", prBoostOneArg<double, double, boost::math::airy_ai_prime<double>>, 1,
1509                     0);
1510     definePrimitive(base, index++, "_AiryBiPrime", prBoostOneArg<double, double, boost::math::airy_bi_prime<double>>, 1,
1511                     0);
1512     definePrimitive(base, index++, "_AiryAiZero", prBoostOneArg<double, int, boost::math::airy_ai_zero<double>>, 1, 0);
1513     definePrimitive(base, index++, "_AiryBiZero", prBoostOneArg<double, int, boost::math::airy_bi_zero<double>>, 1, 0);
1514 
1515     //  Elliptic Integrals:
1516     // Elliptic Integrals - Carlson Form
1517     definePrimitive(base, index++, "_EllintRf",
1518                     prBoostThreeArg<double, double, double, double, boost::math::ellint_rf<double>>, 3, 0);
1519     definePrimitive(base, index++, "_EllintRd",
1520                     prBoostThreeArg<double, double, double, double, boost::math::ellint_rd<double>>, 3, 0);
1521     definePrimitive(base, index++, "_EllintRj",
1522                     prBoostFourArg<double, double, double, double, double, boost::math::ellint_rj<double>>, 4, 0);
1523     definePrimitive(base, index++, "_EllintRc", prBoostTwoArg<double, double, double, boost::math::ellint_rc<double>>,
1524                     2, 0);
1525     definePrimitive(base, index++, "_EllintRg",
1526                     prBoostThreeArg<double, double, double, double, boost::math::ellint_rg<double>>, 3, 0);
1527     // Elliptic Integrals of the First, Second, Third Kind, D - Legendre Form
1528     definePrimitive(base, index++, "_Ellint1", prBoostTwoArg<double, double, double, boost::math::ellint_1<double>>, 2,
1529                     0);
1530     definePrimitive(base, index++, "_Ellint1C", prBoostOneArg<double, double, boost::math::ellint_1<double>>, 1, 0);
1531     definePrimitive(base, index++, "_Ellint2", prBoostTwoArg<double, double, double, boost::math::ellint_2<double>>, 2,
1532                     0);
1533     definePrimitive(base, index++, "_Ellint2C", prBoostOneArg<double, double, boost::math::ellint_2<double>>, 1, 0);
1534     definePrimitive(base, index++, "_Ellint3",
1535                     prBoostThreeArg<double, double, double, double, boost::math::ellint_3<double>>, 3, 0);
1536     definePrimitive(base, index++, "_Ellint3C", prBoostTwoArg<double, double, double, boost::math::ellint_3<double>>, 2,
1537                     0);
1538     definePrimitive(base, index++, "_EllintD", prBoostTwoArg<double, double, double, boost::math::ellint_d<double>>, 2,
1539                     0);
1540     definePrimitive(base, index++, "_EllintDC", prBoostOneArg<double, double, boost::math::ellint_d<double>>, 1, 0);
1541     // Jacobi Zeta, Heuman Lambda Function
1542     definePrimitive(base, index++, "_JacobiZeta",
1543                     prBoostTwoArg<double, double, double, boost::math::jacobi_zeta<double>>, 2, 0);
1544     definePrimitive(base, index++, "_HeumanLambda",
1545                     prBoostTwoArg<double, double, double, boost::math::heuman_lambda<double>>, 2, 0);
1546 
1547 
1548     //  Jacobi Elliptic Functions:
1549     definePrimitive(base, index++, "_JacobiCd", prBoostTwoArg<double, double, double, boost::math::jacobi_cd<double>>,
1550                     2, 0);
1551     definePrimitive(base, index++, "_JacobiCn", prBoostTwoArg<double, double, double, boost::math::jacobi_cn<double>>,
1552                     2, 0);
1553     definePrimitive(base, index++, "_JacobiCs", prBoostTwoArg<double, double, double, boost::math::jacobi_cs<double>>,
1554                     2, 0);
1555     definePrimitive(base, index++, "_JacobiDc", prBoostTwoArg<double, double, double, boost::math::jacobi_dc<double>>,
1556                     2, 0);
1557     definePrimitive(base, index++, "_JacobiDn", prBoostTwoArg<double, double, double, boost::math::jacobi_dn<double>>,
1558                     2, 0);
1559     definePrimitive(base, index++, "_JacobiDs", prBoostTwoArg<double, double, double, boost::math::jacobi_ds<double>>,
1560                     2, 0);
1561     definePrimitive(base, index++, "_JacobiNc", prBoostTwoArg<double, double, double, boost::math::jacobi_nc<double>>,
1562                     2, 0);
1563     definePrimitive(base, index++, "_JacobiNd", prBoostTwoArg<double, double, double, boost::math::jacobi_nd<double>>,
1564                     2, 0);
1565     definePrimitive(base, index++, "_JacobiNs", prBoostTwoArg<double, double, double, boost::math::jacobi_ns<double>>,
1566                     2, 0);
1567     definePrimitive(base, index++, "_JacobiSc", prBoostTwoArg<double, double, double, boost::math::jacobi_sc<double>>,
1568                     2, 0);
1569     definePrimitive(base, index++, "_JacobiSd", prBoostTwoArg<double, double, double, boost::math::jacobi_sd<double>>,
1570                     2, 0);
1571     definePrimitive(base, index++, "_JacobiSn", prBoostTwoArg<double, double, double, boost::math::jacobi_sn<double>>,
1572                     2, 0);
1573 
1574     //  Zeta Function:
1575     definePrimitive(base, index++, "_Zeta", prBoostOneArg<double, double, boost::math::zeta<double>>, 1, 0);
1576 
1577     //  Exponential Integrals:
1578     definePrimitive(base, index++, "_ExpintEn", prBoostTwoArg<double, double, double, boost::math::expint<double>>, 2,
1579                     0);
1580     definePrimitive(base, index++, "_ExpintEi", prBoostOneArg<double, double, boost::math::expint<double>>, 1, 0);
1581 
1582     //  Basic functions:
1583     definePrimitive(base, index++, "_SinPi", prBoostOneArg<double, double, boost::math::sin_pi<double>>, 1, 0);
1584     definePrimitive(base, index++, "_CosPi", prBoostOneArg<double, double, boost::math::cos_pi<double>>, 1, 0);
1585     definePrimitive(base, index++, "_Log1p", prBoostOneArg<double, double, boost::math::log1p<double>>, 1, 0);
1586     definePrimitive(base, index++, "_ExpM1", prBoostOneArg<double, double, boost::math::expm1<double>>, 1, 0);
1587     definePrimitive(base, index++, "_Cbrt", prBoostOneArg<double, double, boost::math::cbrt<double>>, 1, 0);
1588     definePrimitive(base, index++, "_Sqrt1pm1", prBoostsqrt1pm1<double, boost::math::sqrt1pm1<double>>, 1, 0);
1589     definePrimitive(base, index++, "_PowM1", prBoostTwoArg<double, double, double, boost::math::powm1<double>>, 2, 0);
1590     definePrimitive(base, index++, "_Hypot2", prBoostTwoArg<double, double, double, boost::math::hypot<double>>, 2, 0);
1591 
1592     // Sinus Cardinal and Hyperbolic Sinus Cardinal Functions:
1593     definePrimitive(base, index++, "_SincPi", prBoostOneArg<double, double, boost::math::sinc_pi<double>>, 1, 0);
1594     definePrimitive(base, index++, "_SinhcPi", prBoostOneArg<double, double, boost::math::sinhc_pi<double>>, 1, 0);
1595 
1596     // Inverse Hyperbolic Functions:
1597     definePrimitive(base, index++, "_Asinh", prBoostOneArg<double, double, boost::math::asinh<double>>, 1, 0);
1598     definePrimitive(base, index++, "_Acosh", prBoostOneArg<double, double, boost::math::acosh<double>>, 1, 0);
1599     definePrimitive(base, index++, "_Atanh", prBoostOneArg<double, double, boost::math::atanh<double>>, 1, 0);
1600 
1601     //	Owen's T function:
1602     definePrimitive(base, index++, "_OwensT", prBoostTwoArg<double, double, double, boost::math::owens_t<double>>, 2,
1603                     0);
1604 }
1605