20090110.mandpoints.c

#include <stdio.h>
#include <stdlib.h>

#include <math.h>
#include <complex.h>

#include <assert.h>

#ifndef MAXITER
#define MAXITER 1000
#endif

#ifndef LEAVERADIUS
#define LEAVERADIUS 1.e5
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef M_PIl
#define M_PIl 3.1415926535897932384626433832795029L
#endif

#if 1
#define udouble long double
#define uexp2(x) exp2l(x)
#define ucexp(z) cexpl(z)
#define uclog(z) clogl(z)
#define ucreal(z) creall(z)
#define ucimag(z) cimagl(z)
#define ufabs(x) fabsl(x)
#define ucabs(z) cabsl(z)
#define PI M_PIl
#define UFFORMAT "Lf"
#define UEFORMAT "Le"
#else
#define udouble double
#define uexp2(x) exp2(x)
#define ucexp(z) cexp(z)
#define uclog(z) clog(z)
#define ucreal(z) creal(z)
#define ucimag(z) cimag(z)
#define ufabs(x) fabs(x)
#define ucabs(z) cabs(z)
#define PI M_PI
#define UFFORMAT "f"
#define UEFORMAT "e"
#endif

#ifndef DEBUG
#define DEBUG 0
#endif

char
examine_point (udouble complex c, long int *rn2,
               udouble complex *ru, udouble complex *rup)
/* Compute Boettcher coordinate (u) and its derivative (up);
   u is defined only mod 2pi/(2^n2). */
{
  long int n = 1;
  long int n2 = 0;
  udouble complex z = c;
  udouble complex u = uclog(z);
  udouble complex up = 1/z;
  while ( n < MAXITER )
    {
      if ( ( ucreal(z)>LEAVERADIUS/2 || ucreal(z)<-LEAVERADIUS/2
             || ucimag(z)>LEAVERADIUS/2 || ucimag(z)<-LEAVERADIUS/2 )
           && ucabs(z)>LEAVERADIUS )
        {
#if DEBUG
          printf ("***\tc = %+.15" UFFORMAT "%+.15" UFFORMAT "*I :\n\tEscaped after %ld iterations\n"
                  "\tu = %+.15" UFFORMAT "%+.15" UFFORMAT "*I+k*%.15" UFFORMAT "*I\n\tu'= %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n",
                  ucreal(c), ucimag(c), n, ucreal(u), ucimag(u), 2*PI*uexp2(-n2),
                  ucreal(up), ucimag(up));
#endif
          if ( rn2 )
            *rn2 = n2;
          if ( ru )
            *ru = u;
          if ( rup )
            *rup = up;
          return 1;
        }
      z = z*z + c;
      u -= uclog(1-c/z)*uexp2(-n);
      up -= (c*up-uexp2(-n))/z;
      if ( ( ucreal(z)<=2 && ucreal(z)>=-2 && ucimag(z)<=2 && ucimag(z)>=-2 )
           && ucabs(z)<=4 )
        n2 = n;
      n++;
    }
  return 0;
}

