43 #include "exception.h"
48 class SBasis :
public std::vector<Linear>{
51 explicit SBasis(
double a) {
55 std::vector<Linear>(a)
62 typedef double output_type;
63 inline bool isZero()
const {
64 if(empty())
return true;
65 for(
unsigned i = 0; i < size(); i++) {
66 if(!(*
this)[i].isZero())
return false;
70 inline bool isConstant()
const {
71 if (empty())
return true;
72 for (
unsigned i = 0; i < size(); i++) {
73 if(!(*
this)[i].isConstant())
return false;
78 bool isFinite()
const;
79 inline double at0()
const {
80 if(empty())
return 0;
else return (*
this)[0][0];
82 inline double at1()
const{
83 if(empty())
return 0;
else return (*
this)[0][1];
86 double valueAt(
double t)
const {
88 double p0 = 0, p1 = 0;
91 for(
unsigned k = 0; k < size(); k++) {
92 p0 += sk*(*this)[k][0];
93 p1 += sk*(*this)[k][1];
96 return (1-t)*p0 + t*p1;
98 double valueAndDerivative(
double t,
double &der)
const {
100 double p0 = 0, p1 = 0;
103 for(
unsigned k = 0; k < size(); k++) {
104 p0 += sk*(*this)[k][0];
105 p1 += sk*(*this)[k][1];
110 return (1-t)*p0 + t*p1;
112 double operator()(
double t)
const {
116 std::vector<double> valueAndDerivatives(
double ,
unsigned )
const {
118 throwNotImplemented(0);
123 double tailError(
unsigned tail)
const;
128 Linear operator[](
unsigned i)
const {
130 return std::vector<Linear>::operator[](i);
134 Linear& operator[](
unsigned i) {
return this->at(i); }
138 while(!empty() && 0 == back()[0] && 0 == back()[1])
141 void truncate(
unsigned k) {
if(k < size()) resize(k); }
145 inline SBasis Linear::toSBasis()
const {
return SBasis(*
this); }
148 Interval bounds_exact(SBasis
const &a);
149 Interval bounds_fast(SBasis
const &a,
int order = 0);
150 Interval bounds_local(SBasis
const &a,
const Interval &t,
int order = 0);
152 inline SBasis reverse(SBasis
const &a) {
154 result.reserve(a.size());
155 for(
unsigned k = 0; k < a.size(); k++)
156 result.push_back(reverse(a[k]));
161 inline SBasis operator-(
const SBasis& p) {
162 if(p.isZero())
return SBasis();
164 result.reserve(p.size());
166 for(
unsigned i = 0; i < p.size(); i++) {
167 result.push_back(-p[i]);
171 SBasis operator*(SBasis
const &a,
double k);
172 inline SBasis operator*(
double k, SBasis
const &a) {
return a*k; }
173 inline SBasis operator/(SBasis
const &a,
double k) {
return a*(1./k); }
174 SBasis& operator*=(SBasis& a,
double b);
175 inline SBasis& operator/=(SBasis& a,
double b) {
return (a*=(1./b)); }
178 SBasis operator+(
const SBasis& a,
const SBasis& b);
179 SBasis operator-(
const SBasis& a,
const SBasis& b);
180 SBasis& operator+=(SBasis& a,
const SBasis& b);
181 SBasis& operator-=(SBasis& a,
const SBasis& b);
184 inline SBasis operator+(
const SBasis & a, Linear
const & b) {
185 if(b.isZero())
return a;
186 if(a.isZero())
return b;
191 inline SBasis operator-(
const SBasis & a, Linear
const & b) {
192 if(b.isZero())
return a;
197 inline SBasis& operator+=(SBasis& a,
const Linear& b) {
204 inline SBasis& operator-=(SBasis& a,
const Linear& b) {
213 inline SBasis operator+(
const SBasis & a,
double b) {
214 if(a.isZero())
return Linear(b, b);
219 inline SBasis operator-(
const SBasis & a,
double b) {
220 if(a.isZero())
return Linear(-b, -b);
225 inline SBasis& operator+=(SBasis& a,
double b) {
227 a.push_back(Linear(b,b));
232 inline SBasis& operator-=(SBasis& a,
double b) {
234 a.push_back(Linear(-b,-b));
240 SBasis shift(SBasis
const &a,
int sh);
241 SBasis shift(Linear
const &a,
int sh);
243 inline SBasis truncate(SBasis
const &a,
unsigned terms) {
245 c.insert(c.begin(), a.begin(), a.begin() + std::min(terms, (
unsigned)a.size()));
249 SBasis multiply(SBasis
const &a, SBasis
const &b);
251 SBasis integral(SBasis
const &c);
252 SBasis derivative(SBasis
const &a);
254 SBasis sqrt(SBasis
const &a,
int k);
257 SBasis reciprocal(Linear
const &a,
int k);
258 SBasis divide(SBasis
const &a, SBasis
const &b,
int k);
260 inline SBasis operator*(SBasis
const & a, SBasis
const & b) {
261 return multiply(a, b);
264 inline SBasis& operator*=(SBasis& a, SBasis
const & b) {
271 valuation(SBasis
const &a,
double tol=0){
273 while( val<a.size() &&
274 fabs(a[val][0])<tol &&
275 fabs(a[val][1])<tol )
281 SBasis compose(SBasis
const &a, SBasis
const &b);
282 SBasis compose(SBasis
const &a, SBasis
const &b,
unsigned k);
283 SBasis inverse(SBasis a,
int k);
286 SBasis compose_inverse(SBasis
const &f, SBasis
const &g,
unsigned order=2,
double tol=1e-3);
288 inline SBasis portion(
const SBasis &t,
double from,
double to) {
return compose(t, Linear(from, to)); }
292 SBasis::operator()(SBasis
const & g)
const {
293 return compose(*
this, g);
296 inline std::ostream &operator<< (std::ostream &out_file,
const Linear &bo) {
297 out_file <<
"{" << bo[0] <<
", " << bo[1] <<
"}";
301 inline std::ostream &operator<< (std::ostream &out_file,
const SBasis & p) {
302 for(
unsigned i = 0; i < p.size(); i++) {
303 out_file << p[i] <<
"s^" << i <<
" + ";
308 SBasis sin(Linear bo,
int k);
309 SBasis cos(Linear bo,
int k);
311 std::vector<double> roots(SBasis
const & s);
312 std::vector<std::vector<double> > multi_roots(SBasis
const &f,
313 std::vector<double>
const &levels,