1 #include <2geom/poly.h>
3 Poly Poly::operator*(const Poly& p) const {
4 Poly result;
5 result.resize(degree() + p.degree()+1);
7 for(unsigned i = 0; i < size(); i++) {
8 for(unsigned j = 0; j < p.size(); j++) {
9 result[i+j] += (*this)[i] * p[j];
10 }
11 }
12 return result;
13 }
14 #define HAVE_GSL
15 #ifdef HAVE_GSL
16 #include <gsl/gsl_poly.h>
17 #endif
19 /*double Poly::eval(double x) const {
20 return gsl_poly_eval(&coeff[0], size(), x);
21 }*/
23 void Poly::normalize() {
24 while(back() == 0)
25 pop_back();
26 }
28 void Poly::monicify() {
29 normalize();
31 double scale = 1./back(); // unitize
33 for(unsigned i = 0; i < size(); i++) {
34 (*this)[i] *= scale;
35 }
36 }
39 #ifdef HAVE_GSL
40 std::vector<std::complex<double> > solve(Poly const & pp) {
41 Poly p(pp);
42 p.normalize();
43 gsl_poly_complex_workspace * w
44 = gsl_poly_complex_workspace_alloc (p.size());
46 gsl_complex_packed_ptr z = new double[p.degree()*2];
47 double* a = new double[p.size()];
48 for(unsigned int i = 0; i < p.size(); i++)
49 a[i] = p[i];
50 std::vector<std::complex<double> > roots;
51 //roots.resize(p.degree());
53 gsl_poly_complex_solve (a, p.size(), w, z);
54 delete[]a;
56 gsl_poly_complex_workspace_free (w);
58 for (unsigned int i = 0; i < p.degree(); i++) {
59 roots.push_back(std::complex<double> (z[2*i] ,z[2*i+1]));
60 //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
61 }
62 delete[] z;
63 return roots;
64 }
66 std::vector<double > solve_reals(Poly const & p) {
67 std::vector<std::complex<double> > roots = solve(p);
68 std::vector<double> real_roots;
70 for(unsigned int i = 0; i < roots.size(); i++) {
71 if(roots[i].imag() == 0) // should be more lenient perhaps
72 real_roots.push_back(roots[i].real());
73 }
74 return real_roots;
75 }
76 #endif
78 double polish_root(Poly const & p, double guess, double tol) {
79 Poly dp = derivative(p);
81 double fn = p(guess);
82 while(fabs(fn) > tol) {
83 guess -= fn/dp(guess);
84 fn = p(guess);
85 }
86 return guess;
87 }
89 Poly integral(Poly const & p) {
90 Poly result;
92 result.reserve(p.size()+1);
93 result.push_back(0); // arbitrary const
94 for(unsigned i = 0; i < p.size(); i++) {
95 result.push_back(p[i]/(i+1));
96 }
97 return result;
99 }
101 Poly derivative(Poly const & p) {
102 Poly result;
104 if(p.size() <= 1)
105 return Poly(0);
106 result.reserve(p.size()-1);
107 for(unsigned i = 1; i < p.size(); i++) {
108 result.push_back(i*p[i]);
109 }
110 return result;
111 }
113 Poly compose(Poly const & a, Poly const & b) {
114 Poly result;
116 for(unsigned i = a.size(); i > 0; i--) {
117 result = Poly(a[i-1]) + result * b;
118 }
119 return result;
121 }
123 /* This version is backwards - dividing taylor terms
124 Poly divide(Poly const &a, Poly const &b, Poly &r) {
125 Poly c;
126 r = a; // remainder
128 const unsigned k = a.size();
129 r.resize(k, 0);
130 c.resize(k, 0);
132 for(unsigned i = 0; i < k; i++) {
133 double ci = r[i]/b[0];
134 c[i] += ci;
135 Poly bb = ci*b;
136 std::cout << ci <<"*" << b << ", r= " << r << std::endl;
137 r -= bb.shifted(i);
138 }
140 return c;
141 }
142 */
144 Poly divide(Poly const &a, Poly const &b, Poly &r) {
145 Poly c;
146 r = a; // remainder
147 assert(b.size() > 0);
149 const unsigned k = a.degree();
150 const unsigned l = b.degree();
151 c.resize(k, 0.);
153 for(unsigned i = k; i >= l; i--) {
154 double ci = r.back()/b.back();
155 c[i-l] += ci;
156 Poly bb = ci*b;
157 //std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
158 // << bb.shifted(i-l) << " r= " << r << std::endl;
159 r -= bb.shifted(i-l);
160 r.pop_back();
161 }
162 //std::cout << "r= " << r << std::endl;
163 r.normalize();
164 c.normalize();
166 return c;
167 }
169 Poly gcd(Poly const &a, Poly const &b, const double /*tol*/) {
170 if(a.size() < b.size())
171 return gcd(b, a);
172 if(b.size() <= 0)
173 return a;
174 if(b.size() == 1)
175 return a;
176 Poly r;
177 divide(a, b, r);
178 return gcd(b, r);
179 }
183 /*Poly divide_out_root(Poly const & p, double x) {
184 assert(1);
185 }*/
188 /*
189 Local Variables:
190 mode:c++
191 c-file-style:"stroustrup"
192 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
193 indent-tabs-mode:nil
194 fill-column:99
195 End:
196 */
197 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :