← Back to branch summary

~valavanisalex/ubuntu/oneiric/inkscape/inkscape_0.48.1-2ubuntu4

« back to all changes in this revision

Viewing changes to src/2geom/sbasis-math.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Kees Cook, Kees Cook, Ted Gould
  • Date: 2008-02-10 14:20:16 UTC
  • mfrom: (1.1.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20080210142016-vcnb2zqyhszu0xvb
Tags: 0.46~pre1-0ubuntu1
https://launchpad.net/bugs/183547
[ Kees Cook ]
* debian/control:
  - add libgtkspell-dev build-dep to gain GtkSpell features (LP: #183547).
  - update Standards version (no changes needed).
  - add Vcs and Homepage fields.
  - switch to new python-lxml dep.
* debian/{control,rules}: switch from dpatch to quilt for more sanity.
* debian/patches/20_fix_glib_and_gxx43_ftbfs.patch:
  - merged against upstream fixes.
  - added additional fixes for newly written code.
* debian/rules: enable parallel building.

[ Ted Gould ]
* Updating POTFILES.in to make it so things build correctly.
* debian/control:
  - add ImageMagick++ and libboost-dev to build-deps

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/2geom/sbasis-math.cpp
 
1
/*
 
2
 *  sbasis-math.cpp - some std functions to work with (pw)s-basis
 
3
 *
 
4
 *  Authors:
 
5
 *   Jean-Francois Barraud
 
6
 *
 
7
 * Copyright (C) 2006-2007 authors
 
8
 *
 
9
 * This library is free software; you can redistribute it and/or
 
10
 * modify it either under the terms of the GNU Lesser General Public
 
11
 * License version 2.1 as published by the Free Software Foundation
 
12
 * (the "LGPL") or, at your option, under the terms of the Mozilla
 
13
 * Public License Version 1.1 (the "MPL"). If you do not alter this
 
14
 * notice, a recipient may use your version of this file under either
 
15
 * the MPL or the LGPL.
 
16
 *
 
17
 * You should have received a copy of the LGPL along with this library
 
18
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
 
19
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
20
 * You should have received a copy of the MPL along with this library
 
21
 * in the file COPYING-MPL-1.1
 
22
 *
 
23
 * The contents of this file are subject to the Mozilla Public License
 
24
 * Version 1.1 (the "License"); you may not use this file except in
 
25
 * compliance with the License. You may obtain a copy of the License at
 
26
 * http://www.mozilla.org/MPL/
 
27
 *
 
28
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
 
29
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
 
30
 * the specific language governing rights and limitations.
 
31
 */
 
32
 
 
33
//this a first try to define sqrt, cos, sin, etc...
 
34
//TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>.
 
35
//TODO: in all these functions, compute 'order' according to 'tol'.
 
36
 
 
37
#include "sbasis-math.h"
 
38
//#define ZERO 1e-3
 
39
 
 
40
 
 
41
namespace Geom {
 
42
 
 
43
 
 
44
#include <stdio.h>
 
45
#include <math.h>
 
46
 
 
47
//-|x|-----------------------------------------------------------------------
 
48
Piecewise<SBasis> abs(SBasis const &f){
 
49
    return abs(Piecewise<SBasis>(f));
 
50
}
 
51
Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
 
52
    Piecewise<SBasis> absf=partition(f,roots(f));
 
53
    for (unsigned i=0; i<absf.size(); i++){
 
54
        if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
 
55
    }
 
56
    return absf;
 
57
}
 
58
 
 
59
//-max(x,y), min(x,y)--------------------------------------------------------
 
60
Piecewise<SBasis> max(          SBasis  const &f,           SBasis  const &g){
 
61
    return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
 
62
}
 
63
Piecewise<SBasis> max(Piecewise<SBasis> const &f,           SBasis  const &g){
 
64
    return max(f,Piecewise<SBasis>(g));
 
65
}
 
66
Piecewise<SBasis> max(          SBasis  const &f, Piecewise<SBasis> const &g){
 
67
    return max(Piecewise<SBasis>(f),g);
 
68
}
 
69
Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
 
70
    Piecewise<SBasis> max=partition(f,roots(f-g));
 
71
    Piecewise<SBasis> gg =partition(g,max.cuts);
 
72
    max = partition(max,gg.cuts);
 
73
    for (unsigned i=0; i<max.size(); i++){
 
74
        if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
 
75
    }
 
76
    return max;
 
77
}
 
78
 
 
79
Piecewise<SBasis> 
 
80
min(          SBasis  const &f,           SBasis  const &g){ return -max(-f,-g); }
 
81
Piecewise<SBasis> 
 
82
min(Piecewise<SBasis> const &f,           SBasis  const &g){ return -max(-f,-g); }
 
83
Piecewise<SBasis> 
 
84
min(          SBasis  const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
 
85
Piecewise<SBasis> 
 
86
min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
 
87
 
 
88
 
 
89
//-sign(x)---------------------------------------------------------------
 
90
Piecewise<SBasis> signSb(SBasis const &f){
 
91
    return signSb(Piecewise<SBasis>(f));
 
92
}
 
93
Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
 
94
    Piecewise<SBasis> sign=partition(f,roots(f));
 
95
    for (unsigned i=0; i<sign.size(); i++){
 
96
        sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
 
97
    }
 
98
    return sign;
 
99
}
 
100
 
 
101
//-Sqrt----------------------------------------------------------
 
102
static Piecewise<SBasis> sqrt_internal(SBasis const &f, 
 
103
                                    double tol, 
 
104
                                    int order){
 
105
    SBasis sqrtf;
 
106
    if(f.isZero() || order == 0){
 
107
        return Piecewise<SBasis>(sqrtf);
 
108
    }
 
109
    if (f.at0()<-tol*tol && f.at1()<-tol*tol){
 
110
        return sqrt_internal(-f,tol,order);
 
111
    }else if (f.at0()>tol*tol && f.at1()>tol*tol){
 
112
        sqrtf.resize(order+1, Linear(0,0));
 
113
        sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1]));
 
