/*
 * Arithmetic on polynomials with (integer) coefficients
 *
 * Copyright © 2002  Bart Massey.
 * All Rights Reserved.  See the file COPYING in this directory
 * for licensing information.
 */

namespace Polynomial {

 /*
  * Change the typedef to rational for rational
  * coefficients, then walk carefully through code below.
  */
 typedef int coefficient;

 /*
  * Polynomials are stored as an array of
  * coefficients indexed by degree.  There
  * is always a constant term, and the leading
  * term of a polynomial is always nonzero (except
  * for the zero polynomial).
  */
 public typedef coefficient[*] polynomial;


 /*
  * Bug: There's currently no Nickle way to normalize the
  * polynomial by shortening the array.  This should be
  * fixed.  In the meantime, we will just make a copy.
  */
 polynomial normalize(polynomial p) {
     int n = dim(p) - 1;
     if (p[n] != 0)
	 return p;
     while (n > 0 && p[n] == 0)
	 --n;
     coefficient[n + 1] q = { [i] = p[i] };
     return q;
 }


 /*
  * Addition is performed using an array comprehension
  * block.  I'm not sure whether this is better or worse
  * than just doing it with a loop.
  */
 public polynomial add(polynomial a, polynomial b) {
     coefficient na = dim(a);
     coefficient nb = dim(b);
     coefficient nc = max(na, nb);
     coefficient[nc] c = { [i] {
	 coefficient get_coeff(coefficient[*] a) {
	     if (i < dim(a))
		 return a[i];
	     return 0;
	 }
	 coefficient ai = get_coeff(a);
         coefficient bi = get_coeff(b);
	 return ai + bi;
     } };
     return normalize(c);
 }


 /*
  * Multiplication efficiency might be minutely improved by
  * avoiding the all-zeros initialization.  Or it might
  * not...
  */
 public polynomial mult(polynomial a, polynomial b) {
     coefficient na = dim(a);
     coefficient nb = dim(b);
     coefficient nc = na + nb - 1;
     coefficient[nc] c = {0 ...};

     for (coefficient i = 0; i < na; i++)
	 for (coefficient j = 0; j < nb; j++)
	     c[i + j] += a[i] * b[j];
     return c;
 }


 /*
  * The quotient and remainder for polynomial division are
  * computed simultaneously, and thus should be returned
  * together.
  */
 public typedef struct {
     polynomial quot;
     polynomial rem;
 } quot_rem;


 /*
  * Return quotient and remainder of a div b.  This code
  * assumes that with integer coefficients, leading
  * coefficient of b must be 1.
  */
 public quot_rem div_mod(polynomial a, polynomial b) {
     coefficient na = dim(a);
     coefficient nb = dim(b);
     if (b[nb - 1] != 1)
	 abort("leading integer coefficient in divisor should be 1");
     coefficient nq = na - nb + 1;
     if (nq <= 0)
	 return (quot_rem) { .quot = (polynomial){0}, .rem = b };
     coefficient[nq] q = {0 ...};

     while (na >= nb) {
         coefficient xn = na - nb;
	 coefficient cx = a[na - 1] / b[nb - 1];
	 for (int i = 0; i < nb; i++)
	     a[i + xn] -= cx * b[i];
	 q[xn] += cx;
	 if (a[na - 1] != 0)
	     abort("leading coefficient error");
	 --na;
	 while (na >= 1 && a[na - 1] == 0)
	     --na;
     }
     return (quot_rem) { .quot = q, .rem = normalize(a) };
 }


 /*
  * Returns b ** e computed using the usual fast algorithm
  * (O(log e) multiplies).
  */
 public polynomial
 power(polynomial b, int e) {
     if (e < 0)
	 abort("negative exponent error");
     if (e == 0)
	 return (polynomial) {1};
     if (e == 1)
	 return b;
     polynomial tmp = power(b, e // 2);
     polynomial tmp2 = mult(tmp, tmp);
     if (e % 2 == 0)
	 return tmp2;
     return mult(tmp2, b);
 }


 /*
  * Returns b ** e (mod m) computed using the usual
  * fast algorithm (O(log e) multiplies, each requiring
  * O(|m|*|b**2|) primitive steps).
  */
 public polynomial
 power_mod(polynomial b, int e, polynomial m) {
     if (e < 0)
	 abort("negative exponent error");
     if (e == 0)
	 return (polynomial) {1};
     if (e == 1)
	 return div_mod(b, m).rem;
     polynomial tmp = power_mod(b, e // 2, m);
     polynomial tmp2 = mult(tmp, tmp);
     if (e % 2 == 0)
	 return div_mod(tmp2, m).rem;
     return div_mod(mult(tmp2, b), m).rem;
 }

 /*
  * Returns a string representation of p suitable for Maple
  * input.  x is the name of the polynomial basis variable
  * (usually just "x").  There's a lot of boundary cases
  * here.
  */
 public string maple(polynomial p, string x) {
     string s = "";
     bool started = false;
     for (int i = 0; i < dim(p); i++) {
	 /* Don't print absent terms. */
	 if (p[i] == 0)
	     continue;
	 /* If a term has already been printed, this one
	    must be added or subtracted */
	 if (started) {
	     if (p[i] > 0) {
		 s += " + ";
	     } else {
		 s += " - ";
		 p[i] = -p[i];
	     }
	 }
	 /* Don't print 1*x or 1*x^n. */
	 string star = "*";
	 if (p[i] == 1 && i > 0)
	     star = "";
	 else
	     s += sprintf("%d", p[i]);
	 /* Don't print x^0 or x^1.  If there's no leading
	    coefficient due to previous logic, don't print
	    the multiply symbol. */
	 if (i == 1)
	     s += sprintf("%s%s", star, x);
	 else if (i > 1)
	     s += sprintf("%s%s^%d", star, x, i);
	 /* A term has now been printed */
	 started = true;
     }
     return s;
 }

}