poly.cpp revision 981b809bc6ed10a21e89444d9447e5475801874f
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm#include "poly.h"
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmPoly Poly::operator*(const Poly& p) const {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm Poly result;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm result.resize(degree() + p.degree()+1);
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned i = 0; i < size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned j = 0; j < p.size(); j++) {
ddc251b3cf95b0097b6a5ee39ea132bd4d7d5cbcjohanengelen result[i+j] += (*this)[i] * p[j];
ddc251b3cf95b0097b6a5ee39ea132bd4d7d5cbcjohanengelen }
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm }
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm return result;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm}
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm#ifdef HAVE_GSL
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm#include <gsl/gsl_poly.h>
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm#endif
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
6656f193fdace606d1b162d6dea0223bc295f0a6cilix/*double Poly::eval(double x) const {
6656f193fdace606d1b162d6dea0223bc295f0a6cilix return gsl_poly_eval(&coeff[0], size(), x);
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm }*/
77a4a003111bd5cfb771d4849801c898aeb889b0cilix
77a4a003111bd5cfb771d4849801c898aeb889b0cilixvoid Poly::normalize() {
04c99c338ffdc6e10cb6f5c18f6f06b3f555e8ebcilix while(back() == 0)
04c99c338ffdc6e10cb6f5c18f6f06b3f555e8ebcilix pop_back();
04c99c338ffdc6e10cb6f5c18f6f06b3f555e8ebcilix}
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmvoid Poly::monicify() {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm normalize();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm double scale = 1./back(); // unitize
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned i = 0; i < size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm (*this)[i] *= scale;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm }
ddc251b3cf95b0097b6a5ee39ea132bd4d7d5cbcjohanengelen}
ddc251b3cf95b0097b6a5ee39ea132bd4d7d5cbcjohanengelen
2e2b8bd323e3693f9d86f545ce049d3f1b45d1c2cilix
ddc251b3cf95b0097b6a5ee39ea132bd4d7d5cbcjohanengelen#ifdef HAVE_GSL
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmstd::vector<std::complex<double> > solve(Poly const & pp) {
8c39cbeab9949a0a7d6ae66b768a7352019e42f8johanengelen Poly p(pp);
072916d0ef7dccd696b59381f50bcf776abccefbjohanengelen p.normalize();
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud gsl_poly_complex_workspace * w
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud = gsl_poly_complex_workspace_alloc (p.size());
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm gsl_complex_packed_ptr z = new double[p.degree()*2];
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm double* a = new double[p.size()];
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(int i = 0; i < p.size(); i++)
8d9f5d586a04809427ce1df284a5720112177991cilix a[i] = p[i];
70eb1fc448cb08acf3468f80fa2296c03b32afd2cilix std::vector<std::complex<double> > roots;
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyak //roots.resize(p.degree());
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen
0cc5b8d2f7b87c4222ee3662071bef1cb1f22b06bgk gsl_poly_complex_solve (a, p.size(), w, z);
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen delete[]a;
f4db63be4e929f4706410914295deccaceea19cdcilix
ab99111a42436818e6902e044c8f3af2b724263bcilix gsl_poly_complex_workspace_free (w);
76db360f5f052775326e6d406b9e1e9e2966e11acilix
3d0482af18ffb591c1d8ddecf516629e1bcd2ae4cilix for (int i = 0; i < p.degree(); i++) {
64aee804a6a47424f7994e60558351b8cf2ea4dbcilix roots.push_back(std::complex<double> (z[2*i] ,z[2*i+1]));
b320a8d186114a5122ddc3afbe95110eb6cb10cecilix //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
044d712d4d03f8354962d54e47cfac2346a69ccccilix }
61cfd957cd023c4f432ea0c7307784a56bf978e9cilix delete[] z;
2f5c0701b333a695eedb1680beb1adf95c0723dacilix return roots;
add2ffae3c4686b50d888775bbdf083a4726a210johanengelen}
0a75b58e47d3de42550c4f7960e253ea3befc09ajohanengelen
4afe3fc6b9c122bc5c02b27aea3845ba41384d2acilixstd::vector<double > solve_reals(Poly const & p) {
ecf048161ae0284d01d34ca5844204775c0d3fdecilix std::vector<std::complex<double> > roots = solve(p);
3515994554d167522343ce57417648b39370ccabcilix std::vector<double> real_roots;
e54ce05030e6aab675331e18f46f029f55ed1bf0cilix
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(int i = 0; i < roots.size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm if(roots[i].imag() == 0) // should be more lenient perhaps
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm real_roots.push_back(roots[i].real());
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm }
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm return real_roots;
72b7b31db250f20b90730d2888e6a554b434a407johanengelen}
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm#endif
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelendouble polish_root(Poly const & p, double guess, double tol) {
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen Poly dp = derivative(p);
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen double fn = p(guess);
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen while(fabs(fn) > tol) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm guess -= fn/dp(guess);
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen fn = p(guess);
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm }
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen return guess;
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen}
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelenPoly integral(Poly const & p) {
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen Poly result;
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
3515994554d167522343ce57417648b39370ccabcilix result.reserve(p.size()+1);
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen result.push_back(0); // arbitrary const
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen for(unsigned i = 0; i < p.size(); i++) {
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen result.push_back(p[i]/(i+1));
ecf048161ae0284d01d34ca5844204775c0d3fdecilix }
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen return result;
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen}
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelenPoly derivative(Poly const & p) {
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen Poly result;
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen if(p.size() <= 1)
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen return Poly(0);
4afe3fc6b9c122bc5c02b27aea3845ba41384d2acilix result.reserve(p.size()-1);
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen for(unsigned i = 1; i < p.size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm result.push_back(i*p[i]);
72b7b31db250f20b90730d2888e6a554b434a407johanengelen }
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm return result;
a45e91fd09273632e9acf849bb72621832156f07cilix}
494c671e141564431d7d05f141c885d9a2789db5cilix
a45e91fd09273632e9acf849bb72621832156f07cilixPoly compose(Poly const & a, Poly const & b) {
a45e91fd09273632e9acf849bb72621832156f07cilix Poly result;
a45e91fd09273632e9acf849bb72621832156f07cilix
a45e91fd09273632e9acf849bb72621832156f07cilix for(unsigned i = a.size()-1; i >=0; i--) {
8791d7447034c34fdedc50f261bf5c89c34e59f5cilix result = Poly(a[i]) + result * b;
8791d7447034c34fdedc50f261bf5c89c34e59f5cilix }
a45e91fd09273632e9acf849bb72621832156f07cilix return result;
a45e91fd09273632e9acf849bb72621832156f07cilix
a45e91fd09273632e9acf849bb72621832156f07cilix}
a45e91fd09273632e9acf849bb72621832156f07cilix
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm/* This version is backwards - dividing taylor terms
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmPoly divide(Poly const &a, Poly const &b, Poly &r) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm Poly c;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm r = a; // remainder
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm
8c39cbeab9949a0a7d6ae66b768a7352019e42f8johanengelen const unsigned k = a.size();
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen r.resize(k, 0);
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm c.resize(k, 0);
29f9623ba77fc735b89765ae3a13e0c06aabafcecilix
29f9623ba77fc735b89765ae3a13e0c06aabafcecilix for(unsigned i = 0; i < k; i++) {
29f9623ba77fc735b89765ae3a13e0c06aabafcecilix double ci = r[i]/b[0];
072916d0ef7dccd696b59381f50bcf776abccefbjohanengelen c[i] += ci;
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen Poly bb = ci*b;
42e99769805c14a5cc01c805faa3c3b03f9dd1c0johanengelen std::cout << ci <<"*" << b << ", r= " << r << std::endl;
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud r -= bb.shifted(i);
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen }
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud return c;
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen}
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud*/
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud
92fe3142613d000eff89db8a983b3b18b14eee79johanengelenPoly divide(Poly const &a, Poly const &b, Poly &r) {
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud Poly c;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm r = a; // remainder
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm assert(b.size() > 0);
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm const unsigned k = a.degree();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm const unsigned l = b.degree();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm c.resize(k, 0.);
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned i = k; i >= l; i--) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm assert(i >= 0);
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen double ci = r.back()/b.back();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm c[i-l] += ci;
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen Poly bb = ci*b;
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen //std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen // << bb.shifted(i-l) << " r= " << r << std::endl;
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen r -= bb.shifted(i-l);
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen r.pop_back();
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen }
8d9f5d586a04809427ce1df284a5720112177991cilix //std::cout << "r= " << r << std::endl;
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen r.normalize();
8d9f5d586a04809427ce1df284a5720112177991cilix c.normalize();
70eb1fc448cb08acf3468f80fa2296c03b32afd2cilix
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen return c;
70eb1fc448cb08acf3468f80fa2296c03b32afd2cilix}
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyak
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyakPoly gcd(Poly const &a, Poly const &b, const double tol) {
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyak if(a.size() < b.size())
6f4a90e526af850ffc36064f58f09c190f3b633fjohanengelen return gcd(b, a);
6f4a90e526af850ffc36064f58f09c190f3b633fjohanengelen if(b.size() <= 0)
6f4a90e526af850ffc36064f58f09c190f3b633fjohanengelen return a;
f4db63be4e929f4706410914295deccaceea19cdcilix if(b.size() == 1)
f4db63be4e929f4706410914295deccaceea19cdcilix return a;
f4db63be4e929f4706410914295deccaceea19cdcilix Poly r;
ab99111a42436818e6902e044c8f3af2b724263bcilix divide(a, b, r);
ab99111a42436818e6902e044c8f3af2b724263bcilix return gcd(b, r);
ab99111a42436818e6902e044c8f3af2b724263bcilix}
76db360f5f052775326e6d406b9e1e9e2966e11acilix
76db360f5f052775326e6d406b9e1e9e2966e11acilix
b0c42c0dfcd02cc05126371948489a5a88b2e4b3cilix
3d0482af18ffb591c1d8ddecf516629e1bcd2ae4cilix/*Poly divide_out_root(Poly const & p, double x) {
3d0482af18ffb591c1d8ddecf516629e1bcd2ae4cilix assert(1);
3d0482af18ffb591c1d8ddecf516629e1bcd2ae4cilix }*/
64aee804a6a47424f7994e60558351b8cf2ea4dbcilix
64aee804a6a47424f7994e60558351b8cf2ea4dbcilix
64aee804a6a47424f7994e60558351b8cf2ea4dbcilix/*
b320a8d186114a5122ddc3afbe95110eb6cb10cecilix Local Variables:
b320a8d186114a5122ddc3afbe95110eb6cb10cecilix mode:c++
b320a8d186114a5122ddc3afbe95110eb6cb10cecilix c-file-style:"stroustrup"
044d712d4d03f8354962d54e47cfac2346a69ccccilix c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
044d712d4d03f8354962d54e47cfac2346a69ccccilix indent-tabs-mode:nil
044d712d4d03f8354962d54e47cfac2346a69ccccilix fill-column:99
61cfd957cd023c4f432ea0c7307784a56bf978e9cilix End:
61cfd957cd023c4f432ea0c7307784a56bf978e9cilix*/
61cfd957cd023c4f432ea0c7307784a56bf978e9cilix// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
2f5c0701b333a695eedb1680beb1adf95c0723dacilix