GMP Mandelbrot

Some code to calculate the Mandelbrot set at arbitrary precision using GMP. It does orbit detection. I release this code under the GNU Public License.

#include <gmp.h>
#include <iostream>
#include <cassert>

using std::cout;
using std::endl;

typedef unsigned int Iter;
#define ORBIT
 

long long total_iterations=0;
long long saved_iterations=0;
void inner_loop(int bits, const mpz_t a, const mpz_t b, const mpz_t threshold, Iter max_iter, 
                Iter *iterations, bool *in_m_set) {
  //NOTE that a,b is expected to be premultiplied by 2^bits, and threshold is 4*4^bits;
  mpz_t x2,y2,add;
  mpz_init(x2);
  mpz_init(y2);
  mpz_init(add);
  mpz_t xfast,yfast;
  mpz_init_set(xfast,a);
  mpz_init_set(yfast,b);

  //fast starts with an advantage already

#ifdef ORBIT
  mpz_t xslow,yslow;
  mpz_init_set_si(xslow,0);
  mpz_init_set_si(yslow,0);
#endif

  *in_m_set=true;
  for(*iterations=0;*iterations<max_iter;++*iterations){
    mpz_pow_ui(x2,xfast,2);
    mpz_pow_ui(y2,yfast,2);
    mpz_add(add,x2,y2);
    if (mpz_cmp(add,threshold)>0) {
      *in_m_set=false;
      break;
    }
    mpz_mul(yfast,yfast,xfast);
    mpz_tdiv_q_2exp(yfast,yfast,bits-1);  //y=2*x*y

    mpz_sub(xfast,x2,y2);
    mpz_tdiv_q_2exp(xfast,xfast,bits);

#ifdef ORBIT
    mpz_pow_ui(x2,xslow,2);
    mpz_pow_ui(y2,yslow,2);
    mpz_add(add,x2,y2);
    mpz_mul(yslow,yslow,xslow);
    mpz_tdiv_q_2exp(yslow,yslow,bits-1);  //y=2*x*y
    mpz_sub(xslow,x2,y2);
    mpz_tdiv_q_2exp(xslow,xslow,bits);

    if (mpz_cmp(xslow,xfast)==0 && mpz_cmp(yslow,yfast)==0) {
      //orbit detection
      saved_iterations+= (max_iter-*iterations);
      break;
      //in_m_set is already true
    }
#endif
    mpz_add(yfast,yfast,b);    
    mpz_add(xfast,xfast,a);

    ++*iterations;
    ++total_iterations;

    mpz_pow_ui(x2,xfast,2);
    mpz_pow_ui(y2,yfast,2);
    mpz_add(add,x2,y2);
    if (mpz_cmp(add,threshold)>0) {
      *in_m_set=false;
      break;
    }
    mpz_mul(yfast,yfast,xfast);
    mpz_tdiv_q_2exp(yfast,yfast,bits-1);  //y=2*x*y

    mpz_sub(xfast,x2,y2);
    mpz_tdiv_q_2exp(xfast,xfast,bits);
#ifdef ORBIT    
    if (mpz_cmp(xslow,xfast)==0 && mpz_cmp(yslow,yfast)==0) {
      //orbit detection
      saved_iterations+= (max_iter-*iterations);
      break;
      //in_m_set is already true
    }

    mpz_add(xslow,xslow,a);
    mpz_add(yslow,yslow,b);

#endif
    mpz_add(yfast,yfast,b);    
    mpz_add(xfast,xfast,a);
    ++total_iterations;

  }
}

const int pixel_bits=10;
/* two experiments saw convergence around 50.  40 is pushing it.  as low as 10 if
   you don't need it to be perfect.  60 or 70 to be safe. */

const int screen_half_width=100;

void calc2(Iter max_iter, mpz_t Cx0, mpz_t Cy0, int Cbits, int bits) {
  cout << "P5" << endl << 2*screen_half_width+1 << " " << 2*screen_half_width+1 << endl << 255 << endl;
  assert(bits>=Cbits);
  mpz_t cx,cy;
  mpz_init(cx);
  mpz_init(cy);
  mpz_mul_2exp(cx,Cx0,bits-Cbits);
  mpz_mul_2exp(cy,Cy0,bits-Cbits);

  mpz_t threshold;
  mpz_init_set_si(threshold,4);
  mpz_mul_2exp(threshold,threshold,2*bits);
  
  mpz_t a,b;
  mpz_init(a);
  mpz_init(b);
  for(int dx=-screen_half_width;dx<=screen_half_width;++dx){
    mpz_set_si(a,dx);
    mpz_mul_2exp(a,a,pixel_bits);
    mpz_add(a,a,cx);
    for(int dy=-screen_half_width;dy<=screen_half_width;++dy){
      mpz_set_si(b,dy);
      mpz_mul_2exp(b,b,pixel_bits);
      mpz_add(b,b,cy);
      Iter i;
      bool in_m_set;
      inner_loop(bits,a,b,threshold,max_iter, &i,&in_m_set);

      if (in_m_set) {
        putchar(0);
        //putchar('M');
      }

      else {
        putchar(1+(i%255));
        //putchar(' ');
      }
    }

    //putchar('\n');

  }
}


int main(int argc, char**argv){
  assert(argc==2);
  int mag=atoi(argv[1]);
  mpz_t cx,cy;
  mpz_init_set_si(cx,-1314468);  // -1
  mpz_init_set_si(cy,-397404);
  calc2(100000,cx,cy,20,mag+pixel_bits);
  std::cerr << total_iterations << endl << saved_iterations << endl;
}