/*
 * Newton's Method
 * 15-745:L3, Todd Phillips
 * Find the zero (integer near a zero) of an  N th degree polynomial
 * Only use rationals
 * CSE should be vital
 */

void print_long(int i) {}
void print_newline(int i) {}
void *alloc(int count, int size) {}

/////////////////////////////////////////////////////////////////////////////////
// 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;
  if (r->denom == 0)
    r->denom = 1;
  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;
  int max_d;

  // Throw out some precision
  // to prevent overflows
  max_d = 1<<13;
  while (a->denom > max_d || a->denom < (-max_d))
    {
      a->num = a->num >> 1;
      a->denom = a->denom >> 1;
    }

  g = gcd (a->num, a -> denom);
  if (g== 0)
    g = 1;
  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_neg(struct Rational *a)
{
  static struct Rational r;
  r.num = -a->num;
  r.denom = a->denom;
  simplify(&r);
  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_divide(struct Rational *a, struct Rational *b)
{
  static struct Rational r;
  r.num = a->num * b->denom;
  r.denom = a->denom * b->num;
  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);
  if (g==0)
    g = 1;
  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;
}

int is_zero(struct Rational *r)
{
     return (r->num == 0);
}

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

// Coefficients of a polynomial form an integer list
// Low Degree coefficients first
// Poly's stored in newton form    P(x) = c + x P2(x)      (base case is 0)
struct Poly {
  int c;
  struct Poly *next;
};

struct Rational *eval(struct Rational *x, struct Poly *p)
{
  if (p == NULL)
    return int2rat(0);
  return r_sum(int2rat(p->c), r_mul(x,eval(x,p->next)));
}

// chain rule for newton form
struct Rational *deriv(struct Rational *x, struct Poly *p)
{
  if (p == NULL)
    return int2rat(0);
  return r_sum(eval(x,p->next),r_mul(x,deriv(x,p->next)));
}

// for convenience
struct Poly *make_poly(int *coff)
{
  struct Poly *p;
  int d;

  p = alloc (size(coff),Poly);

  for (d=0;d < size(coff); d +=1)
    {
      p[d].c = coff[d];
      if ( d < size(coff) -1)
	p[d].next = &(p[d+1]);
      else
	p[d].next = NULL;
    }

  return p;
}

void print_long(struct Rational *r)
{
  print_int(r->num); 
  print_space(5); 
  print_int(r->denom);
  print_space(15);
  print_int(rat2int(r));
  print_newline();
}

struct Rational *newton(struct Poly *p, struct Rational *g, int its_left)
{
  struct Rational *newg;

  newg = r_sum(g, r_neg( r_divide(eval(g,p),deriv(g,p))));

  print_long(g);
  if (is_zero(eval(newg,p)) || its_left <= 0)
    {
      if (its_left <= 0)
	print_int(1010101);
      else
	print_int(its_left);
      print_newline();
      return newg;
    }
  else
    return newton(p,newg, its_left - 1);
}

void main(void) {
  int *coff;
  struct Poly *cub;
  struct Rational *x, *y;

  coff = alloc(12,4);
  coff[0]= 30; 
  coff[1] = 10;
  coff[2] = 412;
  coff[3] = 1;
  coff[4] = 3;
  coff[5] = 21;
  coff[6] = 10;
  coff[7] = 3;
  coff[8] = 2;
  coff[9] = 7;
  coff[10] = 10;
  coff[11] = 10;

  cub = make_poly(coff);

  x = int2rat(21);
  y = newton(cub,x,20000);

  print_long(y);
  print_long(eval(y,cub));
  print_newline();
}