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

Question

Steven Stadnicki suggested in a comment that I post the following as a question.

The golden ration $\phi$ is given by $$\phi = \frac{1+\sqrt{5}}{2} \approx 1.618033988.$$

A rational approximation is given by $$\frac{987}{610} \approx 1.618032787.$$

This yields the following pandigital formulas: $$\frac{4 * 3 + 975}{\left(8 + 2\right) * 61}$$ $$\frac{985 + 2}{1 + 7 * \left(6 + {3}^{4}\right)}$$ $$\frac{5 + 984 - 2}{\left(7 + 3\right) * 61}$$ $$\frac{987}{2 + 1 + \left({5}^{4}\right) - \left(6 * 3\right)}$$ $$\frac{7 * \left(98 + 43\right)}{\left(61 * 5\right) * 2}$$

Is there another pandigital rational representation that yields a higher number of accurate digits?

The same question may be asked for the base of the natural logarithm $e$: $$e\approx 2.718281828.$$ A rational approximation is $$\frac{1457}{536}\approx 2.718283582.$$ It produces the following pandigital formulas: $$\frac{31 * \left(5 + 42\right)}{\left(76 - 9\right) * 8}$$ $$\frac{\left(\left({3}^{6}\right) * 2\right) - 1}{\left({5}^{4}\right) - \left(97 - 8\right)}$$ $$\frac{\left(6 * \left({3}^{5}\right)\right) - 1}{7 + {\left(\left(8 * 4\right) - 9\right)}^{2}}$$ $$\frac{1 + \left(6 * \left({3}^{5}\right)\right) - 2}{\left(9 * 7 + 4\right) * 8}$$ $$\frac{31 * \left(8 * 5 + 7\right)}{{2}^{9} + 6 * 4}$$ $$\frac{5 + \left(\left({9}^{3}\right) * 2\right) - 6}{8 * \left(71 - 4\right)}$$ $$\frac{31 * \left(\left(8 * 7\right) - 9\right)}{542 - 6}$$ As before we ask for a better rational approximation that produces a higher precision.

Furthermore (see this post) we ask for a more sophisticated non-trivial algorithm to compute these. I was able to obtain an improvement of the running time by a factor of thirty by re-coding the Maple algorithm from the other post in C, including quite a few ideas that are of use to the combinatorics programmmer, like the representation of sets by bit strings or the enumeration of the permutations of a list using the factorial number system. Happy computing! (Another observation is that with the five operations and concatenation at our disposal there are probably threshold values for numerator and denominator of the rational approximations from the continued fraction where we can perhaps give a proof that they cannot be attained.)

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

int blog(int 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);
}

struct {
  int min, max;
} searchints[]  = {
  {1, 50},
  {100, 120},
  {340, 360},
  {-1, -1}
};

inline int access(int val){
  int interval = 0, offs = 0;

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

    if(val >= min && val <= max){
      return offs + val-min;
    }

    offs += max-min+1;
    interval++;
  }

  return -1;
}

int entcount = -1;

typedef struct {
  short leaf;
  long val;
  char plain[128];
  char latex[256];
} expr;

expr *memo[1<<9];

expr *makememo(void)
{
  expr *ents = (expr *)malloc(entcount*sizeof(expr));

  int ind;
  for(ind=0; ind<entcount; ind++){
    ents[ind].val = -1;
  }

  return ents;
}

