continewz: include/continewz/homotopy.hpp Source File

00001 
00002 /******************************************************************************
00003 * MODULE     : homotopy.hpp
00004 * DESCRIPTION: Homotopies
00005 * COPYRIGHT  : (C) 2009  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX_HOMOTOPY_HPP
00014 #define __MMX_HOMOTOPY_HPP
00015 #include <multimix/slp_multivariate.hpp>
00016 
00017 namespace mmx {
00018 #define Homotopy_variant(C) typename homotopy_variant_helper<C>::HV
00019 #define TMPL template<typename C, typename V>
00020 #define Stepper homotopy_stepper<C,V>
00021 #define Stepper_rep homotopy_stepper_rep<C,V>
00022 #define Function slp_tangent<C>
00023 
00024 struct homotopy_abstract {};
00025 template<typename C> struct homotopy_variant_helper;
00026 template<typename C> nat start_index (const vector<C>& v, const C& z);
00027 
00028 /******************************************************************************
00029 * Prototype for homotopy steppers
00030 ******************************************************************************/
00031 
00032 TMPL
00033 struct homotopy_stepper_rep REP_STRUCT {
00034   Function f;
00035 
00036 public:
00037   inline homotopy_stepper_rep (const Function& f2): f (f2) {}
00038 
00039   virtual vector<C>
00040   hom (const vector<C>& init, nat k, const C& t1) = 0;
00041 
00042   virtual vector<C>
00043   homp (const vector<C>& init, nat k, const vector<C>& p) {
00044     //mmerr << "From x = " << init << " via x[" << k << "] = " << p
00045     //<< " [base stepper]" << "\n";
00046     typedef Abs_type(C) R;
00047     R eps= sqrt (sqrt (Accuracy (R)));
00048     vector<C> v= init;
00049     for (nat i= start_index (p, init[k]); i<N(p); i++) {
00050       v= hom (v, k, p[i]);
00051       if (abs (v[k] - p[i]) >= eps) break;
00052     }
00053     return v;
00054   }
00055 
00056   virtual vector<vector<C> >
00057   homv (const vector<vector<C> >& init, nat k, const C& t1) {
00058     //mmerr << "From x = " << init << " to x[" << k << "] = " << t1
00059     //<< " [base stepper]" << "\n";
00060     nat n= N(init);
00061     if (n == 0) return init;
00062     for (nat i=1; i<n; i++)
00063       if (init[i][k] != init[0][k]) {
00064         vector<vector<C> > fin= fill<vector<C> > (n, CF(init));
00065         for (nat j=0; j<n; j++)
00066           fin[j]= hom (init[j], k, t1);
00067         return fin;
00068       }
00069     return homvp (init, k, vec<C> (init[0][k], t1));
00070   }
00071 
00072   virtual vector<vector<C> >
00073   homvp (const vector<vector<C> >& init, nat k, const vector<C>& p) {
00074     //mmerr << "From x = " << init << " via x[" << k << "] = " << p
00075     //<< " [base stepper]" << "\n";
00076     nat n= N(init);
00077     vector<vector<C> > fin= fill<vector<C> > (n, CF(init));
00078     for (nat i=0; i<n; i++)
00079       fin[i]= homp (init[i], k, p);
00080     return fin;
00081   }
00082 };
00083 
00084 TMPL
00085 class homotopy_stepper {
00086 INDIRECT_PROTO_2 (homotopy_stepper, homotopy_stepper_rep, C, V)
00087 public:
00088   inline homotopy_stepper (const Function& f):
00089     rep (new Stepper_rep (f)) {}
00090 };
00091 INDIRECT_IMPL_2 (homotopy_stepper, homotopy_stepper_rep,
00092                  typename C, C, typename V, V)
00093 
00094 /******************************************************************************
00095 * Internal interface and useful subroutines
00096 ******************************************************************************/
00097 
00098 template<typename C, typename V> vector<C>
00099 homx (const Stepper& stepper,
00100       const vector<C>& x, nat k, const C& t1) {
00101   return const_cast<Stepper_rep*> (inside (stepper)) -> hom (x, k, t1);
00102 }
00103 
00104 template<typename C, typename V> vector<C>
00105 homx (const Stepper& stepper,
00106       const vector<C>& x, nat k, const vector<C>& p) {
00107   return const_cast<Stepper_rep*> (inside (stepper)) -> homp (x, k, p);
00108 }
00109 
00110 template<typename C, typename V> vector<vector<C> >
00111 homx (const Stepper& stepper,
00112       const vector<vector<C> >& x, nat k, const C& t1)
00113 {
00114   return const_cast<Stepper_rep*> (inside (stepper)) -> homv (x, k, t1);
00115 }
00116 
00117 template<typename C, typename V> vector<vector<C> >
00118 homx (const Stepper& stepper,
00119       const vector<vector<C> >& x, nat k, const vector<C>& p)
00120 {
00121   return const_cast<Stepper_rep*> (inside (stepper)) -> homvp (x, k, p);
00122 }
00123 
00124 template<typename C> nat
00125 start_index (const vector<C>& v, const C& z) {
00126   typedef Real_type(C) R;
00127   R eps= sqrt (sqrt (Accuracy (R)));
00128   // NOTE: take fourth root in order to compensate precision loss
00129   // when converting from half precision to full precision in
00130   // recursive methods
00131   if (N(v) >= 2)
00132     for (nat i=N(v)-1; i>0; i--) {
00133       C z1= v[i-1], z2= v[i];
00134       if (abs (z1 - z2) >= eps) {
00135         C t = (z - z1) / (z2 - z1);
00136         if (abs (Im (t)) < eps && Re (t) >= -eps &&
00137             Re (t) <= promote (1, eps) + eps)
00138           return i;
00139       }
00140       else if (abs (z2 - z) < eps) return i;
00141     }
00142   ERROR ("point not on path");
00143 }
00144 
00145 /******************************************************************************
00146 * External interface for numeric homotopies
00147 ******************************************************************************/
00148 
00149 template<typename C, typename Path> vector<C>
00150 homotopy (const vector<multivariate<sparse_polynomial<rational> > >& f,
00151           const vector<multivariate_coordinate<> >& vars,
00152           const vector<C>& x,
00153           const multivariate_coordinate<>& var,
00154           const Path& p) {
00155   typedef Homotopy_variant(C) V;
00156   Function ff= as_slp_tangent (f, vars, CF (x));
00157   Stepper stepper (ff);
00158   nat k= find (vars, var);
00159   ASSERT (k < N(vars), "coordinate not found");
00160   return homx (stepper, x, k, p);
00161 }
00162 
00163 template<typename C, typename Path> vector<vector<C> >
00164 homotopy (const vector<multivariate<sparse_polynomial<rational> > >& f,
00165           const vector<multivariate_coordinate<> >& vars,
00166           const vector<vector<C> >& x,
00167           const multivariate_coordinate<>& var,
00168           const Path& p) {
00169   if (N(x) == 0) return x;
00170   typedef Homotopy_variant(C) V;
00171   Function ff= as_slp_tangent (f, vars, CF (x[0]));
00172   Stepper stepper (ff);
00173   nat k= find (vars, var);
00174   ASSERT (k < N(vars), "coordinate not found");
00175   return homx (stepper, x, k, p);
00176 }
00177 
00178 template<typename C, typename Path> matrix<C>
00179 homotopy (const vector<multivariate<sparse_polynomial<rational> > >& f,
00180           const vector<multivariate_coordinate<> >& vars,
00181           const matrix<C>& x,
00182           const multivariate_coordinate<>& var,
00183           const Path& p) {
00184   return row_matrix (homotopy (f, vars, row_vectors (x), var, p));
00185 }
00186 
00187 } // namespace mmx
00188 #undef TMPL
00189 #undef Stepper
00190 #undef Stepper_rep
00191 #undef Function
00192 #endif // __MMX_HOMOTOPY_HPP