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