Fractals/mandelbrot-numerics

< Fractals

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)



Installation

dependencies


#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

git clone http://code.mathr.co.uk/mandelbrot-numerics.git


and

 make -C mandelbrot-numerics/c/lib prefix=${HOME}/opt install
 make -C mandelbrot-numerics/c/bin prefix=${HOME}/opt install

hen 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[3]

how to use it

binaries

 m-box-period

The result : 

 usage: m-box-period precision center-re center-im radius maxperiod

code

These are one file program using code extracted from Claude's library.

shape

Versions

Description[4]

size

size estimate formula from 

 http://ibiblio.org/e-notes/MSet/windows.htm <-- where I got the formula from

Files:

box-period

Computing period under complex quadratic polynomial

/*


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_d_interior.c

Description:

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

/*

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;
}

Images

  1. ↑ atom_domain_size_estimation
  2. ↑ interior_coordinates_in_the_mandelbrot_set
  3. ↑ stackoverflow question how-to-permanently-set-path-on-linux
  4. ↑ math.stackexchange question: perfect-circles-in-the-mandelbrot-set
  5. ↑ math.stackexchange question test-for-membership-in-mandelbrot-bulb-of-period-n/1151953#1151953
  6. ↑ practical_interior_distance_rendering
This article is issued from Wikibooks. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.