//
// Copyright (c) 2005 Mooffie <mooffie@typo.co.il>
//
// This script may not be installed on any server
// other than that contracted by Mooffie.

SOLVE_QUADRATIC = true;
SOLVE_IMAG      = true;

function Polynomial(L)
{
    this.is_whole = Polynomial_is_whole;
    this.toString = Polynomial_toString;
    this.set = Polynomial_set;
    this.evaluate = Polynomial_evaluate;
    this.is_root = Polynomial_is_root;
    this.find_root = Polynomial_find_root;
    this.first = Polynomial_first;
    this.last = Polynomial_last;
    this.factor = Polynomial_factor;
    this.get_sorted_factors = Polynomial_get_sorted_factors;
    this.reduce_by_root = Polynomial_reduce_by_root;
    this.copy = Polynomial_copy;
    this.reduce_by_num = Polynomial_reduce_by_num;
    this.find_common_factor = Polynomial_find_common_factor;
    this.is_solvable_quadratic = Polynomial_is_solvable_quadratic;
    this.factor_quadratic = Polynomial_factor_quadratic;
    this.is_quadratic = Polynomial_is_quadratic;
    this.allocate = Polynomial_allocate;
    this.equals = Polynomial_equals;
    this.is_constant = Polynomial_is_constant;
    this.factor_constant = Polynomial_factor_constant;

    this.p = array_dup(L);
}

function Polynomial_is_constant() {
    return this.p.length == 1;
}

function Polynomial_equals(other) {
    if (this.p.length == other.p.length) {
	for (var i = 0; i < this.p.length; i++) {
	    if (this.p[i] != other.p[i])
		return false;
	}
	return true;
    } else {
	return false;
    }
}

function Polynomial_allocate(len) {
    this.p = [];
    while (len-- > 0)
	this.p[this.p.length] = 0;
}

function Polynomial_set(L) {
    this.p = array_dup(L);
}

function Polynomial_is_whole() {
    for (var i = 0; i < this.p.length; i++)
	if (!is_whole(this.p[i]))
	    return false;
    return true;
}

function Polynomial_evaluate(x) {
    var result = 0;
    for (var i = 0; i < this.p.length; i++) {
	var coef = this.p[i];		    // coefficient
	var deg  = this.p.length - i - 1;   // degree
	result += coef * Math.pow(x, deg);
    }
    return result;
}

function Polynomial_is_root(x) {
    return (this.p.length > 1) && this.evaluate(x) == 0;
}

function Polynomial_first() {
    return this.p[0];
}

function Polynomial_last() {
    return this.p[this.p.length - 1];
}

function Polynomial_copy() {
    return new Polynomial(this.p);
}

// find_common_factor() finds the common factor of all the coefficients.
function Polynomial_find_common_factor() {
    if (this.is_whole() && this.p.length > 1) {
	return gcd_mult(this.p);
    } else {
	return 1;
    }
}

// reduce_by_num(n) - divides the polynomial by a constant
function Polynomial_reduce_by_num(n) {
    for (var i = 0; i < this.p.length; i++) {
	this.p[i] = get_whole(this.p[i]) / n;
    }
}

// reduce_by_root(r) - divides the polynomial by (x-r), using
// synthetic division.
function Polynomial_reduce_by_root(r) {
    if (ISRAT(r))
	r = r.float_val();
    var new_p = [];
    var carry = 0;
    var down;
    for (var i = 0; i < this.p.length; i++) {
	var coef = this.p[i];		    // coefficient
	var deg  = this.p.length - i - 1;   // degree
	down = coef + carry;
	carry = down * r;
	array_append(new_p, down);
    }
    if (down == 0) {
	array_pop(new_p);
	this.set(new_p);
	return true;
    } else {
	return false;
    }
}

// solve_quadratic(a, b, c) is a helper function: it gets a, b and c of
// a quadratic and solves it using the quadratic formula.
function solve_quadratic(a, b, c) {
    a = FVAL(a);
    b = FVAL(b);
    c = FVAL(c);
    
    var r = new Irrational();
    r.a = -b;
    r.b = 1;
    r.c = b*b - 4*a*c;
    r.d = 2*a;
    r.normalize();

    var solution = [];
    solution[0] = solution.irrat = r;
    solution[1] = solution.x1 = (r.a + r.b*Math.sqrt(r.c)) / r.d;
    solution[2] = solution.x2 = (r.a - r.b*Math.sqrt(r.c)) / r.d;
    return solution;
}

function Polynomial_is_quadratic() {
    return this.p.length == 3;
}

function Polynomial_is_solvable_quadratic() {
    if (!SOLVE_QUADRATIC || !this.is_quadratic())
	return false;
    var sol = solve_quadratic(this.p[0], this.p[1], this.p[2]);
    return !sol.irrat.imag || SOLVE_IMAG;
}