expr *repr(int subsind)
{
  if(memo[subsind] != NULL){ return memo[subsind]; };

  int elcount = 0; int elements[9+1];

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

  expr *res = makememo();
  int idx;

  if(elcount == 1){
    idx = access(elements[0]);
    res[idx].leaf = 1;
    res[idx].val = elements[0];
    sprintf(res[idx].plain, "%d", elements[0]);
    sprintf(res[idx].latex, "%d", elements[0]);

    memo[subsind] = res; return res;
  }

  int count = 0;

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

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

    int upper = fact(elcount);
    for(pos=0; pos<upper; pos++){
      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], *digits = "123456789";
    for(idx=0; idx<elcount; idx++){
      buf[idx] = digits[perm[idx]-1];
    }
    buf[idx] = 0;

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

    idx = access(numval);
    if(idx != -1){
      res[idx].leaf = 1;
      res[idx].val = numval;
      sprintf(res[idx].plain, "%d", numval);
      sprintf(res[idx].latex, "%d", numval);

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

    expr *lres = repr(lidx), *rres = repr(ridx);

    if(count>=entcount){ continue; }

    int lpos, rpos;
    for(lpos=0; lpos<entcount; lpos++){
      for(rpos=0; rpos<entcount; rpos++){
      expr exl = lres[lpos], exr = rres[rpos];
      int newval, newind;

      if(exl.val != -1 && exr.val != -1){
        char plainl[256], plainr[256];
        sprintf(plainl, (exl.leaf ? "%s" : "(%s)"), exl.plain);
        sprintf(plainr, (exr.leaf ? "%s" : "(%s)"), exr.plain);

        char latexl[256], latexr[256];
        sprintf(latexl, (exl.leaf ? "%s" : "\\left(%s\\right)"), exl.latex);
        sprintf(latexr, (exr.leaf ? "%s" : "\\left(%s\\right)"), exr.latex);

        newval = exl.val + exr.val;
        if((newind = access(newval))!=-1){
          res[newind].leaf = 0;
          res[newind].val = newval;

          sprintf(res[newind].plain, "%s + %s", exl.plain, exr.plain);
          sprintf(res[newind].latex, "%s + %s", exl.latex, exr.latex);

          count++;
        }

        newval = exl.val - exr.val;
        if(newval>0 && (newind = access(newval))!=-1){
          res[newind].leaf = 0;
          res[newind].val = newval;

          sprintf(res[newind].plain, "%s - %s", plainl, plainr);
          sprintf(res[newind].latex, "%s - %s", latexl, latexr);

          count++;
        }

        if(exl.val>1 && exr.val>1){
          newval = exl.val * exr.val;
          if((newind = access(newval))!=-1){
            res[newind].leaf = 0;
            res[newind].val = newval;

            sprintf(res[newind].plain, "%s * %s", plainl, plainr);
            sprintf(res[newind].latex, "%s * %s", latexl, latexr);

            count++;
          }
        }

        if(exr.val>1 && (exl.val % exr.val) == 0){
          newval = exl.val / exr.val;
          if((newind = access(newval))!=-1){
            res[newind].leaf = 0;
            res[newind].val = newval;

            sprintf(res[newind].plain, "%s / %s", plainl, plainr);
            sprintf(res[newind].latex, "\\frac{%s}{%s}", exl.latex, exr.latex);

            count++;
          }
        }

        if(exl.val>1 && exr.val>1 && exr.val*blog(exl.val)<32){
          newval = ipow(exl.val, exr.val);
          if((newind = access(newval))!=-1){
            res[newind].leaf = 0;
            res[newind].val = newval;

            sprintf(res[newind].plain, "%s ^ %s", plainl, plainr);
            sprintf(res[newind].latex, "{%s}^{%s}", latexl, exr.latex);

            count++;
          }
        }
      }
      }
    }
  }

  memo[subsind] = res; return res;
}


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 = 0, ent;

  entcount = 0;
  while(searchints[ind].min != -1){
    entcount += searchints[ind].max-searchints[ind].min+1;
    ind++;
  }

  for(ind=0; ind<(1<<9); ind++){
    memo[ind] = NULL;
  }

  for(ind=1; ind<(1<<9)-1; ind++){
    expr *leftres = repr(ind), *rightres = repr((1<<9)-1-ind);
    int idxleft = access(p), idxright = access(q);

    if(idxleft != -1 && idxright != -1){
      expr left = leftres[idxleft];
      expr right = rightres[idxright];

      if(left.val != -1 && right.val != -1){
      printf("(%s) / (%s); %d / %d; %.15lf\n", 
             left.plain, right.plain, p, q,
             (double)left.val/(double)right.val);
      printf("\\frac{%s}{%s}\n", left.latex, right.latex);
      }
    }
  }

  return 0;
}

Regarding the golden ratio, isn't it equal to $(1+5^{4/8})/2$? — Lord Soth, Jul 23 '13 at 00:18
Yes it is, but I understood pandigital to mean that one uses all digits. — Marko Riedel, Jul 23 '13 at 00:20
But if you are looking for something that is strictly fractional (whatever that is) then of course my "spoiling the fun" line of thought may not work. — Lord Soth, Jul 23 '13 at 00:24
Yes indeed the post asks for rational approximations, so the $\sqrt{5}$ term will not be examined. — Marko Riedel, Jul 23 '13 at 00:25
What operations/functions are allowed? $+, -, \times, /, \sqrt{},\sqrt[n]{}, \log_{n}(...), \ln(...), !, $ decimal point, power, ... ??? — Oleg567, Jul 23 '13 at 00:54
I learned about this problem from the post mentioned above. It uses addition, subtraction, multiplication, division and exponentiation. Furthermore we are allowed to concatenate single digits into one value. — Marko Riedel, Jul 23 '13 at 00:56
Another feature of the problem specification that you have to decide on is whether you permit irrational values as intermediate values that reduce to a rational number in the end by cancellation. The above algorithm does not examine these. By including them you obviously extend the set of expressions that can be produced and I'm sure those types of answers would also be interesting to see. That would mean using hash tables instead of arrays in the algorithm. — Marko Riedel, Jul 23 '13 at 01:16
The choice of intervals to track in the above code is dictated by the $355/113$ approximation to $\pi$ that was discussed in the other post. — Marko Riedel, Jul 23 '13 at 01:20
The above code is not all that efficient, in particular, it should use expression trees and thereby do away with those text fields that store formulas, as these can be computed on demand by walking the trees. I will have an optimized version soon. — Marko Riedel, Jul 24 '13 at 06:04
Pandigital approximations galore at http://www2.stetson.edu/~efriedma/mathmagic/0804.html — Gerry Myerson, Sep 28 '16 at 13:54
Just by the by... Here is found by me the Pi approximation expressed via a ratio of two pandigital numbers (particular pandigital numbers in which every digit occurs exactly once): 6897253140 / 2195463870 = 3.141592642105287754063563797112270401425462765643235112... — Alex, Aug 29 '22 at 17:31

