/* Program to compute terms of OEIS sequence A88995
   by Robert Munafo, 25th March 2025

    cc a88995.c -lm -o seq88995 && ./seq88995
*/
#include <math.h>
#include <stdio.h>

int main(int argc, char **argv)
{
  long double k,
              l2, l5, /* log of 2 and of 5 */
              p2, p5, /* log of 2^k and of 5^k */
              p2f, p5f, /* fractional part of p2, p5 */
              p2e, p5e, /* base-10 exponent for printing 2^k and 5^k */
              thres;  /* threshold for quick check */
  int i,
      recd,  /* current record number of matched digits */
      elen,  /* exponent length (for formatting and knowing when to stop) */
      going;
  char d2s[50]; char d5s[50]; char fmt1[20];
  recd = 0;
  /* For two numbers to agree in the first 2 digits, the fractional part
     of their log base 10 definitely needs to be less than log_10(1.1).
     For example, this requirement is met for the numbers 10000 and 10999,
     whose logs base 10 differ by log_10(10999)-log_10(10000) = log_10(1.0999)
     However we want a margin of error to avoid roundoff problems, so
     we'll use log_10(1.2). We will also tighten the requirement as
     the current "record" increases. */
  thres = log10l(1.2) / powl(10.0, (long double) (recd-2));
  l2 = log10l(2.0); l5 = log10l(5.0);
  /* Try all exponents up to 9 digits long */
  for(going=1, k=1.0; going && (k<=999999999.0); k+=1.0) {
    /* Using multiplication, "floor", and subtraction we can quickly
       filter out most of the failing exponents */
    p2 = k*l2; p5 = k*l5;
    p2e = floorl(p2); p5e = floorl(p5);
    p2f = p2-p2e; p5f = p5-p5e;
    if ((p2f-p5f < thres) && (p5f-p2f < thres)) {
      /* The logarithms are within the margin, but the numbers might not
         actually agree in all the digits we need to set a new record;
         so we now look at the actual digits */
      p2f = powl(10.0, p2f); p5f = powl(10.0, p5f);
      /* Count the number of digits in the exponent */
      elen = (int) floorl(log10l(p5e) + 1.0);
      /* If the exponent plus needed digit record is now too big for the
         precision, we'll stop. */
      if (elen + recd >= 17) {
        printf("At k=%d, 5^k > 10^%d : insufficient precision to continue\n",
                                                         (int) k, (int) p5e);
        going = 0;
      } else {
        /* Get the mantissas in decimal form as a string */
        snprintf(d2s, 49, "%.18Lf", p2f);
        snprintf(d5s, 49, "%.18Lf", p5f);
        /* Now count how many identical characters we see */
        for(i=0; d2s[i] && d5s[i] && (d2s[i] == d5s[i]); i++) {   }
        if (i>1) { i--; } /* don't count the '.' */
        if (i > recd) {
          recd++;
          while (recd < i) {
            printf(" %12s %2d (same as k=%d)\n", "", recd, (int) k);
            recd++;
          }
          /* We'll pretty-print the numbers in a way that shows only as
             many mantissa digits as can be reliably extracted from our
             'long double' data type */
          snprintf(fmt1, 19, "%%.%dLf...e%%d", 18-elen);
          snprintf(d2s, 49, fmt1, p2f, (int)p2e);
          snprintf(d5s, 49, fmt1, p5f, (int)p5e);
          printf(" k=%10d %2d 2^k=%-25s 5^k=%s\n", (int) k, recd, d2s, d5s);
          /* A new record means we can use a tighter threshold for the
             quick test */
          thres = log10l(1.2) / powl(10.0, (long double) (recd-2));
        }
      }
    }
  }
  return 0;
}