114
        SBasis r = f - multiply(sqrtf, sqrtf); // remainder    
 
115
        for(unsigned i = 1; int(i) <= order and i<r.size(); i++) {
 
116
            Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
 
117
            SBasis cisi = shift(ci, i);
 
118
            r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
 
119
            r.truncate(order+1);
 
120
            sqrtf[i] = ci;
 
121
            if(r.tailError(i) == 0) // if exact
 
122
                break;
 
123
        }
 
124
    }else{
 
125
        sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
 
126
    }
 
127
 
 
128
    double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
 
129
    if (err<tol){
 
130
        return Piecewise<SBasis>(sqrtf);
 
131
    }
 
132
 
 
133
    Piecewise<SBasis> sqrtf0,sqrtf1;
 
134
    sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order);
 
135
    sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order);
 
136
    sqrtf0.setDomain(Interval(0.,.5));
 
137
    sqrtf1.setDomain(Interval(.5,1.));
 
138
    sqrtf0.concat(sqrtf1);
 
139
    return sqrtf0;
 
140
}
 
141
 
 
142
Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
 
143
    return sqrt(max(f,Linear(tol*tol)),tol,order);
 
144
}
 
145
 
 
146
Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
 
147
    Piecewise<SBasis> result;
 
148
    Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol));
 
149
    zero.setDomain(f.domain());
 
150
    Piecewise<SBasis> ff=max(f,zero);
 
151
 
 
152
    for (unsigned i=0; i<ff.size(); i++){
 
153
        Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
 
154
        sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
 
155
        result.concat(sqrtfi);
 
156
    }
 
157
    return result;
 
158
}
 
159
 
 
160
//-Yet another sin/cos--------------------------------------------------------------
 
161
 
 
162
Piecewise<SBasis> sin(          SBasis  const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
 
163
Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
 
164
 
 
165
Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
 
166
    Piecewise<SBasis> result;
 
167
    for (unsigned i=0; i<f.size(); i++){
 
168
        Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
 
169
        cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
 
170
        result.concat(cosfi);
 
171
    }
 
172
    return result;
 
173
}
 
