sbasis-math.cpp revision 29684a16b6c92bee28a94fdc2607bcc143950fa8
/*
* sbasis-math.cpp - some std functions to work with (pw)s-basis
*
* Authors:
* Jean-Francois Barraud
*
* Copyright (C) 2006-2007 authors
*
* This library is free software; you can redistribute it and/or
* modify it either under the terms of the GNU Lesser General Public
* License version 2.1 as published by the Free Software Foundation
* (the "LGPL") or, at your option, under the terms of the Mozilla
* Public License Version 1.1 (the "MPL"). If you do not alter this
* notice, a recipient may use your version of this file under either
* the MPL or the LGPL.
*
* You should have received a copy of the LGPL along with this library
* in the file COPYING-LGPL-2.1; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
* You should have received a copy of the MPL along with this library
* in the file COPYING-MPL-1.1
*
* The contents of this file are subject to the Mozilla Public License
* Version 1.1 (the "License"); you may not use this file except in
* compliance with the License. You may obtain a copy of the License at
* http://www.mozilla.org/MPL/
*
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
* OF ANY KIND, either express or implied. See the LGPL or the MPL for
* the specific language governing rights and limitations.
*/
//this a first try to define sqrt, cos, sin, etc...
//TODO: define a truncated compose(sb,sb, order) and extend it to pw<sb>.
//TODO: in all these functions, compute 'order' according to 'tol'.
#include "sbasis-math.h"
//#define ZERO 1e-3
namespace Geom {
#include <stdio.h>
#include <math.h>
//-|x|-----------------------------------------------------------------------
Piecewise<SBasis> abs(SBasis const &f){
return abs(Piecewise<SBasis>(f));
}
Piecewise<SBasis> abs(Piecewise<SBasis> const &f){
Piecewise<SBasis> absf=partition(f,roots(f));
for (unsigned i=0; i<absf.size(); i++){
if (absf.segs[i](.5)<0) absf.segs[i]*=-1;
}
return absf;
}
//-maxSb(x,y), minSb(x,y)--------------------------------------------------------
Piecewise<SBasis> maxSb( SBasis const &f, SBasis const &g){
return maxSb(Piecewise<SBasis>(f),Piecewise<SBasis>(g));
}
Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, SBasis const &g){
return maxSb(f,Piecewise<SBasis>(g));
}
Piecewise<SBasis> maxSb( SBasis const &f, Piecewise<SBasis> const &g){
return maxSb(Piecewise<SBasis>(f),g);
}
Piecewise<SBasis> maxSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){
Piecewise<SBasis> maxSb=partition(f,roots(f-g));
Piecewise<SBasis> gg =partition(g,maxSb.cuts);
maxSb = partition(maxSb,gg.cuts);
for (unsigned i=0; i<maxSb.size(); i++){
if (maxSb.segs[i](.5)<gg.segs[i](.5)) maxSb.segs[i]=gg.segs[i];
}
return maxSb;
}
Piecewise<SBasis>
minSb( SBasis const &f, SBasis const &g){ return -maxSb(-f,-g); }
Piecewise<SBasis>
minSb(Piecewise<SBasis> const &f, SBasis const &g){ return -maxSb(-f,-g); }
Piecewise<SBasis>
minSb( SBasis const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
Piecewise<SBasis>
minSb(Piecewise<SBasis> const &f, Piecewise<SBasis> const &g){ return -maxSb(-f,-g); }
//-sign(x)---------------------------------------------------------------
Piecewise<SBasis> signSb(SBasis const &f){
return signSb(Piecewise<SBasis>(f));
}
Piecewise<SBasis> signSb(Piecewise<SBasis> const &f){
Piecewise<SBasis> sign=partition(f,roots(f));
for (unsigned i=0; i<sign.size(); i++){
sign.segs[i] = (sign.segs[i](.5)<0)? Linear(-1.):Linear(1.);
}
return sign;
}
//-Sqrt----------------------------------------------------------
static Piecewise<SBasis> sqrt_internal(SBasis const &f,
double tol,
int order){
SBasis sqrtf;
if(f.isZero() || order == 0){
return Piecewise<SBasis>(sqrtf);
}
if (f.at0()<-tol*tol && f.at1()<-tol*tol){
return sqrt_internal(-f,tol,order);
}else if (f.at0()>tol*tol && f.at1()>tol*tol){
sqrtf.resize(order+1, Linear(0,0));
sqrtf[0] = Linear(std::sqrt(f[0][0]), std::sqrt(f[0][1]));
SBasis r = f - multiply(sqrtf, sqrtf); // remainder
for(unsigned i = 1; int(i) <= order and i<r.size(); i++) {
Linear ci(r[i][0]/(2*sqrtf[0][0]), r[i][1]/(2*sqrtf[0][1]));
SBasis cisi = shift(ci, i);
r -= multiply(shift((sqrtf*2 + cisi), i), SBasis(ci));
r.truncate(order+1);
sqrtf[i] = ci;
if(r.tailError(i) == 0) // if exact
break;
}
}else{
sqrtf = Linear(std::sqrt(fabs(f.at0())), std::sqrt(fabs(f.at1())));
}
double err = (f - multiply(sqrtf, sqrtf)).tailError(0);
if (err<tol){
return Piecewise<SBasis>(sqrtf);
}
Piecewise<SBasis> sqrtf0,sqrtf1;
sqrtf0 = sqrt_internal(compose(f,Linear(0.,.5)),tol,order);
sqrtf1 = sqrt_internal(compose(f,Linear(.5,1.)),tol,order);
sqrtf0.setDomain(Interval(0.,.5));
sqrtf1.setDomain(Interval(.5,1.));
sqrtf0.concat(sqrtf1);
return sqrtf0;
}
Piecewise<SBasis> sqrt(SBasis const &f, double tol, int order){
return sqrt(maxSb(f,Linear(tol*tol)),tol,order);
}
Piecewise<SBasis> sqrt(Piecewise<SBasis> const &f, double tol, int order){
Piecewise<SBasis> result;
Piecewise<SBasis> ff=maxSb(f,Linear(tol*tol));
for (unsigned i=0; i<ff.size(); i++){
Piecewise<SBasis> sqrtfi = sqrt_internal(ff.segs[i],tol,order);
sqrtfi.setDomain(Interval(ff.cuts[i],ff.cuts[i+1]));
result.concat(sqrtfi);
}
return result;
}
//-Yet another sin/cos--------------------------------------------------------------
Piecewise<SBasis> sin( SBasis const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
Piecewise<SBasis> sin(Piecewise<SBasis> const &f, double tol, int order){return(cos(-f+M_PI/2,tol,order));}
Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
Piecewise<SBasis> result;
for (unsigned i=0; i<f.size(); i++){
Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
result.concat(cosfi);
}
return result;
}
Piecewise<SBasis> cos( SBasis const &f, double tol, int order){
double alpha = (f.at0()+f.at1())/2.;
SBasis x = f-alpha;
double d = x.tailError(0),err=1;
//estimate cos(x)-sum_0^order (-1)^k x^2k/2k! by the first neglicted term
for (int i=1; i<=2*order; i++) err*=d/i;
if (err<tol){
SBasis xk=Linear(1), c=Linear(1), s=Linear(0);
for (int k=1; k<=2*order; k+=2){
xk*=x/k;
//take also truncature errors into account...
err+=xk.tailError(order);
xk.truncate(order);
s+=xk;
xk*=-x/(k+1);
//take also truncature errors into account...
err+=xk.tailError(order);
xk.truncate(order);
c+=xk;
}
if (err<tol){
return Piecewise<SBasis>(std::cos(alpha)*c-std::sin(alpha)*s);
}
}
Piecewise<SBasis> c0,c1;
c0 = cos(compose(f,Linear(0.,.5)),tol,order);
c1 = cos(compose(f,Linear(.5,1.)),tol,order);
c0.setDomain(Interval(0.,.5));
c1.setDomain(Interval(.5,1.));
c0.concat(c1);
return c0;
}
//--1/x------------------------------------------------------------
//TODO: this implementation is just wrong. Remove or redo!
void truncateResult(Piecewise<SBasis> &f, int order){
if (order>=0){
for (unsigned k=0; k<f.segs.size(); k++){
f.segs[k].truncate(order);
}
}
}
Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
Piecewise<SBasis> reciprocal_fn;
//TODO: deduce R from tol...
double R=2.;
SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
double a=range.min(), b=range.max();
if (a*b<0){
b=std::max(fabs(a),fabs(b));
a=0;
}else if (b<0){
a=-range.max();
b=-range.min();
}
if (a<=tol){
reciprocal_fn.push_cut(0);
int i0=(int) floor(std::log(tol)/std::log(R));
a=pow(R,i0);
reciprocal_fn.push(Linear(1/a),a);
}else{
int i0=(int) floor(std::log(a)/std::log(R));
a=pow(R,i0);
reciprocal_fn.cuts.push_back(a);
}
while (a<b){
reciprocal_fn.push(reciprocal1_R/a,R*a);
a*=R;
}
if (range.min()<0 || range.max()<0){
Piecewise<SBasis>reciprocal_fn_neg;
//TODO: define reverse(pw<sb>);
reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
for (unsigned i=0; i<reciprocal_fn.size(); i++){
int idx=reciprocal_fn.segs.size()-1-i;
reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
}
if (range.max()>0){
reciprocal_fn_neg.concat(reciprocal_fn);
}
reciprocal_fn=reciprocal_fn_neg;
}
return(reciprocal_fn);
}
Piecewise<SBasis> reciprocal(SBasis const &f, double tol, int order){
Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
Piecewise<SBasis> result=compose(reciprocal_fn,f);
truncateResult(result,order);
return(result);
}
Piecewise<SBasis> reciprocal(Piecewise<SBasis> const &f, double tol, int order){
Piecewise<SBasis> reciprocal_fn=reciprocalOnDomain(bounds_fast(f), tol);
Piecewise<SBasis> result=compose(reciprocal_fn,f);
truncateResult(result,order);
return(result);
}
}
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :