summaryrefslogtreecommitdiff
path: root/graphics/asymptote/runmath.in
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /graphics/asymptote/runmath.in
Initial commit
Diffstat (limited to 'graphics/asymptote/runmath.in')
-rw-r--r--graphics/asymptote/runmath.in451
1 files changed, 451 insertions, 0 deletions
diff --git a/graphics/asymptote/runmath.in b/graphics/asymptote/runmath.in
new file mode 100644
index 0000000000..fa9dbf3e69
--- /dev/null
+++ b/graphics/asymptote/runmath.in
@@ -0,0 +1,451 @@
+/*****
+ * runmath.in
+ *
+ * Runtime functions for math operations.
+ *
+ *****/
+
+pair => primPair()
+realarray* => realArray()
+pairarray* => pairArray()
+
+#include <inttypes.h>
+
+#include "mathop.h"
+#include "path.h"
+
+#ifdef __CYGWIN__
+extern "C" double yn(int, double);
+extern "C" double jn(int, double);
+extern "C" int __signgam;
+#define signgam __signgam
+#endif
+
+using namespace camp;
+
+typedef array realarray;
+typedef array pairarray;
+
+using types::realArray;
+using types::pairArray;
+
+using run::integeroverflow;
+using vm::frame;
+
+const char *invalidargument="invalid argument";
+
+extern uint32_t CLZ(uint32_t a);
+
+inline unsigned intbits() {
+ static unsigned count=0;
+ if(count > 0) return count;
+ while((1ULL << count) < Int_MAX)
+ ++count;
+ ++count;
+ return count;
+}
+
+static const unsigned char BitReverseTable8[256]=
+{
+#define R2(n) n, n+2*64, n+1*64, n+3*64
+#define R4(n) R2(n),R2(n+2*16),R2(n+1*16),R2(n+3*16)
+#define R6(n) R4(n),R4(n+2*4 ),R4(n+1*4 ),R4(n+3*4 )
+R6(0),R6(2),R6(1),R6(3)
+};
+#undef R2
+#undef R4
+#undef R6
+
+unsigned long long bitreverse8(unsigned long long a)
+{
+ return
+ (unsigned long long) BitReverseTable8[a];
+}
+
+unsigned long long bitreverse16(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 8)]);
+}
+
+unsigned long long bitreverse24(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 16)]);
+}
+
+unsigned long long bitreverse32(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 24) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 24)]);
+}
+
+unsigned long long bitreverse40(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 32) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 24) |
+ ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 32)]);
+}
+
+unsigned long long bitreverse48(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 40) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 32) |
+ ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 24) |
+ ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 40)]);
+}
+
+unsigned long long bitreverse56(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 48) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 40) |
+ ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 32) |
+ ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 24) |
+ ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 40) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 48)]);
+}
+
+unsigned long long bitreverse64(unsigned long long a)
+{
+ return
+ ((unsigned long long) BitReverseTable8[a & 0xff] << 56) |
+ ((unsigned long long) BitReverseTable8[(a >> 8) & 0xff] << 48) |
+ ((unsigned long long) BitReverseTable8[(a >> 16) & 0xff] << 40) |
+ ((unsigned long long) BitReverseTable8[(a >> 24) & 0xff] << 32) |
+ ((unsigned long long) BitReverseTable8[(a >> 32) & 0xff] << 24) |
+ ((unsigned long long) BitReverseTable8[(a >> 40) & 0xff] << 16) |
+ ((unsigned long long) BitReverseTable8[(a >> 48) & 0xff] << 8) |
+ ((unsigned long long) BitReverseTable8[(a >> 56)]);
+}
+
+// From Warren, Jr., Henry S. (2013) [2002]. Hacker's Delight (2 ed.).
+// Addison Wesley - Pearson Education, Inc. pp. 81-96.
+Int popcount(Int a)
+{
+ const uint64_t m1 = 0x5555555555555555; //binary: 0101...
+ const uint64_t m2 = 0x3333333333333333; //binary: 00110011..
+ const uint64_t m4 = 0x0f0f0f0f0f0f0f0f; //binary: 4 zeros, 4 ones ...
+ const uint64_t h01 = 0x0101010101010101; //the sum of 256 to the power of 0,1,2,3...
+
+// This algorithm uses 12 arithmetic operations, one of which is a multiply.
+ a -= (a >> 1) & m1; //put count of each 2 bits into those 2 bits
+ a=(a & m2)+((a >> 2) & m2); //put count of each 4 bits into those 4 bits
+ a=(a+(a >> 4)) & m4; //put count of each 8 bits into those 8 bits
+ return (a * h01) >> 56; //returns left 8 bits of a+(a << 8)+(a << 16)+(a << 24)+...
+}
+
+// Return the factorial of a non-negative integer using a lookup table.
+Int factorial(Int n)
+{
+ static Int *table;
+ static Int size=0;
+ if(size == 0) {
+ Int f=1;
+ size=2;
+ while(f <= Int_MAX/size)
+ f *= (size++);
+ table=new Int[size];
+ table[0]=f=1;
+ for(Int i=1; i < size; ++i) {
+ f *= i;
+ table[i]=f;
+ }
+ }
+ if(n >= size) integeroverflow(0);
+ return table[n];
+}
+
+static inline Int Round(double x)
+{
+ return Int(x+((x >= 0) ? 0.5 : -0.5));
+}
+
+inline Int sgn(double x)
+{
+ return (x > 0.0 ? 1 : (x < 0.0 ? -1 : 0));
+}
+
+static bool initializeRandom=true;
+
+void Srand(Int seed)
+{
+ initializeRandom=false;
+ const int n=256;
+ static char state[n];
+ initstate(intcast(seed),state,n);
+}
+
+// Autogenerated routines:
+
+
+real ^(real x, Int y)
+{
+ return pow(x,y);
+}
+
+pair ^(pair z, Int y)
+{
+ return pow(z,y);
+}
+
+Int quotient(Int x, Int y)
+{
+ return quotient<Int>()(x,y);
+}
+
+Int abs(Int x)
+{
+ return Abs(x);
+}
+
+Int sgn(real x)
+{
+ return sgn(x);
+}
+
+Int rand()
+{
+ if(initializeRandom)
+ Srand(1);
+ return random();
+}
+
+void srand(Int seed)
+{
+ Srand(seed);
+}
+
+// a random number uniformly distributed in the interval [0,1]
+real unitrand()
+{
+ return ((real) random())/RANDOM_MAX;
+}
+
+Int ceil(real x)
+{
+ return Intcast(ceil(x));
+}
+
+Int floor(real x)
+{
+ return Intcast(floor(x));
+}
+
+Int round(real x)
+{
+ if(validInt(x)) return Round(x);
+ integeroverflow(0);
+}
+
+Int Ceil(real x)
+{
+ return Ceil(x);
+}
+
+Int Floor(real x)
+{
+ return Floor(x);
+}
+
+Int Round(real x)
+{
+ return Round(Intcap(x));
+}
+
+real fmod(real x, real y)
+{
+ if (y == 0.0) dividebyzero();
+ return fmod(x,y);
+}
+
+real atan2(real y, real x)
+{
+ return atan2(y,x);
+}
+
+real hypot(real x, real y)
+{
+ return hypot(x,y);
+}
+
+real remainder(real x, real y)
+{
+ return remainder(x,y);
+}
+
+real Jn(Int n, real x)
+{
+ return jn(n,x);
+}
+
+real Yn(Int n, real x)
+{
+ return yn(n,x);
+}
+
+real erf(real x)
+{
+ return erf(x);
+}
+
+real erfc(real x)
+{
+ return erfc(x);
+}
+
+Int factorial(Int n) {
+ if(n < 0) error(invalidargument);
+ return factorial(n);
+}
+
+Int choose(Int n, Int k) {
+ if(n < 0 || k < 0 || k > n) error(invalidargument);
+ Int f=1;
+ Int r=n-k;
+ for(Int i=n; i > r; --i) {
+ if(f > Int_MAX/i) integeroverflow(0);
+ f=(f*i)/(n-i+1);
+ }
+ return f;
+}
+
+real gamma(real x)
+{
+#ifdef HAVE_TGAMMA
+ return tgamma(x);
+#else
+ real lg = lgamma(x);
+ return signgam*exp(lg);
+#endif
+}
+
+realarray *quadraticroots(real a, real b, real c)
+{
+ quadraticroots q(a,b,c);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.t1;
+ if(q.roots == 2) (*roots)[1]=q.t2;
+ return roots;
+}
+
+pairarray *quadraticroots(explicit pair a, explicit pair b, explicit pair c)
+{
+ Quadraticroots q(a,b,c);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.z1;
+ if(q.roots == 2) (*roots)[1]=q.z2;
+ return roots;
+}
+
+realarray *cubicroots(real a, real b, real c, real d)
+{
+ cubicroots q(a,b,c,d);
+ array *roots=new array(q.roots);
+ if(q.roots >= 1) (*roots)[0]=q.t1;
+ if(q.roots >= 2) (*roots)[1]=q.t2;
+ if(q.roots == 3) (*roots)[2]=q.t3;
+ return roots;
+}
+
+
+// Logical operations
+
+bool !(bool b)
+{
+ return !b;
+}
+
+bool :boolMemEq(frame *a, frame *b)
+{
+ return a == b;
+}
+
+bool :boolMemNeq(frame *a, frame *b)
+{
+ return a != b;
+}
+
+bool :boolFuncEq(callable *a, callable *b)
+{
+ return a->compare(b);
+}
+
+bool :boolFuncNeq(callable *a, callable *b)
+{
+ return !(a->compare(b));
+}
+
+
+// Bit operations
+
+Int AND(Int a, Int b)
+{
+ return a & b;
+}
+
+Int OR(Int a, Int b)
+{
+ return a | b;
+}
+
+Int XOR(Int a, Int b)
+{
+ return a ^ b;
+}
+
+Int NOT(Int a)
+{
+ return ~a;
+}
+
+Int CLZ(Int a)
+{
+ if((unsigned long long) a > 0xFFFFFFFF)
+ return CLZ((uint32_t) ((unsigned long long) a >> 32));
+ else {
+ int bits=intbits();
+ if(a != 0) return bits-32+CLZ((uint32_t) a);
+ return bits;
+ }
+}
+
+Int popcount(Int a)
+{
+ return popcount(a);
+}
+
+Int CTZ(Int a)
+{
+ return popcount((a&-a)-1);
+}
+
+// bitreverse a within a word of length bits.
+Int bitreverse(Int a, Int bits)
+{
+ typedef unsigned long long Bitreverse(unsigned long long a);
+ static Bitreverse *B[]={bitreverse8,bitreverse16,bitreverse24,bitreverse32,
+ bitreverse40,bitreverse48,bitreverse56,bitreverse64};
+ int maxbits=intbits()-1; // Drop sign bit
+#if Int_MAX2 >= 0x7fffffffffffffffLL
+ --maxbits; // Drop extra bit for reserved values
+#endif
+ if(bits <= 0 || bits > maxbits || a < 0 ||
+ (unsigned long long) a >= (1ULL << bits))
+ return -1;
+ unsigned int bytes=(bits+7)/8;
+ return B[bytes-1]((unsigned long long) a) >> (8*bytes-bits);
+}