Real algebraic numbers

00001 /**********************************************************************
00002  *   This file is part of the source code of the SYNAPS library.      *
00003  *   Author(s): B. Mourrain, GALAAD, INRIA                            *
00004  *   $Id: Dichotomy.h,v 1.1 2005/07/11 10:03:50 jppavone Exp $
00005  **********************************************************************/
00006 #ifndef SYNAPS_SOLVE_DICHOTOMY_H
00007 #define SYNAPS_SOLVE_DICHOTOMY_H
00008 
00009 #include <synaps/init.h>
00010 #include <synaps/arithm/sign.h>
00011 
00012 
00013 __BEGIN_NAMESPACE_SYNAPS
00014 //----------------------------------------------------------------------
00018 struct Dichotomy {
00019   double eps;
00020   Dichotomy():eps(1.e-6){}
00021   Dichotomy(double x):eps(x){}
00022 };
00023 //----------------------------------------------------------------------
00028 template<class T,class C>
00029 C RootOf(const T &p, const C & a, const C & b, const Dichotomy & Mth)
00030 {
00031         register C u0=a, u1=b, um;
00032         int s0=sign(p(a)), s1=sign(p(b));
00033         while((u1-u0)>Mth.eps && s0 != s1)
00034         {
00035                 um = C((u0+u1)/2);
00036                 int sm= sign(p(um));
00037                 if(sm != s0)
00038                 { s1=sm; u1=um; }
00039                 else if(sm != s1)
00040                 { s0=sm; u0=um; }
00041                 else if(sm!=0)
00042                 {
00043                         std::cerr <<"Error in sign evaluation"<< std::endl;
00044                         return 100000;//          exit(1);
00045                 }
00046         }
00047         if ((s0 == s1) && (a != b))
00048                 std::cerr <<"No opposite sign at "<<a<<" and "<<b<<std::endl;
00049         return (u0+u1)/2;
00050 }
00051 
00052 __END_NAMESPACE_SYNAPS
00053 
00054 /********************************************************************/
00055 #endif // SYNAPS_SOLVE_DICHOTOMY_H
synaps/usolve/Dichotomy.h