174
 
 
175
Piecewise<SBasis> cos(          SBasis  const &f, double tol, int order){
 
176
    double alpha = (f.at0()+f.at1())/2.;
 
177
    SBasis x = f-alpha;
 
178
    double d = x.tailError(0),err=1;
 
179
    //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
 
180
    for (int i=1; i<=2*order; i++) err*=d/i;
 
181
    
 
182
    if (err<tol){
 
183
        SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
 
184
        for (int k=1; k<=2*order; k+=2){
 
185
            xk*=x/k;
 
186
            //take also truncature errors into account...
 
187
            err+=xk.tailError(order);
 
188
            xk.truncate(order);
 
189
            s+=xk;
 
190
            xk*=-x/(k+1);
 
191
            //take also truncature errors into account...
 
192
            err+=xk.tailError(order);
 
193
            xk.truncate(order);
 
194
            c+=xk;
 
195
        }
 
196
        if (err<tol){
 
197
            return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
 
198
        }
 
199
    }
 
200
    Piecewise<SBasis> c0,c1;
 
201
    c0 = cos(compose(f,Linear(0.,.5)),tol,order);
 
202
    c1 = cos(compose(f,Linear(.5,1.)),tol,order);
 
203
    c0.setDomain(Interval(0.,.5));
 
204
    c1.setDomain(Interval(.5,1.));
 
205
    c0.concat(c1);
 
206
    return c0;
 
207
}
 
208
 
 
209
//--1/x------------------------------------------------------------
 
210
//TODO: this implementation is just wrong. Remove or redo!
 
211
 
 
212
void truncateResult(Piecewise<SBasis> &f, int order){
 
213
    if (order>=0){
 
214
        for (unsigned k=0; k<f.segs.size(); k++){
 
215
            f.segs[k].truncate(order);
 
216
        }
 
217
    }
 
218
}
 
219
 
 
220
Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
 
221
    Piecewise<SBasis> reciprocal_fn;
 
222
    //TODO: deduce R from tol...
 
223
    double R=2.;
 
224
    SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
 
225
    double a=range.min(), b=range.max();
 
226
    if (a*b<0){
 
227
        b=std::max(fabs(a),fabs(b));
 
228
        a=0;
 
229
    }else if (b<0){
 
230
        a=-range.max();
 
231
        b=-range.min();
 
232
    }
 
233
 
 
234
    if (a<=tol){
 
235
        reciprocal_fn.push_cut(0);
 
236
        int i0=(int) floor(std::log(tol)/std::log(R));
 
237
        a=pow(R,i0);
 
238
        reciprocal_fn.push(Linear(1/a),a);
 
239
    }else{
 
240
        int i0=(int) floor(std::log(a)/std::log(R));
 
241
        a=pow(R,i0);
 
242
        reciprocal_fn.cuts.push_back(a);
 
243
    }  
 
244
 
 
245
    while (a<b){
 
246
        reciprocal_fn.push(reciprocal1_R/a,R*a);
 
247
        a*=R;
 
248
    }
 
249
    if (range.min()<0 || range.max()<0){
 
250
        Piecewise<SBasis>reciprocal_fn_neg;
 
251
        //TODO: define reverse(pw<sb>);
 
252
        reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
 
253
        for (unsigned i=0; i<reciprocal_fn.size(); i++){
 
254
            int idx=reciprocal_fn.segs.size()-1-i;
 
255
            reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
 
256
            reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
 
257
        }
 
258
        if (range.max()>0){
 
259
            reciprocal_fn_neg.concat(reciprocal_fn);
 
260
        }
 
261
        reciprocal_fn=reciprocal_fn_neg;
 
262
    }
 
263
 
 
264
    return(reciprocal_fn);
 
265
}
 
266
 
 
267
Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
 
268
    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
 
269
    Piecewise<SBasis> result=compose(reciprocal_fn,f);
 
270
    truncateResult(result,order);
 
271
    return(result);
 
272
}
 
273
Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
 
274
    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
 
275
    Piecewise<SBasis> result=compose(reciprocal_fn,f);
 
276
    truncateResult(result,order);
 
277
    return(result);
 
278
}
 
279
 
 
280
}
 
281
 
 
282
/*
 
283
  Local Variables:
 
284
  mode:c++
 
285
  c-file-style:"stroustrup"
 
286
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 
287
  indent-tabs-mode:nil
 
288
  fill-column:99
 
289
  End:
 
290
*/
 
291
// vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :