/*****
 * runpair.in
 *
 * Runtime functions for pair operations.
 *
 *****/

pair   => primPair()

#include "pair.h"

using namespace camp;

namespace run {
extern pair zero;
}

pair sin(pair z)
{
  return pair(sin(z.getx())*cosh(z.gety()),cos(z.getx())*sinh(z.gety()));
}

pair exp(pair z)
{
  return exp(z.getx())*expi(z.gety());
}

pair gamma(pair z)
{
  static double p[]={0.99999999999980993,676.5203681218851,-1259.1392167224028,
                     771.32342877765313,-176.61502916214059,12.507343278686905,
                     -0.13857109526572012,9.9843695780195716e-6,
                     1.5056327351493116e-7};
  static int n=sizeof(p)/sizeof(double);
  static double root2pi=sqrt(2*PI);
  if(z.getx() < 0.5)
    return PI/(sin(PI*z)*gamma(1.0-z));
  z -= 1.0;
  pair x=p[0];
  for(int i=1; i < n; ++i)
    x += p[i]/(z+i);
  pair t=n-1.5+z;
  return root2pi*pow(t,z+0.5)*exp(-t)*x;
}

// Autogenerated routines:


pair :pairZero()
{
  return zero;
}

pair :realRealToPair(real x, real y)
{
  return pair(x,y);
}

pair :pairNegate(pair z)
{
  return -z;
}

real xpart:pairXPart(pair z)
{
  return z.getx();
}

real ypart:pairYPart(pair z)
{
  return z.gety();
}

real length(pair z)
{
  return z.length();
}

real abs(pair z)
{
  return z.length();
}

real abs2(pair z)
{
  return z.abs2();
}

pair sqrt(explicit pair z)
{
  return Sqrt(z);
}

// Return the angle of z in radians.
real angle(pair z, bool warn=true)
{
  return z.angle(warn);
}

// Return the angle of z in degrees in the interval [0,360).
real degrees(pair z, bool warn=true)
{
  return principalBranch(degrees(z.angle(warn)));
}

// Convert degrees to radians.
real radians(real degrees)
{
  return radians(degrees);
}

// Convert radians to degrees.
real degrees(real radians)
{
  return degrees(radians);
}

// Convert radians to degrees in [0,360).
real Degrees(real radians)
{
  return principalBranch(degrees(radians));
}

real Sin(real deg)
{
  int n=(int) (deg/90.0);
  if(deg == n*90.0) {
    int m=n % 4;
    if(m < 0) m += 4;
    if(m == 1) return 1;
    if(m == 3) return -1;
    return 0.0;
  }
  return sin(radians(deg));
}

real Cos(real deg)
{
  int n=(int) (deg/90.0);
  if(deg == n*90.0) {
    int m=n % 4;
    if(m < 0) m += 4;
    if(m == 0) return 1;
    if(m == 2) return -1;
    return 0.0;
  }
  return cos(radians(deg));
}

real Tan(real deg)
{
  int n=(int) (deg/90.0);
  if(deg == n*90.0) {
    int m=n % 4;
    if(m < 0) m += 4;
    if(m == 1) return HUGE_VAL;
    if(m == 3) return -HUGE_VAL;
    return 0.0;
  }
  return tan(radians(deg));
}

real aSin(real x)
{
  return degrees(asin(x));
}

real aCos(real x)
{
  return degrees(acos(x));
}

real aTan(real x)
{
  return degrees(atan(x));
}

pair unit(pair z)
{
  return unit(z);
}

pair dir(real degrees)
{
  return expi(radians(degrees));
}

pair dir(explicit pair z)
{
  return unit(z);
}

pair expi(real angle)
{
  return expi(angle);
}

pair exp(explicit pair z)
{
  return exp(z);
}

pair log(explicit pair z)
{
  return pair(log(z.length()),z.angle());
}

pair sin(explicit pair z)
{
  return sin(z);
}

pair cos(explicit pair z)
{
  return pair(cos(z.getx())*cosh(z.gety()),-sin(z.getx())*sinh(z.gety()));
}

// Complex Gamma function
pair gamma(explicit pair z)
{
  return gamma(z);
}

pair conj(pair z)
{
  return conj(z);
}

pair realmult(pair z, pair w)
{
  return pair(z.getx()*w.getx(),z.gety()*w.gety());
}

// To avoid confusion, a dot product requires explicit pair arguments.
real dot(explicit pair z, explicit pair w)
{
  return dot(z,w);
}

// Return the 2D scalar cross product z.x*w.y-z.y*w.x.
real cross(explicit pair z, explicit pair w)
{
  return cross(z,w);
}

pair bezier(pair a, pair b, pair c, pair d, real t)
{
  real onemt=1-t;
  real onemt2=onemt*onemt;
  return onemt2*onemt*a+t*(3.0*(onemt2*b+t*onemt*c)+t*t*d);
}

pair bezierP(pair a, pair b, pair c, pair d, real t)
{
  return 3.0*(t*t*(d-a+3.0*(b-c))+t*(2.0*(a+c)-4.0*b)+b-a);
}

pair bezierPP(pair a, pair b, pair c, pair d, real t)
{
  return 6.0*(t*(d-a+3.0*(b-c))+a+c)-12.0*b;
}

pair bezierPPP(pair a, pair b, pair c, pair d)
{
  return 6.0*(d-a)+18.0*(b-c);
}