/***** * runmath.in * * Runtime functions for math operations. * *****/ pair => primPair() realarray* => realArray() pairarray* => pairArray() #include #include "mathop.h" #include "path.h" #ifdef __CYGWIN__ extern "C" double yn(int, double); extern "C" double jn(int, double); #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()(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) (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); }