Fractals/mandelbrot-numerics

mandelbrot-numerics is a library for advanced numeric calculations related with Mandelbrot set by Claude Heiland-Allen. This is unofficial wiki about it, containing most of maximus/book 'lib' and 'bin' (but no rendering) See: dictionary notation fork of the library Libraries by Claude Heiland-Allen exrtact various tools to manipulate EXR images compare with OpenExr kf-extras programs for manipulating output from Kalles Fraktaler 2 mandelbrot-symbolics - symbolic algorithms related to the Mandelbrot set mandelbrot-graphics - CPU-based visualisation of the Mandelbrot set mandelbrot-text - parsing and pretty printing related to the Mandelbrot set Notation directories bin for programs lib for functions with 2 versions d for double precision r for arbitrary precision include precision Precision here means precision of floating point numbers in the numerical computing. It is a positive integer. It is a number of bits of binary number ( Function names ( notation) double precision: m_d_*() functions taking pointers to double _Complex use them for output arbitrary precision: m_r_*() how to find precision ? start with minimal precision ( double = 53 bits ) compute new precision = 53 + f(parameter) where parameter is a feature of the procedure ( for example minimal size of the component for m-interior ). In other words new precision is grater then 53 Installation " you need to "make install" mandelbrot-symbolics lib before trying anything with mandelbrot-numerics" dependencies pkg-config gmp (libgmp) mpfr (libmpfr) mpc (libmpc) haskell - ghc sndfile using : sudo apt-get install libsndfile1-devaskubuntu question: package-was-not-found-in-the-pkg-config-search-path #include <complex.h> #include <math.h> #include <stdbool.h> #include <stdint.h> #include <stdlib.h> #include <gmp.h> #include <mpfr.h> #include <mpc.h> Gcc flags : gcc -lmpc -lmpfr -lgmp -lm From console : sudo aptitude install libmpc-dev git clone From directory in which you want to clone git repo: git clone https://code.mathr.co.uk/mandelbrot-numerics.git compile and from the same directory : make -C mandelbrot-numerics/c/lib prefix=${HOME}/opt install make -C mandelbrot-numerics/c/bin prefix=${HOME}/opt install then to run do:stackoverflow question: linux-error-while-loading-shared-libraries-cannot-open-shared-object-file-no-s export LD_LIBRARY_PATH=${HOME}/opt/lib check : echo $LD_LIBRARY_PATH result : /home/a/opt/lib or export PATH=${HOME}/opt/bin:${PATH} and check : echo $PATH To set it permanently change file .profilestackoverflow question how-to-permanently-set-path-on-linux sudo gedit ~/.profile test m-describe 53 20 100 0 0 4 LD_LIBRARY_PATH=${HOME}/opt/lib m-primary-separators update git From console opened in the mandelbrot-numerics directory : git pull If you made some local changes you can undu them : git checkout -f then git pull Now install again structures //mandelbrot-numerics.h struct m_d_mat2 { double _Complex a, b, c, d; }; typedef struct m_d_mat2 m_d_mat2; how to use it to include library c source should *only* have #include <mandelbrot-numerics.h> compile and link with pkg-config: see mandelbrot-numerics/c/bin/Makefile for an example quickest way to get started is to just put your file in mandelbrot-numerics/c/bin and run make Misiurewicz point m-feature-database m-misiurewicz m-misiurewicz usage: m-misiurewicz precision guess-re guess-im preperiod period maxsteps Description uses preperiod of critical point Examples using double precision m-misiurewicz double -2.1 0 2 1 100 -2.0000000000000000e+00 0.0000000000000000e+00 m-misiurewicz double -1.6 0 3 1 100 -1.5436890126920764e+00 0.0000000000000000e+00 m-misiurewicz double -1.5 0 5 2 100 -1.4303576324513072e+00 0.0000000000000000e+00 m-misiurewicz double -1.4 0 9 4 100 -1.4074051181647020e+00 0.0000000000000000e+00 Using bits : ./m-misiurewicz 200 -1.4 0 9 4 100 -1.4074051181647020225078282291990509777838059260945479350908287e+00 0.0000000000000000000000000000000000000000000000000000000000000e+00 m-feature-database m-feature-database usage: ./m-feature-database preperiod period m-feature-database 1 1 .1(0) 1/2 -2.00 0.00 m-describe m-describe official page this procedure gives description of point c from the parameter plane. Use without parameters: m-describe gives usage description: usage: m-describe precision maxperiod maxiters re im ncpus Check: point c = 0 = 0 +0*I= re+ im*I with precision = 53 bits max period= 20 max iters = 100 m-describe 53 20 100 0 0 4 output is : the input point was +0e+00 + +0e+00 i the point didn't escape after 100 iterations nearby hyperbolic components to the input point: - a period 1 cardioid with nucleus at +0e+00 + +0e+00 i the component has size 1.00000e+00 and is pointing west the atom domain has size 0.00000e+00 the atom domain coordinates of the input point are -nan + -nan i the atom domain coordinates in polar form are nan to the east the nucleus is 0.00000e+00 to the east of the input point the input point is interior to this component at radius 0.00000e+00 and angle 0.000000000000000000 (in turns) the multiplier is +0.00000e+00 + +0.00000e+00 i a point in the attractor is +0e+00 + +0e+00 i external angles of this component are: .(0) .(1) Another example: m-describe 1000 1000 10000000 -1.749512 0.0 8 the input point was -1.74951199999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999e+00 + +0e+00 i the point didn't escape after 10000000 iterations nearby hyperbolic components to the input point: - a period 1 cardioid with nucleus at +0e+00 + +0e+00 i the component has size 1.00000e+00 and is pointing west the atom domain has size 0.00000e+00 the atom domain coordinates of the input point are -nan + -nan i the atom domain coordinates in polar form are nan to the east the nucleus is 1.74951e+00 to the east of the input point the input point is exterior to this component at radius 1.82808e+00 and angle 0.500000000000000000 (in turns) the multiplier is -1.82808e+00 + +0.00000e+00 i a point in the attractor is -9.14032e-01 + +0e+00 i external angles of this component are: .(0) .(1) - a period 2 circle with nucleus at -1e+00 + +0e+00 i the component has size 5.00000e-01 and is pointing west the atom domain has size 1.00000e+00 the atom domain coordinates of the input point are -0.74951 + -0 i the atom domain coordinates in polar form are 0.74951 to the west the nucleus is 7.49512e-01 to the east of the input point the input point is exterior to this component at radius 2.99805e+00 and angle 0.500000000000000000 (in turns) the multiplier is -2.99805e+00 + +0.00000e+00 i a point in the attractor is -1.49976e+00 + +0e+00 i - a period 3 cardioid with nucleus at -1.754878e+00 + +0e+00 i the component has size 1.90355e-02 and is pointing west the atom domain has size 2.34487e-01 the atom domain coordinates of the input point are -0.022921 + +0 i the atom domain coordinates in polar form are 0.022921 to the west the nucleus is 5.36567e-03 to the west of the input point the input point is exterior to this component at radius 1.24010e+00 and angle 0.000000000000000000 (in turns) the multiplier is +1.24010e+00 + +0.00000e+00 i a point in the attractor is -6.8613231e-02 + +0e+00 i external angles of this component are: .(011) .(100) - a period 237 cardioid with nucleus at -1.7495120000000000000000000000000116053893024686343718029506670414e+00 + +0e+00 i the component has size 5.10303e-60 and is pointing west the atom domain has size 1.34480e-32 the atom domain coordinates of the input point are -0.86149 + -0 i the atom domain coordinates in polar form are 0.86149 to the west the nucleus is 1.16054e-32 to the west of the input point the input point is exterior to this component at radius 5.87695e+33 and angle 0.500000000000000000 (in turns) the multiplier is -5.87695e+33 + +0.00000e+00 i a point in the attractor is +2.5892947214253619092848904520135803680047966975848204640792969525e-02 + +0e+00 i - a period 756 cardioid with nucleus at -1.74951200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000437443253498601125336414441820009561885311706078736706424485680024027767522121183394207260374919267049e+00 + +0e+00 i the component has size 2.81459e-198 and is pointing west the atom domain has size 8.60400e-102 the atom domain coordinates of the input point are -0.51039 + +0 i the atom domain coordinates in polar form are 0.51039 to the west the nucleus is 4.37443e-102 to the west of the input point the input point is exterior to this component at radius 1.99523e+68 and angle 0.500000000000000000 (in turns) the multiplier is -1.99523e+68 + +0.00000e+00 i a point in the attractor is -1.32154302717027851153258720563184483110950266513429522776342823541752294814605666085090894143487158886716259121852162061940587472967643029955949630872182624723698360235112807868903937707485681031295828661e-02 + +0e+00 i shape Versions m_d.shape.c m_r_shape.c Descriptionmath.stackexchange question: perfect-circles-in-the-mandelbrot-set size atom domain Atom domain size estimateatom_domain_size_estimation by Claude Heiland-Allen deriving_the_size_estimate by Claude Heiland-Allen formula from http://ibiblio.org/e-notes/MSet/windows.htm <-- where I got the formula from component size Files: m-size.c m_d_size.c output: complex double size m_r_size.c output : mpc_t size // mandelbrot-numerics -- numerical algorithms related to the Mandelbrot set // Copyright (C) 2015-2018 Claude Heiland-Allen // License GPL3+ http://www.gnu.org/licenses/gpl.html #include <mandelbrot-numerics.h> extern double _Complex m_d_size(double _Complex nucleus, int period) { double _Complex l = 1; double _Complex b = 1; double _Complex z = 0; for (int i = 1; i < period; ++i) { z = z * z + nucleus; l = 2 * z * l; b = b + 1 / l; } return 1 / (b * l * l); } Use without parameters: m-size gives usage description : usage: ./m-size precision nucleus-re nucleus-im period Parameters: precisionon_the_precision_required_for_size_estimates by Claude Heiland-Allen nucleus ( complex number c ) period positive integer string "double" Results are 2 numbers is: cabs(size) = size and carg(size) = an orientation estimate mpc_abs(r, size, MPFR_RNDN); mpc_arg(t, size, MPFR_RNDN); It gives the size as a complex number, the magnitude is the size relative to the period 1 continent and the phase / argument is the angle of rotation. Code: $ m-size 53 -8.6915874972342078e-01 2.5565708568620021e-01 48 1.1864737161932281e-08 -5.9691937849660148e-01 my program splits it into magnitude and phase (in radians) for convenience, so the minibrot is 1.1865e-8 times the size of the period 1 continent, so multiply the view size of your initial view by this number to get the minibrot to appear about the same size. Example usage: m-size double 0 0 1 output : 1.0000000000000000e+00 0.0000000000000000e+00 In other words: the component has size 1.00000e+00 and is pointing west Minibrot period 48fractalforums.com : deepest-e10000 m-size 53 -8.6915874972342078e-01 2.5565708568620021e-01 48 1.1864737161932281e-08 -5.9691937849660148e-01 period box-period Computing period under complex quadratic polynomial "But, the box period method detects the period of a *nucleus* within the box - if there is no nucleus there it will fail." Find the period of the central minibrot. You can do this by iterating the 4 c corners of a box around the view, and seeing at which iteration count the 4 z points surround the origin. See http://www.mrob.com/pub/muency/period.html for a more detailed description of the method, I have an implementation which you can read here: http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/c/lib/m_d_box_period.c (there's also an arbitrary precision version in the same directory) See: http://www.mrob.com/pub/muency/period.html binaries m-box-period The result : usage: m-box-period precision center-re center-im radius maxperiod Example: m-box-period 53 -0.8691524744 0.2556487868 1.25e-5 1000 48 code This is one file program using code extracted from Claude's library. /* http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/c/bin/m-box-period.c numerical algorithms related to the Mandelbrot set by Claude Heiland-Allen https://mathr.co.uk/blog/ GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 http://code.mathr.co.uk/mandelbrot-numerics/blob/HEAD:/COPYING mandelbrot-numerics Numerical algorithms related to the Mandelbrot set: ray tracing, nucleus location, bond points, etc. --------------------- Dependencies sudo aptitude install libmpc-dev --------------------------------- https://gitlab.com/adammajewski/mandelbrot-numerics-box-periood Existing folder or Git repository cd existing_folder git init git remote add origin git@gitlab.com:adammajewski/mandelbrot-numerics-box-periood.git git add p.c git commit -m "sasa" git push -u origin master ------------------- gcc p.c -lmpc -lmpfr -lgmp -Wall ./a.out precision center-re center-im radius maxperiod ./a.out 100 0.37496784 0.21687214 1.0 20 */ #include <stdio.h> #include <stdint.h> #include <stdlib.h> #include <errno.h> #include <string.h> #include <math.h> #include <complex.h> #include <stdbool.h> #include <gmp.h> #include <mpfr.h> #include <mpc.h> // static const double twopi = 6.283185307179586; /* --------------- http://code.mathr.co.uk/mandelbrot-numerics/blob_plain/HEAD:/c/lib/m_d_util.h ---- */ static inline int sgn(double z) { if (z > 0) { return 1; } if (z < 0) { return -1; } return 0; } static inline bool odd(int a) { return a & 1; } static inline double cabs2(double _Complex z) { return creal(z) * creal(z) + cimag(z) * cimag(z); } static inline bool cisfinite(double _Complex z) { return isfinite(creal(z)) && isfinite(cimag(z)); } static const double pi = 3.141592653589793; static const double twopi = 6.283185307179586; // last . takeWhile (\x -> 2 /= 2 + x) . iterate (/2) $ 1 :: Double static const double epsilon = 4.440892098500626e-16; // epsilon^2 static const double epsilon2 = 1.9721522630525295e-31; /* --------------- utils.h -------------------------------------- */ static inline bool arg_precision(const char *arg, bool *native, int *bits) { if (0 == strcmp("double", arg)) { *native = true; *bits = 53; return true; } else { char *check = 0; errno = 0; long int li = strtol(arg, &check, 10); bool valid = ! errno && arg != check && ! *check; int i = li; if (valid && i > 1) { *native = false; *bits = i; return true; } } return false; } static inline bool arg_double(const char *arg, double *x) { char *check = 0; errno = 0; double d = strtod(arg, &check); if (! errno && arg != check && ! *check) { *x = d; return true; } return false; } static inline bool arg_int(const char *arg, int *x) { char *check = 0; errno = 0; long int li = strtol(arg, &check, 10); if (! errno && arg != check && ! *check) { *x = li; return true; } return false; } static inline bool arg_rational(const char *arg, mpq_t x) { int ok = mpq_set_str(x, arg, 10); mpq_canonicalize(x); return ok == 0; } static inline bool arg_mpfr(const char *arg, mpfr_t x) { return 0 == mpfr_set_str(x, arg, 10, MPFR_RNDN); } static inline bool arg_mpc(const char *re, const char *im, mpc_t x) { int ok = mpfr_set_str(mpc_realref(x), re, 10, MPFR_RNDN) + mpfr_set_str(mpc_imagref(x), im, 10, MPFR_RNDN); return ok == 0; } /* -----------c/lib/m_d_box_period.c -----------------------*/ static double cross(double _Complex a, double _Complex b) { return cimag(a) * creal(b) - creal(a) * cimag(b); } static bool crosses_positive_real_axis(double _Complex a, double _Complex b) { if (sgn(cimag(a)) != sgn(cimag(b))) { double _Complex d = b - a; int s = sgn(cimag(d)); int t = sgn(cross(d, a)); return s == t; } return false; } static bool surrounds_origin(double _Complex a, double _Complex b, double _Complex c, double _Complex d) { return odd ( crosses_positive_real_axis(a, b) + crosses_positive_real_axis(b, c) + crosses_positive_real_axis(c, d) + crosses_positive_real_axis(d, a) ); } typedef struct m_d_box_period { double _Complex c[4]; double _Complex z[4]; int p; } m_d_box_period_t ; m_d_box_period_t *m_d_box_period_new(double _Complex center, double radius) { m_d_box_period_t *box = (m_d_box_period_t *) malloc(sizeof(*box)); if (! box) { return 0; } box->z[0] = box->c[0] = center + ((-radius) + I * (-radius)); box->z[1] = box->c[1] = center + (( radius) + I * (-radius)); box->z[2] = box->c[2] = center + (( radius) + I * ( radius)); box->z[3] = box->c[3] = center + ((-radius) + I * ( radius)); box->p = 1; return box; } void m_d_box_period_delete(m_d_box_period_t *box) { if (box) { free(box); } } extern bool m_d_box_period_step(m_d_box_period_t *box) { if (! box) { return false; } bool ok = true; for (int i = 0; i < 4; ++i) { box->z[i] = box->z[i] * box->z[i] + box->c[i]; ok = ok && cisfinite(box->z[i]); } box->p = box->p + 1; return ok; } extern bool m_d_box_period_have_period(const m_d_box_period_t *box) { if (! box) { return true; } return surrounds_origin(box->z[0], box->z[1], box->z[2], box->z[3]); } extern int m_d_box_period_get_period(const m_d_box_period_t *box) { if (! box) { return 0; } return box->p; } extern int m_d_box_period_do(double _Complex center, double radius, int maxperiod) { m_d_box_period_t *box = m_d_box_period_new(center, radius); if (! box) { return 0; } int period = 0; for (int i = 0; i < maxperiod; ++i) { if (m_d_box_period_have_period(box)) { period = m_d_box_period_get_period(box); break; } if (! m_d_box_period_step(box)) { break; } } m_d_box_period_delete(box); return period; } /* -----------http://code.mathr.co.uk/mandelbrot-numerics/blob_plain/HEAD:/c/lib/m_r_box_period.c -----------------------*/ static void cross_r(mpfr_t out, const mpc_t a, const mpc_t b, mpfr_t t) { mpfr_mul(out, mpc_imagref(a), mpc_realref(b), MPFR_RNDN); mpfr_mul(t, mpc_realref(a), mpc_imagref(b), MPFR_RNDN); mpfr_sub(out, out, t, MPFR_RNDN); } static bool crosses_positive_real_axis_r(const mpc_t a, const mpc_t b, mpc_t d, mpfr_t t1, mpfr_t t2) { if (mpfr_sgn(mpc_imagref(a)) != mpfr_sgn(mpc_imagref(b))) { // d = b - a; mpc_sub(d, b, a, MPC_RNDNN); // s = sgn(cimag(d)); int s = mpfr_sgn(mpc_imagref(d)); // t = sgn(cross(d, a)); cross_r(t1, d, a, t2); int t = mpfr_sgn(t1); return s == t; } return false; } static bool surrounds_origin_r(const mpc_t a, const mpc_t b, const mpc_t c, const mpc_t d, mpc_t e, mpfr_t t1, mpfr_t t2) { return odd ( crosses_positive_real_axis_r(a, b, e, t1, t2) + crosses_positive_real_axis_r(b, c, e, t1, t2) + crosses_positive_real_axis_r(c, d, e, t1, t2) + crosses_positive_real_axis_r(d, a, e, t1, t2) ); } typedef struct m_r_box_period { mpc_t c[4]; mpc_t z[4]; int p; mpc_t e; mpfr_t t1; mpfr_t t2; } m_r_box_period_t; extern m_r_box_period_t *m_r_box_period_new(const mpc_t center, const mpfr_t radius) { m_r_box_period_t *box = (m_r_box_period_t *) malloc(sizeof(*box)); if (! box) { return 0; } // prec mpfr_prec_t precr, preci, prec; mpc_get_prec2(&precr, &preci, center); prec = precr > preci ? precr : preci; // init for (int i = 0; i < 4; ++i) { mpc_init2(box->c[i], prec); mpc_init2(box->z[i], prec); } mpc_init2(box->e, prec); mpfr_init2(box->t1, prec); mpfr_init2(box->t2, prec); // box mpfr_sub(mpc_realref(box->c[0]), mpc_realref(center), radius, MPFR_RNDN); mpfr_sub(mpc_imagref(box->c[0]), mpc_imagref(center), radius, MPFR_RNDN); mpfr_add(mpc_realref(box->c[1]), mpc_realref(center), radius, MPFR_RNDN); mpfr_sub(mpc_imagref(box->c[1]), mpc_imagref(center), radius, MPFR_RNDN); mpfr_add(mpc_realref(box->c[2]), mpc_realref(center), radius, MPFR_RNDN); mpfr_add(mpc_imagref(box->c[2]), mpc_imagref(center), radius, MPFR_RNDN); mpfr_sub(mpc_realref(box->c[3]), mpc_realref(center), radius, MPFR_RNDN); mpfr_add(mpc_imagref(box->c[3]), mpc_imagref(center), radius, MPFR_RNDN); mpc_set(box->z[0], box->c[0], MPC_RNDNN); mpc_set(box->z[1], box->c[1], MPC_RNDNN); mpc_set(box->z[2], box->c[2], MPC_RNDNN); mpc_set(box->z[3], box->c[3], MPC_RNDNN); box->p = 1; return box; } extern void m_r_box_period_delete(m_r_box_period_t *box) { if (box) { for (int i = 0; i < 4; ++i) { mpc_clear(box->c[i]); mpc_clear(box->z[i]); } mpc_clear(box->e); mpfr_clear(box->t1); mpfr_clear(box->t2); free(box); } } extern bool m_r_box_period_step(m_r_box_period_t *box) { if (! box) { return false; } bool ok = true; for (int i = 0; i < 4; ++i) { // box->z[i] = box->z[i] * box->z[i] + box->c[i]; mpc_sqr(box->z[i], box->z[i], MPC_RNDNN); mpc_add(box->z[i], box->z[i], box->c[i], MPC_RNDNN); ok = ok && mpfr_number_p(mpc_realref(box->z[i])) && mpfr_number_p(mpc_imagref(box->z[i])); } box->p = box->p + 1; return ok; } extern bool m_r_box_period_have_period(const m_r_box_period_t *box) { if (! box) { return true; } m_r_box_period_t *ubox = (m_r_box_period_t *) box; // const cast return surrounds_origin_r(box->z[0], box->z[1], box->z[2], box->z[3], ubox->e, ubox->t1, ubox->t2); } extern int m_r_box_period_get_period(const m_r_box_period_t *box) { if (! box) { return 0; } return box->p; } extern int m_r_box_period_do(const mpc_t center, const mpfr_t radius, int maxperiod) { m_r_box_period_t *box = m_r_box_period_new(center, radius); if (! box) { return 0; } int period = 0; for (int i = 0; i < maxperiod; ++i) { if (m_r_box_period_have_period(box)) { period = m_r_box_period_get_period(box); break; } if (! m_r_box_period_step(box)) { break; } } m_r_box_period_delete(box); return period; } /*--------------- */ /* c:0.37496784+%i*0.21687214; ./a.out 100 0.37496784 0.21687214 1.0 20 ./a.out 200 -1.121550281113895 +0.265176187855967 0.005 40 so radius = domain size */ static void usage(const char *progname) { fprintf ( stderr , "usage: %s precision center-re center-im radius maxperiod\n" , progname ); } /* ========================= mAIN =========================== */ extern int main(int argc, char **argv) { if (argc != 6) { usage(argv[0]); return 1; } bool native = true; int bits = 0; if (! arg_precision(argv[1], &native, &bits)) { return 1; } if (native) { double cre = 0; double cim = 0; double radius = 0; int maxperiod = 0; if (! arg_double(argv[2], &cre)) { return 1; } if (! arg_double(argv[3], &cim)) { return 1; } if (! arg_double(argv[4], &radius)) { return 1; } if (! arg_int(argv[5], &maxperiod)) { return 1; } int period = m_d_box_period_do(cre + I * cim, radius, maxperiod); if (period > 0) { printf("%d\n", period); return 0; } } else { mpc_t center; mpfr_t radius; int maxperiod = 0; mpc_init2(center, bits); mpfr_init2(radius, bits); if (! arg_mpc(argv[2], argv[3], center)) { return 1; } if (! arg_mpfr(argv[4], radius)) { return 1; } if (! arg_int(argv[5], &maxperiod)) { return 1; } int period = m_r_box_period_do(center, radius, maxperiod); if (period > 0) { printf("%d\n", period); } mpc_clear(center); mpfr_clear(radius); return period <= 0; } return 1; } m-interior Interior coordinates in the Mandelbrot set Steps: C is fixed throughout this process - it determines the pixel we want to colour. The first step is to find a special z (maybe call it w) such that p iterations of z to z^2 + c starting from w, gets you back to w. This is the limit cycle attractor of c. This is done by m_d_attractor, using Newton's method to solve f_c^p(z) - z = 0. Then to find the interior coordinate for c, take your w, and iterate it p times, calculating the derivative w.r.t. z: z = w; dz = 1; for (int i = 0; i < p; ++i) dz = 2 * z * dz z = z * z + c then the interior coordinate is simply dz after the loop. if |dz| < 1 then c is interior to a component of period p you can colour pixels based on dz (eg angle -> hue, radius -> saturation) for highlighting the edge you can optionally use the w to find the interior distance estimate (m_d_interior_de). Description from math.stackexchangemath.stackexchange question test-for-membership-in-mandelbrot-bulb-of-period-n/1151953#1151953math.stackexchange question: what-is-my-method-of-rendering-the-mandelbrot-set-called blog interior_coordinates_in_the_mandelbrot_set "traced boundaries using Newton's method in two variables" theory an algorithm to find points on the boundary of the Mandelbrot set, given a particular hyperbolic component and the desired internal angle. It involves Newton's method in two complex variables to solve system of 2 equations: where: F is function ( complex quadratic polynomial) is iterated polynomial p is the period of the target component with |r|≤1. θ is the desired internal angle. First equation describes periodic point second equation describes attractive periodic point ( magnitude r of multiplier b is not grater then 1 ) The resulting c is the coordinates of the point on the boundary. It can also be modified to find points in the interior, simply set b=re2πiθ with |r|≤1. The algorithm for finding numerically interior coordinate b by Claude: When c is outside the Mandelbrot set, give up now; For each period p, starting from 1 and increasing: Find z0 such that Fp(z0,c)=z0 using Newton's method in one complex variable; Find b by evaluating first derivative with respect to z at z0; If |b|≤1 then return b, otherwise continue with the next p. usage complex double z = 0; complex double c = 0; m_d_interior(&z, &c, zre + I * zim, cre + I * cim, ir * cexp(I * twopi * it), period, maxsteps); printf("z = %.16e%+.16e\t c = %.16e%+.16e\n", creal(z), cimag(z), creal(c), cimag(c)); Use without parameters: m-interior gives usage usage: ./m-interior precision z-guess-re z-guess-im c-guess-re c-guess-im interior-r interior-t period maxsteps + Input parameter variable description precision int bits number of binary digits = width of significant = precision. Here one can use : double or 53 ( for double numbers) positive integer = number of bits for MPFR z-guess-re double zre real part of the initial guess z_0 for Newton method, use z = 0 z-guess-im double zim c-guess-re double cre real part of the initial guess c_0 for Newton method, use nucleus of given hyperbolic component c-guess-im double cim imaginary part of c_0 interior-r double ir Example interior-t double it Example period int period Period of given hyperbolic component maxsteps int maxsteps Example z-guess = z-guess-re+z-guess-im*I c-guess = c-guess-re + c-guess=im*I Output: point c with given internal coordinate from hyperbolic componenet described by period and nucleus periodic point z (wucleus) Example : ./a.out double 0 0 0 0 1 0 1 100 5.0000000000000000e-01 0.0000000000000000e+00 2.5000000000000000e-01 0.0000000000000000e+00 How to read it : z = 0.5 ( fixed point for c = 1/4) c = 0.25 General algorithm: choose period compute / choose center c from period = initial aproximation compute size of component with center c. Compute precision from size choose t = p/q and r compute code Program m-interior from bin directory Source code : mandelbrot-numerics/c/bin/ m-interior.c mandelbrot-numerics/c/lib/m_d_interior.c mandelbrot-numerics/c/lib/m_r_interior.c Functions: m_d_interior(&z, &c, zre + I * zim, cre + I * cim, ir * cexp(I * twopi * it), period, maxsteps); m_r_interior(z, c, z, c, interior, period, maxsteps); #include <mandelbrot-numerics.h> #include "m_d_util.h" extern m_newton m_d_interior_step(double _Complex *z_out, double _Complex *c_out, double _Complex z_guess, double _Complex c_guess, double _Complex interior, int period) { double _Complex c = c_guess; double _Complex z = z_guess; double _Complex dz = 1; double _Complex dc = 0; double _Complex dzdz = 0; double _Complex dcdz = 0; for (int p = 0; p < period; ++p) { dcdz = 2 * (z * dcdz + dc * dz); dzdz = 2 * (z * dzdz + dz * dz); dc = 2 * z * dc + 1; dz = 2 * z * dz; z = z * z + c; } double _Complex det = (dz - 1) * dcdz - dc * dzdz; double _Complex z_new = z_guess - (dcdz * (z - z_guess) - dc * (dz - interior)) / det; double _Complex c_new = c_guess - ((dz - 1) * (dz - interior) - dzdz * (z - z_guess)) / det; if (cisfinite(z_new) && cisfinite(c_new)) { *z_out = z_new; *c_out = c_new; if (cabs2(z_new - z_guess) <= epsilon2 && cabs2(c_new - c_guess) <= epsilon2) { return m_converged; } else { return m_stepped; } } else { *z_out = z_guess; *c_out = c_guess; return m_failed; } } extern m_newton m_d_interior(double _Complex *z_out, double _Complex *c_out, double _Complex z_guess, double _Complex c_guess, double _Complex interior, int period, int maxsteps) { m_newton result = m_failed; double _Complex z = z_guess; double _Complex c = c_guess; for (int i = 0; i < maxsteps; ++i) { if (m_stepped != (result = m_d_interior_step(&z, &c, z, c, interior, period))) { break; } } *z_out = z; *c_out = c; return result; } m-interior-de descriptionpractical_interior_distance_rendering m-attractor to find a special z (maybe call it w) such that p iterations of z to z^2 + c starting from w, gets you back to w. This is the limit cycle attractor of c. It is done using Newton's method to solve f_c^p(z) - z = 0. usage: %s precision z-guess-re z-guess-im c-re c-im period maxsteps Finds one periodic z point ( point on dynamical plane ) near z-guess for : period c= c-re + c-im* i using Newton method with : z-guess = z-guess-re + z-guess-im* i ( initial aproximation of the root) maxima number of steps = maxsteps using number with width of significant = precision. Here one can use : double or 53 positive integer = number of bits for MPFR ./m-attractor double 0 0 -0.965 0.085 2 100 or ./m-attractor 53 0 0 -0.965 0.085 2 100 result : -2.7669310528133408e-02 -8.9979332165608591e-02 The precision is set manually in m-attractor to let user experiment. m-exray-in Parameter external ray:mathoverflow question: algorithm-for-computing-external-angles-for-the-mandelbrot-set traced inward ( from infinity to the boundary ) using Newton method When tracing inwards, one peels off the most-significant bit (aka angle doubling) each time the ray crosses a dwell band (integer part of normalized iteration count increases by 1). binaries m-exray-in usage: m-exray-in precision angle sharpness maxsteps m-exray-in double 1/3 4 200 m-exray-in double 0 4 200 2.5405429038323935e-01 0.0000000000000000e+00 Code Source code: c/bin/m-exray-in.c = program showing usage of the library function c/lib/m_d_exray_in.c = double ( machine double precision ) version of procedure c/lib/m_r_exray_in.c = mpfr version of procedure = (dynamically changed as necessary) precision variables : native ( if true then use double= native precision) precision (can be double , see arg_precision ) double version : m_newton m_d_exray_in_step(m_d_exray_in *ray, int maxsteps) { if (ray->j >= ray->sharpness) { mpq_mul_2exp(ray->angle, ray->angle, 1); if (mpq_cmp_ui(ray->angle, 1, 1) >= 0) { mpq_sub(ray->angle, ray->angle, ray->one); } ray->k = ray->k + 1; ray->j = 0; } double r = pow(ray->er, pow(0.5, (ray->j + 0.5) / ray->sharpness)); double a = twopi * mpq_get_d(ray->angle); double _Complex target = r * (cos(a) + I * sin(a)); double _Complex c = ray->c; for (int i = 0; i < maxsteps; ++i) { double _Complex z = 0; double _Complex dc = 0; for (int p = 0; p <= ray->k; ++p) { dc = 2 * z * dc + 1; z = z * z + c; } double _Complex c_new = c - (z - target) / dc; double d2 = cabs2(c_new - c); if (cisfinite(c_new)) { c = c_new; } else { break; } if (d2 <= epsilon2) { break; } } ray->j = ray->j + 1; double d2 = cabs2(c - ray->c); if (d2 <= epsilon2) { ray->c = c; return m_converged; } if (cisfinite(c)) { ray->c = c; return m_stepped; } else { return m_failed; } } mpfr version : /* it is a c console program which - computes external parameter ray - for complex quadratic polynomial fc(z) = z^2+c - uses arbitrary precision ( mpfr) with dynamic precision adjustment - uses Newton method ( described by T Kawahira) http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf gcc -o rayin e.c -std=c99 -Wall -Wextra -pedantic -lmpfr -lgmp -lm ./rayin 1/3 25 ./rayin 1/3 25 >ray13.txt ./rayin 1/3 25 | ./rescale 53 53 -0.75 0 1.5 0 >rray13.txt program uses the code by Claude Heiland-Allen from https://gitorious.org/maximus/book/ " ray_in computes 4 points per dwell band, and ray_out currently computes 16. Moreover ray_in has dynamic precision adjustment and ray_out doesn't (so you need a high precision to get through the thin gaps) - fixing both of these is on my mental todo list. ray_out will always be slower though, as it needs to compute extra points when crossing dwell bands to accumulate angle bits (and to work correctly)." What algorithm do you use for drawing external ray ? "essentially the Newton's method one in http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf precision is automatic, effectively it checks how close dwell bands are together and uses higher precision when they are narrow and lower precision when they are wide. in general precision will increase along the ray as it gets closer to the boundary. " "periodic rays stop near the landing point at the current zoom level, manually retrace if you zoom in and find the gap annoying. tracing periodic rays is a lot slower than tracing preperiodic rays, especially if you zoom in to the cusp - suggestions welcome for improvements here (perhaps a larger radius than a single pixel for landing point nearness threshold?). " */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <gmp.h> #include <mpfr.h> // arbitrary precision const float pif = 3.14159265358979323846264338327950288419716939937510f; const double pi = 3.14159265358979323846264338327950288419716939937510; const long double pil = 3.14159265358979323846264338327950288419716939937510l; struct mandelbrot_external_ray_in { mpq_t angle; mpq_t one; unsigned int sharpness; // number of steps to take within each dwell band unsigned int precision; // delta needs this many bits of effective precision unsigned int accuracy; // epsilon is scaled relative to magnitude of delta double escape_radius; mpfr_t epsilon2; mpfr_t cx; mpfr_t cy; unsigned int j; unsigned int k; // temporaries mpfr_t zx, zy, dzx, dzy, ux, uy, vx, vy, ccx, ccy, cccx, cccy; }; // new : allocates memory, inits and sets variables // input angle = external angle in turns extern struct mandelbrot_external_ray_in *mandelbrot_external_ray_in_new(mpq_t angle) { struct mandelbrot_external_ray_in *r = malloc(sizeof(struct mandelbrot_external_ray_in)); mpq_init(r->angle); mpq_set(r->angle, angle); mpq_init(r->one); mpq_set_ui(r->one, 1, 1); r->sharpness = 4; r->precision = 4; r->accuracy = 4; r->escape_radius = 65536.0; mpfr_init2(r->epsilon2, 53); mpfr_set_ui(r->epsilon2, 1, GMP_RNDN); // external angle in radions double a = 2.0 * pi * mpq_get_d(r->angle); // initial value c = mpfr_init2(r->cx, 53); mpfr_init2(r->cy, 53); mpfr_set_d(r->cx, r->escape_radius * cos(a), GMP_RNDN); mpfr_set_d(r->cy, r->escape_radius * sin(a), GMP_RNDN); r->k = 0; r->j = 0; // initialize temporaries mpfr_inits2(53, r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0); return r; } // delete extern void mandelbrot_external_ray_in_delete(struct mandelbrot_external_ray_in *r) { mpfr_clear(r->epsilon2); mpfr_clear(r->cx); mpfr_clear(r->cy); mpq_clear(r->angle); mpq_clear(r->one); mpfr_clears(r->zx, r->zy, r->dzx, r->dzy, r->ux, r->uy, r->vx, r->vy, r->ccx, r->ccy, r->cccx, r->cccy, (mpfr_ptr) 0); free(r); } // step extern int mandelbrot_external_ray_in_step(struct mandelbrot_external_ray_in *r) { if (r->j >= r->sharpness) { mpq_mul_2exp(r->angle, r->angle, 1); if (mpq_cmp_ui(r->angle, 1, 1) >= 0) { mpq_sub(r->angle, r->angle, r->one); } r->k++; r->j = 0; } // r0 <- er ** ((1/2) ** ((j + 0.5)/sharpness)) // t0 <- r0 cis(2 * pi * angle) double r0 = pow(r->escape_radius, pow(0.5, (r->j + 0.5) / r->sharpness)); double a0 = 2.0 * pi * mpq_get_d(r->angle); double t0x = r0 * cos(a0); double t0y = r0 * sin(a0); // c <- r->c mpfr_set(r->ccx, r->cx, GMP_RNDN); mpfr_set(r->ccy, r->cy, GMP_RNDN); for (unsigned int i = 0; i < 64; ++i) { // FIXME arbitrary limit // z <- 0 // dz <- 0 mpfr_set_ui(r->zx, 0, GMP_RNDN); mpfr_set_ui(r->zy, 0, GMP_RNDN); mpfr_set_ui(r->dzx, 0, GMP_RNDN); mpfr_set_ui(r->dzy, 0, GMP_RNDN); // iterate for (unsigned int p = 0; p <= r->k; ++p) { // dz <- 2 z dz + 1 mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN); mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN); mpfr_mul(r->vx, r->zx, r->dzy, GMP_RNDN); mpfr_mul(r->vy, r->zy, r->dzx, GMP_RNDN); mpfr_sub(r->dzx, r->ux, r->uy, GMP_RNDN); mpfr_add(r->dzy, r->vx, r->vy, GMP_RNDN); mpfr_mul_2ui(r->dzx, r->dzx, 1, GMP_RNDN); mpfr_mul_2ui(r->dzy, r->dzy, 1, GMP_RNDN); mpfr_add_ui(r->dzx, r->dzx, 1, GMP_RNDN); // z <- z^2 + c mpfr_sqr(r->ux, r->zx, GMP_RNDN); mpfr_sqr(r->uy, r->zy, GMP_RNDN); mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN); mpfr_mul(r->vy, r->zx, r->zy, GMP_RNDN); mpfr_mul_2ui(r->vy, r->vy, 1, GMP_RNDN); mpfr_add(r->zx, r->vx, r->ccx, GMP_RNDN); mpfr_add(r->zy, r->vy, r->ccy, GMP_RNDN); } // c' <- c - (z - t0) / dz mpfr_sqr(r->ux, r->dzx, GMP_RNDN); mpfr_sqr(r->uy, r->dzy, GMP_RNDN); mpfr_add(r->vy, r->ux, r->uy, GMP_RNDN); mpfr_sub_d(r->zx, r->zx, t0x, GMP_RNDN); mpfr_sub_d(r->zy, r->zy, t0y, GMP_RNDN); mpfr_mul(r->ux, r->zx, r->dzx, GMP_RNDN); mpfr_mul(r->uy, r->zy, r->dzy, GMP_RNDN); mpfr_add(r->vx, r->ux, r->uy, GMP_RNDN); mpfr_div(r->ux, r->vx, r->vy, GMP_RNDN); mpfr_sub(r->cccx, r->ccx, r->ux, GMP_RNDN); mpfr_mul(r->ux, r->zy, r->dzx, GMP_RNDN); mpfr_mul(r->uy, r->zx, r->dzy, GMP_RNDN); mpfr_sub(r->vx, r->ux, r->uy, GMP_RNDN); mpfr_div(r->uy, r->vx, r->vy, GMP_RNDN); mpfr_sub(r->cccy, r->ccy, r->uy, GMP_RNDN); // delta^2 = |c' - c|^2 mpfr_sub(r->ux, r->cccx, r->ccx, GMP_RNDN); mpfr_sub(r->uy, r->cccy, r->ccy, GMP_RNDN); mpfr_sqr(r->vx, r->ux, GMP_RNDN); mpfr_sqr(r->vy, r->uy, GMP_RNDN); mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN); // enough_bits = 0 < 2 * prec(eps^2) + exponent eps^2 - 2 * precision int enough_bits = 0 < 2 * (mpfr_get_prec(r->epsilon2) - r->precision) + mpfr_get_exp(r->epsilon2); if (enough_bits) { // converged = delta^2 < eps^2 int converged = mpfr_less_p(r->ux, r->epsilon2); if (converged) { // eps^2 <- |c' - r->c|^2 >> (2 * accuracy) mpfr_sub(r->ux, r->cccx, r->cx, GMP_RNDN); mpfr_sub(r->uy, r->cccy, r->cy, GMP_RNDN); mpfr_sqr(r->vx, r->ux, GMP_RNDN); mpfr_sqr(r->vy, r->uy, GMP_RNDN); mpfr_add(r->ux, r->vx, r->vy, GMP_RNDN); mpfr_div_2ui(r->epsilon2, r->ux, 2 * r->accuracy, GMP_RNDN); // j <- j + 1 r->j = r->j + 1; // r->c <- c' mpfr_set(r->cx, r->cccx, GMP_RNDN); mpfr_set(r->cy, r->cccy, GMP_RNDN); return 1; } } else { // bump precision mpfr_prec_t prec = mpfr_get_prec(r->cx) + 32; mpfr_prec_round(r->cx, prec, GMP_RNDN); mpfr_prec_round(r->cy, prec, GMP_RNDN); mpfr_prec_round(r->epsilon2, prec, GMP_RNDN); mpfr_set_prec(r->ccx, prec); mpfr_set_prec(r->ccy, prec); mpfr_set_prec(r->cccx, prec); mpfr_set_prec(r->cccy, prec); mpfr_set_prec(r->zx, prec); mpfr_set_prec(r->zy, prec); mpfr_set_prec(r->dzx, prec); mpfr_set_prec(r->dzy, prec); mpfr_set_prec(r->ux, prec); mpfr_set_prec(r->uy, prec); mpfr_set_prec(r->vx, prec); mpfr_set_prec(r->vy, prec); i = 0; } mpfr_set(r->ccx, r->cccx, GMP_RNDN); mpfr_set(r->ccy, r->cccy, GMP_RNDN); } return 0; } // get extern void mandelbrot_external_ray_in_get(struct mandelbrot_external_ray_in *r, mpfr_t x, mpfr_t y) { mpfr_set_prec(x, mpfr_get_prec(r->cx)); mpfr_set(x, r->cx, GMP_RNDN); mpfr_set_prec(y, mpfr_get_prec(r->cy)); mpfr_set(y, r->cy, GMP_RNDN); } void usage(const char *progname) { fprintf(stderr, "usage: %s angle depth\n" "outputs space-separated complex numbers on stdout\n" "example : \n %s 1/3 25\n" "or \n %s 1/3 25 > r.txt\n", progname, progname, progname); } // needs input int main(int argc, char **argv) { // read and check input if (argc < 3) { usage(argv[0]); return 1; } mpq_t angle; mpq_init(angle); mpq_set_str(angle, argv[1], 0); mpq_canonicalize(angle); int depth = atoi(argv[2]); struct mandelbrot_external_ray_in *ray = mandelbrot_external_ray_in_new(angle); mpfr_t cre, cim; mpfr_init2(cre, 53); mpfr_init2(cim, 53); // compute and send output for (int i = 0; i < depth * 4; ++i) { mandelbrot_external_ray_in_step(ray); mandelbrot_external_ray_in_get(ray, cre, cim); // send new c to output mpfr_out_str(0, 10, 0, cre, GMP_RNDN); putchar(' '); mpfr_out_str(0, 10, 0, cim, GMP_RNDN); putchar('\n'); } // clear mpfr_clear(cre); mpfr_clear(cim); mandelbrot_external_ray_in_delete(ray); mpq_clear(angle); return 0; } m-exray-out thumb|right|Collecting bits when tracing outwards The trick when tracing outwards is to prepend bits when crossing dwell bands, depending if the outer cell was entered from its left or right inner cell. A picture may make it clearer: Parameter external ray: traced outward ( from point c ( from exterior of Mandelbrot set ) to the infinity ) using Newton method "there are m_r_exray_out_* functions that can be used to get the external angle by tracing an external ray. but it is O(dwell^2), too slow to be practical beyond a few 1000, and it does not yet have adaptive precision so it cam get stuck in tight gaps between spirals..." checkin what ray crosses given c point land on that c point arguments: precision: input string is converted to integer number, in bits "replace double with 1000 to use 1000 bits instead of native 53". It can be "double" string = native precision. c-re = real part of complex c c-im = imaginary part of complex c int sharpness ( replace 8 with 4 to be less sharp, but don't go too low or it may break ) int maxdwell used only for formatting output (" if you have a desired (pre)period, use that to get the angle formatted nicely" ) preperiod ( int) period (int) usage : m-exray-out precision c-re c-im sharpness maxdwell preperiod period Example m-exray-out double -0.1 0.651 8 10000 100000 1 | tail -n 1 output ( only binary expansion of externala angle): .00100100100100100100100100100100100100100100100100100100100100100100100100() or full output: m-exray-out double -1.7904997268969142e-01 1.0856197533070304e+00 4 20 100 1 -1.7904997268969142e-01 1.0856197533070304e+00 -1.7935443511056054e-01 1.0854961235734493e+00 -1.7969015435633140e-01 1.0853555826931407e+00 -1.8051917357961708e-01 1.0850119514784207e+00 1 -1.8269681168554189e-01 1.0840986080730519e+00 -1.8930978756749961e-01 1.0809967205891908e+00 0 -1.9360134634853554e-01 1.0737664210007076e+00 1 -2.0801776265667046e-01 1.0704315884719302e+00 0 -2.4370145596832851e-01 1.0627226337492655e+00 1 -2.7278498707172455e-01 1.0400463032627758e+00 1 -3.6142246174842857e-01 9.5627327660306660e-01 -4.0875942404913329e-01 9.4318016425218687e-01 -4.5809426255241681e-01 9.9407738931502432e-01 1 -6.0223759082842288e-01 9.3674887784190841e-01 -8.0664626534435835e-01 8.8455840042031386e-01 0 -8.3743553612789834e-01 9.3674304615806725e-01 -8.7287981833135664e-01 9.9952074265609137e-01 -9.1381642842183941e-01 1.0750773818725967e+00 -9.6138941844723436e-01 1.1665050307228815e+00 0 -1.0172147464760770e+00 1.2782611560409487e+00 -1.0836516873794186e+00 1.4168346080502183e+00 -1.1642543320399381e+00 1.5917680433570660e+00 -1.2645279518271537e+00 1.8173272090846146e+00 0 -1.3932089082214374e+00 2.1153654197070240e+00 -1.5644798765083927e+00 2.5204286897782238e+00 -1.8019543103478257e+00 3.0891768094384502e+00 -2.1462162170666716e+00 3.9184694463323830e+00 1 -2.6699272827383060e+00 5.1817748128576184e+00 -3.5099957053533952e+00 7.2066815445348764e+00 -4.9407140183321250e+00 1.0650884836446950e+01 -7.5526171722332034e+00 1.6932133257514206e+01 0 -1.2727799761990893e+01 2.9370322256332429e+01 -2.4030360674257757e+01 5.6528056606204473e+01 -5.1753810887543764e+01 1.2313571132413860e+02 -1.2986129881788275e+02 3.1079020626913712e+02 1 -3.8951035577778663e+02 9.3459826885312373e+02 -1.4412515103402823e+03 3.4614095793299703e+03 -6.8365543576278487e+03 1.6423640211807055e+04 -4.3547867632624519e+04 1.0462387704049867e+05 0 -3.9380510085478675e+05 9.4611788566395245e+05 -5.4017847803638447e+06 1.2977803446720989e+07 -1.2161023621086851e+08 2.9216894171553040e+08 .01010001110101() Input : ./ray_in ".(0101101011010110101101011010110101100)" 1000 Example: m-exray-out 100 -0.7432918908524301 0.1312405523087976 8 1000 24 4 result : .010101010101010101010100(1010) which is .01010101010101010101010(01) external ray tracing implementation, that switches between: * perturbation methods (including acceleration with bilinear approximation) * arbitrary precision iterations depending on how close neighbouring points on the ray are. Trying to be correct in all cases but faster than just using arbitrary precision for the whole thing. m-feigenbaum thumb|right|Part of the Mandelbrot set: 12 hyperbolic components period doubling cascade from period 2^0 = 1 to period 2^11 = 2 048, under exponnetial mapping "Estimate location of period 2P atom using size of period P island, use Newton's method to find its nucleus and size more precisely, repeat" computes c centers of hyperbolic components for period doubling cascade on the real axis ( so imaginary part is 0): center for period 1 : c = 0 ( omitted) center for period 2 : c = -1 center for period 2^2= 4 : c = -1.310702641336833008 ... center for period 2^n ./m-feigenbaum -1.000000000000000000e+00 -1.310702641336833008e+00 -1.381547484432061657e+00 -1.396945359704560685e+00 -1.400253081214782869e+00 -1.400961962944840655e+00 -1.401113804939776664e+00 -1.401146325826946537e+00 -1.401153290849922906e+00 -1.401154782546613964e+00 -1.401155102022467736e+00 -1.401155170444417619e+00 -1.401155185098302614e+00 -1.401155188236698823e+00 -1.401155188908866478e+00 -1.401155189052811778e+00 -1.401155189083657104e+00 -1.401155189090263598e+00 -1.401155189091694675e+00 -1.401155189093319819e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 -1.401155189093314712e+00 Small modification: printf("period = %d \t\t\t nucleus = %.18e \tsize = %.18e\n", period, creal(nucleus), size); m-feigenbaum period = 4 nucleus = -1.310702641336833008e+00 size = 1.179602276312223114e-01 period = 8 nucleus = -1.381547484432061657e+00 size = 2.590331788620404280e-02 period = 16 nucleus = -1.396945359704560685e+00 size = 5.574904503817590395e-03 period = 32 nucleus = -1.400253081214783091e+00 size = 1.195288587163239984e-03 period = 64 nucleus = -1.400961962944842210e+00 size = 2.560528805309172169e-04 period = 128 nucleus = -1.401113804939776664e+00 size = 5.484142131458955168e-05 period = 256 nucleus = -1.401146325826946537e+00 size = 1.174547734653405196e-05 period = 512 nucleus = -1.401153290849925570e+00 size = 2.515527294624666474e-06 period = 1024 nucleus = -1.401154782546613964e+00 size = 5.387491809938387059e-07 period = 2048 nucleus = -1.401155102022462628e+00 size = 1.153835913118543914e-07 period = 4096 nucleus = -1.401155170444413400e+00 size = 2.471163334462234791e-08 period = 8192 nucleus = -1.401155185098302614e+00 size = 5.292465190668938379e-09 period = 16384 nucleus = -1.401155188236713700e+00 size = 1.133481216591179550e-09 period = 32768 nucleus = -1.401155188908850047e+00 size = 2.427895317585669274e-10 period = 65536 nucleus = -1.401155189052770256e+00 size = 5.204848444681121846e-11 period = 131072 nucleus = -1.401155189083645558e+00 size = 1.116488091172383448e-11 period = 262144 nucleus = -1.401155189090176112e+00 size = 2.401481531760724511e-12 period = 524288 nucleus = -1.401155189091663589e+00 size = 6.463020532652650290e-13 period = 1048576 nucleus = -1.401155189092033737e+00 size = 3.635753686080000682e-14 period = 2097152 nucleus = -1.401155189093066022e+00 size = 9.366255511607845319e-19 period = 4194304 nucleus = -1.401155189093066022e+00 size = 2.943937947459073103e-26 period = 8388608 nucleus = -1.401155189093066022e+00 size = 5.438985544182134629e-40 period = 16777216 nucleus = -1.401155189093066022e+00 size = 5.967915741322236566e-68 period = 33554432 nucleus = -1.401155189093066022e+00 size = 8.736600707368986815e-124 period = 67108864 nucleus = -1.401155189093066022e+00 size = 1.222084981000085251e-235 period = 134217728 nucleus = -1.401155189093066022e+00 size = 0.000000000000000000e+00 period = 268435456 nucleus = -1.401155189093066022e+00 size = 0.000000000000000000e+00 It seems that double precision is not enough m-feigenbaum3 "estimate location of tip of antenna (renormalized from -2+0i) of period P island using size estimate, from there find nucleus of period 3P island using Newton's method, repeat" fractalforums.org : re-feigenbaum-stretch usage: m-feigenbaum3 re(location) im(location) periodFactor m-feigenbaum3 -2 0 3 #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.659992713573920042e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.492879663528381684e+01,0.000000000000000000e+00 T = 43 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.526412200825988208e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.524573293871185342e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.524447443629173904e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 5.510191709584042741e+01,0.000000000000000000e+00 T = 44 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 4.825684575389948350e+01,0.000000000000000000e+00 T = 42 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = 6.901913875598086001e+00,0.000000000000000000e+00 T = 21 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = inf,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset #preset C = -1.786440255563804813e+00,0.000000000000000000e+00 Z = -nan,-nan T = -2147483648 #endpreset m-feigenbaum3 -2 0 3 period = 3 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.659992713573920042e+01 +0.000000000000000000e+00 T = 44 size = 1.000000000000000000e+00 period = 9 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.492879663528381684e+01 +0.000000000000000000e+00 T = 43 size = 1.903551591313245098e-02 period = 27 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.526412200825988208e+01 +0.000000000000000000e+00 T = 44 size = 3.464641937910663597e-04 period = 81 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.524573293871185342e+01 +0.000000000000000000e+00 T = 44 size = 6.269877899222824167e-06 period = 243 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.524447443629173904e+01 +0.000000000000000000e+00 T = 44 size = 1.134900148763894982e-07 period = 729 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 5.510191709584042741e+01 +0.000000000000000000e+00 T = 44 size = 2.054226124560898465e-09 period = 2187 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 4.825684575389948350e+01 +0.000000000000000000e+00 T = 42 size = 3.718249148787694584e-11 period = 6561 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = 6.901913875598086001e+00 +0.000000000000000000e+00 T = 21 size = 6.727899706332398247e-13 period = 19683 C = -1.786440255563804813e+00 +0.000000000000000000e+00 Z = inf -nan T = -2147483648 size = 1.225497697338603077e-14 m-nearest-roots prints: period + 1 distance m-nearest-roots 5 -6.0077592739339081e-01 9 -1.2345690574077126e+00 17 -1.5948261888043576e+00 33 -1.7912740981919455e+00 65 -1.8940314652359109e+00 129 -1.9466050018630559e+00 257 -1.9731986196121953e+00 513 -1.9865731876029906e+00 1025 -1.9932800441369114e+00 2049 -1.9966383822549216e+00 4097 -1.9983187808738916e+00 8193 -1.9991592878359983e+00 nucleus complex double *c_out; const complex double c_guess = 0.0; m_d_nucleus(c_out, c_guess, period, maxsteps); printf("period = %d \t c = %.16f using double\n", period, creal(*c_out)); // Example m-nucleus 53 -0.8691524744 0.2556487868 48 64 -8.6915874972342078e-01 2.5565708568620021e-01 Using mandelbrot-numerics code to find a period 20000 minibrot near c=i: time ( m-size 50000 `m-nucleus 50000 0 1 20000 32 | tee /dev/stderr` 20000 ) result: Real: -7.9437163502071985489173816041278335565557165593422097986415443718666526638183975269482741283979864796150039328820600171824357494796904100157145276457381327038677304389065791544759044068226298284798485074795465139753011290574940547565632295268309965890736856044225816065002060784840401957436001592831073232116289224835376766024837019729052593614702740064116283928435871912023279827933594666780735810543887516779374159354732526673506932730255777860864410617962251986612293528344539369634355763385689632183660361344255846109415319987535013048216370406305683033823649556367343903980453528078550807179972805086726177286676432146600418321802012705947183289186915513177714879362378374485500397855616627594598759238285491434600163661656097093569019464247435079968347293144424227930353208249468172472673021579306852412165125553964864288143879390896150333653457265029413832025627226640151187411034335007580482643346655743497049637725677688986423224436940553132460564725885956480187539465814934618297444604422558707577928894218074677873014006311238296593575712979281521376321522028222864880290047385477309492930859941138933523459312865570596181676682065664641625195711959850159891958399655264891695837255722367048386224829361833770848853678382189241483407917363077785688704033915945145445141382032098503113386121141821618976435909826878012810311315238838927942886145795725258557067351714858363660176807474230411907332988496361646952820878046993327965167592505974342069923145251125405511298210735382002351444475111195883871537375357019225324966837329585172700415064929059719568669415835536305907605195015067551532717806048959958720723351540772045136711267008091434857778994767951971042965880533276128160757945671341024842631669739934842190421628745823057485321876759312698540929225608881625608367714997633793882137116252066659684867523320525998235559519150039569822804344024248346045328260662350506766959441224921002078134780987834208399572026751630341064076128585432804822899645121908649438079186104626223225086835026727973207053106933545605807581086173838972807582278495224764921308243782567040001796469511628492096380473864043528112882638850008499650397617471312161946511543108954047732214732604192518761590429256338997424173451213591713726558272002992188667406903054868988361837689051267067909007940396165686808795106868838064855005108089621233950746964278729824349462021919452828912238487948591999838204600089889473079741079883341554696124986558349029858769442566619743815502763844097701527868328322323759228519682883781772486449895023296886890395156702488270397371139695710695990682206272068513184347868985875355437735986097344366110220146098454669911676872422507466726262772844396687992023980393180412252620962951600694487648989414046125141359107741643310781724540269445598327469901679273898223648390876129705793738640048570537095417043836345968852339629450002567700189777577964846790537104618393106879939835732394904152494946398727231384913787249745823001738289311645727103810197887877515539458354319263104285589267937979420741297746874900374889388549539632838441649971129535778931268439396997823639463581253024656692066840729143014425564553465530544851212792187013660966565334102233995602996284005853152619476167135763137548181086648741615317739115246457134590756369220988612182423106888272870894941766951250687466714807759040424513113776919266305284693482786000012318355370203000510649927436937799090573931965327514602174194247869031450849888578824607292360533650605659626633451928857198284911861895518252866768510849271811271674466215339007494425799589266330368037494808523322072015355702363872561150182139487288949855852914835191151647214953580577992765768923589452140092207774849659839859152689402490049014646780450855312634212041907597487788211223125953582387474055558632251459696857778410866837388535466178076110887271500918602073599110008782828349828601931359617357488907006692730574882684291152057556083434843531655449794230421218496348063440957859520929997115553713223272459893888222599714500315885758226251226383907100867373946949159210357809409066172026448365305726879690385903084965357293868883959241495118228357461132946108953197235701605108487932988139914383636638205336622907254521211337215354959781371499838940575460882923252535688983896279577141026183725803250025941366898551995317381343295602432854880680913105417171403478847572564220151501945641257407256340713758082071177905993565684176799173836494504063303062708493883917533694417688017619038230941678760662664057011965910800231539615509203568524423393735430913683843408798930423887383733707108315296851537552773795174571699872575675465307644486366568477860016882096508306899350495745941389774331465508906343167441651902671583925484007761590126001875408698776037548173000537310514627103673719106723875477565011787438594421891881293376341836817387959998884620172811186225741103273771975399857174630178915209374630124944654809434747663998641540089249451044747343793912738684530410065341666811576503875144080041376635588445419247720870201032961075580784574643803834224040349515934580424790949301603070840039881116983632470720266625847077403948483754828340322244591513248736220612802960652020455915703123428617077016498405385864432461997065871572070257516341957020978947252938690548377219964733666078050655650666746125819484651038400639111767173754889123755628492910863471339403636445116498685130083472801719762230149643141939154039001471145670460001527706331210349655036747526294372592292737874421182976105207938866601322001196843303227985801688904265826419497214272904851901075941514972117647200278204061021869787050146340891622184720663076877881386068868416366064471354402886662288552667078706482937945747653944295880393480095834773836964538912825459369034140530990210203217740376593495298910340327466503746890581057555405215105328622050072523451061799089554563726997933998099494022843643356810171935875956083213474450964753871742597596389159807474860548937449082552764665987353842664468415487133534606684963715041757839071452744227881047305656416036812546802774367092104370634445506994254300152537889204524138962819925474418837039981003771901730655105777080171387232652844063902421538914319305968703705473681853824215543410439200761681559696856513162963772239637111929231912298561552761476051768494603150247970255274087139475417437666525318837699314477624013690658187766683186987513666300409400279403768729919447855449000797444258300988278145284834950271127363489449639121339411432241394689045849625215128792858854683514551608992163329455011784996662465048627222172830045774045434488939393408845069630045746752624798539349705914385735569651595250218651352938215626071449610093033419595771219480026325620655464072741468432773768618336398345601286278453384744487239861535583484271197660008924457670571819830930013272839428183593854035761689365351460247606170808699221660877760009044931352005199906138059024299151054829073698920465120899137076026463287302396273561452735416778600019151509061871192996649944544635251255611180939657319086749495750406650703523104631376966634580854744528028348644766661861687803874476503066257968690726504793327378255680115059130846355163503033322685363135207492701334721066801514843289941471860778761567160868619351770348032409538305377970923902100885197699029395920849667846810993094612702293706672531075824284278002289430678551378565535722768390905768325693280507976229296441994336567584568603213159337614179963913967140245071426406010126568343361152854104386871950971193764271019726152905089469444204032522564569281201767167121304137166303415015387588952606612266240324729892259808131866274224382895297824105879377663513842630692659814953190082184611299068317174185652027168995389944464563646470703089248951504008599661061822888980138794720478053691503842822034170421882089481794863846177177647112938213372331349756322623026418328473098549683499447951740645257409574863864331936413558808652101756422536878243221845141405047163092445446309652337506812857818627712672119535182802674487218494516703692576818047705180818990227687611349865697646700417084858530111203684655604862529459824809506713267533309643092879524182664685875269816703424729371134632023211821682751089970776093301423791969260953108431784555846379656677296878357053153268291167052511154697766977051318490201590272694053975247144309260438403419516059314927198140682502647753867140477723365586992936626937070645112420562978420494969180857648414886732128746833394793338951136536812507919552105463278901514836838501627409427945166072127389034634662362535768301094003301978784609662103275295887490983814104476638090203027305898662646192580113323630177677789960362015187337059766586565991561230871896736278580930538113756251722826451445081763396147220450655330971575966542768165629827950058414694457612107002519813422137935933949051802958001466676705218971585089638867365108649794271967367583738948555076277154455873640422096572950582246453668223767071601874641607213366687118498245978253390109444418995997750392780107313255139840820773405754956577553794289770409116409571420228693380230025899437603870877054542784281317586453244622402660818643831821462226891660685072457573117889485902077797412004329040165095691873500314917343659203512085778821361204452523870368136469226149736015867804245577929236392814077339316439392508508206520959727928573800885094153800553314652906939293080711662049154470341556580077893261447671240157808246582347200333618553406085052554826594036992297841269796755362033887196089006252179737733109866573340227539862034141198839388081519872826063624800513941456162727625429496413744959884317434642624473125087427423717708382301927586834616549672826513360395335642137272724139545454739463104449628322251582845457933561149241779746747083923968435596911958159633854978255431018163676494053663975309474153510434284245769525836241401382297602993004553331413623115085862521547918233116144303815147644592323844409534033312591272415940304330316480139161832153930361866639526736202903837917442316305999918735106124476264198373838671169480263093073548771272893575807108321917398785294490102026462289081419098788347663282918233719927616971290489000448605479370675888150101653436770635409683915707260833378384329978382632001535926336764826666274559700141470555541357684564795628678577384149810990074344641965294243434863683046415340191004534007747861984920997073615420680497487884717648552467181792474987367112740072255272190588800538041606909185987209170159028084428459572622258803142739029327530812085898685496500033506607514853966450810131101320181576443944528119520678262542516665667772589225658367409437310308191692681701747929941825717813914760436490027296973273696885927024209476169849577627774609821735288543579823873933746165418041525243007331797265820048090493269547699807443692298819880238271902411894022717673853039050368938062710759904516779434177910027164576717375460649435897796076715236692576082399547459593785193303729528593113309230666964408518409037311263890279582964898612945710648622089308353146591558950556457774096875906864050442683109026765888050991130358931886565546257805826690073777858095011429817715273547968084945928031301599621650496157700528460843729211007404968844131661389160468042607031052067033392469344850696050409683463254404597341777664046659525995421309405891229261642956098092899701872432918414448769927311569703420232182442026732210734846569104330914950450386817272059918946706401776070037996180381660685319130717417939927179286304468696802295571418406036768424662411467135246405133085173144139264971120263053400301163186239059701230700612814579728270543021200439105129200019275455983470874139432086875339191481441412783871068993677566173792832145720122415779284407893352668293527601354952001547147784537072825711418553102839482249087891314635317734225217087825443501818465384174890328302318165740210720929441444953628896415212119842220277809228686748800353964121521758377906665150446944597690876685717728001020923030903467397889577625987702929819042793216861981297432410751862996742939821733552376022465063980433476877095489742083667438482524358479920093362671073499261536442509635978938254572171219604223534931048801135996840036883391975783540196339746427425612440311802930181139735185695582618989988117958811316075569114339755674685677959928036900764099458526085419821675058081439927368204872326934818289759660177557064039914505754639463377998991196112023961998610668733702044390628198372939453058274311754372582831610070101869681001548640766315884308521449756188946074250519766336654408201104386777792941840610188681466688436569953569101322404299961245271858446078822593644050595237420612302915602324865912944449105260539387579372432085618031718618467399412498198037020264352404192035385525718637738445737173126887883204422338914094077117405266369732555437529835445785310643028012476150339747654769477837975324497094712742252873395275094111602092150983608511685837417691370283134112558707615144700568498489604300808338456202632290195011163561765905032989492777753130863693458023240541703789339666010447835045257252715361960269883076313236873064974259132787196196424346621374024293770156313978590441642103302299007768492108804953992601544447761538374942979300239461633293490859599063452938858711250267104956507397978474403844807309379849329624027364559546136877836071647522487335680653875008575928199095693277072357337411294063778987393546727728838736254506361668708240985190446787175682032349692908267897330716513026626394037531214887715636133426679029411503800223468119703782419347164191227385529597982634049525824249489381004584297761575054940323722350162865941897918756430316221819391580856455929590012212693971767718111156134625119611537902990743877562663928969788178546328357672482203921172030547393621262720602994243945777523017683019828452189053409346073165497827641541546661811362000069860698473847235166780178734237150200969007414767921753122493459975351206764098974709508735766435838682683431012232653877433197670329155560765634097241957053112237676513811748995099145494510366979493271174202393093949872617145305822631758071249330647970105326491208165384679106654911636291369219591774072389480094363848652230439114251230301727437540976950510736579784395802523606453713668228567047296997827776205239775215773191504628314406040463975203778623513875837323952073230041221445555570861593964350888214523314641085730437948149961122240579538109762211134983888224068986964032346290276113268596102816682853305085149030325552748215495078408622679995786080249473690737026626873405306052299932424760848269117260311258313117827223249824645230833114931474245579229315184388294163598032419249401202852156302446910810247963489869968887090519584584819820604231857986183254977904943759345535508926487273426498495065483968097609125809575651864291878638656751118948252709080366966631434561807538043142623374192000538765039906244673936922990047910893164884788344826422647704033143076607040449945061462449827844843237262351978322868811378865263527055247746845940099739111104279244801339852071272925084450912408790993841951019053768851872131163025230275e-7527 Imag: 9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999991734492077330513689380886733519291674946671408282583644209407630575650409140342995770119063424352210271488566807058196756290031171701011701529100242419560679125380260046815348728362732098903518927842731522442015459899781490157079648835177627682202584659756122982363439113383649152007085924164482218692663011713033118422900677796130976701520211513973179804073856412336912323471275216141504900024638942128288839288824544726257598363037946954747139286140996713614519984757554081285137303793792145483088100948602288971711699818733272109930278128457054114530565415061361089457611033887291209920519084835617420429501562494329154552253147974514089623762774850187343594170152609219811762528104878460467096558776750623330946583444771677244364552796775144576896199277889280199535802388272476453809146195137520963756883307724481370820965867991787274442477494687690041794672045778026327225231447931749253744034412377634436209688157234219689149647284294317994921943258044571538072204447829486195849511596206224271784639410492783745179840157550028751542998761457025592735390256235705578527400511578961509524795765723846769216602534827598544122595055537226394226769431805638838745187514550034115461736580630942772164850138959148037936542189861498291521024384584123266494676747874024658839560362669842111747120028879375611008493691777571650981915670778877856645267853315285564690838429624046285747718787113149379509342042014222183841341852480309391774998684638396165511282456072444392493337309576584765407098415136785336850803964219190455543773565500777469684197751437542701373860103750789485898916431227179802485712595100542939658261974899323094584451988998658155488568870156711079534034726043188161080239900290517603670651552229213022394434669013398782592380368779183995497512609086992015378492374681966407092817332227547134424773876378525948648281148884775757472839432230631294074294857393651417886848443276428721134947690017433025325153248913554312348436870298461202014417691273225926050595655721157069521076214438687068333072136368413085150575362642380849505084969356878401879793263084179759550747507234523232524796967394646964932368873283962438332014394396863216097188067789543751355227029409828459507781736661768248386270817788738777462668450623755227855820114090319893003480371488432847626428399909880285573597465643610948439912464793498238988585822403705164724699464120093570367708183907545235654250619728658435398562469329842236925494721158612870658140940335567711385774552936065430987224064165463886334630367918901941265514184193460454762359263731196410400313787736101813913548190929105563969609759961286142020359681386163860459140037407728122634484712324766777888892902726372897663970591495119930657196797922969697210887944061663587643330053324801023593783956257456402354277193486606699249716660922975324918446780356978648595624510216079738304712240699383022296050307095326386181436021280821946461416911955375685450958489937715061871379614514046821562054907141306926601353784074703466872196180068651947552564031583073723325119005346348879674247302910116126291796234873627140147503137007974533422374524109527394238215360495340869525049828516165693842590343161532628841468410993746576931920685810272338948639553304457274026775879365952258661615059382775614181775292223518160756734117381397998310327294281664957921427253102145159956251792002855401711344356095059093958391981346288299586778567037745773029382822174296514329430014003230736855824042705301686236442231943132404408336068221178832509202932258219292166235157083978711447788498538296242489741089788080914827695688672104190838035163401216262600615305503117497776805939741146793074105899920727912025006743184132116016580499916303874955586731672402804758464494865562667287513471323587649568052936006244205760085680463430770593090850578873410241227234316091614508294678455447144047019836067353592195113713166048538989096099123016992319315487579497656388673505433776336703189690278955717239710851820593750575802469723477164557613094688400569544664083906632328980489050231322562719313784223293757762537765674860550223296757517936309514311364098955067987810622990290044610233314787617884818714076827914250567437281973993353197482098813435409755235361075969985798441337126312593960627842039335713762009624169971072389098442582559772999481646704068278833015939333008141841109592846547157670464817674715015636259928604821071600800873091706767125393539557966484701666419162259684521882415687361386486795039924743902847684455653672009117814002731651676792734130082374726001233740571183340055316038466019397098875483608598573963531197795231533784669488383851034898599484418365666793002675472937237274092221819708169896011572109417029618685719634982095269652618039495655688706701247477025085922221279463320915913152801169322907305498072372506045260384015734591714087786003188152532841620806596086485229059207306028656858812580376023757671759911398575408116826175902909012761714255020785822061641547239324759162257787133638766091640926969579726988971447262832031549969840158881281750736965695337486321518286113798675255459838184121812112863705855160188444420663952500796910005915012793169847550253471470237546273447099368433714408055787441435236639758429439479732597655485321756562005631462458721857009810050979024000068037899966702594286779178167104283058516755533365765054701540652974134109847005430449255509769055026949291803646099149879515878433840836779489981856269066079459746553702885573001041905761796108832155227593049875211636581275432414690028646002788224021719021326741635129313266431353578862049903277695541296438178355751015037527905950478175463466493661217977735637802129256184758010305201810374991157129887897715190473625904450982408862004860880041855444963715403816816697913347488258213796231921246982456767604532357185085595619225419494972418775913218700576674855723280963208713608759050096965340865736024194082384061468544443620372024356238984943165014000830348667282990315062294958017713685112256010523333894116609975573376955381051903582439853973475874565611435725163035826136598817082772994463193067795074551604854663708173556403995799133023516315205958185369423580248522691815920365525381967586987223820592252750246790439201396825488793939257364115218262946083474860523784830924246085630585814123698134711292541643275834774769358550038707208687561681985220208081109264672682079000578528676995067870581728389481309632471022423883273871806201922953495909101020921632302048824742636087442222271981512160368757719560642803804335208325783581324630948528372742921249668136599699430864754216512734431140871256723937412438287901118364475021875556572692586266456274813336044845327973318471952277585776779170860292001768956960996018482974850927059488179484838649054397646499203792764009751267114995747305588264268827308358967751669508793420818356792715912995146476118789747260518386453554450611500509015784210629209568851076896828558589881187034131236514543108275118274492680820722023356279609283225501210503057769959987375674460963328162312089681988114274578205642753188886259905616611476198660210570725766467627170954257896651227694959336299713781484347595473261014179723849074993539787736555059799178437730200032338748364417885122128183167598798597900535307364189904009133503543999817322306032952189385928613024146079985906321984855056816088225874686809752502871137380525052265954021704331783968465977328015884627311680559023242904016721675131323566011657154289720702412494820964921312304644491537963385820077869566295857565120487324475146126593516710180032457589644029021496e-01 Size: 2.395769e-15051 /* console program in c language based on the code by Claude Heiland-Allen http://code.mathr.co.uk/mandelbrot-numerics it computes center ( nucleus) of hyperbolic componnet of Mandelbrot set using Newton method and double precision https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers to compile gcc m.c -Wall -std=c99 -lm to run ./a.out ------------- Newton method converged after 6 steps and near c= ( 0.000000 ;1.000000 ) found nucleus = center = ( -0.1225611668766536, 0.7448617666197442 ) for period = 3 Newton method converged after 5 steps and near c= ( 0.250000 ;0.500000 ) found nucleus = center = ( 0.2822713907669139, 0.5300606175785253 ) for period = 4 Newton method converged after 5 steps and near c= ( 0.250000 ;-0.500000 ) found nucleus = center = ( 0.2822713907669139, -0.5300606175785253 ) for period = 4 Newton method converged after 8 steps and near c= ( 0.300000 ;0.300000 ) found nucleus = center = ( 0.3795135880159238, 0.3349323055974976 ) for period = 5 Newton method converged after 5 steps and near c= ( -0.122561 ;0.844862 ) found nucleus = center = ( -0.1134186559494366, 0.8605694725015731 ) for period = 6 Newton method converged after 7 steps and near c= ( -1.250000 ;0.250000 ) found nucleus = center = ( -1.1380006666509646, 0.2403324012620983 ) for period = 6 Newton method converged after 35 steps and near c= ( -3.250000 ;0.250000 ) found nucleus = center = ( -1.9667732163929286, -0.0000000000000000 ) for period = 6 compare https://gitlab.com/adammajewski/cpp-mandelbrot-center-by-knighty https://gitlab.com/adammajewski/lawrence https://gitlab.com/adammajewski/mandelbrot-numerics-nucleus */ #include <complex.h> #include <math.h> #include <stdbool.h> #include <stdint.h> #include <stdlib.h> #include <stdio.h> // // from mandelbrot-numerics.h enum m_newton { m_failed, m_stepped, m_converged }; typedef enum m_newton m_newton; // from m_d_util.h static inline bool cisfinite(double _Complex z) { return isfinite(creal(z)) && isfinite(cimag(z)); } static inline double cabs2(double _Complex z) { return creal(z) * creal(z) + cimag(z) * cimag(z); } // last . takeWhile (\x -> 2 /= 2 + x) . iterate (/2) $ 1 :: Double static const double epsilon = 4.440892098500626e-16; // epsilon^2 static const double epsilon2 = 1.9721522630525295e-31; int PrintResult(int iResult, double _Complex c_guess, double _Complex nucleus, int period, int steps) { static char* sResult; switch (iResult){ case 0 : sResult = "failed"; break; case 1 : sResult = "stepped"; break; case 2 : sResult = "converged"; break; default : sResult = "unknown"; } printf("Newton method %s after %d steps and near c= ( %f ;%f ) found nucleus = center = ( %.16f, %.16f ) for period = %d \n ", sResult, steps, creal(c_guess), cimag(c_guess), creal(nucleus), cimag(nucleus), period ); return 0; } // ----------------------- from m_d_nucleus.c ------------------------------------------------- // nucleus = center of hyperbolic componnet of Mandelbrot set // https://en.wikibooks.org/wiki/Fractals/Mathematics/Newton_method#center // extern m_newton m_d_nucleus_step(double _Complex *c_out, double _Complex c_guess, int period) { double _Complex z = 0; double _Complex dc = 0; for (int i = 0; i < period; ++i) { dc = 2 * z * dc + 1; z = z * z + c_guess; } if (cabs2(dc) <= epsilon2) { *c_out = c_guess; return m_converged; } double _Complex c_new = c_guess - z / dc; double _Complex d = c_new - c_guess; if (cabs2(d) <= epsilon2) { *c_out = c_new; return m_converged; } if (cisfinite(d)) { *c_out = c_new; return m_stepped; } else { *c_out = c_guess; return m_failed; } } // compute nucleus using double precision and Newton method extern m_newton m_d_nucleus(double _Complex *c_out, double _Complex c_guess, int period, int maxsteps) { m_newton result = m_failed; double _Complex c = c_guess; int i; for (i = 0; i < maxsteps; ++i) { if (m_stepped != (result = m_d_nucleus_step(&c, c, period))) { break; } } *c_out = c; PrintResult(result,c_guess, c , period, i); return result; } /* ==================== main ================================================================================================*/ int main() { double _Complex c3, c4a, c4b, c5, c3c2, c2c3; m_d_nucleus(&c3, 0 + I * 1, 3, 64); m_d_nucleus(&c4a, 0.25 + 0.5 * I, 4, 64); m_d_nucleus(&c4b, 0.25 - 0.5 * I, 4, 64); m_d_nucleus(&c5, 0.3 + 0.3 * I, 5, 64); m_d_nucleus(&c3c2, c3 + I * 0.1, 6, 64); m_d_nucleus(&c2c3, -1 - 0.25 + 0.25 * I, 6, 64); m_d_nucleus(&c2c3, -3 - 0.25 + 0.25 * I, 6, 64); return 0; } m-render-wakes Draw wakes filled with solid color, uses mpq type from gmp library m-render-wakes > m-render-wakes.ppm Images Problems error while loading shared libraries cd ~/mandelbrot-numerics/c/bin ./m-interior ./m-interior: error while loading shared libraries: libmandelbrot-numerics.so: cannot open shared object file: No such file or directory export LD_LIBRARY_PATH=${HOME}/opt/lib export PATH=${HOME}/opt/bin:${PATH} ./m-interior usage: ./m-interior precision z-guess-re z-guess-im c-guess-re c-guess-im interior-r interior-t period maxsteps References