udouble complex
trace_ray (udouble l)
/* Attempt to locate a point with external argument l (that's 2pi*l)
   which is close enough to the Mandelbrot set. */
{
  udouble ltwopi = 2*PI*l;
  udouble complex z = 4.*ucexp(ltwopi*I);
  long int n2;
  udouble complex u, up;
#if DEBUG
  printf ("Target arg:                   %+.15" UFFORMAT "\n", ltwopi);
  printf ("Startw/ z = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(z), ucimag(z));
#endif
  examine_point (z, &n2, &u, &up);
  if ( fabs(ucimag(u)-ltwopi) > PI )
    ltwopi += round((ucimag(u)-ltwopi)/(2*PI))*(2*PI);
#if DEBUG
  printf ("Doing initial Newton's method...\n");
#endif
  for ( int i=0 ; i<5 ; i++ )
    {
      udouble diff = ltwopi-ucimag(u);
      udouble modulus = 2*PI*uexp2(-n2);
      assert (n2 == 0);
#if DEBUG
      printf ("diff to fix -> %.5" UEFORMAT "\n", diff);
#endif
      diff -= round(diff/modulus)*modulus;
      z += diff*I/up;
      examine_point (z, &n2, &u, &up);
    }
  long int stepcnt = 0;
  do
    {
#if DEBUG
      printf ("MAJOR STEP\n");
      printf ("Current z = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(z), ucimag(z));
      printf ("Target arg:                   %+.15" UFFORMAT "\n", ltwopi);
      printf ("Doing Runge-Kutta for descent...\n");
#endif
      udouble h = ucreal(u)*0.1;
      udouble complex k1, k2, k3, k4;
      examine_point (z, NULL, NULL, &k1);
      examine_point (z-0.5*h/k1, NULL, NULL, &k2);
      examine_point (z-0.5*h/k2, NULL, NULL, &k3);
      examine_point (z-h/k3, NULL, NULL, &k4);
      z -= h*(1./k1+2./k2+2./k3+1./k4)/6;
#if DEBUG
      printf ("RK: now z = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(z), ucimag(z));
#endif
      examine_point (z, &n2, &u, &up);
      if ( n2 < 25 )
        {
#if DEBUG
          printf ("Doing corrective Newton's method...\n");
#endif
          for ( int i=0 ; i<2 ; i++ )
            {
              udouble diff = ltwopi-ucimag(u);
              udouble modulus = 2*PI*uexp2(-n2);
              diff -= round(diff/modulus)*modulus;
#if DEBUG
              printf ("diff to fix -> %.5" UEFORMAT " (modulo %.5" UEFORMAT ")\n", diff, modulus);
#endif
              z += diff*I/up;
              examine_point (z, &n2, &u, &up);
            }
#if DEBUG
          printf ("Nt: now z = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(z), ucimag(z));
#endif
        }
      stepcnt++;
    }
  while ( ucreal(u) > 3.e-10 && stepcnt < 1000 );
#if DEBUG
  printf ("%+.17" UFFORMAT "\t%+.17" UFFORMAT "\n", ucreal(z), ucimag(z));
#endif
  return z;
}

void
doubling_period (long int p, long int q,
                 long int *rper, long int *rpreper)
/* Determine period and preperiod for point with angle p/q (i.e., just
   the doubling period of p/q mod 1). */
{
  assert ( ! ( p%2 == 0 && q%2 == 0 ) );
  long int n = 0;
  while ( q%2 == 0 )
    {
      q /= 2;
      n++;
    }
  if ( rpreper )
    *rpreper = n;
  n = 0;
  long int tp = p%q;
  do
    {
      p = (p*2)%q;
      n++;
    }
  while ( p != tp );
  if ( rper )
    *rper = n;
}

long int *c_sucks;

static int
c_sucks_compare (const void *a, const void *b)
{
  long int na = *(long int *)a;
  long int nb = *(long int *)b;
  return (c_sucks[na]<c_sucks[nb] ? -1 :
          (c_sucks[na]>c_sucks[nb] ? 1 : 0));
}

long int
check_subperiod (long int p, long int q,
                 long int per, long int subper)
/* Attempt to test the parent period (here known as "subperiod") for
   the component with angle p/q (q odd, hence periodic), period per
   having already been computed.  If this succeeds (component is a bud
   and not a cardioid), i.e., if subper is the desired subperiod, then
   it returns the associated rotation angle, expressed in units of
   (per/subper). */
{
#if DEBUG
  printf ("Checking subperiod %ld...\n", subper);
#endif
  assert ( q%2 != 0 );
  assert ( per%subper == 0 );
  long int valence = per/subper;
  long int ptab[per];
  ptab[0] = p;
  for ( int i=1 ; i<per ; i++ )
    ptab[i] = (ptab[i-1]*2)%q;
  assert ( (ptab[per-1]*2)%q == p );
#if DEBUG
  printf ("Doubling structure mod %6ld:", q);
  for ( long int i=0 ; i<per ; i++ )
    printf("%6ld", ptab[i]);
  printf ("\n");
#endif
  c_sucks = ptab;
  long int cind[subper][valence];  // Doubling indexes arranged cyclically in each set
  long int icyc[per]; // Cyclic indexes in the respective set
  for ( long int k=0 ; k<subper ; k++ )
    {
      for ( long int j=0 ; j<valence ; j++ )
        cind[k][j] = j*subper + k;
      qsort (&(cind[k]), valence, sizeof(long int), c_sucks_compare);
#if DEBUG
      printf ("Cyclic order for set %3ld: ", k);
      for ( long int j=0 ; j<valence ; j++ )
        printf("%3ld", cind[k][j]);
      printf ("\n");
#endif
      for ( long int j=0 ; j<valence ; j++ )
        icyc[cind[k][j]] = j;
    }
#if DEBUG
  printf ("Cyclic indexes:               ");
  for ( long int i=0 ; i<per ; i++ )
    printf("%6ld", icyc[i]);
  printf ("\n");
#endif
  // Check whether sets are not intertwined:
  for ( long int k0=0 ; k0<subper ; k0++ )
    for ( long int k1=k0+1 ; k1<subper ; k1++ )
      {
        long int j0;
        for ( j0=0 ; j0<valence-1 ; j0++ )
          if ( ptab[cind[k1][0]] > ptab[cind[k0][j0]]
               && ptab[cind[k1][0]] < ptab[cind[k0][j0+1]] )
            break;
        for ( long int j1=0 ; j1<valence ; j1++ )
          if ( ( j0 < valence-1 )
               ? ( ptab[cind[k1][j1]] < ptab[cind[k0][j0]]
                   || ptab[cind[k1][j1]] > ptab[cind[k0][j0+1]] )
               : ( ptab[cind[k1][j1]] < ptab[cind[k0][j0]]
                   && ptab[cind[k1][j1]] > ptab[cind[k0][0]] ) )
            {
#if DEBUG
              printf ("Sets %ld and %ld are intertwined (index %ld of %ld does not occur between indices %ld and %ld of set %ld as index 0 does)\n",
                      k0, k1, j1, k1, j0, (j0+1)%valence, k0);
#endif
              return 0;
            }
      }
  // Check whether doubling preserves cyclic order between sets:
  for ( long int k=0 ; k<subper ; k++ )
    {
      long int t = (icyc[(cind[k][0]+1)%per]-icyc[cind[k][0]]+valence)%valence;
      for ( long int j=0 ; j<valence ; j++ )
        if ( icyc[(cind[k][j]+1)%per] != (icyc[cind[k][j]]+t)%valence )
          {
#if DEBUG
            printf ("Cyclic order not preserved between sets %ld and %ld\n",
                    k, (k+1)%subper);
#endif
            return 0;
          }
    }
  // Sanity check...
  long int rot = (icyc[subper]-icyc[0]+valence)%valence;
  for ( long int i=0 ; i<per ; i++ )
    assert (icyc[(i+subper)%per] == (icyc[i]+rot)%valence);
  // Determine smallest angle in a set:
  long int crit_angle = q;
  long int crit_count = 1;
  long int crit_k = 0;
  long int crit_j = icyc[0];
  for ( long int k=0 ; k<subper ; k++ )
    for ( long int j=0 ; j<valence ; j++ )
      {
        long int angle = ( (j<valence-1)
                           ?(ptab[cind[k][j+1]]-ptab[cind[k][j]])
                           :(ptab[cind[k][0]]-ptab[cind[k][j]]+q) );
        if ( angle < crit_angle )
          {
            crit_angle = angle;
            crit_count = 1;
            crit_k = k;
            crit_j = j;
          }
        else if ( angle == crit_angle )
          crit_count++;
      }
  if ( crit_count > 1 || crit_k != 0 || ( crit_j != icyc[0]
                                          && (crit_j+1)%valence != icyc[0] ) )
    {
#if DEBUG
      printf ("Doubling structure is coherent, but critical angle in bad place.\n");
#endif
      return 0;
    }
#if DEBUG
  printf ("This (%ld) appears to be the correct subperiod!, with rotation %ld/%ld.\n",
          subper, rot, valence);
#endif
  return rot;
}

void
p_and_pp (udouble complex c, long int n,
          udouble complex *rp, udouble complex *rpp)
/* Compute n-th iterate of c under z->z^2+c, and its derivative wrt c */
{
  long int i = 0;
  udouble complex p = c;
  udouble complex pp = 1;
  for ( i=0 ; i<n ; i++ )
    {
      pp = 2*p*pp + 1;
      p = p*p + c;
    }
  if ( rp )
    *rp = p;
  if ( rpp )
    *rpp = pp;
}

void
p_and_derivs (udouble complex z, udouble complex c, long int n,
              udouble complex *rp, udouble complex *rpz, udouble complex *rpz2,
              udouble complex *rpc, udouble complex *rpzc)
/* Compute n-th iterate of z under z->z^2+c, and its derivatives wrt z
   and c and second derivatives wrt z twice and z-and-c */
{
  long int i = 0;
  udouble complex p = z;
  udouble complex pz = 1;
  udouble complex pz2 = 0;
  udouble complex pc = 0;
  udouble complex pzc = 0;
  for ( i=0 ; i<n ; i++ )
    {
      pzc = 2*(pz*pc+p*pzc);
      pc = 2*p*pc + 1;
      pz2 = 2*(pz*pz+p*pz2);
      pz = 2*p*pz;
      p = p*p + c;
    }
  if ( rp )
    *rp = p;
  if ( rpz )
    *rpz = pz;
  if ( rpz2 )
    *rpz2 = pz2;
  if ( rpc )
    *rpc = pc;
  if ( rpzc )
    *rpzc = pzc;
}

void
process_angle (long int p, long int q)
{
  while ( p%2 == 0 && q%2 == 0 )
    {
      p /= 2;
      q /= 2;
    }
#if DEBUG
  printf ("========= ANGLE %ld/%ld =========\n", p, q);
  printf ("****** FIRST PART: RAY TRACING ******\n");
#endif
  udouble complex c0 = trace_ray ((udouble)p/q);
  long int preper, per;
#if DEBUG
  printf ("****** COMBINATORIC INTERLUDE: DOUBLING PERIOD ******\n");
#endif
  doubling_period (p, q, &per, &preper);
#if DEBUG
  printf ("%ld/%ld has doubling period %ld after preperiod %ld\n",
          p, q, per, preper);
#endif
  long int subper = per, valence = 1, rot = 0;
  if ( preper == 0 )
    {
      char have_subper = 0;
      for ( long int sb=1 ; sb<per ; sb++ )
        if ( per%sb == 0 )
          {
            long int rt = check_subperiod (p, q, per, sb);
            if ( rt )
              {
                assert (! have_subper);
                have_subper = 1;
                subper = sb;
                valence = per/sb;
                rot = rt;
              }
          }
    }
#if DEBUG
  printf ("****** SECOND PART: NEWTON AGAIN ******\n");
#endif
  if ( preper > 0 )
    {
      // It's a Misiurewicz point...
      udouble complex c = c0;
      udouble complex p0, pp0, p1, pp1;
      for ( int i=0 ; i<10 ; i++ )
        {
#if DEBUG
          printf ("Current c = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(c), ucimag(c));
#endif
          p_and_pp (c, preper, &p0, &pp0);
          p_and_pp (c, preper+per, &p1, &pp1);
          udouble complex diff = p1-p0;
#if DEBUG
          printf ("|diff| to fix -> %.5" UEFORMAT "\n", ucabs(diff));
#endif
          c -= diff/(pp1-pp0);
        }
#if DEBUG
      printf ("%+.17" UFFORMAT "\t%+.17" UFFORMAT "\n", ucreal(c), ucimag(c));
#endif
        udouble ci_cheat = ucimag(c);
        if ( ufabs(ci_cheat) < 1.e-7 )
          ci_cheat = 0;
      printf ("%ld/%ld\t%+.17" UFFORMAT "\t%+.17" UFFORMAT "\t"
              "Misiurewicz\t%ld(%ld)\n",
              p, q, ucreal(c), ci_cheat, per, preper);
      fflush (stdout);
    }
  else
    {
      // It's a component...
      {
#if DEBUG
        printf ("Looking for center...\n");
#endif
        udouble complex c = c0;
        udouble complex pt, pp;
        for ( int i=0 ; i<10 ; i++ )
          {
#if DEBUG
            printf ("Current c = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(c), ucimag(c));
#endif
            p_and_pp (c, per-1, &pt, &pp);
            udouble complex diff = pt;
#if DEBUG
            printf ("|diff| to fix -> %.5" UEFORMAT "\n", ucabs(diff));
#endif
            c -= pt/pp;
          }
#if DEBUG
        printf ("Component center:\n");
        printf ("%+.17" UFFORMAT "\t%+.17" UFFORMAT "\n", ucreal(c), ucimag(c));
#endif
        udouble ci_cheat = ucimag(c);
        if ( ufabs(ci_cheat) < 1.e-7 )
          ci_cheat = 0;
        printf ("%ld/%ld\t%+.17" UFFORMAT "\t%+.17" UFFORMAT "\t"
                "CompCenter\t%ld\n",
                p, q, ucreal(c), ci_cheat, per);
        fflush (stdout);
      }
      {
#if DEBUG
        printf ("Looking for bud root...\n");
#endif
        udouble complex c = c0;
        udouble complex z = 0;
        udouble complex pt, pz, pz2, pc, pzc;
        udouble complex zeta = ucexp(2*PI*I*((udouble)rot/valence));
#if DEBUG
        printf ("subper: %ld\n", subper);
        printf ("   rot: %ld/%ld\n", rot, valence);
        printf ("     zeta = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(zeta), ucimag(zeta));
#endif
        for ( int i=0 ; i<10 ; i++ )
          {
#if DEBUG
            printf ("Current c = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(c), ucimag(c));
            printf ("Current z = %+.15" UFFORMAT "%+.15" UFFORMAT "*I\n", ucreal(z), ucimag(z));
#endif
            p_and_derivs (z, c, subper, &pt, &pz, &pz2, &pc, &pzc);
            udouble complex diff0 = pt - z;
            udouble complex diff1 = pz - zeta;
#if DEBUG
            printf ("|diff| to fix -> %.5" UEFORMAT "\n", ucabs(diff0)+ucabs(diff1));
#endif
            udouble complex a00=pz-1, a01=pc, a10=pz2, a11=pzc;
            udouble complex det = a00*a11-a01*a10;
            z -= (a11*diff0-a01*diff1)/det;
            c -= (-a10*diff0+a00*diff1)/det;
          }
#if DEBUG
        printf ("Bud root:\n");
        printf ("%+.17" UFFORMAT "\t%+.17" UFFORMAT "\n", ucreal(c), ucimag(c));
#endif
        udouble ci_cheat = ucimag(c);
        if ( ufabs(ci_cheat) < 1.e-7 )
          ci_cheat = 0;
        printf ("%ld/%ld\t%+.17" UFFORMAT "\t%+.17" UFFORMAT "\t"
                "CompRoot\t%ld@%ld/%ld\n",
                p, q, ucreal(c), ci_cheat, subper, rot, valence);
        fflush (stdout);
      }
    }
}

int
main (void)
{
  for ( long int r=1 ; r<=12 ; r++ )
    for ( long int r1=1 ; r1<=r ; r1++ )
      {
        long int r2 = r-r1;
        long int q0 = ((1<<r1)-1)<<r2;
        for ( long int p0=0 ; p0<=q0/2 ; p0++ )
          {
            long int per, preper;
            if ( p0%2==0 && q0%2 == 0 )
              continue;
            doubling_period (p0, q0, &per, &preper);
            assert (preper == r2);
            if ( per < r1 )
              continue;
            long int u = p0;
            long int v = q0;
            while ( u )
              {
                long int t = v%u;
                v = u;
                u = t;
              }
            long int p = p0/v;
            long int q = q0/v;
            process_angle (p, q);
          }
      }
  return 0;
}

See https://twitter.com/gro_tsen/status/1222934904511639552 for discussion.

https://gitlab.com/adammajewski/parameter_external_angle