Subversion Repositories Scribus

Rev

Rev 15347 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
12065 fschmid 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
18122 mrdocs 19
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
12065 fschmid 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
 
15347 fschmid 37
#include "sbasis-math.h"
12065 fschmid 38
//#define ZERO 1e-3
39
 
15347 fschmid 40
#include <cstdio>
41
#include <cmath>
42
#include "angle.h"
12665 jghali 43
 
12065 fschmid 44
namespace Geom {
45
 
46
//-|x|-----------------------------------------------------------------------
47
Piecewise<SBasis> abs(SBasis const &f){
48
    return abs(Piecewise<SBasis>(f));
49
}
50
Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
51
    Piecewise<SBasis> absf=partition(f,roots(f));
52
    for (unsigned i=0; i<absf.size(); i++){
53
        if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
54
    }
55
    return absf;
56
}
57
 
58
//-max(x,y), min(x,y)--------------------------------------------------------
59
Piecewise<SBasis> max(          SBasis  const &f,           SBasis  const &g){
60
    return max(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
61
}
62
Piecewise<SBasis> max(Piecewise<SBasis> const &f,           SBasis  const &g){
63
    return max(f,Piecewise<SBasis>(g));
64
}
65
Piecewise<SBasis> max(          SBasis  const &f, Piecewise<SBasis> const &g){
66
    return max(Piecewise<SBasis>(f),g);
67
}
68
Piecewise<SBasis> max(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
69
    Piecewise<SBasis> max=partition(f,roots(f-g));
70
    Piecewise<SBasis> gg =partition(g,max.cuts);
71
    max = partition(max,gg.cuts);
72
    for (unsigned i=0; i<max.size(); i++){
73
        if (max.segs[i](.5)<gg.segs[i](.5)) max.segs[i]=gg.segs[i];
74
    }
75
    return max;
76
}
77
 
78
Piecewise<SBasis>
79
min(          SBasis  const &f,           SBasis  const &g){ return -max(-f,-g); }
80
Piecewise<SBasis>
81
min(Piecewise<SBasis> const &f,           SBasis  const &g){ return -max(-f,-g); }
82
Piecewise<SBasis>
83
min(          SBasis  const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
84
Piecewise<SBasis>
85
min(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -max(-f,-g); }
86
 
87
 
88
//-sign(x)---------------------------------------------------------------
89
Piecewise<SBasis> signSb(SBasis const &f){
90
    return signSb(Piecewise<SBasis>(f));
91
}
92
Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
93
    Piecewise<SBasis> sign=partition(f,roots(f));
94
    for (unsigned i=0; i<sign.size(); i++){
95
        sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
96
    }
97
    return sign;
98
}
99
 
100
//-Sqrt----------------------------------------------------------
101
static Piecewise<SBasis> sqrt_internal(SBasis const &f,
102
                                    double tol,
103
                                    int order){
104
    SBasis sqrtf;
105
    if(f.isZero() || order == 0){
106
        return Piecewise<SBasis>(sqrtf);
107
    }
108
    if (f.at0()<-tol*tol && f.at1()<-tol*tol){
109
        return sqrt_internal(-f,tol,order);
110
    }else if (f.at0()>tol*tol && f.at1()>tol*tol){
111
        sqrtf.resize(order+1, Linear(0,0));
112
        sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1]));
113
        SBasis r = f - multiply(sqrtf, sqrtf); // remainder    
15347 fschmid 114
        for(unsigned i = 1; int(i) <= order && i<r.size(); ++i) {
12065 fschmid 115
            Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
116
            SBasis cisi = shift(ci, i);
117
            r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
118
            r.truncate(order+1);
119
            sqrtf[i] = ci;
120
            if(r.tailError(i) == 0) // if exact
121
                break;
122
        }
123
    }else{
124
        sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
125
    }
126
 
127
    double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
128
    if (err<tol){
129
        return Piecewise<SBasis>(sqrtf);
130
    }
131
 
132
    Piecewise<SBasis> sqrtf0,sqrtf1;
133
    sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order);
134
    sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order);
135
    sqrtf0.setDomain(Interval(0.,.5));
136
    sqrtf1.setDomain(Interval(.5,1.));
137
    sqrtf0.concat(sqrtf1);
138
    return sqrtf0;
139
}
140
 
141
Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
142
    return sqrt(max(f,Linear(tol*tol)),tol,order);
143
}
144
 
145
Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
146
    Piecewise<SBasis> result;
147
    Piecewise<SBasis> zero = Piecewise<SBasis>(Linear(tol*tol));
148
    zero.setDomain(f.domain());
149
    Piecewise<SBasis> ff=max(f,zero);
150
 
151
    for (unsigned i=0; i<ff.size(); i++){
152
        Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
153
        sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
154
        result.concat(sqrtfi);
155
    }
156
    return result;
157
}
158
 
159
//-Yet another sin/cos--------------------------------------------------------------
160
 
15347 fschmid 161
Piecewise<SBasis> sin(          SBasis  const &f, double tol, int order){return cos(-f+M_PI_2,tol,order);}
162
Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return cos(-f+M_PI_2,tol,order);}
12065 fschmid 163
 
164
Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
165
    Piecewise<SBasis> result;
166
    for (unsigned i=0; i<f.size(); i++){
167
        Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
168
        cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
169
        result.concat(cosfi);
170
    }
171
    return result;
172
}
173
 
174
Piecewise<SBasis> cos(          SBasis  const &f, double tol, int order){
175
    double alpha = (f.at0()+f.at1())/2.;
176
    SBasis x = f-alpha;
177
    double d = x.tailError(0),err=1;
178
    //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
179
    for (int i=1; i<=2*order; i++) err*=d/i;
180
 
181
    if (err<tol){
182
        SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
183
        for (int k=1; k<=2*order; k+=2){
184
            xk*=x/k;
185
            //take also truncature errors into account...
186
            err+=xk.tailError(order);
187
            xk.truncate(order);
188
            s+=xk;
189
            xk*=-x/(k+1);
190
            //take also truncature errors into account...
191
            err+=xk.tailError(order);
192
            xk.truncate(order);
193
            c+=xk;
194
        }
195
        if (err<tol){
196
            return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
197
        }
198
    }
199
    Piecewise<SBasis> c0,c1;
200
    c0 = cos(compose(f,Linear(0.,.5)),tol,order);
201
    c1 = cos(compose(f,Linear(.5,1.)),tol,order);
202
    c0.setDomain(Interval(0.,.5));
203
    c1.setDomain(Interval(.5,1.));
204
    c0.concat(c1);
205
    return c0;
206
}
207
 
15347 fschmid 208
 
12065 fschmid 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){
15347 fschmid 268
    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
12065 fschmid 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){
15347 fschmid 274
    Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
12065 fschmid 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
*/
15347 fschmid 291
// vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :