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