/*
 * 3D Cubic Bezier Curve (rational coordinates)
 * 15-745:L3, Todd Phillips
 * Some basic 3d bezier curve predicates for rational coordinate points
 * Should benefit greatly from propagating the "resolution" constant deep
 */

// Basic 2D Point


/////////////////////////////////////////////////////////////////////////////////
// very basic rational class 
// doesn't do equality checking
// doesn't bother checking for num/denominator overflow (so watch out)
struct Rational {
  int num;
  int denom;
};

// Round to the nearest
int rat2int(struct Rational *r)
{
  int q,remain;
  q = r->num / r->denom;
  remain = r->num % r->denom;
  if (remain >= ((r->denom >> 1) + (r->denom % 2)))
    q = q + 1;
  return q;
}

int gcd(int a, int b)
{
  if(b == 0)
    return a;
  else
    return gcd(b, a % b);
}

void simplify(struct Rational *a)
{
  int g;
  g = gcd (a->num, a -> denom);
  a->num = a->num / g;
  a->denom = a->denom / g;
  return;
}

struct Rational *int2rat(int i)
{
  static struct Rational r;
  r.num = i;
  r.denom = 1;
  return &r;
}

struct Rational *r_mul(struct Rational *a, struct Rational *b)
{
  static struct Rational r;
  r.num = a->num * b->num;
  r.denom = a->denom * b->denom;
  simplify(&r);
  return &r;
}

struct Rational *r_scale(struct Rational *a, int i)
{
  static struct Rational r;
  r.num = a->num * i;
  r.denom = a->denom;
  simplify(&r);
  return &r;
}

// Really just a quick constructor
struct Rational *r_quot(int i, int j)
{
  static struct Rational r;
  r.num = i;
  r.denom = j;
  simplify(&r);
  return &r;
}

struct Rational *r_sum(struct Rational *a, struct Rational *b)
{
  static struct Rational r;
  int g, aog, bog;
  simplify(a);
  simplify(b);
  g = gcd (a->denom, b->denom);
  aog = a->denom / g;
  bog = b->denom / g;
  r.num = a->num * bog  + b->num * aog;
  r.denom = g * aog * bog;
  simplify(&r);
  return &r;
}

////////////////////////////////////////////////////////////////////

struct Rational *r_neg(struct Rational *a)
{
  static struct Rational r;
  r.num = - a->num;
  r.denom = a->denom;
  simplify(&r);
  return &r;
}

void print_rational(struct Rational *a)
{
  // nearest int
  print_int(rat2int(a));
}

struct Point3D {
  struct Rational *x;
  struct Rational *y;
  struct Rational *z;
};

// negation
struct Point3D *minus(struct Point3D a)
{
  static struct Point3D  b;
  b.x = r_neg(a.x);
  b.y = r_neg(a.y);
  b.z = r_neg(a.z);
  return &b;
}


// add two vectors                                               
struct Point3D *sum(struct Point3D a, struct Point3D b)
{
  static struct Point3D c;
  c.x = r_sum(a.x,b.x);
  c.y = r_sum(a.y,b.y); 
  c.z = r_sum(a.z,b.z); 
  return &c;                                                  
}                                                       

struct Point3D *dilate(struct Point3D a, struct Rational *scale)
{
  static struct Point3D b;
  b.x = r_mul(a.x,scale);
  b.y = r_mul(a.y,scale);
  b.z = r_mul(a.z,scale);
  return &b;
}


// a and b should be in scale multiplied coordinates
// a_weight/resolution should be the contribution of a
// b_weight/resolution ( = 1- a_weight/resolution) is the contribution of b 
struct Point3D *affine_combine(struct Point3D a, struct Point3D b, int a_weight, int resolution)
{
  int b_weight;
  struct Rational *a_frac, *b_frac;
  struct Point3D *c;

  b_weight = resolution - a_weight;

  a_frac = r_quot(b_weight, resolution);
  b_frac = r_quot(a_weight, resolution);

  c = sum( *(dilate(a,a_frac)) , *(dilate(b,b_frac)));

  return c;
}

struct CubicCurve {
  struct Point3D a;
  struct Point3D b;
  struct Point3D c;
  struct Point3D c;
};

void print_point(struct Point3D a)
{
  print_rational(a.x);
  print_space(1);
  print_rational(a.y);
  print_space(1);
  print_rational(a.z);
  print_space(1);
}

// Performs deCasteljou's algo on a cubic bezier curve
struct Point3D *curve_point (struct CubicCurve curve, int a_weight, int resolution)
{
  struct Point3D *qa,*qb,*qc;   // The first (quadratic) polar
  struct Point3D *la,*lb;      // The second (linear) polar
  struct Point3D *r;          // The third (point) polar, i.e. the answer

  qa = affine_combine(curve.a,curve.b,a_weight,resolution);
  qb = affine_combine(curve.b,curve.c,a_weight,resolution);
  qc = affine_combine(curve.c,curve.d,a_weight,resolution);
  
  la = affine_combine ( *qa, *qb, a_weight, resolution);
  lb = affine_combine ( *qb, *qc, a_weight, resolution);

  r = affine_combine ( *la, *lb, a_weight, resolution);

  return r;
}
     
void main(void) {
  int i,j,w;
  Point3D *p;
  struct CubicCurve curve;
  int resolution;
  int steps;   // This just makes it do a lot of work
  struct Rational *t;

  steps = 35;

  // Hopefully this constant makes it deep into the calculation
  // Make sure resolution**4 << 2**32 otherwise we have overflow (not a bother, just not right answers)
  resolution = 100;

  // Initialize a basic space curve
  // hopefully these constants make it deep
  curve.a.x = int2rat(0);
  curve.a.y = int2rat(0);
  curve.a.z = int2rat(0);

  curve.b.x = int2rat(0);
  curve.b.y = int2rat(0);
  curve.b.z = int2rat(12);

  curve.c.x = int2rat(0);
  curve.c.y = int2rat(13);
  curve.c.z = int2rat(14);

  curve.d.z = int2rat(16);

  for (i = 0; i < steps; i+= 1)
    {
      curve.d.x = int2rat(17 + i);
      for(j = 0; j < steps; j += 1)
  	{
	  curve.d.y = int2rat(15 - i );

	  for (w = 0; w <= resolution; w += 1)
	    {
	      p = curve_point(curve, w, resolution);
	      print_point((*p));
	      print_newline();
	    }
	}
    }
}