mightymandel v16

GPU-based Mandelbrot set explorer
atom.c
Go to the documentation of this file.
1 // mightymandel -- GPU-based Mandelbrot Set explorer
2 // Copyright (C) 2012,2013,2014,2015 Claude Heiland-Allen
3 // License GPL3+ http://www.gnu.org/licenses/gpl.html
4 
10 #include <mpfr.h>
11 
12 #include "atom.h"
13 
23 int crosses_positive_real_axis(const mpfr_t ax, const mpfr_t ay, const mpfr_t bx, const mpfr_t by) {
24  if (mpfr_sgn(ay) != mpfr_sgn(by)) {
25  mpfr_t dx, dy;
26  mpfr_inits2(mpfr_get_prec(ax), dx, dy, (mpfr_ptr) 0);
27  // d = b - a
28  mpfr_sub(dx, bx, ax, MPFR_RNDN);
29  mpfr_sub(dy, by, ay, MPFR_RNDN);
30  int s = mpfr_sgn(dy);
31  // dy = d `cross` a
32  mpfr_mul(dy, dy, ax, MPFR_RNDN);
33  mpfr_mul(dx, dx, ay, MPFR_RNDN);
34  mpfr_sub(dy, dy, dx, MPFR_RNDN);
35  int t = mpfr_sgn(dy);
36  return s == t;
37  }
38  return 0;
39 }
40 
54 int surrounds_origin(const mpfr_t ax, const mpfr_t ay, const mpfr_t bx, const mpfr_t by, const mpfr_t cx, const mpfr_t cy, const mpfr_t dx, const mpfr_t dy) {
55  return 1 & (
56  crosses_positive_real_axis(ax, ay, bx, by) +
57  crosses_positive_real_axis(bx, by, cx, cy) +
58  crosses_positive_real_axis(cx, cy, dx, dy) +
59  crosses_positive_real_axis(dx, dy, ax, ay) );
60 }
61 
71 int did_escaped(const mpfr_t x, const mpfr_t y) {
72  return (mpfr_regular_p(x) && mpfr_get_exp(x) > 10) || (mpfr_regular_p(y) && mpfr_get_exp(y) > 10) || ! mpfr_number_p(x) || ! mpfr_number_p(y);
73 }
74 
87 unsigned int boxperiod(const mpfr_t cx, const mpfr_t cy, const mpfr_t r, unsigned int maxperiod, void *abort_data, abort_t abort_fn) {
88  // allocate box
89  mpfr_prec_t prec = mpfr_get_prec(cx);
90  mpfr_t v[18];
91  for (int i = 0; i < 18; ++i) {
92  mpfr_init2(v[i], prec);
93  }
94  // initialize box
95  for (int i = 0; i < 8; ++i) {
96  mpfr_set_ui(v[i], 0, MPFR_RNDN);
97  }
98  mpfr_sub(v[ 8], cx, r, MPFR_RNDN);
99  mpfr_sub(v[ 9], cy, r, MPFR_RNDN);
100  mpfr_add(v[10], cx, r, MPFR_RNDN);
101  mpfr_sub(v[11], cy, r, MPFR_RNDN);
102  mpfr_add(v[12], cx, r, MPFR_RNDN);
103  mpfr_add(v[13], cy, r, MPFR_RNDN);
104  mpfr_sub(v[14], cx, r, MPFR_RNDN);
105  mpfr_add(v[15], cy, r, MPFR_RNDN);
106  // iterate
107  unsigned int period = 0;
108  for (unsigned int p = 1; p <= maxperiod; ++p) {
109  for (unsigned int i = 0; i < 8; i += 2) {
110  // z = z^2 + c
111  mpfr_sqr(v[16], v[i], MPFR_RNDN);
112  mpfr_sqr(v[17], v[i+1], MPFR_RNDN);
113  mpfr_mul(v[i+1], v[i], v[i+1], MPFR_RNDN);
114  mpfr_mul_2ui(v[i+1], v[i+1], 1, MPFR_RNDN);
115  mpfr_sub(v[i], v[16], v[17], MPFR_RNDN);
116  mpfr_add(v[i], v[i], v[i+8], MPFR_RNDN);
117  mpfr_add(v[i+1], v[i+1], v[i+9], MPFR_RNDN);
118  // check escaped
119  if (did_escaped(v[i], v[i+1])) {
120  goto done;
121  }
122  }
123  if (surrounds_origin(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7])) {
124  period = p;
125  goto done;
126  }
127  if (abort_fn(abort_data)) { break; }
128  }
129 done:
130  for (int i = 0; i < 18; ++i) {
131  mpfr_clear(v[i]);
132  }
133  return period;
134 }
135 
151 int muatom(int period, mpfr_t x, mpfr_t y, mpfr_t z, void *abort_data, abort_t abort_fn) {
152  const mpfr_prec_t accuracy = 32;
153  mpfr_prec_t bits = mpfr_get_prec(x);
154 #define VARS nx, ny, bx, by, wx, wy, zz, Ax, Ay, Bx, By, Cx, Cy, Dx, Dy, Ex, Ey, t1, t2, t3, t4, t5, t6, t7, t8, t9, t0, t01, t02, t03, t04
155  mpfr_t VARS;
156  mpfr_inits2(bits, VARS, (mpfr_ptr) 0);
157  mpfr_set(nx, x, MPFR_RNDN);
158  mpfr_set(ny, y, MPFR_RNDN);
159  mpfr_set_prec(zz, mpfr_get_prec(z));
160  int n_ok, w_ok, b_ok;
161  int r = 0;
162  do {
163  mpfr_set_prec(t0, bits - accuracy);
164  mpfr_set_prec(t01, bits - accuracy);
165  mpfr_set_prec(t02, bits - accuracy);
166  mpfr_set_prec(t03, bits - accuracy);
167  mpfr_set_prec(t04, bits - accuracy);
168  int limit = 100;
169  do {
170  // refine nucleus
171  mpfr_set_ui(Ax, 0, MPFR_RNDN);
172  mpfr_set_ui(Ay, 0, MPFR_RNDN);
173  mpfr_set_ui(Dx, 0, MPFR_RNDN);
174  mpfr_set_ui(Dy, 0, MPFR_RNDN);
175  for (int i = 0; i < period; ++i) {
176  if (abort_fn(abort_data)) { goto cleanup; }
177  // D <- 2AD + 1 : Dx <- 2 (Ax Dx - Ay Dy) + 1 ; Dy = 2 (Ax Dy + Ay Dx)
178  mpfr_mul(t1, Ax, Dx, MPFR_RNDN);
179  mpfr_mul(t2, Ay, Dy, MPFR_RNDN);
180  mpfr_mul(t3, Ax, Dy, MPFR_RNDN);
181  mpfr_mul(t4, Ay, Dx, MPFR_RNDN);
182  mpfr_sub(Dx, t1, t2, MPFR_RNDN);
183  mpfr_mul_2ui(Dx, Dx, 1, MPFR_RNDN);
184  mpfr_add_ui(Dx, Dx, 1, MPFR_RNDN);
185  mpfr_add(Dy, t3, t4, MPFR_RNDN);
186  mpfr_mul_2ui(Dy, Dy, 1, MPFR_RNDN);
187  // A <- A^2 + n : Ax <- Ax^2 - Ay^2 + nx ; Ay <- 2 Ax Ay + ny
188  mpfr_sqr(t1, Ax, MPFR_RNDN);
189  mpfr_sqr(t2, Ay, MPFR_RNDN);
190  mpfr_mul(t3, Ax, Ay, MPFR_RNDN);
191  mpfr_sub(Ax, t1, t2, MPFR_RNDN);
192  mpfr_add(Ax, Ax, nx, MPFR_RNDN);
193  mpfr_mul_2ui(t3, t3, 1, MPFR_RNDN);
194  mpfr_add(Ay, t3, ny,MPFR_RNDN);
195  }
196  // n <- n - A / D = (nx + ny i) - ((Ax Dx + Ay Dy) + (Ay Dx - Ax Dy)i) / (Dx^2 + Dy^2)
197  mpfr_sqr(t1, Dx, MPFR_RNDN);
198  mpfr_sqr(t2, Dy, MPFR_RNDN);
199  mpfr_add(t1, t1, t2, MPFR_RNDN);
200 
201  mpfr_mul(t2, Ax, Dx, MPFR_RNDN);
202  mpfr_mul(t3, Ay, Dy, MPFR_RNDN);
203  mpfr_add(t2, t2, t3, MPFR_RNDN);
204  mpfr_div(t2, t2, t1, MPFR_RNDN);
205  mpfr_sub(t2, nx, t2, MPFR_RNDN);
206 
207  mpfr_mul(t3, Ay, Dx, MPFR_RNDN);
208  mpfr_mul(t4, Ax, Dy, MPFR_RNDN);
209  mpfr_sub(t3, t3, t4, MPFR_RNDN);
210  mpfr_div(t3, t3, t1, MPFR_RNDN);
211  mpfr_sub(t3, ny, t3, MPFR_RNDN);
212 
213  mpfr_set(t01, t2, MPFR_RNDN);
214  mpfr_set(t02, t3, MPFR_RNDN);
215  mpfr_set(t03, nx, MPFR_RNDN);
216  mpfr_set(t04, ny, MPFR_RNDN);
217  mpfr_sub(t01, t01, t03, MPFR_RNDN);
218  mpfr_sub(t02, t02, t04, MPFR_RNDN);
219  mpfr_add_ui(t01, t01, 1, MPFR_RNDN);
220  mpfr_sub_ui(t01, t01, 1, MPFR_RNDN);
221  mpfr_add_ui(t02, t02, 1, MPFR_RNDN);
222  mpfr_sub_ui(t02, t02, 1, MPFR_RNDN);
223  n_ok = mpfr_zero_p(t01) && mpfr_zero_p(t02);
224 
225  mpfr_set(nx, t2, MPFR_RNDN);
226  mpfr_set(ny, t3, MPFR_RNDN);
227  } while (! n_ok && --limit);
228  if (! limit) goto cleanup;
229  mpfr_set(wx, nx, MPFR_RNDN);
230  mpfr_set(wy, ny, MPFR_RNDN);
231  mpfr_set(bx, nx, MPFR_RNDN);
232  mpfr_set(by, ny, MPFR_RNDN);
233  limit = 100;
234  do {
235  // refine bond
236  mpfr_set(Ax, wx, MPFR_RNDN);
237  mpfr_set(Ay, wy, MPFR_RNDN);
238  mpfr_set_ui(Bx, 1, MPFR_RNDN);
239  mpfr_set_ui(By, 0, MPFR_RNDN);
240  mpfr_set_ui(Cx, 0, MPFR_RNDN);
241  mpfr_set_ui(Cy, 0, MPFR_RNDN);
242  mpfr_set_ui(Dx, 0, MPFR_RNDN);
243  mpfr_set_ui(Dy, 0, MPFR_RNDN);
244  mpfr_set_ui(Ex, 0, MPFR_RNDN);
245  mpfr_set_ui(Ey, 0, MPFR_RNDN);
246  for (int i = 0; i < period; ++i) {
247  if (abort_fn(abort_data)) { goto cleanup; }
248  // E <- 2 ( A E + B D )
249  mpfr_mul(t1, Ax, Ex, MPFR_RNDN);
250  mpfr_mul(t2, Ay, Ey, MPFR_RNDN);
251  mpfr_sub(t1, t1, t2, MPFR_RNDN);
252  mpfr_mul(t2, Ax, Ey, MPFR_RNDN);
253  mpfr_mul(t3, Ay, Ex, MPFR_RNDN);
254  mpfr_add(t2, t2, t3, MPFR_RNDN);
255  mpfr_mul(t3, Bx, Dx, MPFR_RNDN);
256  mpfr_mul(t4, By, Dy, MPFR_RNDN);
257  mpfr_sub(t3, t3, t4, MPFR_RNDN);
258  mpfr_add(Ex, t1, t3, MPFR_RNDN);
259  mpfr_mul(t3, Bx, Dy, MPFR_RNDN);
260  mpfr_mul(t4, By, Dx, MPFR_RNDN);
261  mpfr_add(t3, t3, t4, MPFR_RNDN);
262  mpfr_add(Ey, t2, t3, MPFR_RNDN);
263  mpfr_mul_2ui(Ex, Ex, 1, MPFR_RNDN);
264  mpfr_mul_2ui(Ey, Ey, 1, MPFR_RNDN);
265  // D <- 2 A D + 1
266  mpfr_mul(t1, Ax, Dx, MPFR_RNDN);
267  mpfr_mul(t2, Ay, Dy, MPFR_RNDN);
268  mpfr_mul(t3, Ax, Dy, MPFR_RNDN);
269  mpfr_mul(t4, Ay, Dx, MPFR_RNDN);
270  mpfr_sub(Dx, t1, t2, MPFR_RNDN);
271  mpfr_mul_2ui(Dx, Dx, 1, MPFR_RNDN);
272  mpfr_add_ui(Dx, Dx, 1, MPFR_RNDN);
273  mpfr_add(Dy, t3, t4, MPFR_RNDN);
274  mpfr_mul_2ui(Dy, Dy, 1, MPFR_RNDN);
275  // C <- 2 ( B^2 + A C )
276  mpfr_sqr(t1, Bx, MPFR_RNDN);
277  mpfr_sqr(t2, By, MPFR_RNDN);
278  mpfr_sub(t1, t1, t2, MPFR_RNDN);
279  mpfr_mul(t2, Bx, By, MPFR_RNDN);
280  mpfr_mul_2ui(t2, t2, 1, MPFR_RNDN);
281  mpfr_mul(t3, Ax, Cx, MPFR_RNDN);
282  mpfr_mul(t4, Ay, Cy, MPFR_RNDN);
283  mpfr_sub(t3, t3, t4, MPFR_RNDN);
284  mpfr_add(t1, t1, t3, MPFR_RNDN);
285  mpfr_mul(t3, Ax, Cy, MPFR_RNDN);
286  mpfr_mul(t4, Ay, Cx, MPFR_RNDN);
287  mpfr_add(t3, t3, t4, MPFR_RNDN);
288  mpfr_add(t2, t2, t3, MPFR_RNDN);
289  mpfr_mul_2ui(Cx, t1, 1, MPFR_RNDN);
290  mpfr_mul_2ui(Cy, t2, 1, MPFR_RNDN);
291  // B <- 2 A B
292  mpfr_mul(t1, Ax, Bx, MPFR_RNDN);
293  mpfr_mul(t2, Ay, By, MPFR_RNDN);
294  mpfr_mul(t3, Ax, By, MPFR_RNDN);
295  mpfr_mul(t4, Ay, Bx, MPFR_RNDN);
296  mpfr_sub(Bx, t1, t2, MPFR_RNDN);
297  mpfr_mul_2ui(Bx, Bx, 1, MPFR_RNDN);
298  mpfr_add(By, t3, t4, MPFR_RNDN);
299  mpfr_mul_2ui(By, By, 1, MPFR_RNDN);
300  // A <- A^2 + b
301  mpfr_sqr(t1, Ax, MPFR_RNDN);
302  mpfr_sqr(t2, Ay, MPFR_RNDN);
303  mpfr_mul(t3, Ax, Ay, MPFR_RNDN);
304  mpfr_sub(Ax, t1, t2, MPFR_RNDN);
305  mpfr_add(Ax, Ax, bx, MPFR_RNDN);
306  mpfr_mul_2ui(t3, t3, 1, MPFR_RNDN);
307  mpfr_add(Ay, t3, by,MPFR_RNDN);
308  }
309  // (w) <- (w) - (B-1 D)^-1 (A - w)
310  // (b) <- (b) ( C E) (B + 1)
311  //
312  // det = (B-1)E - CD
313  // inv = (E -D)
314  // (-C (B-1) / det
315  //
316  // w - (E(A-w) - D(B+1))/det
317  // b - (-C(A-w) + (B-1)(B+1))/det ; B^2 - 1
318  //
319  // A -= w
320  mpfr_sub(Ax, Ax, wx, MPFR_RNDN);
321  mpfr_sub(Ay, Ay, wy, MPFR_RNDN);
322  // (t1,t2) = B^2 - 1
323  mpfr_mul(t1, Bx, Bx, MPFR_RNDN);
324  mpfr_mul(t2, By, By, MPFR_RNDN);
325  mpfr_sub(t1, t1, t2, MPFR_RNDN);
326  mpfr_sub_ui(t1, t1, 1, MPFR_RNDN);
327  mpfr_mul(t2, Bx, By, MPFR_RNDN);
328  mpfr_mul_2ui(t2, t2, 1, MPFR_RNDN);
329  // (t1,t2) -= C(A-w)
330  mpfr_mul(t3, Cx, Ax, MPFR_RNDN);
331  mpfr_mul(t4, Cy, Ay, MPFR_RNDN);
332  mpfr_sub(t3, t3, t4, MPFR_RNDN);
333  mpfr_sub(t1, t1, t3, MPFR_RNDN);
334  mpfr_mul(t3, Cx, Ay, MPFR_RNDN);
335  mpfr_mul(t4, Cy, Ax, MPFR_RNDN);
336  mpfr_add(t3, t3, t4, MPFR_RNDN);
337  mpfr_sub(t2, t2, t3, MPFR_RNDN);
338  // (t3, t4) = (B-1)E - CD
339  mpfr_sub_ui(t3, Bx, 1, MPFR_RNDN);
340  mpfr_mul(t4, t3, Ex, MPFR_RNDN);
341  mpfr_mul(t5, By, Ey, MPFR_RNDN);
342  mpfr_sub(t4, t4, t5, MPFR_RNDN);
343  mpfr_mul(t5, t3, Ey, MPFR_RNDN);
344  mpfr_mul(t6, By, Ex, MPFR_RNDN);
345  mpfr_sub(t5, t5, t6, MPFR_RNDN);
346  mpfr_mul(t6, Cx, Dx, MPFR_RNDN);
347  mpfr_mul(t7, Cy, Dy, MPFR_RNDN);
348  mpfr_sub(t6, t6, t7, MPFR_RNDN);
349  mpfr_sub(t3, t4, t6, MPFR_RNDN);
350  mpfr_mul(t6, Cx, Dy, MPFR_RNDN);
351  mpfr_mul(t7, Cy, Dx, MPFR_RNDN);
352  mpfr_add(t6, t6, t7, MPFR_RNDN);
353  mpfr_sub(t4, t5, t6, MPFR_RNDN);
354  // (t3, t4) = 1/(t3, t4) ; z^-1 = z* / (z z*)
355  mpfr_sqr(t5, t3, MPFR_RNDN);
356  mpfr_sqr(t6, t4, MPFR_RNDN);
357  mpfr_add(t5, t5, t6, MPFR_RNDN);
358  mpfr_div(t3, t3, t5, MPFR_RNDN);
359  mpfr_div(t4, t4, t5, MPFR_RNDN);
360  mpfr_neg(t4, t4, MPFR_RNDN);
361  // (t1, t2) *= (t3, t4)
362  mpfr_mul(t5, t1, t3, MPFR_RNDN);
363  mpfr_mul(t6, t2, t4, MPFR_RNDN);
364  mpfr_mul(t7, t1, t4, MPFR_RNDN);
365  mpfr_mul(t8, t2, t3, MPFR_RNDN);
366  mpfr_sub(t1, t5, t6, MPFR_RNDN);
367  mpfr_add(t2, t7, t8, MPFR_RNDN);
368 
369  // (t1, t2) = b - (t1, t2)
370  mpfr_sub(t1, bx, t1, MPFR_RNDN);
371  mpfr_sub(t2, by, t2, MPFR_RNDN);
372 
373  mpfr_set(t01, t1, MPFR_RNDN);
374  mpfr_set(t02, t2, MPFR_RNDN);
375  mpfr_set(t03, bx, MPFR_RNDN);
376  mpfr_set(t04, by, MPFR_RNDN);
377  mpfr_sub(t01, t01, t03, MPFR_RNDN);
378  mpfr_sub(t02, t02, t04, MPFR_RNDN);
379  mpfr_add_ui(t01, t01, 1, MPFR_RNDN);
380  mpfr_sub_ui(t01, t01, 1, MPFR_RNDN);
381  mpfr_add_ui(t02, t02, 1, MPFR_RNDN);
382  mpfr_sub_ui(t02, t02, 1, MPFR_RNDN);
383  b_ok = mpfr_zero_p(t01) && mpfr_zero_p(t02);
384 
385  // (t5, t6) = E (A-w)
386  mpfr_mul(t5, Ex, Ax, MPFR_RNDN);
387  mpfr_mul(t6, Ey, Ay, MPFR_RNDN);
388  mpfr_sub(t5, t5, t6, MPFR_RNDN);
389  mpfr_mul(t6, Ex, Ay, MPFR_RNDN);
390  mpfr_mul(t7, Ey, Ax, MPFR_RNDN);
391  mpfr_add(t6, t6, t7, MPFR_RNDN);
392  // B += 1
393  mpfr_add_ui(Bx, Bx, 1, MPFR_RNDN);
394  // (t7, t8) = D (B+1)
395  mpfr_mul(t7, Dx, Bx, MPFR_RNDN);
396  mpfr_mul(t8, Dy, By, MPFR_RNDN);
397  mpfr_sub(t7, t7, t8, MPFR_RNDN);
398  mpfr_mul(t8, Dx, By, MPFR_RNDN);
399  mpfr_mul(t9, Dy, Bx, MPFR_RNDN);
400  mpfr_add(t8, t8, t9, MPFR_RNDN);
401  // (t5, t6) -= (t7, t8)
402  mpfr_sub(t5, t5, t7, MPFR_RNDN);
403  mpfr_sub(t6, t6, t8, MPFR_RNDN);
404  // (t7, t8) = (t5, t6) * (t3, t4)
405  mpfr_mul(t7, t3, t5, MPFR_RNDN);
406  mpfr_mul(t8, t4, t6, MPFR_RNDN);
407  mpfr_sub(t7, t7, t8, MPFR_RNDN);
408  mpfr_mul(t8, t3, t6, MPFR_RNDN);
409  mpfr_mul(t9, t4, t5, MPFR_RNDN);
410  mpfr_add(t8, t8, t9, MPFR_RNDN);
411 
412  // (t3, t4) = w - (t7, t8)
413  mpfr_sub(t3, wx, t7, MPFR_RNDN);
414  mpfr_sub(t4, wy, t8, MPFR_RNDN);
415 
416  mpfr_set(t01, t3, MPFR_RNDN);
417  mpfr_set(t02, t4, MPFR_RNDN);
418  mpfr_set(t03, wx, MPFR_RNDN);
419  mpfr_set(t04, wy, MPFR_RNDN);
420  mpfr_sub(t01, t01, t03, MPFR_RNDN);
421  mpfr_sub(t02, t02, t04, MPFR_RNDN);
422  mpfr_add_ui(t01, t01, 1, MPFR_RNDN);
423  mpfr_sub_ui(t01, t01, 1, MPFR_RNDN);
424  mpfr_add_ui(t02, t02, 1, MPFR_RNDN);
425  mpfr_sub_ui(t02, t02, 1, MPFR_RNDN);
426  w_ok = mpfr_zero_p(t01) && mpfr_zero_p(t02);
427 
428  mpfr_set(bx, t1, MPFR_RNDN);
429  mpfr_set(by, t2, MPFR_RNDN);
430  mpfr_set(wx, t3, MPFR_RNDN);
431  mpfr_set(wy, t4, MPFR_RNDN);
432  } while (!(w_ok && b_ok) && --limit);
433  if (! limit) goto cleanup;
434  // t1 = |n - b|
435  mpfr_sub(t1, nx, bx, MPFR_RNDN);
436  mpfr_sqr(t1, t1, MPFR_RNDN);
437  mpfr_sub(t2, ny, by, MPFR_RNDN);
438  mpfr_sqr(t2, t2, MPFR_RNDN);
439  mpfr_add(t1, t1, t2, MPFR_RNDN);
440  mpfr_sqrt(t1, t1, MPFR_RNDN);
441  bits = 2 * bits;
442  mpfr_prec_round(nx, bits, MPFR_RNDN);
443  mpfr_prec_round(ny, bits, MPFR_RNDN);
444  mpfr_prec_round(wx, bits, MPFR_RNDN);
445  mpfr_prec_round(wy, bits, MPFR_RNDN);
446  mpfr_prec_round(bx, bits, MPFR_RNDN);
447  mpfr_prec_round(by, bits, MPFR_RNDN);
448  mpfr_set(zz, t1, MPFR_RNDN);
449  mpfr_set_prec(Ax, bits);
450  mpfr_set_prec(Ay, bits);
451  mpfr_set_prec(Bx, bits);
452  mpfr_set_prec(By, bits);
453  mpfr_set_prec(Cx, bits);
454  mpfr_set_prec(Cy, bits);
455  mpfr_set_prec(Dx, bits);
456  mpfr_set_prec(Dy, bits);
457  mpfr_set_prec(Ex, bits);
458  mpfr_set_prec(Ey, bits);
459  mpfr_set_prec(t1, bits);
460  mpfr_set_prec(t2, bits);
461  mpfr_set_prec(t3, bits);
462  mpfr_set_prec(t4, bits);
463  mpfr_set_prec(t5, bits);
464  mpfr_set_prec(t6, bits);
465  mpfr_set_prec(t7, bits);
466  mpfr_set_prec(t8, bits);
467  mpfr_set_prec(t9, bits);
468  } while (mpfr_zero_p(zz) || bits + mpfr_get_exp(zz) < mpfr_get_prec(zz) + accuracy);
469 
470  // copy to output
471  r = 1;
472 cleanup:
473  bits = mpfr_get_prec(nx);
474  mpfr_set_prec(x, bits);
475  mpfr_set_prec(y, bits);
476  mpfr_set(x, nx, MPFR_RNDN);
477  mpfr_set(y, ny, MPFR_RNDN);
478  mpfr_set(z, zz, MPFR_RNDN);
479  mpfr_clears(VARS, (mpfr_ptr) 0);
480 #undef VARS
481  return r;
482 }