db4cf5e08dee315c62eee2d354dbfe0a6f09d5f9
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 */
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'.
37 #include "sbasis-math.h"
38 //#define ZERO 1e-3
41 namespace Geom {
44 #include <stdio.h>
45 #include <math.h>
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 }
59 //-maxSb(x,y), minSb(x,y)--------------------------------------------------------
60 Piecewise<SBasis> maxSb( SBasis const &f, SBasis const &g){
61 return maxSb(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
62 }
63 Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, SBasis const &g){
64 return maxSb(f,Piecewise<SBasis>(g));
65 }
66 Piecewise<SBasis> maxSb( SBasis const &f, Piecewise<SBasis> const &g){
67 return maxSb(Piecewise<SBasis>(f),g);
68 }
69 Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
70 Piecewise<SBasis> maxSb=partition(f,roots(f-g));
71 Piecewise<SBasis> gg =partition(g,maxSb.cuts);
72 maxSb = partition(maxSb,gg.cuts);
73 for (unsigned i=0; i<maxSb.size(); i++){
74 if (maxSb.segs[i](.5)<gg.segs[i](.5)) maxSb.segs[i]=gg.segs[i];
75 }
76 return maxSb;
77 }
79 Piecewise<SBasis>
80 minSb( SBasis const &f, SBasis const &g){ return -maxSb(-f,-g); }
81 Piecewise<SBasis>
82 minSb(Piecewise<SBasis> const &f, SBasis const &g){ return -maxSb(-f,-g); }
83 Piecewise<SBasis>
84 minSb( SBasis const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
85 Piecewise<SBasis>
86 minSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
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 }
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 }
128 double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
129 if (err<tol){
130 return Piecewise<SBasis>(sqrtf);
131 }
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 }
142 Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
143 return sqrt(maxSb(f,Linear(tol*tol)),tol,order);
144 }
146 Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
147 Piecewise<SBasis> result;
148 Piecewise<SBasis> ff=maxSb(f,Linear(tol*tol));
150 for (unsigned i=0; i<ff.size(); i++){
151 Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
152 sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
153 result.concat(sqrtfi);
154 }
155 return result;
156 }
158 //-Yet another sin/cos--------------------------------------------------------------
160 Piecewise<SBasis> sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
161 Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
163 Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
164 Piecewise<SBasis> result;
165 for (unsigned i=0; i<f.size(); i++){
166 Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
167 cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
168 result.concat(cosfi);
169 }
170 return result;
171 }
173 Piecewise<SBasis> cos( SBasis const &f, double tol, int order){
174 double alpha = (f.at0()+f.at1())/2.;
175 SBasis x = f-alpha;
176 double d = x.tailError(0),err=1;
177 //estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
178 for (int i=1; i<=2*order; i++) err*=d/i;
180 if (err<tol){
181 SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
182 for (int k=1; k<=2*order; k+=2){
183 xk*=x/k;
184 //take also truncature errors into account...
185 err+=xk.tailError(order);
186 xk.truncate(order);
187 s+=xk;
188 xk*=-x/(k+1);
189 //take also truncature errors into account...
190 err+=xk.tailError(order);
191 xk.truncate(order);
192 c+=xk;
193 }
194 if (err<tol){
195 return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
196 }
197 }
198 Piecewise<SBasis> c0,c1;
199 c0 = cos(compose(f,Linear(0.,.5)),tol,order);
200 c1 = cos(compose(f,Linear(.5,1.)),tol,order);
201 c0.setDomain(Interval(0.,.5));
202 c1.setDomain(Interval(.5,1.));
203 c0.concat(c1);
204 return c0;
205 }
207 //--1/x------------------------------------------------------------
208 //TODO: this implementation is just wrong. Remove or redo!
210 void truncateResult(Piecewise<SBasis> &f, int order){
211 if (order>=0){
212 for (unsigned k=0; k<f.segs.size(); k++){
213 f.segs[k].truncate(order);
214 }
215 }
216 }
218 Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
219 Piecewise<SBasis> reciprocal_fn;
220 //TODO: deduce R from tol...
221 double R=2.;
222 SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
223 double a=range.min(), b=range.max();
224 if (a*b<0){
225 b=std::max(fabs(a),fabs(b));
226 a=0;
227 }else if (b<0){
228 a=-range.max();
229 b=-range.min();
230 }
232 if (a<=tol){
233 reciprocal_fn.push_cut(0);
234 int i0=(int) floor(std::log(tol)/std::log(R));
235 a=pow(R,i0);
236 reciprocal_fn.push(Linear(1/a),a);
237 }else{
238 int i0=(int) floor(std::log(a)/std::log(R));
239 a=pow(R,i0);
240 reciprocal_fn.cuts.push_back(a);
241 }
243 while (a<b){
244 reciprocal_fn.push(reciprocal1_R/a,R*a);
245 a*=R;
246 }
247 if (range.min()<0 || range.max()<0){
248 Piecewise<SBasis>reciprocal_fn_neg;
249 //TODO: define reverse(pw<sb>);
250 reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
251 for (unsigned i=0; i<reciprocal_fn.size(); i++){
252 int idx=reciprocal_fn.segs.size()-1-i;
253 reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
254 reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
255 }
256 if (range.max()>0){
257 reciprocal_fn_neg.concat(reciprocal_fn);
258 }
259 reciprocal_fn=reciprocal_fn_neg;
260 }
262 return(reciprocal_fn);
263 }
265 Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
266 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
267 Piecewise<SBasis> result=compose(reciprocal_fn,f);
268 truncateResult(result,order);
269 return(result);
270 }
271 Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
272 Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
273 Piecewise<SBasis> result=compose(reciprocal_fn,f);
274 truncateResult(result,order);
275 return(result);
276 }
278 }
280 /*
281 Local Variables:
282 mode:c++
283 c-file-style:"stroustrup"
284 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
285 indent-tabs-mode:nil
286 fill-column:99
287 End:
288 */
289 // vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :