\\ Kevin Ryde, November 2021
\\
\\ Usage; gp <a014121.gp
\\
\\ This is some PARI/GP code for calculating terms of two sequences
\\
\\      A014121 = numbers of the form abs(2^i-3^j)
\\      A128760 = number of ways to write n in the form abs(2^i-3^j)
\\
\\ The strategy is a search of i,j seeking differences <= a target limit.
\\
\\ Is there some theory on how big i and j could be and get abs(2^i-3^j)
\\ small compared to 2^i and 3^j ?  The search range (of j) chosen is hoped
\\ to be a very generous over-estimate, but without proof.
\\

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

\\ all_differences(limit) returns a vector of differences abs(2^i-3^j) which
\\ are <= limit, for i>=0,j>=0, including duplicates from different i,j
\\ pairs, and in no particular order.
\\
all_differences(limit) =
{
  my(ret=List([]),
     three_j=1);  \\ 3^j
  for(j=0, if(limit,2*logint(limit,3)) + 15,
     my(i = logint(three_j,2),  \\ so 2^i just below 3^j
        two_i = 1<<i);          \\ 2^i

     \\ 3^j - 2^i
     while(two_i > 0,
           my(t = three_j - two_i);  \\ increasing differences
           if(t>limit,break);
           listput(ret,t);
           two_i>>=1);

     \\ 2^i - 3^j
     two_i = 1<<i;
     while(1,
           two_i<<=1;
           my(t = two_i - three_j);  \\ increasing differences
           if(t>limit,break);
           listput(ret,t));

     three_j*=3);

  Vec(ret);
}

\\ Return a vector of all terms of A014121 which are <= "limit".
\\
A014121_distinct_differences(limit) = Set(all_differences(limit));

\\ Return a vector of terms A128760(n) for n=0 through to n="limit"
\\ inclusive.  So v[1..limit+1] = A128760(0..n).
\\
A128760_num_ways(limit) =
{
  my(diffs=all_differences(limit),
     v=vector(limit+1));
  for(i=1,#diffs, v[diffs[i]+1]++);
  v;
}

{
  print("A014121 differences abs(2^i-3^j)");
  my(v=A014121_distinct_differences(70));
  print1("  = ");
  for(t=1,#v, print1(v[t],","));
  print("...");
}{
  print("A128760 number of ways n = abs(2^i-3^j), for n>=0");
  my(v=A128760_num_ways(30));
  print1("  = ");
  for(t=1,#v, print1(v[t],","));
  print("...");
}

\\-----------------------------------------------------------------------------
\\ The strategy in all_differences() is to take each power 3^j and consider
\\ powers 2^i near it.  i = logint(3^j,2) gives the power 2^i equal to or
\\ closest below 3^j.  This is i=A056576(j) and 2^i=A098232(j).
\\
\\ Working i downwards gives differences 3^j - 2^i for the return.
\\ Differences increase as i decreases so stop on exceeding "limit" or at
\\ i=0.
\\
\\ Working i upwards gives differences 2^i - 3^j for the return.
\\ Differences increase as i increases so stop on exceeding "limit".
\\
\\ A014121_distinct_differences() un-duplicates all_differences().
\\ Duplicates include, after the first appearance of each,
\\
\\     1, 1, 1, 5, 5, 7, 7, 13, 23
\\     
\\ These are where A128760(n) >= 2.  Reinhard Zumkeller in A128760
\\ conjectures there may exist some n beyond which no more duplicates.
\\ Trying terms to limit=2^1024 finds no more.  On that basis there's no
\\ progressive un-duplicating here, just leave it to the end.


\\ Differences increasing away from 3^j can be illustrated,
{
  my(j=6);
  3^j == 729           || error();
  my(i=logint(3^j,2));
  i == 9               || error();
  2^i == 512           || error();

  2^12 - 3^j  == 3367  || error();
  2^11 - 3^j  == 1319  || error();
  2^10 - 3^j  == 295   || error();  \\ growing differences upwards

  3^j - 2^9   == 217   || error();  \\ growing differences downwards
  3^j - 2^8   == 473   || error();
  3^j - 2^7   == 601   || error();
}


\\ The maximum j for the outer 3^j loop is based on what ought to be a
\\ generous over-estimate.
\\
\\     j <= 2*logint(limit,3) + 15
\\
\\ This means 3^j double the length in digits (or bits) of "limit", plus a
\\ fixed 15 extra head-room.
\\
\\ For the 3^j - 2^i side, the problem is when 3^j has a lot of high 0-bits.
\\ For example,
{
  my(j=12, i=19);
  3^j     == 531441  || error();
  2^i     == 524288  || error();
  3^j-2^i ==   7153  || error();
  binary(3^j)     == [1, 0,0,0,0,0,0, 1,1,0,1,1,1,1,1,1,0,0,0,1]  || error();
  binary(2^i)     == [1, 0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0]  || error();
  binary(3^j-2^i) == [                1,1,0,1,1,1,1,1,1,0,0,0,1]  || error();
}
\\ Here 3^12 has 6 high 0-bits so that subtracting 2^i leaves a relatively
\\ small difference 7153 from what was a larger power 531441.  Or
\\ equivalently, it happened 3^j is only a little bigger than a 2^i.
\\
\\ For the 2^i - 3^j side, the problem is instead 3^j having a lot of high
\\ 1-bits.  A small example illustrates,
{
  my(j=5, i=8);
  3^j     == 243  || error();
  2^i-3^j ==  13  || error();
  binary(2^i)     == [1, 0,0,0,0, 0,0,0,0]  || error();
  binary(3^j)     == [   1,1,1,1, 0,0,1,1]  || error();
  binary(2^i-3^j) == [            1,1,0,1]  || error();
}
\\ 2^i in this case is immediately above the bits of 3^j so the subtract
\\ bit-flip the high 1s of 3^j, clearing them to 0s and so relatively small
\\ difference 13 = binary 1101.


\\ The "2*length + 15" limit chosen means an unexamined 3^j would have to be
\\ a bit pattern like
\\
\\     3^j = binary  10000000000000  000000000000   [.....any......]
\\     length:          >= "limit"    15*log2(3)       "limit"
\\
\\ This is 3^j comprising at least half plus 15 ternary digits worth of high
\\ 0-bits (or similarly 1-bits).  When j is bigger, the proportion of 0s out
\\ of total length becomes yet greater to make the "any" part <= limit.
\\
\\ Inspecting some bit patterns of 3^j suggests the number of high 0-bits
\\ (or 1-bits) grows without bound, but slowly, and seemingly much slower
\\ than total bit length and hence decreasing proportion of total length --
\\ as opposed to increasing (somewhere) which would have to happen to evade
\\ this search range.

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