COCONUT API: prointerval.hpp Source File

00001 // Projective Interval Class -*- C++ -*-
00002 
00003 // $Id: prointerval.h 827 2008-11-26 23:28:47Z hermannb $
00004 // Copyright (C) 2001-2004 Hermann Schichl and Xuan-Ha VU
00005 //
00006 // This file is part of the COCONUT API.  This library
00007 // is free software; you can redistribute it and/or modify it under the
00008 // terms of the Library GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // Library GNU General Public License for more details.
00016 
00017 // As a special exception, you may use this file as part of a free software
00018 // library without restriction.  Specifically, if other files instantiate
00019 // templates or use macros or inline functions from this file, or you compile
00020 // this file and link it with other files to produce an executable, this
00021 // file does not by itself cause the resulting executable to be covered by
00022 // the Library GNU General Public License.  This exception does not however
00023 // invalidate any other reasons why the executable file might be covered by
00024 // the Library GNU General Public License.
00025 
00028 #ifndef _PROINTERVAL_HPP_
00029 #define _PROINTERVAL_HPP_
00030 
00031 #include <api_exception.h>
00032 #include <iostream>
00033 
00034 namespace coco {
00035 
00036 inline
00037 unsigned int proj_rational::scm(unsigned int& n, unsigned int& m) const
00038 {
00039   unsigned int g(gcd(n, m));
00040   unsigned int r(m);
00041   if(g != 1)
00042   {
00043     n/=g;
00044     m/=g;
00045   }
00046   return n*r;
00047 }
00048 
00049 inline
00050 unsigned int proj_rational::gcd(int n, unsigned int m) const
00051 {
00052   unsigned int a(abs(n));
00053   while(m != 0)
00054   {
00055     unsigned int h(a%m);
00056     a=m;
00057     m=h;
00058   }
00059   return a;
00060 }
00061 
00062 inline
00063 void proj_rational::cancel()
00064 {
00065   unsigned int g(gcd(nom, denom));
00066   if(g != 1)
00067   {
00068     nom /= g;
00069     denom /= g;
00070   }
00071 }
00072 
00073 inline
00074 proj_rational proj_rational::operator+(const proj_rational& R) const
00075 {
00076   unsigned int rR(R.denom);
00077   unsigned int r(denom);
00078   proj_rational ret;
00079   ret.denom = scm(r, rR);
00080   ret.nom = nom*(int)rR+R.nom*(int)r;
00081   return ret;
00082 }
00083 
00084 inline
00085 proj_rational proj_rational::operator-(const proj_rational& R) const
00086 {
00087   unsigned int rR(R.denom);
00088   unsigned int r(denom);
00089   proj_rational ret;
00090   ret.denom = scm(r, rR);
00091   ret.nom = nom*rR-R.nom*r;
00092   return ret;
00093 }
00094 
00095 inline
00096 proj_rational proj_rational::operator*(const proj_rational& R) const
00097 {
00098   int g1=gcd(nom, R.denom);
00099   int g2=gcd(R.nom, denom);
00100   proj_rational ret;
00101   ret.nom = nom/g1;
00102   ret.nom *= R.nom/g2;
00103   ret.denom = denom/g2;
00104   ret.denom *= R.denom/g1;
00105   return ret;
00106 }
00107 
00108 inline
00109 proj_rational proj_rational::operator/(const proj_rational& R) const
00110 {
00111   if(R.nom == 0)
00112     throw api_exception(apiee_internal,
00113         std::string("division by zero in proj_rational"));
00114   int g1=gcd(nom, abs(R.nom));
00115   int g2=gcd(R.denom, denom);
00116   proj_rational ret;
00117   ret.nom = nom/g1;
00118   ret.nom *= R.denom/g2;
00119   ret.denom = denom/g2;
00120   ret.denom *= abs(R.nom)/g1;
00121   if(R.nom<0) ret.nom = -ret.nom;
00122   return ret;
00123 }
00124 
00125 inline
00126 proj_rational& proj_rational::operator+=(const proj_rational& R)
00127 {
00128   unsigned int rR(R.denom);
00129   unsigned int r(denom);
00130   denom = scm(r, rR);
00131   nom *= rR;
00132   nom += R.nom*r;
00133   return *this;
00134 }
00135 
00136 inline
00137 proj_rational& proj_rational::operator-=(const proj_rational& R)
00138 {
00139   unsigned int rR(R.denom);
00140   unsigned int r(denom);
00141   denom = scm(r, rR);
00142   nom *= rR;
00143   nom -= R.nom*r;
00144   return *this;
00145 }
00146 
00147 inline
00148 proj_rational& proj_rational::operator*=(const proj_rational& R)
00149 {
00150   int g1=gcd(nom, R.denom);
00151   int g2=gcd(R.nom, denom);
00152   proj_rational ret;
00153   nom /= g1;
00154   nom *= R.nom/g2;
00155   denom /= g2;
00156   denom *= R.denom/g1;
00157   return *this;
00158 }
00159 
00160 inline
00161 proj_rational& proj_rational::operator/=(const proj_rational& R)
00162 {
00163   if(R.nom == 0)
00164     throw api_exception(apiee_internal,
00165         std::string("division by zero in proj_rational"));
00166   int g1=gcd(nom, abs(R.nom));
00167   int g2=gcd(R.denom, denom);
00168   nom /= g1;
00169   nom *= R.denom/g2;
00170   denom /= g2;
00171   denom *= abs(R.nom)/g1;
00172   if(R.nom<0) nom = -nom;
00173   return *this;
00174 }
00175 
00176 } // namespace coco
00177 
00178 namespace std {
00179 //inline
00180 ostream& operator<<(ostream& o, const coco::proj_rational& x)
00181 {
00182   o << x.nominator();
00183   if(x.denominator() != 1)
00184     o << '/' << x.denominator();
00185   return o;
00186 }
00187 
00188 template <class I>
00189 ostream& operator<<(ostream& o, const coco::projective_interval<I>& a)
00190 {
00191   o << '[' << a.pi();
00192   if(a.pr() != 0)
00193     o << '/' << a.pt() << "^(" << a.pr() << ")";
00194   return o << "]";
00195 }
00196 
00197 } // namespace std
00198 
00199 namespace coco {
00200 
00202 template <class I>
00203 void projective_interval<I>::shift(const proj_rational& s)
00204 {
00205   if(s != this->r)
00206   {
00207     this->i *= pow(*this->t,(s-this->r).ival());
00208     this->r = s;
00209   }
00210 }
00211   
00212 template <class I>
00213 typename projective_interval<I>::interval
00214         projective_interval<I>::shiftr(const proj_rational& s) const
00215 {
00216   interval ret(this->i);
00217   if(s != this->r)
00218     ret *= pow(*this->t,(s-this->r).ival());
00219   return ret;
00220 }
00221   
00222 template <class I>
00223 typename projective_interval<I>::interval
00224         projective_interval<I>::shiftr(const proj_interval& x,
00225                                                 const proj_rational& s) const
00226 {
00227   interval ret(x.i);
00228   if(s != x.r)
00229     ret *= pow(*this->t,(s-x.r).ival());
00230   return ret;
00231 }
00232   
00233 template <class I>
00234 typename projective_interval<I>::interval
00235         projective_interval<I>::shiftintv(const proj_rational& r,
00236                                            const proj_rational& s) const
00237 {
00238   if(s != r)
00239     return pow(*this->t,(r-s).ival());
00240   else
00241     return interval(1.);
00242 }
00243   
00244 template <class I>
00245 //i.set(a,b) sets i to [a, b]
00246 void projective_interval<I>::set(double lo, double up, const proj_interval& _p)
00247 {
00248   this->i = interval(lo, up);
00249   this->r = 0;
00250   this->t = _p.t;
00251 }
00252 
00253 template <class I>
00254 //i.set(a,b) sets i to [a, b]
00255 void projective_interval<I>::set_i(double lo, double up)
00256 {
00257   this->i.set(lo, up);
00258 }
00259 
00260 template <class I>
00261 //i.set(a,b) sets i to [a, b]
00262 void projective_interval<I>::set_i(const interval& j)
00263 {
00264   this->i = j;
00265 }
00266 
00267 template <class I>
00268 //i.set(a,b) sets i to [a, b]
00269 void projective_interval<I>::set_t(double lo, double up)
00270 {
00271   this->t->set(lo, up);
00272 }
00273 
00274 template <class I>
00275 //i.set(a,b) sets i to [a, b]
00276 void projective_interval<I>::set_t(const interval& u)
00277 {
00278   *this->t = u;
00279 }
00280 
00281 template <class I>
00282 //i.set(a,b) sets i to [a, b]
00283 void projective_interval<I>::set_r(int n, unsigned int m)
00284 {
00285   this->r.set(n,m);
00286 }
00287 
00288 template <class I>
00289 //i.set(a,b) sets i to [a, b]
00290 void projective_interval<I>::set_r(const proj_rational& s)
00291 {
00292   this->r = s;
00293 }
00294 
00295 template <class I>
00296 bool projective_interval<I>::certain_subset (const proj_interval& x) const    // certain subset
00297 {
00298   if(this->r == x.r)                // t powers coincide
00299     return this->i.subset(x.i);
00300   else
00301     return shiftr(x.r).subset(x.i);
00302 }
00303 
00304 template <class I>
00305 bool projective_interval<I>::possible_subset (const proj_interval& x) const   // possible subset
00306 {
00307   if(this->r == x.r)                // t powers coincide
00308     return this->i.subset(x.i);
00309   else
00310     return this->i.subset(shiftr(x,this->r));
00311 }
00312 
00313 template <class I>
00314 bool projective_interval<I>::certain_proper_subset (const proj_interval& x) const //certain proper subset
00315 {
00316   if(this->r == x.r)                // t powers coincide
00317     return this->i.proper_subset(x.i);
00318   else
00319     return shiftr(x.r).proper_subset(x.i);
00320 }
00321 
00322 template <class I>
00323 bool projective_interval<I>::possible_proper_subset (const proj_interval& x) const//possible proper subset
00324 {
00325   if(this->r == x.r)                // t powers coincide
00326     return this->i.proper_subset(x.i);
00327   else
00328     return this->i.proper_subset(shiftr(x,this->r));
00329 }
00330 
00331 template <class I>
00332 bool projective_interval<I>::disjoint (const proj_interval& x) const 
00333 {
00334   if(this->r == x.r)                // t powers coincide
00335     return this->i.disjoint(x.i);
00336   else
00337   {
00338     proj_interval hlpr(this->r-x.r);
00339     if(this->r < x.r)               // x has higher t power
00340     {
00341       shift(x.r);
00342       return this->i.disjoint(x.i);
00343     }
00344     else                       // x has lower t power
00345     {
00346       return this->i.disjoint(shiftr(x,this->r));
00347     }
00348   }
00349 }
00350 
00351 // y.intersect(x) intersects y with x
00352 template <class I>
00353 projective_interval<I>&
00354         projective_interval<I>::intersectwith(const proj_interval& x)
00355 {
00356   if(this->r == x.r)                // t powers coincide
00357     this->i.intersect(x.i);
00358   else
00359   {
00360     if(this->r < x.r)               // x has higher t power
00361     {
00362       shift(x.r);
00363       this->i.intersect(x.i);
00364     }
00365     else                       // x has lower t power
00366       this->i.intersect(shiftr(x,this->r));
00367   }
00368   return *this;
00369 }
00370 
00371 // y.hull(x) builds the hull of y with x
00372 template <class I>
00373 projective_interval<I>& projective_interval<I>::hull(const proj_interval& x)
00374 {
00375   if(this->r == x.r)                // t powers coincide
00376     this->i.hull(x.i);
00377   else
00378   {
00379     if(this->r < x.r)               // x has higher t power
00380     {
00381       shift(x.r);
00382       this->i.hull(x.i);
00383     }
00384     else                       // x has lower t power
00385       this->i.hull(shiftr(x,this->r));
00386   }
00387   return *this;
00388 }
00389 
00390 template <class I>
00391 projective_interval<I>& projective_interval<I>::hulldiff(const proj_interval& x)
00392 {
00393   if(this->r == x.r)                // t powers coincide
00394     this->i.hulldiff(x.i);
00395   else
00396   {
00397     this->i*=shiftintv(x.r,this->r);
00398     this->i.hulldiff(x.i);
00399   }
00400   return *this;
00401 }
00402 
00404 template <class I>
00405 projective_interval<I>& projective_interval<I>::ipow(int n)
00406 {
00407   this->i.ipow(n);
00408   this->r *= n;
00409   return *this;
00410 }
00412 template <class I>
00413 projective_interval<I>& projective_interval<I>::imin(const proj_interval& x)
00414 {
00415   if(this->r == x.r)                // t powers coincide
00416     this->i.imin(x.i);
00417   else
00418   {
00419     if(this->r < x.r)               // x has higher t power
00420     {
00421       shift(x.r);
00422       this->i.imin(x.i);
00423     }
00424     else                       // x has lower t power
00425       this->i.imin(shiftr(x,this->r));
00426   }
00427   return *this;
00428 }
00429 
00431 template <class I>
00432 projective_interval<I>& projective_interval<I>::imax(const proj_interval& x)
00433 {
00434   if(this->r == x.r)                // t powers coincide
00435     this->i.imax(x.i);
00436   else
00437   {
00438     if(this->r < x.r)               // x has higher t power
00439     {
00440       shift(x.r);
00441       this->i.imax(x.i);
00442     }
00443     else                       // x has lower t power
00444       this->i.imax(shiftr(x,this->r));
00445   }
00446   return *this;
00447 }
00448 
00450 template <class I>
00451 projective_interval<I>& projective_interval<I>::round_to_integer()
00452 {
00453   // this is impossible for projective intervals, except for r==0;
00454   if(this->r == 0)
00455     this->i.round_to_integer();
00456   return *this;
00457 }
00458 
00459 template <class I>
00460 projective_interval<I>& projective_interval<I>::intersect_div(const proj_interval& __i,
00461                                                      const proj_interval& __j)
00462 {
00463   projective_interval<I> hlp(__i);
00464   hlp.shift(this->r+__j.r);
00465   this->i.intersect_div(hlp.i, __j.i);
00466   return *this;
00467 }
00468 
00469 template <class I>
00470 projective_interval<I>& projective_interval<I>::intersect_power(const proj_interval& __i,
00471                                                        int __n)
00472 {
00473   interval zi(this->i);
00474   proj_rational ir(__i.r*__n);
00475   if(ir != this->r);
00476   zi *= shiftintv(ir, this->r);
00477   zi.intersect_power(__i.i, __n);
00478   this->i = zi;
00479   return *this;
00480 }
00481 
00482 template <class I>
00483 projective_interval<I>& projective_interval<I>::intersect_invsqr_wc(const proj_interval& __i,
00484                                                            double __l, double __d)
00485 {
00486   //TODO
00487   api_exception(apiee_nyi, "projective_interval intersect_invsqr_wc is not yet implemented!");
00488   return *this;
00489 }
00490 
00491 template <class I>
00492 projective_interval<I>& projective_interval<I>::intersect_invpower_wc(const proj_interval& __i,
00493                                                              double __l, int __n)
00494 {
00495   //TODO
00496   api_exception(apiee_nyi, "projective_interval intersect_invpower_wc is not yet implemented!");
00497   return *this;
00498 }
00499 
00500 template <class I>
00501 projective_interval<I>& projective_interval<I>::intersect_invsin_wc(const proj_interval& __i,
00502                                                            double __l, double __d)
00503 {
00504   // TODO: This is suboptimal!
00505   this->t.intersectwith(__i.t);
00506   projective_interval<I> ret(__i);
00507   ret.t = this->t;
00508   ret.shift(0);
00509   proj_rational rs(this->r);
00510   shift(0);
00511   this->i.intersect_invsin_wc(ret.i, __l, __d);
00512   shift(rs);
00513   return *this;
00514 }
00515 
00516 template <class I>
00517 projective_interval<I>& projective_interval<I>::intersect_invcos_wc(const proj_interval& __i,
00518                                                            double __l, double __d)
00519 {
00520   // TODO: This is suboptimal!
00521   projective_interval<I> ret(__i);
00522   ret.shift(0);
00523   proj_rational rs(this->r);
00524   shift(0);
00525   this->i.intersect_invcos_wc(ret.i, __l, __d);
00526   shift(rs);
00527   return *this;
00528 }
00529 
00530 template <class I>
00531 projective_interval<I>& projective_interval<I>::intersect_invgauss_wc(const proj_interval& __i,
00532                                                 double __l, double __m, double __s)
00533 {
00534   // TODO: This is suboptimal!
00535   projective_interval<I> ret(__i);
00536   ret.shift(0);
00537   proj_rational rs(this->r);
00538   shift(0);
00539   this->i.intersect_invgauss_wc(ret.i, __l, __d, __s);
00540   shift(rs);
00541   return *this;
00542 }
00543 
00544 template <class I>
00545 double projective_interval<I>::project(double __d) const
00546 {
00547   if(this->r > 0)
00548     return __d;
00549   else if(this->r < 0)
00550   {
00551     interval hlp(this->i*pow(this->t,(-this->r).ival()));
00552     return hlp.project(__d);
00553   }
00554   else
00555     return this->i.project(__d);
00556 }
00557 
00558 template <class I>
00559 projective_interval<I>& projective_interval<I>::operator+=(const proj_interval& a)
00560 {
00561   if(this->r == a.r)
00562     this->i += a.i;
00563   else
00564   {
00565     if(this->r < a.r)
00566     {
00567       shift(a.r);
00568       this->i += a.i;
00569     }
00570     else
00571       this->i += a.i*shiftintv(this->r,a.r);
00572   }
00573   return *this;
00574 }
00575 
00576 template <class I>
00577 projective_interval<I>& projective_interval<I>::operator+=(double a)
00578 {
00579   if(this->r == 0)
00580     this->i += a;
00581   else
00582   {
00583     if(this->r < 0)
00584     {
00585       shift(0);
00586       this->i += a;
00587     }
00588     else
00589       this->i += a*shiftintv(this->r,0.);
00590   }
00591   return *this;
00592 }
00593 
00594 template <class I>
00595 projective_interval<I>& projective_interval<I>::operator-=(const proj_interval& a)
00596 {
00597   if(this->r == a.r)
00598     this->i -= a.i;
00599   else
00600   {
00601     if(this->r < a.r)
00602     {
00603       shift(a.r);
00604       this->i -= a.i;
00605     }
00606     else
00607       this->i -= a.i*shiftintv(this->r,a.r);
00608   }
00609   return *this;
00610 }
00611 
00612 template <class I>
00613 projective_interval<I>& projective_interval<I>::operator-=(double a)
00614 {
00615   if(this->r == 0)
00616     this->i -= a;
00617   else
00618   {
00619     if(this->r < 0)
00620     {
00621       shift(0);
00622       this->i -= a;
00623     }
00624     else
00625       this->i -= a*shiftintv(this->r,0);
00626   }
00627   return *this;
00628 }
00629 
00630 template <class I>
00631 projective_interval<I>& projective_interval<I>::operator*=(const proj_interval& a)
00632 {
00633   this->i *= a.i;
00634   this->r += a.r;
00635   return *this;
00636 }
00637 
00638 template <class I>
00639 projective_interval<I>& projective_interval<I>::operator*=(double a)
00640 {
00641   this->i *= a;
00642   return *this;
00643 }
00644 
00645 template <class I>
00646 double projective_interval<I>::mid() const
00647 {
00648   return this->i.mid();
00649 }
00650 
00651 template <class I>
00652 double projective_interval<I>::diam() const
00653 {
00654   return this->i.diam();
00655 
00656 }
00657 
00658 template <class I>
00659 double projective_interval<I>::width() const
00660 {
00661   return this->i.width();
00662 }
00663 
00664 template <class I>
00665 double projective_interval<I>::relDiam() const
00666 {
00667   return this->i.relDiam();
00668 }
00669 
00670 template <class I>
00671 double projective_interval<I>::rad() const
00672 {
00673   return this->i.rad();
00674 }
00675 
00676 template <class I>
00677 double projective_interval<I>::mig() const
00678 {
00679   return this->i.mig();
00680 }
00681 
00682 template <class I>
00683 double projective_interval<I>::mag() const
00684 {
00685   return this->i.mag();
00686 }
00687 
00688 template <class I>
00689 double projective_interval<I>::dist(const proj_interval& x) const
00690 {
00691   // TODO
00692   api_exception(apiee_nyi, "projective_interval dist is not yet implemented!");
00693   return 0;
00694 }
00695 
00696 template <class I>
00697 double projective_interval<I>::safeguarded_mid() const
00698 {
00699   return this->i.safeguarded_mid();
00700 }
00701 
00702 template <class I>
00703 projective_interval<I> operator+(const projective_interval<I>& a,
00704                                  const projective_interval<I>& b)
00705 {
00706   projective_interval<I> ret(a);
00707   if(ret.r == b.r)
00708     ret.i += b.i;
00709   else
00710   {
00711     if(ret.r < b.r)
00712     {
00713       ret.shift(b.r);
00714       ret.i += b.i;
00715     }
00716     else
00717       ret.i += b.i*ret.shiftintv(ret.r,b.r);
00718   }
00719   return ret;
00720 }
00721 
00722 template <class I>
00723 projective_interval<I> operator+(const projective_interval<I>& a, double b)
00724 {
00725   projective_interval<I> ret(a);
00726   if(ret.r == 0)
00727     ret.i += b;
00728   else
00729   {
00730     if(ret.r < 0)
00731     {
00732       ret.shift(0);
00733       ret.i += b;
00734     }
00735     else
00736       ret.i += b*ret.shiftintv(ret.r,0);
00737   }
00738   return ret;
00739 }
00740 
00741 template <class I>
00742 projective_interval<I> operator+(double b, const projective_interval<I>& a)
00743 {
00744   projective_interval<I> ret(a);
00745   if(ret.r == 0)
00746     ret.i += b;
00747   else
00748   {
00749     if(ret.r < 0)
00750     {
00751       ret.shift(0);
00752       ret.i += b;
00753     }
00754     else
00755       ret.i += b*ret.shiftintv(ret.r,0);
00756   }
00757   return ret;
00758 }
00759 
00760 template <class I>
00761 projective_interval<I> operator-(const projective_interval<I>& a,
00762                                  const projective_interval<I>& b)
00763 {
00764   projective_interval<I> ret(a);
00765   if(ret.r == b.r)
00766     ret.i -= b.i;
00767   else
00768   {
00769     if(ret.r < b.r)
00770     {
00771       ret.shift(b.r);
00772       ret.i -= b.i;
00773     }
00774     else
00775       ret.i -= b.i*ret.shiftintv(ret.r,b.r);
00776   }
00777   return ret;
00778 }
00779 
00780 template <class I>
00781 projective_interval<I> operator-(const projective_interval<I>& a, double b)
00782 {
00783   projective_interval<I> ret(a);
00784   if(ret.r == 0)
00785     ret.i -= b;
00786   else
00787   {
00788     if(ret.r < 0)
00789     {
00790       shift(0);
00791       ret.i -= b;
00792     }
00793     else
00794       ret.i -= b*ret.shiftintv(ret.r,0);
00795   }
00796   return ret;
00797 }
00798 
00799 template <class I>
00800 projective_interval<I> operator-(double b, const projective_interval<I>& a)
00801 {
00802   projective_interval<I> ret(a);
00803   if(ret.r == 0)
00804     ret.i -= b;
00805   else
00806   {
00807     if(ret.r < 0)
00808     {
00809       shift(0);
00810       ret.i = b-ret.i;
00811     }
00812     else
00813       ret.i = b*ret.shiftintv(ret.r,0) - ret.i;
00814   }
00815   return ret;
00816 }
00817 
00818 template <class I>
00819 projective_interval<I> operator*(const projective_interval<I>& a,
00820                                  const projective_interval<I>& b)
00821 {
00822   projective_interval<I> ret(a);
00823   ret.i *= b.i;
00824   ret.r += b.r;
00825   return ret;
00826 }
00827 
00828 template <class I>
00829 projective_interval<I> operator*(const projective_interval<I>& a, double b)
00830 {
00831   projective_interval<I> ret(a);
00832   ret.i *= b;
00833   return ret;
00834 }
00835 
00836 template <class I>
00837 projective_interval<I> operator*(double b, const projective_interval<I>& a)
00838 {
00839   projective_interval<I> ret(a);
00840   ret.i *= b;
00841   return ret;
00842 }
00843 
00844 template <class I>
00845 projective_interval<I> operator/(const projective_interval<I>& a,
00846                                  const projective_interval<I>& b)
00847 {
00848   projective_interval<I> ret(a);
00849   ret.i /= b.i;
00850   ret.r -= b.r;
00851   return ret;
00852 }
00853 
00854 template <class I>
00855 projective_interval<I> operator/(const projective_interval<I>& a, double b)
00856 {
00857   projective_interval<I> ret(a);
00858   ret.i /= b;
00859   return ret;
00860 }
00861 
00862 template <class I>
00863 projective_interval<I> operator/(double b, const projective_interval<I>& a)
00864 {
00865   projective_interval<I> ret(a);
00866   ret.i /= b;
00867   return ret;
00868 }
00869 
00870 template <class I>
00871 projective_interval<I> acos(const projective_interval<I>& x)
00872 {
00873   projective_interval<I> ret(x);
00874   ret.shift(0);
00875   ret.i = acos(ret.i);
00876   return ret;
00877 }
00878 
00879 template <class I>
00880 projective_interval<I> abs(const projective_interval<I>& x)
00881 {
00882   projective_interval<I> ret(x);
00883   ret.i = abs(ret.i);
00884   return ret;
00885 }
00886 
00887 template <class I>
00888 projective_interval<I> acosh(const projective_interval<I>& x)
00889 {
00890   projective_interval<I> ret(x);
00891   ret.shift(0);
00892   ret.i = acosh(ret.i);
00893   return ret;
00894 }
00895 
00896 template <class I>
00897 projective_interval<I> acot(const projective_interval<I>& x)
00898 {
00899   projective_interval<I> ret(x);
00900   ret.shift(0);
00901   ret.i = acot(ret.i);
00902   return ret;
00903 }
00904 
00905 template <class I>
00906 projective_interval<I> acoth(const projective_interval<I>& x)
00907 {
00908   projective_interval<I> ret(x);
00909   ret.shift(0);
00910   ret.i = acoth(ret.i);
00911   return ret;
00912 }
00913 
00914 template <class I>
00915 projective_interval<I> asin(const projective_interval<I>& x)
00916 {
00917   projective_interval<I> ret(x);
00918   ret.shift(0);
00919   ret.i = asin(ret.i);
00920   return ret;
00921 }
00922 
00923 template <class I>
00924 projective_interval<I> asinh(const projective_interval<I>& x)
00925 {
00926   projective_interval<I> ret(x);
00927   ret.shift(0);
00928   ret.i = asinh(ret.i);
00929   return ret;
00930 }
00931 
00932 template <class I>
00933 projective_interval<I> atan(const projective_interval<I>& x)
00934 {
00935   projective_interval<I> ret(x);
00936   ret.shift(0);
00937   ret.i = atan(ret.i);
00938   return ret;
00939 }
00940 
00941 template <class I>
00942 projective_interval<I> atanh(const projective_interval<I>& x)
00943 {
00944   projective_interval<I> ret(x);
00945   ret.shift(0);
00946   ret.i = atanh(ret.i);
00947   return ret;
00948 }
00949 
00950 template <class I>
00951 projective_interval<I> cos(const projective_interval<I>& x)
00952 {
00953   projective_interval<I> ret(x);
00954   ret.shift(0);
00955   ret.i = cos(ret.i);
00956   return ret;
00957 }
00958 
00959 template <class I>
00960 projective_interval<I> cosh(const projective_interval<I>& x)
00961 {
00962   projective_interval<I> ret(x);
00963   ret.shift(0);
00964   ret.i = cosh(ret.i);
00965   return ret;
00966 }
00967 
00968 template <class I>
00969 projective_interval<I> cot(const projective_interval<I>& x)
00970 {
00971   projective_interval<I> ret(x);
00972   ret.shift(0);
00973   ret.i = cot(ret.i);
00974   return ret;
00975 }
00976 
00977 template <class I>
00978 projective_interval<I> coth(const projective_interval<I>& x)
00979 {
00980   projective_interval<I> ret(x);
00981   ret.shift(0);
00982   ret.i = coth(ret.i);
00983   return ret;
00984 }
00985 
00986 template <class I>
00987 projective_interval<I> exp(const projective_interval<I>& x)
00988 {
00989   projective_interval<I> ret(x);
00990   ret.shift(0);
00991   ret.i = exp(ret.i);
00992   return ret;
00993 }
00994 
00995 template <class I>
00996 projective_interval<I> exp10(const projective_interval<I>& x)
00997 {
00998   projective_interval<I> ret(x);
00999   ret.shift(0);
01000   ret.i = exp10(ret.i);
01001   return ret;
01002 }
01003 
01004 template <class I>
01005 projective_interval<I> exp2(const projective_interval<I>& x)
01006 {
01007   projective_interval<I> ret(x);
01008   ret.shift(0);
01009   ret.i = exp2(ret.i);
01010   return ret;
01011 }
01012 
01013 template <class I>
01014 projective_interval<I> expm1(const projective_interval<I>& x)
01015 {
01016   projective_interval<I> ret(x);
01017   ret.shift(0);
01018   ret.i = expm1(ret.i);
01019   return ret;
01020 }
01021 
01022 template <class I>
01023 projective_interval<I> log(const projective_interval<I>& x)
01024 {
01025   projective_interval<I> ret(x);
01026   ret.shift(0);
01027   ret.i = log(ret.i);
01028   return ret;
01029 }
01030 
01031 template <class I>
01032 projective_interval<I> log10(const projective_interval<I>& x)
01033 {
01034   projective_interval<I> ret(x);
01035   ret.shift(0);
01036   ret.i = log10(ret.i);
01037   return ret;
01038 }
01039 
01040 template <class I>
01041 projective_interval<I> log1p(const projective_interval<I>& x)
01042 {
01043   projective_interval<I> ret(x);
01044   ret.shift(0);
01045   ret.i = log1p(ret.i);
01046   return ret;
01047 }
01048 
01049 template <class I>
01050 projective_interval<I> log2(const projective_interval<I>& x)
01051 {
01052   projective_interval<I> ret(x);
01053   ret.shift(0);
01054   ret.i = log2(ret.i);
01055   return ret;
01056 }
01057 
01058 template <class I>
01059 projective_interval<I> power(const projective_interval<I>& x, int n)
01060 {
01061   projective_interval<I> ret(x);
01062   ret.i.ipow(n);
01063   ret.r *= n;
01064   return ret;
01065 }
01066 
01067 template <class I>
01068 projective_interval<I> pow(const projective_interval<I>& x,
01069                            const projective_interval<I>& y)
01070 {
01071   projective_interval<I> ret(x);
01072   projective_interval<I> hlp(y);
01073   ret.shift(0);
01074   hlp.shift(0);
01075   ret.i = pow(ret.i,hlp.i);
01076   return ret;
01077 }
01078 
01079 template <class I>
01080 projective_interval<I> sin(const projective_interval<I>& x)
01081 {
01082   projective_interval<I> ret(x);
01083   ret.shift(0);
01084   ret.i = sin(ret.i);
01085   return ret;
01086 }
01087 
01088 template <class I>
01089 projective_interval<I> sinh(const projective_interval<I>& x)
01090 {
01091   projective_interval<I> ret(x);
01092   ret.shift(0);
01093   ret.i = sinh(ret.i);
01094   return ret;
01095 }
01096 
01097 template <class I>
01098 projective_interval<I> sqr(const projective_interval<I>& x)
01099 {
01100   projective_interval<I> ret(x);
01101   ret.i = sqr(ret.i);
01102   ret.r *= 2;
01103   return ret;
01104 }
01105 
01106 template <class I>
01107 projective_interval<I> sqrt(const projective_interval<I>& x)
01108 {
01109   projective_interval<I> ret(x);
01110   ret.i = sqrt(ret.i);
01111   ret.r /= 2;
01112   return ret;
01113 }
01114 
01115 template <class I>
01116 projective_interval<I> tan(const projective_interval<I>& x)
01117 {
01118   projective_interval<I> ret(x);
01119   ret.shift(0);
01120   ret.i = tan(ret.i);
01121   return ret;
01122 }
01123 
01124 template <class I>
01125 projective_interval<I> tanh(const projective_interval<I>& x)
01126 {
01127   projective_interval<I> ret(x);
01128   ret.shift(0);
01129   ret.i = tanh(ret.i);
01130   return ret;
01131 }
01132 
01133 template <class I>
01134 projective_interval<I> imax(const projective_interval<I>& x,
01135                             const projective_interval<I>& y)
01136 {
01137   projective_interval<I> ret(x);
01138   ret.imax(y);
01139   return ret;
01140 }
01141 
01142 template <class I>
01143 projective_interval<I> imin(const projective_interval<I>& x,
01144                             const projective_interval<I>& y)
01145 {
01146   projective_interval<I> ret(x);
01147   ret.imin(y);
01148   return ret;
01149 }
01150 
01151 template <class I>
01152 projective_interval<I> division_part1(const projective_interval<I>& x,
01153       const projective_interval<I>& y, bool& b)
01154 {
01155   projective_interval<I> ret(x);
01156   division_part1(ret.i, y.i, b);
01157   ret.r -= b.r;
01158   return ret;
01159 }
01160 
01161 template <class I>
01162 projective_interval<I> division_part2(const projective_interval<I>& x,
01163       const projective_interval<I>& y, bool b)
01164 {
01165   projective_interval<I> ret(x);
01166   division_part2(ret.i, y.i, b);
01167   ret.r -= b.r;
01168   return ret;
01169 }
01170 
01171 } // namespace coco
01172 
01173 #endif // _PROINTERVAL_HPP_