Lord Soth · Answer 1 · 2013-07-23T18:39:00.970

12

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

edited Jul 23 '13 at 18:39

answered Jul 23 '13 at 00:12

Lord Soth

7,750
20
37

Very nice and upvoted. Should have looked at the page before posting. – Marko Riedel Jul 23 '13 at 00:15
+1. Since $z$ is integer, then $(1+1/z)$ is rational, and $(1+1/z)^z$ is rational too. – Oleg567 Jul 23 '13 at 01:24
2

Since $\left(1+\frac1x\right)^x=e\left(1-\frac1{2x}+\dots\right)$, we get that the accuracy is $1$ part in $2\cdot3^{2^{85}}$ which is $$\lfloor\begin{align}2^{85}\log_{10}(3)+\log_{10}(2)\rfloor= 18457734525360901453873570\text{ digits}\end{align}$$ – robjohn Jul 31 '13 at 13:33
Of course, this rational number has $$\lfloor3^{2^{85}}2^{85}\log_{10}(3)+1\rfloor\text{ digits}$$ in the denominator. – robjohn Jul 31 '13 at 15:11

Oleg567 · Answer 2 · 2013-07-31T12:42:02.913

5

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.

edited Jul 31 '13 at 12:42

answered Jul 23 '13 at 02:03

Oleg567

17,295

Excellent. 1+ You obviously have a very good algorithm that implements fractions as intermediate results, which the above C code does not do. I might implement it later without posting it and just post some results. What programming language are you using? – Marko Riedel Jul 23 '13 at 02:38
Thanks a lot. I use C++ language. – Oleg567 Jul 23 '13 at 02:40

Marko Riedel · Answer 3 · 2013-08-01T23:41:01.920

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

For the approximation $$e\approx\frac{2721}{1001}\approx 2.718281718,$$ the above program found $$\frac{{7}^{4} + \left({2}^{6}\right) * 5}{91 * \left(8 + 3\right)}$$ after 30 minutes of computation time. — Marko Riedel, Jul 26 '13 at 03:21
The approximation $$\frac{2721}{1001}$$ can also be expressed as $$\frac{7 + 46 * 59}{1 + {\left(2 + 8\right)}^{3}}$$ and $$3 - \left(\frac{2}{7}\right) + \frac{4}{15 + 986}.$$ — Marko Riedel, Jul 27 '13 at 21:36

Marko Riedel · Answer 4 · 2013-08-11T22:51:54.303

I found the next rational approximation in the series for $\phi$, it is $$\phi \approx\frac{4181}{2584} =1.618034.$$ The pandigital formula is $$ \left(\frac{13}{8}\right) - \left(\frac{2 + 7}{5 + \left({6}^{4}\right) - 9}\right).$$ Another one is $$\left(\frac{13}{8}\right) - \left(\frac{{2}^{- 4}}{9 - \left({6}^{- \left(7 - 5\right)}\right)}\right).$$ EDIT Aug 6 Even better, using the setting


struct {
  int min, max;
} numerints[]  = {
  {1, 550},
  {800, 1000},
  {1580, 1600},
  {2570, 2590},
  {4150, 4200},
  {6720, 6820},
  {-1, -1}
}, denomints[] = {
  {1, 136},
  {2550, 2600},
  {4150, 4200},
  {-1, -1}
};

I computed the rational approximation $$\phi\approx\frac{6765}{4181}\approx 1.618033963,$$ with pandigital formula $$\frac{5 * \left(2 + 9\right)}{34 - \left(\frac{7}{861}\right)} .$$ Even though $861$ is not listed as a possible denominator it was nonetheless admitted because $7/861$ is equal to $1/123,$ which does fit the value mask.

Edit Aug 12 The following setting

struct {
  int min, max;
} numerints[]  = {
  {1, 850},
  {10920, 10980},
  {-1, -1}
}, denomints[] = {
  {1, 210},
  {6755, 6780},
  {-1, -1}
};

will produce $$\frac{{3}^{6} + \frac{2 + 9}{7 + 8}}{451},$$ which was also found by others.

score 0 · Answer 5 · answered Aug 13 '13 at 00:26

0

The approximation $$\frac{23225}{8544} \approx 2.718281835205993$$ has the pandigital representation $$\displaystyle\frac{\left({3}^{5}\right) - \left(1 + \frac{{2}^{- 4}}{\frac{6}{7}}\right)}{89}.$$

answered Aug 13 '13 at 00:26

Marko Riedel

61,317

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

5 Answers5

Linked