poly.cpp revision 981b809bc6ed10a21e89444d9447e5475801874f
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];
6656f193fdace606d1b162d6dea0223bc295f0a6cilix/*double Poly::eval(double x) const {
6656f193fdace606d1b162d6dea0223bc295f0a6cilix return gsl_poly_eval(&coeff[0], size(), x);
04c99c338ffdc6e10cb6f5c18f6f06b3f555e8ebcilix while(back() == 0)
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned i = 0; i < size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm (*this)[i] *= scale;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmstd::vector<std::complex<double> > solve(Poly const & pp) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm double* a = new double[p.size()];
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(int i = 0; i < p.size(); i++)
8d9f5d586a04809427ce1df284a5720112177991cilix a[i] = p[i];
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyak //roots.resize(p.degree());
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]);
61cfd957cd023c4f432ea0c7307784a56bf978e9cilix delete[] z;
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm if(roots[i].imag() == 0) // should be more lenient perhaps
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelendouble polish_root(Poly const & p, double guess, double tol) {
9de5f3bfa8e35ad5e6adcd6eaad4d0d49e7043a3johanengelen for(unsigned i = 0; i < p.size(); i++) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm/* This version is backwards - dividing taylor terms
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrmPoly divide(Poly const &a, Poly const &b, Poly &r) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm r = a; // remainder
8c39cbeab9949a0a7d6ae66b768a7352019e42f8johanengelen const unsigned k = a.size();
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen r.resize(k, 0);
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm c.resize(k, 0);
29f9623ba77fc735b89765ae3a13e0c06aabafcecilix for(unsigned i = 0; i < k; i++) {
29f9623ba77fc735b89765ae3a13e0c06aabafcecilix double ci = r[i]/b[0];
92fe3142613d000eff89db8a983b3b18b14eee79johanengelen Poly bb = ci*b;
42e99769805c14a5cc01c805faa3c3b03f9dd1c0johanengelen std::cout << ci <<"*" << b << ", r= " << r << std::endl;
dc98accfae7a38326b92d74fa4330ac8ccb5b778jfbarraud r -= bb.shifted(i);
92fe3142613d000eff89db8a983b3b18b14eee79johanengelenPoly divide(Poly const &a, Poly const &b, Poly &r) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm r = a; // remainder
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm const unsigned k = a.degree();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm const unsigned l = b.degree();
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm for(unsigned i = k; i >= l; i--) {
f07bfd5a05d43a6d11f7cd442f085149092dea88pjrm c[i-l] += ci;
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen //std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
0563fd55cbad59e8a878e6d4cbbdd8e47f74488djohanengelen // << bb.shifted(i-l) << " r= " << r << std::endl;
8d9f5d586a04809427ce1df284a5720112177991cilix //std::cout << "r= " << r << std::endl;
c169f6cddd2da06cfb761339f445bbd8866f72a8buliabyakPoly gcd(Poly const &a, Poly const &b, const double tol) {
6f4a90e526af850ffc36064f58f09c190f3b633fjohanengelen return gcd(b, a);
6f4a90e526af850ffc36064f58f09c190f3b633fjohanengelen if(b.size() <= 0)
ab99111a42436818e6902e044c8f3af2b724263bcilix return gcd(b, r);
3d0482af18ffb591c1d8ddecf516629e1bcd2ae4cilix/*Poly divide_out_root(Poly const & p, double x) {
b320a8d186114a5122ddc3afbe95110eb6cb10cecilix Local Variables:
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// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :