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-dev

#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:

 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 .profile

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

Description

size

atom domain

Atom domain size estimate 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:

• precision
• 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 48

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.stackexchange
• blog
• "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:

${\begin{cases}F^{p}(z,c)=z\\{\dfrac {\partial }{\partial z}}F^{p}(z,c)=b\end{cases}}$

where:

• F is function ( complex quadratic polynomial)
• $F^{p}$ is iterated polynomial
• p is the period of the target component
• $b=re^{2\pi i\theta }$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

description

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:

• 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

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

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"

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

• Wakes
• 
  Wakes near the period 1 continent
• 
  Wakes near the period 3 island
• 
  Wakes along the main antenna

• 
  Escape route 1/2

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