\\ Kevin Ryde, June 2022

\\ This is some code finding terms of A354466, which is numbers whose sum of
\\ reciprocals of digits,
\\
\\     term = d[1] d[2] ... d[t]  integer decimals
\\     f = 1/d[1] + 1/d[2] + ... + 1/d[t]
\\
\\ has decimal expansion of f starting with those same digits of the term.
\\ The decimal point can be after d[1] or anywhere later (so f>=1),
\\
\\     f = d[1] d[2] . d[3] ... d[t] ...
\\
\\ The strategy here is to try all possible fractions f.  Terms found this
\\ way are not in ascending order.  a_vector() finds terms up to a target
\\ maximum and sorts.
\\
\\ For use in PARI/GP 2.9.x and higher.


default(strictargs,1);
default(recover,0);


\\ cden = LCM of digits 1..9, so a common denominator for any sum of 1/d.
cden = lcm([1..9]);

\\ f_periodic(f) takes a fraction f ranging 0 <= f < 1.
\\ If the decimal digits of f are purely periodic starting from the digit
\\ immediately below the decimal point then return a vector of those digits.
\\ If not then return symbol 'not.
\\     f_periodic(1/7) == [1,4,2,8,5,7]
\\
f_periodic(f) =
{
  \\ This is a simple implementation watching for a repeat of a remaining
  \\ digits part and is enough for the uses here.

  my(seen=Map(),
     ret=List([]),
     prev);
  while(!mapisdefined(seen,f,&prev),
        mapput(seen,f,#ret);
        my(d=floor(f*=10));
        listput(ret,d);
        f-=d);
  if(prev==0, Vec(ret), 'not);
}

\\ periodic_recips_and_digs is a Map() of
\\     f  =>  [recips, digs]
\\ which is when fraction f is in the range 1/10 <= f < 1, has purely
\\ periodic decimal digits starting immediately below the decimal point,
\\ and no 0 digits among those.
\\
\\ recips is sum(1/d) of those digits.
\\ digs is how many digits.
\\
periodic_recips_and_digs = Map();

{
  \\ initialise periodic_recips_and_digs
  my(verbose=0);
  if(verbose,
     print("periodic_recips_and_digs /",cden));
  for(r=1,cden-1,
     my(f=r/cden,
        v=f_periodic(f));
     if(v=='not || vecmin(v)==0, next);
     my(sd=[sum(i=1,#v, 1/v[i]),   \\ recips
            #v]);                  \\ digs
     mapput(periodic_recips_and_digs,f,sd);
     if(verbose,
        printf("  r=%-4d f=%-5s  %-20s  sd=%s\n",
               r,f,v,sd)));

  if(verbose,
     print(#periodic_recips_and_digs," many periodic r\n"));
}

\\ f_matching() takes a fraction f and if it equals the sum of reciprocals
\\ of some number of its own initial decimal digits, then roturn those
\\ digits as an integer.
\\ Return symbol 'not if not or if more than maxdigs.
\\
\\ The implemenation demands f >= 1/10000, which suffices for the uses here
\\ and anything smaller would be 'not anyway.
\\
f_matching(f,maxdigs=oo) =
{
  \\ "scale" is enough factors of 10 to reach 1/10 <= f*10^scale < 1
  \\ r is the residue after taking some decimal digits from the start of f.
  \\ r ranges 1/10 <= r < 1 and the next digit is immediately below the
  \\ decimal point in r.
  \\
  \\ s = f - sum(1/d) for the digits seen so far.
  \\ If s==0 is reached then f is indeed its own initial etc.
  \\
  \\ if r is in periodic_recips_and_digs then it has a perodic part starting
  \\ immediately below the deicmal point.
  \\ Divide s so as to subtract q copies of those, and r is unchanged.
  \\ If q==0 then no copies as s is already small and will have only a few
  \\ more one-by-one digits until an answer.
  \\
  \\ if r is not in periodic_recips_and_digs then have not reached the
  \\ periodic part of f yet, or have reached but it has a 0 digit.
  \\ In both cases a few more one-by-one steps reach the periodic part or
  \\ the disallowed 0.

  my(digs=0, \\ number of digits for return
     scale = if(f,3-logint(floor(f*10000),10)),
     s = f,
     r = f*10^scale);

  while(s>0 && digs < maxdigs,
        my(sd);
        if(mapisdefined(periodic_recips_and_digs,r,&sd),
           my(q); [q,s]=divrem(s,sd[1]);
           digs += q*sd[2];
           if(q,next));

        my(d=floor(r*=10));  \\ digit 0..9
        if(d==0, return('not));
        r -= d;
        s -= 1/d;
        digs++);

  if(s, 'not,
     floor(f*10^(digs+scale)));
}

\\ Return a vector of all terms of A354466 which are <= maxterm.
\\
a_vector(maxterm) =
{
  my(maxdigs=if(maxterm,1+logint(maxterm,10)),
     ret=List([]));
  forstep(f=1, maxdigs, 1/cden,
     my(m=f_matching(f));
     if(m!='not && m <= maxterm,
        listput(ret,m)));

  vecsort(Vec(ret));
}

{
  my(v=a_vector(10^5));
  print1("A354466 = ");
  for(i=1,#v,
     print1(v[i],", "));
  print("...");
}

\\ Fractions f tried have denominator cden = 2520 since that is the common
\\ denominator for any sum of digit reciprocals, so nothing else could
\\ succeed.
\\
\\ maxdigs is the number of decimal digits of maxterm.
\\ The largest possible reciprocal for one digit is 1/1 and if all of
\\ maxdigs are 1 then that would be f = 1/1 * maxdigs.
\\ So forstep() doesn't need to go beyond f = maxdigs.
\\
\\ This bound is approached by a fraction like f = 911..11 + 1/9 which has
\\ digit length digs = Ceil(f), ie. f only a little smaller than digs.
\\
\\ About 2520*maxdigs fractions are considered.  The speed in f_matching()
\\ is reasonable, and suffices for a b-file of maxterm = 10^500.
\\
\\ Speed in f_matching() is helped mainly by noticing the periodic part and
\\ taking repeats of it by a division instead of individually.
\\ 
\\ An option not taken for a little more speed could be stay in integers by
\\ r = f*cden and digits by divrem cden*10^highpos.  But there should be far
\\ better possible by considering the integer part of f.  A given integer
\\ part looks like it has only a few options for repeating part (far fewer
\\ than all 2520 fraction parts).


\\-----------------------------------------------------------------------------
\\ Predicate
\\
\\ By way of comparison, a simple predicate function is as follows.
\\ f is multiplied by a suitable power of 10 so that its highest non-0 digit
\\ is aligned to the high digit of m, ready to compare against m.

\\ Return 1 if m is a term of A354466.
is(m) =
{
  my(v=digits(m));
  vecmin(v) || return(0);  \\ contains 0 digit
  my(f=vecsum(apply(d->1/d, v)));
  if(f<1, return(0));
  m == floor(f*10^(#v-1 - logint(floor(f),10)));
}

{
  \\ is() predicate vs a_vector()
  my(maxterm = 10^5,
     by_pred = List([]));
  for(m=1,maxterm, if(is(m), listput(by_pred,m)));
  by_pred = Vec(by_pred);
  my(by_vec = a_vector(maxterm));
  by_pred == by_vec  || error();
}

\\-----------------------------------------------------------------------------
print("end");