src/m_cpoly.h

   1 /*$Id: m_cpoly.h,v 1.2 2009-12-10 14:34:44 felix Exp $ -*- C++ -*-
   2  * Copyright (C) 2001 Albert Davis
   3  * Author: Albert Davis <aldavis@gnu.org>
   4  *
   5  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  20  * 02110-1301, USA.
  21  *------------------------------------------------------------------
  22  * structs for fixed order polynomials, in 2 different forms
  23  * FPOLY is by function -- f0 = eval of function
  24  *                         f1 = 1st derivative
  25  *                         f2 = 2nd derivative
  26  *                         etc.
  27  *                      f(x) = f0
  28  *                      f(t) = f0 + f1*(t-x) + f2*(t-x)^2 + ...
  29  * CPOLY is by series -- c0 = coeff of x^0 term (constant)
  30  *                       c1 = coeff of x^1 term (1st derivative)
  31  *                       c2 = coeff of x^2 term
  32  *                       etc.
  33  *                      f(x) = c0 + f1*x + f2*x^2 + ...
  34  *                      f(t) = c0 + f1*t + f2*t^2 + ...
  35  */
  36 //testing=script,sparse 2006.07.13
  37 #ifndef M_CPOLY_H
  38 #define M_CPOLY_H
  39 #include "constant.h"
  40 #include "io_.h"
  41 /*--------------------------------------------------------------------------*/
  42 // fixme. move (where?)
  43 inline double operator*(const DPAIR& t, const DPAIR& a)
  44 { itested();
  45   double ret = t.first * a.first;
  46   ret += t.second * a.second;
  47   return ret;
  48 }
  49 /*--------------------------------------------------------------------------*/
  50 struct FPOLY1;
  51 struct CPOLY1;
  52 /*--------------------------------------------------------------------------*/
  53 struct FPOLY1{          /* first order polynomial       */
  54   hp_float_t   x;               /* the argument                 */
  55   hp_float_t   f0;              /* the function (c0 + x*f1)     */
  56   hp_float_t   f1;              /* the first derivative         */
  57   explicit FPOLY1() : x(0), f0(0), f1(0) {}
  58            FPOLY1(const FPOLY1& p) : x(p.x), f0(p.f0), f1(p.f1) {}
  59   explicit FPOLY1(hp_float_t X,hp_float_t F0,hp_float_t F1) : x(X), f0(F0), f1(F1) {}
  60   explicit FPOLY1(const CPOLY1& p);
  61   explicit FPOLY1(const DPAIR& l, const DPAIR& r);
  62   ~FPOLY1() {}
  63
  64   bool     operator==(const FPOLY1& p)const
  65                                 {return (f1==p.f1 && f0==p.f0 && x==p.x);}
  66   FPOLY1&  operator*=(hp_float_t s)     {f0*=s; f1*=s; return *this;}
  67   FPOLY1&  operator*=(const FPOLY1& s);
  68   FPOLY1&  operator+=(hp_float_t f) {untested(); f0+=f; return *this;}
  69   FPOLY1&  operator+=(const FPOLY1& s)
  70                 {untested(); assert(x==s.x); f0+=s.f0; f1+=s.f1; return *this;}
  71   FPOLY1   operator-()const     {untested(); return FPOLY1(x, -f0, -f1);}
  72   hp_float_t   c1()const                {assert(f1 == f1); return f1;}
  73   hp_float_t   c0()const
  74   {assert(f0==f0); assert(f1==f1); assert(x==x); assert(f0!=LINEAR); return (f0 - x * f1);}
  75   double f(double xx)const {return (f0 + (xx-x)*f1);}
  76   double operator()(double xx)const {return f(xx);}
  77 };
  78 /*--------------------------------------------------------------------------*/
  79 struct CPOLY1{          /* first order polynomial       */
  80   hp_float_t   x;               /* the argument                 */
  81   hp_float_t   c0;              /* f(x) - x*f'(x), or f0 - x*f1 */
  82   hp_float_t   c1;              /* the first derivative         */
  83   explicit CPOLY1() : x(0), c0(0), c1(0) {}
  84            CPOLY1(const CPOLY1& p) : x(p.x), c0(p.c0), c1(p.c1) {}
  85   explicit CPOLY1(hp_float_t X,hp_float_t C0,hp_float_t C1) : x(X), c0(C0), c1(C1) {}
  86   explicit CPOLY1(const FPOLY1& p);
  87   ~CPOLY1() {}
  88
  89   bool     operator==(const CPOLY1& p)const
  90                                 {return (c1==p.c1 && c0==p.c0 && x==p.x);}
  91   CPOLY1&  operator*=(hp_float_t s)     {c0*=s; c1*=s; return *this;}
  92   CPOLY1&  operator+=(const CPOLY1& s)  {c0+=s.c0; c1+=s.c1; return *this;}
  93   CPOLY1&  operator-=(const CPOLY1& s)  {c0-=s.c0; c1-=s.c1; return *this;}
  94   hp_float_t   f1()const                {return c1;}
  95   hp_float_t   f0()const                {return (c0 + x * c1);}
  96   hp_float_t   f(double xx)const                {return (c0 + xx * c1);}
  97   DPAIR intersect(const CPOLY1&) const;
  98   double zero()const;
  99   double zero(const FPOLY1&)const;
 100 };
 101 /*--------------------------------------------------------------------------*/
 102 inline FPOLY1::FPOLY1(const CPOLY1& p)
 103   :x(p.x),
 104    f0(p.f0()),
 105    f1(p.f1())
 106 {
 107   assert(p == p);
 108   assert(*this == *this);
 109 }
 110 /*--------------------------------------------------------------------------*/
 111 inline FPOLY1::FPOLY1(const DPAIR& l, const DPAIR& r):
 112   x(l.first),
 113   f0(l.second),
 114   f1((r.second - l.second ) / (r.first - l.first))
 115 {
 116 }
 117 /*--------------------------------------------------------------------------*/
 118   inline CPOLY1::CPOLY1(const FPOLY1& p)
 119 :x(p.x),
 120   c0(p.c0()),
 121   c1(p.c1())
 122 {
 123   assert(p == p);
 124   assert(x == x);
 125   assert(c1 == c1);
 126   assert(c0 == c0);
 127   assert(*this == *this);
 128 }
 129 /*--------------------------------------------------------------------------*/
 130 inline FPOLY1& FPOLY1::operator*=(const FPOLY1& s)
 131 {
 132   untested();
 133   assert(x == s.x);
 134   *this *= s.f0;
 135   f1 += f0 * s.f1;
 136   return *this;
 137 }
 138 /*--------------------------------------------------------------------------*/
 139 inline FPOLY1 operator*(FPOLY1 a, const FPOLY1& b)
 140 {
 141   untested();
 142   a *= b;
 143   return a;
 144 }
 145 /*--------------------------------------------------------------------------*/
 146 inline FPOLY1 operator+(FPOLY1 a, const FPOLY1& b)
 147 {
 148   untested();
 149   a += b;
 150   return a;
 151 }
 152 /*--------------------------------------------------------------------------*/
 153 inline FPOLY1 operator+(FPOLY1 a, hp_float_t b)
 154 {
 155   untested();
 156   a += b;
 157   return a;
 158 }
 159 /*--------------------------------------------------------------------------*/
 160 inline FPOLY1 operator-(hp_float_t a, const FPOLY1& b)
 161 {
 162   untested();
 163   return -b + a;
 164 }
 165 /*--------------------------------------------------------------------------*/
 166 // compute a zero. return inf if there is none.
 167 inline double CPOLY1::zero() const
 168 {
 169   if(c0==0){
 170     return 0;
 171   }else{
 172     return -c0/c1;
 173   }
 174 }
 175 /*--------------------------------------------------------------------------*/
 176   template<class T>
 177 inline T& operator<<(T& o, const CPOLY1& x)
 178 {
 179   return o << x.c0 << ((x.c1<0)?"":"+") << x.c1 << "*x at x=" << x.x;
 180 }
 181 /*--------------------------------------------------------------------------*/
 182   template<class T>
 183 inline T& operator<<(T& o, const FPOLY1& x)
 184 {
 185   return o << "F " << x.x << " " << x.f0 << " " << x.f1;
 186 }
 187 /*--------------------------------------------------------------------------*/
 188 // compute a zero wrt linear function. return inf if there is none.
 189 inline double CPOLY1::zero(const FPOLY1& ff) const
 190 {
 191   double p1 = ff.x;
 192   double p2 =   ff.f0;
 193   DPAIR az(1., ff.f1);
 194
 195   double t = (p1 - ff.f1*(c0-p2)) / (1.+ff.f1*c1);
 196
 197   trace5("zero", *this, ff, t, p1, p2);
 198   return t;
 199 }
 200 /*--------------------------------------------------------------------------*/
 201 inline DPAIR CPOLY1::intersect(const CPOLY1& that) const
 202 {
 203   CPOLY1 tmp(*this);
 204   tmp -= that;
 205   double z = tmp.zero();
 206   trace4("found zero", *this, that, z, tmp);
 207
 208   return DPAIR(z,f(z));
 209 }
 210 /*--------------------------------------------------------------------------*/
 211 #endif
 212 // vim:ts=8:sw=2:noet: