/*****
 * 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);
#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) (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);
}