00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_RIEMANN_SQUARE_HPP
00014 #define __MMX_RIEMANN_SQUARE_HPP
00015 #include <basix/table.hpp>
00016 #include <numerix/complex_double.hpp>
00017
00018 namespace mmx {
00019 class std_riemann {};
00020 #define TMPL_DEF template<typename R,typename V= std_riemann>
00021 #define TMPL template<typename R,typename V>
00022 #define C complex<R>
00023 #define Square riemann_square<R,V>
00024 #define Square_rep riemann_square_rep<R,V>
00025 #define Squares table<bool,Square >
00026
00027
00028
00029
00030
00031 TMPL_DEF
00032 class riemann_square_rep REP_STRUCT {
00033 private:
00034 static double next_serial;
00035 public:
00036 C c;
00037 R r;
00038 double serial;
00039 public:
00040 inline riemann_square_rep (const C& c2, const R& r2):
00041 c (c2), r (r2), serial (next_serial++) {}
00042 };
00043
00044 TMPL_DEF
00045 class riemann_square {
00046 INDIRECT_PROTO_2 (riemann_square, riemann_square_rep, R, V)
00047 public:
00048 inline riemann_square ():
00049 rep (new Square_rep (C (), R ())) {}
00050 inline riemann_square (const format<R>& fm):
00051 rep (new Square_rep (C (promote (0, fm)), promote (0, fm))) {}
00052 inline riemann_square (const C& c, const R& r):
00053 rep (new Square_rep (c, r)) {}
00054 };
00055 INDIRECT_IMPL_2 (riemann_square, riemann_square_rep,
00056 typename R, R, typename V, V)
00057 DEFINE_UNARY_FORMAT_2 (riemann_square)
00058
00059 TMPL double Square_rep::next_serial= 0;
00060 TMPL inline C center (const Square& s) { return s->c; }
00061 TMPL inline R radius (const Square& s) { return s->r; }
00062 TMPL inline R diameter (const Square& s) { return s->r + s->r; }
00063 TMPL inline format<R> CF (const Square& s) { return get_format (radius (s)); }
00064
00065 TMPL vector<C >
00066 corners (const Square& s) {
00067 R r= radius (s);
00068 return center (s) + vec<C > (C (r, r), C (r, -r), C (-r, r), C (-r, -r));
00069 }
00070
00071
00072
00073
00074
00075 template<typename Op, typename R, typename V> nat
00076 unary_hash (const Square& s) {
00077 return hard_hash (s);
00078 }
00079
00080 template<typename Op, typename R, typename V> bool
00081 binary_test (const Square& s1, const Square& s2) {
00082 return hard_eq (s1, s2);
00083 }
00084
00085 TRUE_IDENTITY_OP_SUGAR(TMPL,Square)
00086 EXACT_IDENTITY_OP_SUGAR(TMPL,Square)
00087
00088 TMPL syntactic
00089 flatten (const Square& z) {
00090 return apply ("square", flatten (center (z)), flatten (radius (z)));
00091 }
00092
00093 TMPL
00094 struct binary_helper<Square >: public void_binary_helper<Square > {
00095 static inline string short_type_name () {
00096 return "Rsqr" * Short_type_name (R); }
00097 static inline generic full_type_name () {
00098 return gen ("Riemann_square", Full_type_name (R)); }
00099 static inline generic disassemble (const Square& s) {
00100 return gen_vec (as<generic> (center (s)),
00101 as<generic> (radius (s))); }
00102 static inline Square assemble (const generic& v) {
00103 return Square (as<C > (vector_access (v, 0)),
00104 as<R > (vector_access (v, 1))); }
00105 static inline void write (const port& out, const Square& s) {
00106 binary_write<C > (out, center (s));
00107 binary_write<R > (out, radius (s)); }
00108 static inline Square read (const port& in) {
00109 C c= binary_read<C > (in);
00110 R r= binary_read<R > (in);
00111 return Square (c, r); }
00112 };
00113
00114
00115
00116
00117
00118 TMPL bool
00119 admissible (const Square& s) {
00120 R r= radius (s);
00121 C z= center (s) / r;
00122 R I= promote (1, r);
00123 R Z= promote (2, r);
00124 R x= Re(z) - Z * floor (Re (z) / Z);
00125 R y= Im(z) - Z * floor (Im (z) / Z);
00126 return x == I && y == I;
00127 }
00128
00129 TMPL bool
00130 intersect (const Square& s1, const Square& s2) {
00131 R d= radius (s1) + radius (s2);
00132 return abs (Re (center (s1)) - Re (center (s2))) <= d &&
00133 abs (Im (center (s1)) - Im (center (s2))) <= d;
00134 }
00135
00136 TMPL bool
00137 included (const Square& s1, const Square& s2) {
00138 R d= radius (s2) - radius (s1);
00139 return radius (s1) <= radius (s2) &&
00140 abs (Re (center (s1)) - Re (center (s2))) <= d &&
00141 abs (Im (center (s1)) - Im (center (s2))) <= d;
00142 }
00143
00144 TMPL bool
00145 adjacent (const Square& s1, const Square& s2) {
00146 R d= radius (s1) + radius (s2);
00147 return
00148 (abs (Re (center (s1)) - Re (center (s2))) == d &&
00149 abs (Im (center (s1)) - Im (center (s2))) < d) ||
00150 (abs (Re (center (s1)) - Re (center (s2))) < d &&
00151 abs (Im (center (s1)) - Im (center (s2))) == d);
00152 }
00153
00154 TMPL bool
00155 inside (const C& z, const Square& s) {
00156 return abs (Re (z) - Re (center (s))) <= radius (s) &&
00157 abs (Im (z) - Im (center (s))) <= radius (s);
00158 }
00159
00160
00161
00162
00163
00164 TMPL void
00165 digital_disk_sub (Squares& ss,
00166 const C& u, const R& d,
00167 const C& z, const R& r, nat prec)
00168 {
00169 if (promote ((int) prec, r) * d >= r) {
00170 R h= d / promote (2, r);
00171 Square s (u + C (h, h), h);
00172 ASSERT (admissible (s), "admissibility hypothesis violated");
00173 vector<C > v= corners (s);
00174 nat count= 0;
00175 for (nat i=0; i<N(v); i++)
00176 if (abs (v[i] - z) < r) count++;
00177 if (count == 0 && !intersect (s, Square (z, r)));
00178 else if (count == 4) ss[s]= true;
00179 else {
00180 digital_disk_sub<R,V> (ss, u, h, z, r, prec);
00181 digital_disk_sub<R,V> (ss, u + C (h, promote (0, h)), h, z, r, prec);
00182 digital_disk_sub<R,V> (ss, u + C (promote (0, h), h), h, z, r, prec);
00183 digital_disk_sub<R,V> (ss, u + C (h, h), h, z, r, prec);
00184 }
00185 }
00186 }
00187
00188 TMPL Squares
00189 digital_disk (const C& z, const R& r, nat prec= 8) {
00190
00191 Squares ss;
00192 R scale= incexp2 (promote (1, r), exponent (r) + 2);
00193 R x= floor ((Re (z) - r) / scale) * scale;
00194 R y= floor ((Im (z) - r) / scale) * scale;
00195 digital_disk_sub<R,V> (ss, C (x, y), scale, z, r, prec);
00196 digital_disk_sub<R,V> (ss, C (x + scale, y), scale, z, r, prec);
00197 digital_disk_sub<R,V> (ss, C (x, y + scale), scale, z, r, prec);
00198 digital_disk_sub<R,V> (ss, C (x + scale, y + scale), scale, z, r, prec);
00199 return ss;
00200 }
00201
00202 #undef TMPL_DEF
00203 #undef TMPL
00204 #undef C
00205 #undef Square
00206 #undef Square_rep
00207 #undef Squares
00208 }
00209 #endif // __MMX_RIEMANN_SQUARE_HPP
