/*****
* runmath.in
*
* Runtime functions for math operations.
*
*****/
pair => primPair()
realarray* => realArray()
pairarray* => pairArray()
#include <inttypes.h>
#include "mathop.h"
#include "path.h"
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);
// 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)
{
if(y == 0) dividebyzero();
if(y == -1) return Negate(x);
// Implementation-independent definition of integer division: round down
return (x-portableMod(x,y))/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())/RAND_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((uint32_t) a > 0xFFFFFFFF) return -1;
return CLZ((uint32_t) a);
}
Int CTZ(Int a)
{
if((uint32_t) a > 0xFFFFFFFF) return -1;
#if __GNUC_PREREQ(3,4)
return __builtin_ctz(a);
#else
// find the number of trailing zeros in a 32-bit number
static const int MultiplyDeBruijnBitPosition[32] = {
0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
};
return MultiplyDeBruijnBitPosition[((uint32_t)((a & -a) * 0x077CB531U))
>> 27];
#endif
}