continewz: include/continewz/riemann

00001 
00002 /******************************************************************************
00003 * MODULE     : riemann_surface.hpp
00004 * DESCRIPTION: digital Riemann surfaces
00005 * COPYRIGHT  : (C) 2008  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_RIEMANN_SURFACE_HPP
00014 #define __MMX_RIEMANN_SURFACE_HPP
00015 #include <continewz/riemann_point.hpp>
00016 
00017 namespace mmx {
00018 #define TMPL_DEF template<typename R,typename V= std_riemann>
00019 #define TMPL template<typename R,typename V>
00020 #define C complex<R>
00021 #define Square riemann_square<R,V>
00022 #define Point riemann_point<R,V>
00023 #define Squares table<bool,Square >
00024 #define Neighbours table<Squares,Square >
00025 #define Upgrades table<Square,Square >
00026 #define Iterator iterator<Square >
00027 #define Surface riemann_surface<R,V>
00028 #define Surface_rep riemann_surface_rep<R,V>
00029 
00030 /******************************************************************************
00031 * Digital Riemann surfaces
00032 ******************************************************************************/
00033 
00034 TMPL_DEF
00035 class riemann_surface_rep REP_STRUCT {
00036 public:
00037   Squares    squares;
00038   Neighbours neighbours;
00039   Squares    invalid;
00040   Squares    replaced;
00041   Upgrades   upgrades;
00042 public:
00043   inline riemann_surface_rep () {}
00044   inline riemann_surface_rep (const format<R>& fm):
00045     squares (false, format<Square > (fm)),
00046     neighbours (Squares (false, format<Square > (fm)), format<Square > (fm)),
00047     invalid (false, format<Square > (fm)),
00048     replaced (false, format<Square > (fm)),
00049     upgrades (format<Square > (fm), format<Square > (fm)) {}
00050   bool add (const Square& s);
00051   bool add (const Squares& ss);
00052   bool connect (const Square& s1, const Square& s2);
00053   bool connect (const Neighbours& t);
00054   bool plane_connect (const Squares& ss);
00055   bool glue (const Square& s1, const Square& s2);
00056   Squares owners (const Square& s, const C& p) const;
00057   void normalize ();
00058   void clean ();
00059   Square upgrade (const Square& s);
00060   Point follow (const Square& orig, const C& from, const C& to) const;
00061   Point move (const Point& p, const C& delta) const;
00062 };
00063 
00064 TMPL_DEF
00065 class riemann_surface {
00066 INDIRECT_PROTO_2 (riemann_surface, riemann_surface_rep, R, V)
00067 public:
00068   inline riemann_surface (): rep (new Surface_rep ()) {}
00069   inline riemann_surface (const format<R>& fm): rep (new Surface_rep (fm)) {}
00070 };
00071 INDIRECT_IMPL_2 (riemann_surface, riemann_surface_rep,
00072                  typename R, R, typename V, V)
00073 DEFINE_UNARY_FORMAT_2 (riemann_surface)
00074 
00075 TMPL inline format<R> CF (const Surface& rs) {
00076   return get_format1 (CF2 (rs->squares)); }
00077 
00078 /******************************************************************************
00079 * Required routines
00080 ******************************************************************************/
00081 
00082 template<typename Op, typename R, typename V> nat
00083 unary_hash (const Surface& rs) {
00084   return hard_hash (rs);
00085 }
00086 
00087 template<typename Op, typename R, typename V> bool
00088 binary_test (const Surface& rs1, const Surface& rs2) {
00089   return hard_eq (rs1, rs2);
00090 }
00091 
00092 TRUE_IDENTITY_OP_SUGAR(TMPL,Surface)
00093 EXACT_IDENTITY_OP_SUGAR(TMPL,Surface)
00094 
00095 TMPL syntactic
00096 flatten (const Surface& rs) {
00097   //return apply ("surface", flatten (rs->squares));
00098   return apply ("surface", flatten (rs->squares), flatten (rs->neighbours));
00099 }
00100 
00101 TMPL
00102 struct binary_helper<Surface >: public void_binary_helper<Surface > {
00103   static inline string short_type_name () {
00104     return "Rsur" * Short_type_name (R); }
00105   static inline generic full_type_name () {
00106     return gen ("Riemann_surface", Full_type_name (R)); }
00107   static inline generic disassemble (const Surface& rs) {
00108     return gen_vec (as<generic> (rs->squares),
00109                     as<generic> (rs->neighbours)); }
00110   static inline Surface assemble (const generic& v) {
00111     Squares sqs= as<Squares > (vector_access (v, 0));
00112     Surface rs (get_format1 (CF2(sqs)));
00113     inside (rs)->add (sqs);
00114     inside (rs)->connect (as<Neighbours > (vector_access (v, 1)));
00115     return rs; }
00116   static inline void write (const port& out, const Surface& rs) {
00117     binary_write<format<R> > (out, CF (rs));
00118     binary_write<Squares > (out, rs->squares);
00119     binary_write<Neighbours > (out, rs->neighbours); }
00120   static inline Surface read (const port& in) {
00121     Surface rs (binary_read<format<R> > (in));
00122     inside (rs)->add (binary_read<Squares > (in));
00123     inside (rs)->connect (binary_read<Neighbours > (in));
00124     return rs; }
00125 };
00126 
00127 /******************************************************************************
00128 * Adding squares and connecting squares
00129 ******************************************************************************/
00130 
00131 TMPL bool
00132 Surface_rep::add (const Square& s) {
00133   // Add s to the Riemann surface
00134   if (squares->contains (s)) return false;
00135   squares[s]= true;
00136   neighbours[s]= Squares (false, get_format (s));
00137   return true;
00138 }
00139 
00140 TMPL bool
00141 Surface_rep::add (const Squares& ss) {
00142   // Add all squares in ss to the Riemann surface
00143   bool r= false;
00144   for (Iterator it= entries (ss); busy (it); ++it)
00145     r= add (*it) || r;
00146   return r;
00147 }
00148 
00149 TMPL bool
00150 Surface_rep::connect (const Square& s1, const Square& s2) {
00151   // Connect s1 and s2, if not already done
00152   // The projections of s1 and s2 are required to be neighbours
00153   if (!(squares->contains (s1) && squares->contains (s2))) return false;
00154   if (neighbours[s1]->contains (s2) && neighbours[s2]->contains (s1))
00155     return false;
00156   //mmout << "    Connect " << s1 << ", " << s2 << "\n";
00157   invalid[s1]= true;
00158   invalid[s2]= true;
00159   neighbours[s1][s2]= true;
00160   neighbours[s2][s1]= true;
00161   return true;
00162 }
00163 
00164 TMPL bool
00165 Surface_rep::connect (const Neighbours& t) {
00166   // Connect all squares in the neighbour table
00167   bool r= false;
00168   for (Iterator it= entries (t); busy (it); ++it)
00169     for (Iterator nb= entries (t[*it]); busy (nb); ++nb)
00170       r= connect (*it, *nb) || r;
00171   return r;
00172 }
00173 
00174 TMPL bool
00175 Surface_rep::plane_connect (const Squares& ss) {
00176   // Connect each two squares in ss whose projections are neighbours
00177   // mmout << "Plane connect " << ss << "\n";
00178   bool r= false;
00179   for (Iterator it1= entries (ss); busy (it1); ++it1)
00180     for (Iterator it2= entries (ss); busy (it2); ++it2)
00181       if (adjacent (*it1, *it2))
00182         r= connect (*it1, *it2) || r;
00183   return r;
00184 }
00185 
00186 TMPL bool
00187 Surface_rep::glue (const Square& s1, const Square& s2) {
00188   // Glue s1 and s2 together
00189   if (!(squares->contains (s1) && squares->contains (s2))) return false;
00190   if (s1 == s2) return false;
00191   if (included (s1, s2) && (!included (s2, s1) || s2->serial < s1->serial)) {
00192     if (!squares->contains (s1)) return false;
00193     //mmout << "  Glue " << s1 << ", " << s2 << "\n";
00194     Squares ns= neighbours[s1];
00195     upgrades[s1]= s2;
00196     replaced [s1]= true;
00197     reset (squares, s1);
00198     reset (neighbours, s1);
00199     for (Iterator nb= entries (ns); busy (nb); ++nb)
00200       reset (neighbours[*nb], s1);
00201     for (Iterator nb= entries (ns); busy (nb); ++nb)
00202       if (adjacent (*nb, s2)) connect (*nb, s2);
00203     for (Iterator nb= entries (ns); busy (nb); ++nb)
00204       if (included (*nb, s2)) glue (*nb, s2);
00205     return true;
00206   }
00207   else {
00208     ASSERT (included (s2, s1), "projections do not match");
00209     return glue (s2, s1);
00210   }
00211 }
00212 
00213 /******************************************************************************
00214 * Normalization
00215 ******************************************************************************/
00216 
00217 TMPL Squares
00218 Surface_rep::owners (const Square& s, const C& p) const {
00219   // Return list of squares which contain the point p (on s)
00220   Squares r (false, get_format (s));
00221   Squares todo (false, get_format (s));
00222   todo[s]= true;
00223   while (N(todo) != 0) {
00224     Squares temp= todo;
00225     todo= Squares (false, get_format (s));
00226     for (Iterator it= entries (temp); busy (it); ++it) {
00227       r[*it]= true;
00228       for (Iterator nb= entries (neighbours[*it]); busy (nb); ++nb)
00229         if (inside (p, *nb) && !r[*nb]) todo[*nb]= true;
00230     }
00231   }
00232   return r;
00233 }
00234 
00235 TMPL inline Squares
00236 owners (const Surface& rs, const Point& p) {
00237   return rs->owners (owner (p), project (p));
00238 }
00239 
00240 TMPL void
00241 Surface_rep::normalize () {
00242   // After making non trivial connections,
00243   // the Riemann surface should be normalized using this routine
00244   while (N(invalid) != 0) {
00245     Squares ss= invalid;
00246     invalid= Squares (false, CF2(ss));
00247     for (Iterator it= entries (ss); busy (it); ++it) {
00248       //mmout << "Normalize " << *it << "\n";
00249       for (Iterator nb1= entries (neighbours[*it]); busy (nb1); ++nb1)
00250         for (Iterator nb2= entries (neighbours[*nb1]); busy (nb2); ++nb2) {
00251           if (included (*it, *nb2)) glue (*it, *nb2);
00252           if (included (*nb2, *it)) glue (*nb2, *it);
00253         }
00254       vector<C > cs= corners (*it);
00255       for (nat i=0; i<N(cs); i++) {
00256         Squares os= owners (*it, cs[i]);
00257         if (plane_connect (os)) continue;
00258         Square first= *entries (os);
00259         R r= diameter (first);
00260         bool same_diameter= true;
00261         for (Iterator it2= entries (os); busy (it2); ++it2)
00262           same_diameter= same_diameter && (r == diameter (*it2));
00263         Square s (cs[i], r);
00264         if (same_diameter && admissible (s) && N(os) == 4) {
00265           //mmout << "  Join " << os << " -> " << s << "\n";
00266           add (s);
00267           glue (s, first);
00268         }
00269       }
00270     }
00271   }
00272 }
00273 
00274 TMPL void
00275 Surface_rep::clean () {
00276   // At the end of normalize, the variable replaced contains
00277   // the list of all squares which have been upgraded during gluing
00278   // The present routine resets this variable and performs some
00279   // appropriate final cleaning
00280   for (Iterator it= entries (replaced); busy (it); ++it)
00281     reset (neighbours, *it);
00282   replaced= Squares (false, CF2(replaced));
00283 }
00284 
00285 TMPL Square
00286 Surface_rep::upgrade (const Square& s) {
00287   if (squares->contains (s)) return s;
00288   ASSERT (upgrades->contains (s), "square not on Riemann surface");
00289   if (!squares->contains (upgrades[s]))
00290     upgrades[s]= upgrade (upgrades[s]);
00291   return upgrades[s];
00292 }
00293 
00294 TMPL inline Point
00295 upgrade (const Surface& rs, const Point& p) {
00296   return Point (inside (rs)->upgrade (owner (p)), project (p));
00297 }
00298 
00299 /******************************************************************************
00300 * Moving on the Riemann surface
00301 ******************************************************************************/
00302 
00303 TMPL Point
00304 Surface_rep::follow (const Square& s1, const C& z1, const C& z2) const {
00305   if (inside (z2, s1)) return Point (s1, z2);
00306   C c    = center (s1);
00307   R r    = radius (s1);
00308   C delta= z2 - z1;
00309   C z3   = center (s1);
00310   int dx = 0;
00311   int dy = 0;
00312   if (Re (z2) > Re (c) + r && Re (delta) != 0) {
00313     R x= Re (c) + r;
00314     R y= Im (z1) + (x - Re (z1)) * Im (delta) / Re (delta);
00315     if (abs (y - Im (c)) <= r) { z3= C (x, y); dx= 1; }
00316   }
00317   if (Re (z2) < Re (c) - r && Re (delta) != 0) {
00318     R x= Re (c) - r;
00319     R y= Im (z1) + (x - Re (z1)) * Im (delta) / Re (delta);
00320     if (abs (y - Im (c)) <= r) { z3= C (x, y); dx= -1; }
00321   }
00322   if (Im (z2) > Im (c) + r && Im (delta) != 0) {
00323     R y= Im (c) + r;
00324     R x= Re (z1) + (y - Im (z1)) * Re (delta) / Im (delta);
00325     if (abs (x - Re (c)) <= r) { z3= C (x, y); dy= 1; }
00326   }
00327   if (Im (z2) < Im (c) - r && Im (delta) != 0) {
00328     R y= Im (c) - r;
00329     R x= Re (z1) + (y - Im (z1)) * Re (delta) / Im (delta);
00330     if (abs (x - Re (c)) <= r) { z3= C (x, y); dy= -1; }
00331   }
00332   ASSERT (z3 != center (s1), "unexpected situation");
00333   for (Iterator it= entries (neighbours[s1]); busy (it); ++it) {
00334     Square s3= *it;
00335     if (inside (z3, s3)) {
00336       //mmout << "Follow to " << s3 << "? " << dx << ", " << dy << "\n";
00337       R d= radius (s1) + radius (s3);
00338       if ((dx > 0 && Re (center (s3) - center (s1)) == d) ||
00339           (dx < 0 && Re (center (s3) - center (s1)) == -d) ||
00340           (dy > 0 && Im (center (s3) - center (s1)) == d) ||
00341           (dy < 0 && Im (center (s3) - center (s1)) == -d))
00342         return follow (s3, z3, z2);
00343     }
00344   }
00345   return Point (s1, z3);
00346 }
00347 
00348 TMPL Point
00349 Surface_rep::move (const Point& p, const C& delta) const {
00350   Square s= owner (p);
00351   C z1= project (p);
00352   C z2= z1 + delta;
00353   return follow (s, z1, z2);
00354 }
00355 
00356 /******************************************************************************
00357 * Public interface
00358 ******************************************************************************/
00359 
00360 TMPL Surface
00361 riemann_disk (const C& z, const R& r, nat prec= 8) {
00362   Surface rs (get_format (r));
00363   Squares ds= digital_disk<R,V> (z, r, prec);
00364   inside (rs)->add (ds);
00365   inside (rs)->plane_connect (ds);
00366   inside (rs)->normalize ();
00367   inside (rs)->clean ();
00368   return rs;
00369 }
00370 
00371 TMPL Squares
00372 glue (Surface& rs, const Point& p1, const Point& p2) {
00373   ASSERT (project (p1) == project (p2), "incompatible glue points");
00374   inside (rs)->glue (owner (upgrade (rs, p1)), owner (upgrade (rs, p2)));
00375   inside (rs)->normalize ();
00376   Squares replaced= rs->replaced;
00377   inside (rs)->clean ();
00378   return replaced;
00379 }
00380 
00381 TMPL Squares
00382 glue (Surface& rs1, const Point& p1, const Surface& rs2, const Point& p2) {
00383   ASSERT (project (p1) == project (p2), "incompatible glue points");
00384   //mmout << "--- Add ---\n";
00385   inside (rs1)->add (rs2->squares);
00386   inside (rs1)->connect (rs2->neighbours);
00387   Squares os1= owners (rs1, upgrade (rs1, p1));
00388   Squares os2= owners (rs2, upgrade (rs2, p2));
00389   for (Iterator it1= entries (os1); busy (it1); ++it1)
00390     for (Iterator it2= entries (os2); busy (it2); ++it2)
00391       if (included (*it1, *it2) || included (*it2, *it1)) {
00392         //mmout << "--- Glue ---\n";
00393         inside (rs1)->glue (*it1, *it2);
00394         //mmout << "--- Normalize ---\n";
00395         inside (rs1)->normalize ();
00396         Squares replaced= rs1->replaced;
00397         inside (rs1)->clean ();
00398         //mmout << "--- Done ---\n";
00399         return replaced;
00400       }
00401   ERROR ("unexpected situation");
00402   return Squares (false, CF2(os1));
00403 }
00404 
00405 TMPL void
00406 glue (Surface& rs1, const Point& p1, const Surface& rs2) {
00407   Point p2= lift (rs2, project (p1));
00408   glue (rs1, p1, rs2, p2);
00409 }
00410 
00411 TMPL bool
00412 is_plane (const Surface& rs) {
00413   for (Iterator it1= entries (rs->squares); busy (it1); ++it1)
00414     for (Iterator it2= entries (rs->squares); busy (it2); ++it2)
00415       if (adjacent (*it1, *it2))
00416         if (!rs->neighbours[*it1]->contains (*it2))
00417           return false;
00418   return true;
00419 }
00420 
00421 TMPL Point
00422 lift (const Surface& rs, const C& z) {
00423   for (Iterator it= entries (rs->squares); busy (it); ++it)
00424     if (inside (z, *it)) return Point (*it, z);
00425   mmout << "rs= " << rs << "\n";
00426   mmout << "z= " << z << "\n";
00427   ERROR ("unable to lift point");
00428   return Point (CF(rs));
00429 }
00430 
00431 TMPL vector<Point >
00432 lifts (const Surface& rs, const C& z) {
00433   vector<Point > v;
00434   for (Iterator it= entries (rs->squares); busy (it); ++it)
00435     if (inside (z, *it)) v << Point (*it, z);
00436   return v;
00437 }
00438 
00439 TMPL inline Point
00440 move (const Surface& rs, const Point& p, const C& delta) {
00441   // Return p + t delta, where t <= 1 is maximal such that
00442   // the segment [p, p + t delta] lies on rs
00443   return rs->move (upgrade (rs, p), delta);
00444 }
00445 
00446 #undef TMPL_DEF
00447 #undef TMPL
00448 #undef C
00449 #undef Square
00450 #undef Point
00451 #undef Squares
00452 #undef Neighbours
00453 #undef Upgrades
00454 #undef Iterator
00455 #undef Surface
00456 #undef Surface_rep
00457 } // namespace mmx
00458 #endif // __MMX_RIEMANN_HPP
include/continewz/riemann_surface.hpp
