← Back to branch summary

~valavanisalex/ubuntu/precise/inkscape/fix-943984

« back to all changes in this revision

Viewing changes to inkscape-0.47pre1/src/2geom/poly.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Bryce Harrington
  • Date: 2009-07-02 17:09:45 UTC
  • mfrom: (1.1.9 upstream)
  • Revision ID: james.westby@ubuntu.com-20090702170945-nn6d6zswovbwju1t
Tags: 0.47~pre1-0ubuntu1
https://launchpad.net/bugs/348842
https://launchpad.net/bugs/391149
* New upstream release.
  - Don't constrain maximization on small resolution devices (pre0)
    (LP: #348842)
  - Fixes segfault on startup (pre0)
    (LP: #391149)

Show diffs side-by-side

added added

removed removed

Lines of Context:
inkscape-0.47pre1/src/2geom/poly.cpp
 
1
#include <2geom/poly.h>
 
2
 
 
3
#define HAVE_GSL
 
4
#ifdef HAVE_GSL
 
5
#include <gsl/gsl_poly.h>
 
6
#endif
 
7
 
 
8
namespace Geom {
 
9
 
 
10
Poly Poly::operator*(const Poly& p) const {
 
11
    Poly result; 
 
12
    result.resize(degree() +  p.degree()+1);
 
13
    
 
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
}
 
21
 
 
22
/*double Poly::eval(double x) const {
 
23
    return gsl_poly_eval(&coeff[0], size(), x);
 
24
    }*/
 
25
 
 
26
void Poly::normalize() {
 
27
    while(back() == 0)
 
28
        pop_back();
 
29
}
 
30
 
 
31
void Poly::monicify() {
 
32
    normalize();
 
33
    
 
34
    double scale = 1./back(); // unitize
 
35
    
 
36
    for(unsigned i = 0; i < size(); i++) {
 
37
        (*this)[i] *= scale;
 
38
    }
 
39
}
 
40
 
 
41
 
 
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());
 
48
       
 
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());
 
55
    
 
56
    gsl_poly_complex_solve (a, p.size(), w, z);
 
57
    delete[]a;
 
58
     
 
59
    gsl_poly_complex_workspace_free (w);
 
60
     
 
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
}
 
68
 
 
69
std::vector<double > solve_reals(Poly const & p) {
 
70
    std::vector<std::complex<double> > roots = solve(p);
 
71
    std::vector<double> real_roots;
 
72
    
 
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
 
80
 
 
81
double polish_root(Poly const & p, double guess, double tol) {
 
82
    Poly dp = derivative(p);
 
83
    
 
84
    double fn = p(guess);
 
85
    while(fabs(fn) > tol) {
 
86
        guess -= fn/dp(guess);
 
87
        fn = p(guess);
 
88
    }
 
89
    return guess;
 
90
}
 
91
 
 
92
Poly integral(Poly const & p) {
 
93
    Poly result;
 
94
    
 
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;
 
101
 
 
102
}
 
103
 
 
104
Poly derivative(Poly const & p) {
 
105
    Poly result;
 
106
    
 
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
}
 
115
 
 
116
Poly compose(Poly const & a, Poly const & b) {
 
117
    Poly result;
 
118
    
 
119
    for(unsigned i = a.size(); i > 0; i--) {
 
120
        result = Poly(a[i-1]) + result * b;
 
121
    }
 
122
    return result;
 
123
    
 
124
}
 
125
 
 
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
 
130
    
 
131
    const unsigned k = a.size();
 
132
    r.resize(k, 0);
 
133
    c.resize(k, 0);
 
134
 
 
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
    }
 
142
    
 
143
    return c;
 
144
}
 
145
*/
 
146
 
 
147
Poly divide(Poly const &a, Poly const &b, Poly &r) {
 
148
    Poly c;
 
149
    r = a; // remainder
 
150
    assert(b.size() > 0);
 
151
    
 
152
    const unsigned k = a.degree();
 
153
    const unsigned l = b.degree();
 
154
    c.resize(k, 0.);
 
155
    
 
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();
 
169
    
 
170
    return c;
 
171
}
 
172
 
 
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
}
 
184
 
 
185
 
 
186
 
 
187
/*Poly divide_out_root(Poly const & p, double x) {
 
188
    assert(1);
 
189
    }*/
 
190
 
 
191
} //namespace Geom
 
192
 
 
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 :