poly.cpp

#include "poly.h"

Poly Poly::operator*(const Poly& p) const {
    Poly result; 
    result.resize(degree() +  p.degree()+1);
    
    for(unsigned i = 0; i < size(); i++) {
        for(unsigned j = 0; j < p.size(); j++) {
            result[i+j] += (*this)[i] * p[j];
        }
    }
    return result;
}
#ifdef HAVE_GSL
#include <gsl/gsl_poly.h>
#endif

/*double Poly::eval(double x) const {
    return gsl_poly_eval(&coeff[0], size(), x);
    }*/

void Poly::normalize() {
    while(back() == 0)
        pop_back();
}

void Poly::monicify() {
    normalize();
    
    double scale = 1./back(); // unitize
    
    for(unsigned i = 0; i < size(); i++) {
        (*this)[i] *= scale;
    }
}


#ifdef HAVE_GSL
std::vector<std::complex<double> > solve(Poly const & pp) {
    Poly p(pp);
    p.normalize();
    gsl_poly_complex_workspace * w 
        = gsl_poly_complex_workspace_alloc (p.size());
       
    gsl_complex_packed_ptr z = new double[p.degree()*2];
    double* a = new double[p.size()];
    for(int i = 0; i < p.size(); i++)
        a[i] = p[i];
    std::vector<std::complex<double> > roots;
    //roots.resize(p.degree());
    
    gsl_poly_complex_solve (a, p.size(), w, z);
    delete[]a;
     
    gsl_poly_complex_workspace_free (w);
     
    for (int i = 0; i < p.degree(); i++) {
        roots.push_back(std::complex<double> (z[2*i] ,z[2*i+1]));
        //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
    }    
    delete[] z;
    return roots;
}

std::vector<double > solve_reals(Poly const & p) {
    std::vector<std::complex<double> > roots = solve(p);
    std::vector<double> real_roots;
    
    for(int i = 0; i < roots.size(); i++) {
        if(roots[i].imag() == 0) // should be more lenient perhaps
            real_roots.push_back(roots[i].real());
    }
    return real_roots;
}
#endif

double polish_root(Poly const & p, double guess, double tol) {
    Poly dp = derivative(p);
    
    double fn = p(guess);
    while(fabs(fn) > tol) {
        guess -= fn/dp(guess);
        fn = p(guess);
    }
    return guess;
}

Poly integral(Poly const & p) {
    Poly result;
    
    result.reserve(p.size()+1);
    result.push_back(0); // arbitrary const
    for(unsigned i = 0; i < p.size(); i++) {
        result.push_back(p[i]/(i+1));
    }
    return result;

}

Poly derivative(Poly const & p) {
    Poly result;
    
    if(p.size() <= 1)
        return Poly(0);
    result.reserve(p.size()-1);
    for(unsigned i = 1; i < p.size(); i++) {
        result.push_back(i*p[i]);
    }
    return result;
}

Poly compose(Poly const & a, Poly const & b) {
    Poly result;
    
    for(unsigned i = a.size(); i > 0; i--) {
        result = Poly(a[i-1]) + result * b;
    }
    return result;
    
}

/* This version is backwards - dividing taylor terms
Poly divide(Poly const &a, Poly const &b, Poly &r) {
    Poly c;
    r = a; // remainder
    
    const unsigned k = a.size();
    r.resize(k, 0);
    c.resize(k, 0);

    for(unsigned i = 0; i < k; i++) {
        double ci = r[i]/b[0];
        c[i] += ci;
        Poly bb = ci*b;
        std::cout << ci <<"*" << b << ", r= " << r << std::endl;
        r -= bb.shifted(i);
    }
    
    return c;
}
*/

Poly divide(Poly const &a, Poly const &b, Poly &r) {
    Poly c;
    r = a; // remainder
    assert(b.size() > 0);
    
    const unsigned k = a.degree();
    const unsigned l = b.degree();
    c.resize(k, 0.);
    
    for(unsigned i = k; i >= l; i--) {
        assert(i >= 0);
        double ci = r.back()/b.back();
        c[i-l] += ci;
        Poly bb = ci*b;
        //std::cout << ci <<"*(" << b.shifted(i-l) << ") = " 
        //          << bb.shifted(i-l) << "     r= " << r << std::endl;
        r -= bb.shifted(i-l);
        r.pop_back();
    }
    //std::cout << "r= " << r << std::endl;
    r.normalize();
    c.normalize();
    
    return c;
}

Poly gcd(Poly const &a, Poly const &b, const double tol) {
    if(a.size() < b.size())
        return gcd(b, a);
    if(b.size() <= 0)
        return a;
    if(b.size() == 1)
        return a;
    Poly r;
    divide(a, b, r);
    return gcd(b, r);
}



/*Poly divide_out_root(Poly const & p, double x) {
    assert(1);
    }*/


/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :