pandigital rational approximations to the golden ratio and the base of the natural logarithm

Solution 1:

Here is one pandigital formula (found by R. Sabey in 2004 according to http://mathworld.wolfram.com/eApproximations.html) that is accurate to $18457734525360901453873570$ digits (so they say): \begin{align} \left(1+9^{-4^{7\cdot 6}}\right)^{3^{2^{85}}}, \end{align} which is just \begin{align} \left(1+1/z\right)^{z} \end{align} with $z=3^{2^{85}}$ in disguise. This relies on the fact that $e = \lim_{x\rightarrow\infty}(1+\frac{1}{x})^x$.

Solution 2:

My best rational approximation of $\large\phi$ $(1.618033988749894848204586834365...)$ is $$ \phi \approx \dfrac{6}{3-\dfrac{2}{7}+\left(1-8 \cdot 9^{-4}\right)^5} = 1.618033988\color{Tan}{913633}...; \tag{1} $$ error is $\Delta = 1.637385... \times 10^{-10}$.

And I'll show other nice $\phi$ approximations: $$ \phi\approx \dfrac{736}{\dfrac{8^4}{9}-\dfrac{5}{21}} = \dfrac{46368}{28657}=\color{blue}{\dfrac{F_{24}}{F_{23}}} = 1.618033988\color{Tan}{205325}... , \tag{2} $$ $$ \phi \approx {\large\dfrac{3^6+\frac{2+9}{8+7}}{451}} = \dfrac{10946}{6765}=\color{blue}{\dfrac{F_{21}}{F_{20}}} = 1.6180339\color{Tan}{985218}... , \tag{3} $$ where $\color{blue}{F_n}$ $-$ Fibonacci numbers.

Solution 3:

I am sending an amazingly fast C program that remedies many of the defects of the intial version that I posted, eliminating those wasteful text fields among others. Using this program I was able to compute pandigital representations of the following approximations to the golden ratio: $$\phi\approx \frac{1597}{987}\approx 1.618034448$$ in matter of minutes (!) on a laptop of modest means (3GB memory, dual 2,16 GHz processor).

The approximations obtained are $$ \frac{4 + \left(\frac{6}{2}\right) * 531}{987}$$ and $$ \frac{\left({\left(8 * 5\right)}^{2}\right) - 3}{\left(\left(9 + 4\right) * 76\right) - 1}.$$

The reader is invited to try this program herself. Where speed is concerned it can search even large ranges very quickly, being limited only by main memory. (Obviously the speed suffers when the machine has to swap data out to disk.)

I suspect you may be able to find the next approximation, again a quotient of Fibonacci numbers, which is $$\frac{2584}{1597}.$$ Then again, maybe not, as the density of pandigital formulas thins out as we increase the integers that are involved. Do try for it. It only takes an adjustment of the ranges for the numerators being searched and the denominators at the top of the program. Enjoy! A pandigital formula for this last approximation would be quite excellent. Compiled with GCC 4.7.1. EDIT as of Sat Jul 27 Uploaded a bug fix where the previous version would miss entries. EDIT II as of Sat Jul 27 Fixed a bug where not all numbers generated by concatenation of single digits where being examined. This code should now be bug-free. EDIT III as of Fri Aug 02 reduced the number of expressions being examined.

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <time.h>

int maxents = -1, numerents = -1, denoments = -1, diffents = -1;


struct {
  int min, max;
} numerints[]  = {
  {1, 400},
  {2550, 2600},
  {4150, 4200},
  {-1, -1}
}, denomints[] = {
  {1, 49},
  {1280, 1300},
  {2550, 2600},
  {-1, -1}
};




int donecount = 0;
time_t starttime, ticktime;

int blog(long n) {
  assert(n > 0);

  int res = -1; 

  while(n>0){
    n >>= 1;
    res++;
  }

  return res;
}

int fact(int n)
{
  assert(n >= 0);

  if(n<2) return 1;

  return n*fact(n-1);
}

long ipow(int m, int n)
{
  assert(n >= 0);

  if(n==0) return 1;
  if(n==1) return m;

  long r = (n%2 ? ipow(m, (n-1)/2) : ipow(m, n/2));
  return (n%2 ? m*r*r : r*r);
}

typedef enum {
  NODE_LEAF = 0,
  NODE_ADD,
  NODE_SUB,
  NODE_MUL,
  NODE_DIV,
  NODE_POW,
  NODE_INVPOW
} NODE_TYPE;

typedef struct _expr {
  NODE_TYPE what;
  long numer, denom;
  struct _expr *left, *right;
} expr;

long gcd(long a, long b)
{
  if(a < b) return gcd(b, a);
  if(b == 0) return a;

  return gcd(b, a%b);
}

void init_expr(expr *ex, long numer, long denom)
{
  long d = gcd(numer, denom);

  ex->numer = numer/d;
  ex->denom = denom/d;
}

void add_expr(expr *a, expr *b, expr *res)
{
  long p, q, d;

  p = a->numer*b->denom + b->numer*a->denom;
  q = a->denom*b->denom;

  d = gcd(abs(p), abs(q));

  res->numer = p/d;
  res->denom = q/d;
}

void sub_expr(expr *a, expr *b, expr *res)
{
  long p, q, d;

  p = a->numer*b->denom - b->numer*a->denom;
  q = a->denom*b->denom;

  d = gcd(abs(p), abs(q));

  res->numer = p/d;
  res->denom = q/d;
}

void mul_expr(expr *a, expr *b, expr *res)
{
  long p, q, d;

  p = a->numer*b->numer;
  q = a->denom*b->denom;

  d = gcd(abs(p), abs(q));

  res->numer = p/d;
  res->denom = q/d;
}


void div_expr(expr *a, expr *b, expr *res)
{
  long p, q, d;

  p = a->numer*b->denom;
  q = a->denom*b->numer;

  d = gcd(abs(p), abs(q));

  res->numer = p/d;
  res->denom = q/d;
}

void ipow_expr(expr *a, int x, expr *res)
{
  if(x==0){
    res->numer = 1;
    res->denom = 1;
    return;
  }

  if(x==1){
    res->numer = a->numer;
    res->denom = a->denom;
    return;
  }

  int est = (blog(a->numer)-blog(a->denom))*x;
  assert(est >= -32 || est <= 32);

  long p, q;

  p = ipow(a->numer, x);
  q = ipow(a->denom, x);

  res->numer = p;
  res->denom = q;
}


typedef struct {
  int count;
  expr **data;
  char *subset;
} resent;
resent memoaux[1<<9];

inline expr *access_frac(int subsind, long numer, long denom)
{
  int interval, numeroffs = 0, denomoffs = 0;

  interval = 0;
  while(numerints[interval].min !=-1){
    int min = numerints[interval].min;
    int max = numerints[interval].max;

    if(numer >= min && numer <= max){
      numeroffs += numer-min;
      break;
    }

    numeroffs += max-min+1;
    interval++;
  }
  if(numerints[interval].min == -1) return NULL;

  interval = 0;
  while(denomints[interval].min !=-1){
    int min = denomints[interval].min;
    int max = denomints[interval].max;

    if(denom >= min && denom <= max){
      denomoffs += denom-min;
      break;
    }

    denomoffs += max-min+1;
    interval++;
  }
  if(denomints[interval].min == -1) return NULL;


  return &(memoaux[subsind].data[numeroffs][denomoffs]);
}

inline expr *access_ind(int subsind, int ind)
{
  long numeroffs = ind/denoments, denomoffs = ind%denoments;

  return &(memoaux[subsind].data[numeroffs][denomoffs]);
}


void repr(int subsind)
{
  if(memoaux[subsind].count != -1) return;

  int elcount = 0; int elements[9+1]; 
  char buf[16], *digits = "123456789"; // for debugging

  int pos, bitpat = subsind;
  for(pos=0; pos<9; pos++){
    if(bitpat & 1 == 1){
      buf[elcount] = digits[pos];
      elements[elcount++] = pos+1;
    }
    bitpat >>= 1;
  }
  elements[elcount] = -1;
  buf[elcount] = 0;

  memoaux[subsind].subset = strdup(buf);

  resent *res = memoaux + subsind;
  int idx;

  if(elcount == 1){
    long numer = elements[0], denom = 1;
    expr *ex;

    if((ex = access_frac(subsind, numer, denom)) != NULL){
      ex->what = NODE_LEAF;

      ex->left = NULL;
      ex->right = NULL;

      ex->numer = numer;
      ex->denom = denom;

      memoaux[subsind].count = 1;
    }
    else{
      memoaux[subsind].count = 0;
    }

    donecount++;
    return;
  }

  memoaux[subsind].count = 0;

  if(elcount >= 2 && elcount <= 4){
    int perm[9+1], el, rpos;

    int upper = fact(elcount);
    for(pos=0; pos<upper; pos++){
      for(rpos=0; rpos<elcount; rpos++){ 
      perm[rpos] = elements[rpos];
      }

      idx = pos;
      int el;
      for(el=elcount; el>0; el--){
      int targ = idx % el;

      int tmp = perm[targ];
      perm[targ] = perm[el-1];
      perm[el-1] = tmp;

      idx /= el;
      }

      char buf[9+1];
      for(idx=0; idx<elcount; idx++){
      buf[idx] = digits[perm[idx]-1];
      }
      buf[idx] = 0;

      int numval; sscanf(buf, "%d", &numval);

      expr *ex;
      if((ex = access_frac(subsind, numval, 1)) != NULL){
      ex->what = NODE_LEAF;

      ex->left = NULL;
      ex->right = NULL;

      ex->numer = numval;
      ex->denom = 1;

      memoaux[subsind].count++;
      }
    }
  }

  int subsubind;
  for(subsubind=1; subsubind<(1<<elcount)-1; subsubind++){
    int leftel[9+1], rightel[9+1], leftpos=0, rightpos=0;

    bitpat = subsubind;
    for(pos=0; pos<elcount; pos++){
      if(bitpat & 1 == 1){
      leftel[leftpos++] = elements[pos];
      }
      else{
      rightel[rightpos++] = elements[pos];
      }

      bitpat >>= 1;
    }

    int lidx = 0, ridx = 0, innerpos;
    for(innerpos=0; innerpos<leftpos; innerpos++){
      lidx += 1<<(leftel[innerpos]-1);
    }
    for(innerpos=0; innerpos<rightpos; innerpos++){
      ridx += 1<<(rightel[innerpos]-1);
    }

    repr(lidx); repr(ridx);

    if(memoaux[subsind].count >= diffents) continue;

    int lpos, rpos, ldone, rdone;
    for(lpos=0, ldone=0; lpos<maxents; lpos++){
      if(ldone >= memoaux[lidx].count) break;

      expr *exl = access_ind(lidx, lpos);

      if(exl == NULL) continue;
      if(exl->numer == -1) continue;
      ldone++;

      for(rpos=0, rdone=0; rpos<maxents; rpos++){
      if(rdone >= memoaux[ridx].count) break;

      expr *exr = access_ind(ridx, rpos);

      if(exr == NULL) continue;
      if(exr->numer == -1) continue;
      rdone++;

      expr cand, *success;

      add_expr(exl, exr, &cand);
      if((success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
         (success->numer==-1 && success->denom==-1)){
        success->numer = cand.numer; 
        success->denom = cand.denom;

        success->what = NODE_ADD;

        success->left = exl;
        success->right = exr;

        memoaux[subsind].count++;
      }

      sub_expr(exl, exr, &cand);
      if(cand.numer>0 && cand.denom>0 && 
         (success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
         (success->numer==-1 && success->denom==-1)){
        success->numer = cand.numer; 
        success->denom = cand.denom;

        success->what = NODE_SUB;

        success->left = exl;
        success->right = exr;

        memoaux[subsind].count++;
      }

      if(!((exl->numer==1 && exl->denom==1) || (exr->numer==1 && exr->denom==1))){
        mul_expr(exl, exr, &cand);
        if((success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
           (success->numer==-1 && success->denom==-1)){
          success->numer = cand.numer; 
          success->denom = cand.denom;

          success->what = NODE_MUL;

          success->left = exl;
          success->right = exr;

          memoaux[subsind].count++;
        }
      }

      if(!(exr->numer==1 && exr->denom==1)){
        div_expr(exl, exr, &cand);
        if((success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
           (success->numer==-1 && success->denom==-1)){
          success->numer = cand.numer; 
          success->denom = cand.denom;

          success->what = NODE_DIV;

          success->left = exl;
          success->right = exr;

          memoaux[subsind].count++;
        }
      }

      if(!((exl->numer==1 && exl->denom==1) || (exr->numer==1 && exr->denom==1))
         && exr->denom==1){
        ipow_expr(exl,exr->numer, &cand); 
        if((success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
           (success->numer==-1 && success->denom==-1)){
          success->numer = cand.numer; 
          success->denom = cand.denom;

          success->what = NODE_POW;

          success->left = exl;
          success->right = exr;

          memoaux[subsind].count++;
        }
      }

      if(!((exl->numer==1 && exl->denom==1) || (exr->numer==1 && exr->denom==1))
         && exr->denom==1){
        expr inv;
        inv.numer = exl->denom;
        inv.denom = exl->numer;

        ipow_expr(&inv,exr->numer, &cand); 
        if((success = access_frac(subsind, cand.numer, cand.denom)) != NULL &&
           (success->numer==-1 && success->denom==-1)){
          success->numer = cand.numer; 
          success->denom = cand.denom;

          success->what = NODE_INVPOW;

          success->left = exl;
          success->right = exr;

          memoaux[subsind].count++;
        }
      }
      }
    }
  }

  donecount++;
  if(donecount % 24 == 0){
    time_t now = time(NULL); 
    double avg = ((double)now-(double)starttime)/(double)donecount;
    int todo = ((1<<9)-1-donecount);
    double rem1 = 
      (double)(now-ticktime)/(double)(24*60)*(double)todo;
    double rem2 =
      (double)(now-starttime)/(double)donecount*(double)todo/(double)60;
    fprintf(stderr, "# %03d subsets compl., "
          "avg. %07.3lf sec. per set, "
          "time rem. %07.3lf / %07.3lf min\n",
          donecount, avg, rem1, rem2);
    ticktime = now;
  }

  return;
}

void print_plain(expr *node, int parens)
{
  if(node->what == NODE_LEAF){
    if(node->denom == 1){
      printf("%ld", node->numer);
    }
    else{
      printf("%ld/%ld", node->numer, node->denom);
    }

    return;
  }

  if(parens && node->what != NODE_LEAF) printf("(");

  switch(node->what){
  case NODE_ADD:
    print_plain(node->left, 0);
    printf(" + ");
    print_plain(node->right, 0);

    break;
  case NODE_SUB:
    print_plain(node->left, 1);
    printf(" - ");
    print_plain(node->right, 1);

    break;
  case NODE_MUL:
    print_plain(node->left, 1);
    printf(" * ");
    print_plain(node->right, 1);

    break;
  case NODE_DIV:
    print_plain(node->left, 1);
    printf(" / ");
    print_plain(node->right, 1);

    break;
  case NODE_POW:
    print_plain(node->left, 1);
    printf(" ** ");
    print_plain(node->right, 1);

    break;
  case NODE_INVPOW:
    print_plain(node->left, 1);
    printf(" ** (- ");
    print_plain(node->right, 1);
    printf(")");

    break;
  default: assert(node->what>=NODE_ADD && node->what<=NODE_INVPOW);
  }

  if(parens && node->what != NODE_LEAF) printf(")");
}


void print_latex(expr *node, int parens)
{
  if(node->what == NODE_LEAF){
    if(node->denom == 1){
      printf("%ld", node->numer);
    }
    else{
      printf("%ld/%ld", node->numer, node->denom);
    }

    return;
  }

  if(parens && node->what != NODE_LEAF) printf("\\left(");

  switch(node->what){
  case NODE_ADD:
    print_latex(node->left, 0);
    printf(" + ");
    print_latex(node->right, 0);

    break;
  case NODE_SUB:
    print_latex(node->left, 1);
    printf(" - ");
    print_latex(node->right, 1);

    break;
  case NODE_MUL:
    print_latex(node->left, 1);
    printf(" * ");
    print_latex(node->right, 1);

    break;
  case NODE_DIV:
    printf("\\frac{");
    print_latex(node->left, 0);
    printf("}{");
    print_latex(node->right, 0);
    printf("}");

    break;
  case NODE_POW:
    printf("{");
    print_latex(node->left, 1);
    printf("}^{");
    print_latex(node->right, 0);
    printf("}");

    break;
  case NODE_INVPOW:
    printf("{");
    print_latex(node->left, 1);
    printf("}^{- ");
    print_latex(node->right, 1);
    printf("}");

    break;
  default: assert(node->what>=NODE_ADD && node->what<=NODE_INVPOW);
  }

  if(parens && node->what != NODE_LEAF) printf("\\right)");
}

int main(int argc, char **argv)
{
  long p, q;

  if(argc!=3){
    fprintf(stderr, "usage: %s <p> <q>\n", argv[0]);
    exit(-1);
  }

  sscanf(argv[1], "%ld", &p);
  sscanf(argv[2], "%ld", &q);

  if(p<1 || q<1){
    fprintf(stderr, 
          "positive integers please, got %ld and %ld\n",
          p, q);
    exit(-2);
  }

  int ind, ind2, offs;

  numerents = 0; ind = 0;
  while(numerints[ind].min != -1){
    int count = numerints[ind].max-numerints[ind].min+1;
    numerents += count; 

    ind++;
  }

  denoments = 0; ind = 0;
  while(denomints[ind].min != -1){
    int count = denomints[ind].max-denomints[ind].min+1;
    denoments += count;

    ind++;
  }

  maxents = numerents*denoments;

  diffents = 0;
  for(ind=0; numerints[ind].min != -1; ind++){
    for(ind2=0; denomints[ind2].min != -1; ind2++){
      int x, y;
      for(x=numerints[ind].min; x<=numerints[ind].max; x++){
      for(y=denomints[ind2].min; y<=denomints[ind2].max; y++){
        if(gcd(x, y) == 1) diffents++;
      }
      }
    }
  }


  fprintf(stderr, "# initializing "
        "(%d numerators and %d denominators, "
        "total %d, unique %d)\n"
        "# ", numerents, denoments, maxents, diffents);

  for(ind=0; ind<(1<<9); ind++){
    resent *res = &(memoaux[ind]);

    res->subset = NULL;

    res->data = (expr **)malloc(numerents*sizeof(expr *));
    assert(res->data != NULL);

    int pos;
    for(pos=0; pos<numerents; pos++){
      res->data[pos] =
      (expr *)malloc(denoments*sizeof(expr));
      assert(res->data[pos] != NULL);
    }

    if(ind>0 && ind%24==0) { fprintf(stderr, "."); fflush(stderr); }

    int dataind = 0;
    for(dataind=0; dataind<maxents; dataind++){
      int numer = dataind/denoments, denom = dataind%denoments;

      res->data[numer][denom].numer = -1;
      res->data[numer][denom].denom = -1;
    }

    res->count = -1;
  }
  fprintf(stderr, "\n");

  starttime = time(NULL); ticktime = starttime;

  for(ind=1; ind<(1<<9)-1; ind++){
    long lidx = ind, ridx = (1<<9)-1-ind;
    repr(lidx); repr(ridx);

    expr *left = access_frac(lidx, p, 1);
    if(left == NULL || left->numer == -1) continue;

    expr *right = access_frac(ridx, q, 1);
    if(right == NULL || right->numer == -1) continue;

    assert(left->denom == 1 && right->denom == 1);

    expr verify; div_expr(left, right, &verify);

    print_plain(left, 1); printf(" / "); print_plain(right, 1);
    printf("; %ld / %ld; ", left->numer, right->numer);
    printf("%.15lf\n", (double)verify.numer/(double)verify.denom);

    printf("\\frac{");
    print_latex(left, 0);
    printf("}{");
    print_latex(right, 0);
    printf("}\n");
  }

  int allidx = (1<<9)-1;
  repr(allidx);

  expr *result = access_frac(allidx, p, q);
  if(result != NULL && (result->numer != -1 && result->denom !=-1)){
    print_plain(result, 0); 
    printf("; %ld / %ld; %.15lf\n",
         result->numer, result->denom,
         (double)result->numer/(double)result->denom);
    print_latex(result, 0); printf("; *\n");
  }

  return 0;
}