// factor_quadratic() factors this polynomial using the quadric formula.
function Polynomial_factor_quadratic() {
    var factors = [];
    var sol = solve_quadratic(this.p[0], this.p[1], this.p[2]);
    if (sol.irrat.c == 0) {
	// This is rational, after all.
	var rat = new Rational(sol.irrat.a, sol.irrat.d);
	array_append(factors, new Polynomial([1, get_neg(rat)]));
	array_append(factors, new Polynomial([1, get_neg(rat)]));
    } else {
	array_append(factors, new Polynomial([1, sol.irrat]));
    }
    if (this.p[0] != 1)
	array_append(factors, new Polynomial([this.p[0]]));
    return factors;
}

// TODO: optimize!
function Polynomial_factor_constant() {
    if (!is_whole(this.p[0]) || this.p[0] < 2)
	return [ new Polynomial([this.p[0]]) ];
    var x = Math.abs(get_whole(this.p[0]));
    var primes = [];
    var i = 2;
    while (i <= x) {
	if (x % i == 0) {
	    array_append(primes, new Polynomial([i]));
	    x /= i;
	} else {
	    i++;
	}
    }
    return primes;
}

// _cmp_factors(a, b) is used for sorting factors.
function _cmp_factors(a, b) {
    if (a.p.length != b.p.length)
	return a.p.length - b.p.length;
    
    if (a.p.length == 1) {
	return (FVAL(a.p[0]) - FVAL(b.p[0]));
    }
    
    if (a.p.length == 2) {
	// Move irrationals to back.
	if (ISIRR(a.p[1])) {
	    return 1;
	}
	if (ISIRR(b.p[1])) {
	    return -1;
	}
	// Move 'x' to front.
	if (FVAL(a.p[1]) == 0) {
	    return -1;
	}
	if (FVAL(b.p[1]) == 0) {
	    return 1;
	}
	// Else, move (x + 2) in front of (x - 1)
	return (FVAL(b.p[1]) - FVAL(a.p[1]));
    }

    // Else, I have nothing interesting to do...
    return (FVAL(a.p[0]) - FVAL(b.p[0]));
}

// get_sorted_factors() is just a wrapper around factor() that sorts
// the factors.
function Polynomial_get_sorted_factors() {
    var fac = this.factor();
    if (this.is_constant())
	return this.factor_constant();
    fac.sort(_cmp_factors);

    // replace the multiplier.
    var mult = 1;
    var mult_count = 0;
    for (var i = 0; (i < fac.length) && (fac[i].p.length == 1); i++) {
	mult *= fac[i].p[0];
	mult_count++;
    }
    dbg(mult_count + " mults total: " + mult);
    // remove the multipliers:
    fac = fac.slice(mult_count);
    // insert the multiplier, if it's not 1:
    if ((mult != 1) || (fac.length == 0))
	fac = [new Polynomial([mult])].concat(fac);

    return fac;
}

function _odd(x) { return x % 2; }
function _even(x) { return !_odd(x); }
function _is_square_wo_sign(x) { return is_whole(Math.sqrt(Math.abs(FVAL(x)))); }
function _sqrt_wo_sign(x) {
    return get_whole(Math.sqrt(Math.abs(FVAL(x)))) * (FVAL(x) < 0 ? -1 : 1);
}

