perturbation based calculations
[maximus:emndl.git] / emndl_calculate.cc
1 /*
2     emndl -- exponentially transformed Mandelbrot Set renderer
3     Copyright (C) 2011,2012  Claude Heiland-Allen <claude@mathr.co.uk>
4
5     This program is free software: you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation, either version 3 of the License, or
8     (at your option) any later version.
9
10     This program is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13     GNU General Public License for more details.
14
15     You should have received a copy of the GNU General Public License
16     along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <sys/types.h>
23 #include <sys/stat.h>
24 #include <pthread.h>
25
26 #include <boost/multiprecision/mpfr.hpp>
27 #include <complex>
28 #include <cmath>
29 #include <vector>
30 #include <list>
31
32 #define likely(x)   __builtin_expect((x),1)
33 #define unlikely(x) __builtin_expect((x),0)
34
35 #define pi 3.141592653589793
36
37 typedef boost::multiprecision::mpfr_float R;
38 typedef std::complex<R> C;
39 typedef double R53;
40 typedef std::complex<R53> C53;
41
42 #define CHANNELS 3
43 static const R53 escapeRadius2 = 65536.0 * 65536.0;
44
45 R53 lower(R x) {
46   return mpfr_get_d(x.backend().data(), MPFR_RNDN);
47 }
48
49 C53 lower(C x) {
50   return C53(lower(x.real()), lower(x.imag()));
51 }
52
53 C higher(C53 x) {
54   return C(R(x.real()), R(x.imag()));
55 }
56
57 /*
58 class doublee {
59   double x;
60   int e;
61   static const double range = 65536.0;
62   static const double rrange = 1.0/65536.0;
63 public:
64   doublee(R xr) {
65     R::default_precision(20);
66     R ranger(range);
67     R rranger(rrange);
68     e = 0;
69     if (abs(xr) > 0) {
70       while (abs(xr) > ranger) {
71         xr *= rranger;
72         e += 1;
73       }
74       while (abs(xr) < rranger) {
75         xr *= ranger;
76         e -= 1;
77       }
78     }
79     x = lower(xr);
80   };
81 };
82 */
83
84 class reference {
85   int prec;
86   C c;
87   R r;
88   int period;
89   std::vector<C53> zs;
90   std::vector<C53> dzs;
91 public:
92   reference(const char *x, const char *y, const char *z0, const char *p) {
93     R::default_precision(20);
94     r = R(z0);
95     prec = int(ceil(20.0 - lower(log10(r))));
96     R::default_precision(prec);
97     c = C(R(x), R(y));
98     period = atoi(p);
99     zs.resize(period);
100     dzs.resize(period);
101     C zero(R(0.0), R(0.0));
102     C one(R(1.0), R(0.0));
103     C two(R(2.0), R(0.0));
104     C z = zero;
105     C dz = zero;
106     for (int n = 0; n < period - 1; ++n) {
107       dz = two * z * dz + one;
108       z = z * z + c;
109       zs[n] = lower(z);
110       dzs[n] = lower(dz);
111     }
112     zs[period - 1] = lower(zero);
113   };
114   C53 getZ(int n) const { return zs[n % period]; };
115   C53 getDZ(int n) const { return dzs[n % period]; };
116   C getC() const { return c; };
117 };
118
119 class context1 {
120   const reference &ref;
121   int i, j;
122   C53 d0;
123   C53 d;
124   C53 a;
125   int n;
126   R53 w2;
127   int m;
128 public:
129   context1(const reference &reff, int ii, int jj, C53 delta)
130   : ref(reff), i(ii), j(jj), d0(delta), d(d0), a(1.0, 0.0), n(0), w2(escapeRadius2), m(0) {
131     // std::cerr << d0 << std::endl;
132   };
133   bool step1() {
134     C53 z = ref.getZ(n);
135     // C53 dz = ref.getDZ(n);
136     C53 w = z + d;
137     R53 zd2 = norm(w);
138     if (zd2 > escapeRadius2) {
139       return true;
140     }
141     a = 2.0 * w * a + 1.0;
142     // l = 2.0 * (dz * d + l * w);
143     d = 2.0 * z * d + d * d + d0;
144     n = n + 1;
145 //    zd2 /= n * n * log2(n);
146     if (zd2 < w2) {
147       w2 = zd2;
148       m = n;
149     }
150     return false;
151   };
152   bool steps(int maxIters) {
153     while (n < maxIters) {
154       bool escaped = step1();
155       if (escaped) {
156         return true;
157       }
158     }
159     return false;
160   };
161   int getI() { return i; };
162   int getJ() { return j; };
163   R53 getN() {
164     C53 w = ref.getZ(n) + d;
165     return n - log2(log(norm(w)) / log(escapeRadius2));
166   };
167   R getD() {
168     R::default_precision(20);
169     R53 w = abs(ref.getZ(n) + d);
170     R dw = abs(higher(a));
171     return R(w * log(w)) / dw;
172   };
173   int getM() { return m; }
174 };
175
176 bool compare(context1 *x, context1 *y) {
177   return x->getJ() < y->getJ();
178 }
179
180 class context {
181   int prec;
182   C c;
183   R r;
184   std::vector<reference> refs;
185   
186   int maxIters;
187   int w, h;
188   R53 x, y;
189   std::list<context1*> active[2];
190   int from, to;
191   FILE *out;
192   int progress;
193   int running;
194   
195 public:
196   context(const char *ofile, const char *width, const char *x0, const char *y0, const char *z0) {
197     R::default_precision(20);
198     r = R(z0);
199     prec = int(ceil(20 - lower(log10(r))));
200     R::default_precision(prec);
201     c = C(R(x0), R(y0));
202     w = atoi(width);
203     h = 0;
204     maxIters = 16;
205     out = fopen(ofile, "wb+");
206     x = 2.0 * pi / w;
207     y = 8.0 / w;
208   };
209   void addReference(reference ref) {
210     refs.push_back(ref);
211   };
212   const reference& getReference(int i, int j) {
213     // TODO pick best reference
214     return refs[refs.size() - 1];
215   };
216   C53 getDelta(const reference& ref, int i, int j) const {
217     double angle = i * x;
218     double radius = 64 * pow(0.5, j * y);
219     R::default_precision(prec);
220     return lower(higher(std::polar(radius, angle)) + c - ref.getC());
221   };
222
223   void todo(int which, context1 *ctx1) {
224     active[which].push_back(ctx1);
225   };
226   void done(int i, int j, double dwell, double distance, double domain) {
227     if (i < 0 || w <= i) {
228       exit(1);
229     }
230     while (j >= h - 1) {
231       if (0 != fseek(out, h * w * CHANNELS * sizeof(float), SEEK_SET)) {
232         fprintf(stderr, "ouch A\n"); exit(1);
233       }
234       int bytes = w * w * CHANNELS * sizeof(float);
235       void *zero = calloc(1, bytes);
236       fwrite(zero, bytes, 1, out);
237       free(zero);
238       h += w;
239     }
240     float result[CHANNELS];
241     if (j >= 0) {
242       progress++;
243       result[0] = dwell;
244       result[1] = distance;
245       result[2] = domain;
246       fseek(out, ((w * j) + i) * CHANNELS * sizeof(float), SEEK_SET);
247       fwrite(result, sizeof(result), 1, out);
248     }
249     result[0] = -1;
250     result[1] =  0;
251     result[2] =  0;
252     for (int dj = 1; dj >= -1; --dj) {
253       int jj = j + dj;
254       if (jj < 0) continue;
255       for (int di = 1; di >= -1; --di) {
256         int ii = (i + di + w) % w;
257         fseek(out, ((w * jj) + ii) * CHANNELS * sizeof(float), SEEK_SET);
258         float current;
259         if (1 != fread(&current, sizeof(current), 1, out)) {
260           fprintf(stderr, "ouch\n");
261           exit(1);
262         }
263         if (current == 0) {
264           current = -1;
265           fseek(out, ((w * jj) + ii) * CHANNELS * sizeof(float), SEEK_SET);
266           fwrite(&current, sizeof(current), 1, out);
267           const reference &ref = getReference(ii, jj);
268           todo(from, new context1(ref, ii, jj, getDelta(ref, ii, jj)));;
269         }
270       }
271     }
272   };
273   void work() {
274     while (! active[from].empty()) {
275       context1 *job = active[from].front();
276       active[from].pop_front();
277       if (job->steps(maxIters)) {
278         progress++;
279         double i = job->getI();
280         double j = job->getJ();
281         double dwell = job->getN();
282         int m = job->getM();
283         R::default_precision(20);
284         R distance(job->getD());
285         R psize(R(4.0) / (pow(R(0.5), R(j * y)) * R(64.0 * 2.0 * pi / w * sqrt(2.0))));
286         delete job;
287         done(i, j, dwell, lower(distance * psize), m);
288 //        std::cerr << "!";
289       } else {
290 //        std::cerr << "?";
291         todo(to, job);
292       }
293     }
294   }
295   void run() {
296     from = 0;
297     to = 1;
298     for (int i = 0; i < w; ++i) {
299       done(i, -1, 0, 0, 0);
300     }
301     running = 1;
302     while (running) {
303       printf("\n%12d\t%12lu\t%12d\t", maxIters, active[from].size(), h);
304       fflush(stdout);
305       running = 0;
306       do {
307         progress = 0;
308         work();
309         running += progress;
310       } while (!active[from].empty() && progress);
311       int tmp = from; from = to; to = tmp;
312       maxIters *= 2;
313     }
314     fclose(out);
315   };
316 };
317
318 int main(int argc, char **argv) {
319   if (argc != 3) {
320     return 1;
321   }
322   context ctx
323     ( argv[1]
324     , argv[2]
325 /*
326     , "-1.76988189381332902217371904140355696905734691301304962233987271036221e+00"
327     , "5.57194337987873531406774258042133639119793043644019875994399531213881e-02"
328     , "0.8e-55"
329 */
330     , "-1.7698818938133290221737190414035569690573469130130496223583710105429654750665986522688236089835043356094242938111219077244089e+00"
331     , "5.57194337987873531406774258042133639119793043644019889893941586623644602512684656153776745982821883497039424741300938397071099e-02"
332     , "1.35980225e-113"
333 /*
334     , "-1.7698818938133290221737190414035569690573469130130496223583710105429654750665986522688236089835043356094242938110655682091641213767416720102434695905599727191518170225623947453566511293626639781277581571799152771558189037167341121079864770695574353965966693948331605429393494869327714330939221107934482464652474399479732131558678897178990667327336397521011366907555946353585138906069842130897974963435e+00"
335     , "5.5719433798787353140677425804213363911979304364401988989394158662364460251268465615377674598282188349703942474155291706021187795454741550027328821580609808859951023647041146769532379567307529103649336344529380903216434731753770444315908351072155367170132486177853100584966076749231455129031539257363431685241625535974645252402968331304863445824203128661550355816354063819082621362277216597030590627685e-02"
336     , "0.5e-386"
337 */
338     );
339   ctx.addReference(reference
340     ( "-1.754877666273387e+00"
341     , "0e+00"
342     , "1.36514865e-02"
343     , "3"
344     ));
345   ctx.addReference(reference
346     ( "-1.76987732307224376882e+00"
347     , "5.572268594288084728183e-02"
348     , "2.08258871e-06"
349     , "10"
350     ));
351   ctx.addReference(reference
352     ( "-1.769881893813268431539866879e+00"
353     , "5.571943379865769811125326774e-02"
354     , "4.09506515e-14"
355     , "41"
356     ));
357   ctx.addReference(reference
358     ( "-1.7698818938133290221737190413207276256351e+00"
359     , "5.5719433798787353140677425767530493278864e-02"
360     , "1.92613713e-29"
361     , "200"
362     ));
363   ctx.addReference(reference
364     ( "-1.76988189381332902217371904140355696905734691301304962233987271036221e+00"
365     , "5.57194337987873531406774258042133639119793043644019875994399531213881e-02"
366     , "2.59580383e-55"
367     , "1051"
368     ));
369   ctx.addReference(reference
370     ( "-1.7698818938133290221737190414035569690573469130130496223583710105429654750665986522688236089835043356094242938111219077244089e+00"
371     , "5.57194337987873531406774258042133639119793043644019889893941586623644602512684656153776745982821883497039424741300938397071099e-02"
372     , "1.35980225e-113"
373     , "6511"
374     ));
375 /*
376   ctx.addReference(reference
377     ( "-1.7698818938133290221737190414035569690573469130130496223583710105429654750665986522688236089835043356094242938110655682091641213767416720102434695905599727191518170225623947453566511293626639781277581571799152771558190585277039489713e+00"
378     , "5.5719433798787353140677425804213363911979304364401988989394158662364460251268465615377674598282188349703942474155291706021187795454741550027328821580609808859951023647041146769532379567307529103649336344529380903216531114065152366358e-02"
379     , "4.16436771e-218"
380     , "47720"
381     ));
382   ctx.addReference(reference
383     ( "-1.7698818938133290221737190414035569690573469130130496223583710105429654750665986522688236089835043356094242938110655682091641213767416720102434695905599727191518170225623947453566511293626639781277581571799152771558189037167341121079864770695574353965966693948331605429393494869327714330939221107934482464652474399479732131558678897178990667327336397521011366907555946353585138906069842130897974963435e+00"
384     , "5.5719433798787353140677425804213363911979304364401988989394158662364460251268465615377674598282188349703942474155291706021187795454741550027328821580609808859951023647041146769532379567307529103649336344529380903216434731753770444315908351072155367170132486177853100584966076749231455129031539257363431685241625535974645252402968331304863445824203128661550355816354063819082621362277216597030590627685e-02"
385     , "1.87572609e-386"
386     , "389322"
387     ));
388 */
389 /*
390   ctx.addReference(reference
391     ( ""
392     , ""
393     , ""
394     , ""
395     ));
396 */
397   ctx.run();
398   return 0;
399 }