// factor() is the central method in the class.
function Polynomial_factor()
{
    var factors = [];

    var terms_idx = [];
    for (var i = 0; i < this.p.length; i++)
	if (this.p[i] != 0)
	    array_append(terms_idx, i);
    dbg(terms_idx)

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

    // Is this a "(a^n*x^n)^2 - n^2" pattern?

    if ((terms_idx.length == 2) &&                          // 1. two terms.
	(_odd(this.p.length) && _even(terms_idx[1])) &&     // 2. even degrees.
	(FVAL(this.p[0]) *
		FVAL(this.p[terms_idx[1]]) < 0) &&	    // 3. different signs.
	(_is_square_wo_sign(this.p[0]) &&		    
		_is_square_wo_sign(this.p[terms_idx[1]])))  // 4. squares of naturals.
    {
	dbg("Pattern found: a^2 - b^2");
	var p = new Polynomial([]);
	p.allocate(_int(this.p.length/2) + 1);
	p.p[0] = _sqrt_wo_sign(this.p[0]);
	p.p[_int(terms_idx[1]/2)] = _sqrt_wo_sign(this.p[terms_idx[1]]);
	dbg("fac1: " + p);
	array_extend(factors, p.factor());

	if (p.p[0] < 0)
	    p.p[0] *= -1;
	else
	    p.p[_int(terms_idx[1]/2)] *= -1;
	dbg("fac2: " + p);
	array_extend(factors, p.factor());
	return factors;
    }
    
    //////////////////////////////////////////////////////////////////////

    // Is this a "a*(x^n)^2 + b*x^n + c" pattern?

    if ((terms_idx.length == 3) &&                          // 1. three terms.
	(_odd(this.p.length) &&				    // 2. x^4n - x^2n
	    (terms_idx[1] == _int(this.p.length/2))) &&
	(terms_idx[2] == this.p.length-1))		    // 3. + c
    {
	dbg("Pattern found: bi-quad");
	var sol = solve_quadratic(this.p[0], this.p[terms_idx[1]], this.p[terms_idx[2]]);
	if (!sol.irrat.imag && is_whole(sol.x1) && is_whole(sol.x2)) {
	    var p = new Polynomial([]);
	    p.allocate(_int(this.p.length/2) + 1);
	    p.p[0] = 1
	    p.p[p.p.length-1] = -sol.x1;
	    dbg("fac1: " + p);
	    array_extend(factors, p.factor());
	    p.p[p.p.length-1] = -sol.x2;
	    dbg("fac2: " + p);
	    array_extend(factors, p.factor());
	    if (this.p[0] != 1)
		array_append(factors, new Polynomial([this.p[0]]));
	    return factors;
	}
    }

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

    // Now, look for roots using find_root(), which uses the rational roots test.
    var root = this.find_root();
    if (ISRAT(root)) {
	dbg("Found root: " + root);
	array_append(factors, new Polynomial([1, get_neg(root)]));
	var p = this.copy();
	// divide by (x-x0):
	p.reduce_by_root(root);
	// and factor the result recursively:
	array_extend(factors, p.factor());
    } else {
	// Didn't find roots.
	var p = this.copy();
	// Perhaps it can be factored using the quadratic formula?
	if (p.is_solvable_quadratic()) {
	    array_extend(factors, p.factor_quadratic());
	} else {
	    // No it can't. Now just try to extract a common factor.
	    var com = p.find_common_factor();
	    if (com != 1) {
		p.reduce_by_num(com);
		array_append(factors, new Polynomial([com]));
	    }
	    array_append(factors, p);
	}
    }

    return factors;
}

// find_root() searches for a rational root using
// the "rational roots test".
function Polynomial_find_root() {
    // As a special case, check for x=0.
    if (this.is_root(0))
	return new Rational(0, 1);
    
    if (!is_whole(this.first()) || !is_whole(this.last())) {
	dbg("!Not whole numbers");
	return false;
    }

    first = Math.abs(get_whole(this.first()));
    last  = Math.abs(get_whole(this.last()));
    for (var last_f = 1; last_f <= last; last_f++) {
	if (last % last_f != 0)
	    continue; // this is not a factor
	for (var first_f = 1; first_f <= first; first_f++) {
	    if (first % first_f != 0)
		continue; // this is not a factor
	    var root = last_f / first_f;
	    //dbg("trying " + last_f + "/" + first_f);
	    if (this.is_root(root))
		return new Rational(last_f, first_f);
	    if (this.is_root(-root))
		return new Rational(-last_f, first_f);
	}
    }
}

function Polynomial_toString() {
    var s = ""
    for (var i = 0; i < this.p.length; i++) {
	var coef = this.p[i];		    // coefficient
	var deg  = this.p.length - i - 1;   // degree
	if ((FVAL(coef) == 0) && (this.p.length != 1))
	    continue;
	if (s) {
	    // output a plus or minus sign in front of the term.
	    if (ISIRR(coef)) {
		if (coef.plus_minus_outside) {
		    s += " +/- ";
		} else {
		    // since this polynomial means (x - x0),
		    // we should use a minus sign.
		    s += " - ";
		}
	    } else if (FVAL(coef) > 0) {
		s += " + ";
	    } else {
		s += " - ";
		coef = get_neg(coef);
	    }
	}
	if (FVAL(coef) == -1 && !s) {
	    // Here we deal with a special case: the leading
	    // coefficient is -1. We want to print "-x4", not "-1x4".
	    s = "-";
	    coef = 1;
	}
	if (FVAL(coef) != 1 || deg == 0) {
	    s += coef;
	}
	if (deg > 0) {
	    s += "x";
	    if (deg > 1)
		s += "^" + deg;
	}
    }
    return s;
}

// The following utility functions deal with arrays.
// While it's true that JavaScript provides many of these,
// not all browsers support all of them, so I prefer to
// implement them myself, 

// array_dup() - duplicates an array.
function array_dup(a) {
    var na = [];
    for (var i = 0; i < a.length; i++)
	na[i] = a[i];
    return na;
}

// We implement array_append() and array_pop() because IE 5.0 doesn't
// support Array.push() and Array.pop().
function array_append(a, v) {
    a[a.length] = v;
}

function array_pop(a) {
    if (a.length > 0) {
	var v = a[a.length - 1];
	a.length--;
	return v;
    }
}

function array_extend(a, L) { // a-la Python
    for (var i = 0; i < L.length; i++)
	array_append(a, L[i]);
}