\\ Kevin Ryde, July 2023
\\
\\ This is some PARI/GP code to calculate terms of A364133,
\\ which are those n where A000127(n) has a new record high
\\ 2-adic valuation (number of low 0 bits).
\\
\\ Terms n < 2^num_bits are found for given "num_bits".
\\ The strategy is a depth-first search through the bits of n
\\ from low to high and pruning the search at valuations found.
\\ Speed is good for a few thousand bits, though no final
\\ results are known until the end of the whole search.
\\
\\ Requires PARI/GP 2.9.x or higher for "oo" infinity.

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


\\ 24*A000127(n)
A364133_poly = x^4 - 6*x^3 + 23*x^2 - 18*x + 24;

\\ Return a vector of all terms of A364133 < 2^num_bits.
A364133_vector_to_bits(num_bits) =
{
  my(best_n = vector(4*num_bits+3,v, oo),
     inc_limit = num_bits - 1,
     n=0); \\ n = candidate term of A364133
  while(1,
        my(t = subst(A364133_poly,'x,n),
           v = valuation(t,2));
        if(n < best_n[v] && n > 0,
           best_n[v] = n);

        v = min(v, inc_limit);
        while(bittest(n,v), v--);
        if(v < 0, break);
        n = bitand(n, bitneg(0,v)) + 1<<v);

  \\ find records
  my(ret=List([]), m=oo);
  forstep(v=#best_n,1,-1,
     if(best_n[v] < m,
        m=best_n[v];
        listput(ret,m)));
  Vecrev(ret);
}

{
  my(v=A364133_vector_to_bits(20));
  print1("A364133 = ");
  for(i=1,#v, print1(v[i],","));
  print("...");
}


\\-----------------------------------------------------------------------------

\\ Implementation Notes
\\ --------------------
\\
\\ P(n) = 24*A000127(n) = subst(A364133_poly,'x,n) ignores the
\\ /24 part of A000127 since
\\
\\     valuation(24*k) = valuation(k) + 3
\\
\\ which is +3 from factor 8, and no change from factor 3.
\\
\\ n values are considered by a depth-first search through
\\ bit values 0 or 1 at bit positions from low to high.
\\ For example, num_bits = 3 would consider
\\
\\     n = 000        n values considered in
\\         100        order of their bit reversal
\\         010
\\         110
\\         001
\\         101
\\         011
\\         111
\\
\\ Bit positions 0..i of n (position 0 the least significant)
\\ determine the corresponding bit positions 0..i in P(n),
\\
\\ v = valuation(P(n)) means P(n) ends ..100..00 in bit
\\ positions 0..v.  n is a residue class in the sense that
\\ any with the same bits 0..v has the same valuation,
\\
\\     valuation(P(m)) = v
\\     for all m in residue class n mod 2^(v+1)
\\
\\ This means bit positions v+1 and above in n can be skipped,
\\ as any bits there don't change the valuation.
\\
\\ The search order and skipping means n < 2^(v+1), ie. n is
\\ always the smallest of the residue class.  This is since
\\ any reduction "n mod 2^(v+1)" would mean clearing some high
\\ 1-bits of n and that would be earlier in depth-first order.
\\
\\ A given valuation v might occur at various n.  The best,
\\ meaning the smallest, is held
\\
\\     best_n[v] = n smallest having valuation(P(n))=v
\\
\\
\\ At a given v, skip bits above 0..v means increment the bit
\\ at position v in n.  If already a 1-bit then it can't be
\\ increased, so drop to try an increment at v-1, or v-2, etc.
\\ Those drops are backtracking, or back up in the search tree.
\\
\\ The code finds the increment position at the bittest() loop.
\\ At a 0-bit, the increment means set that bit to 1 and clear
\\ to 0s all bits above (bitand() of a mask).
\\
\\ If n has no 0 bits among 0..v then the loop reaches v = -1
\\ (by bittest(n,-1) = 0), and there are no more n to consider.
\\
\\
\\ The final step is to pick record highs from best_n[],
\\
\\     n = best_n[v] is index of a record high valuation
\\     when best_n[v'] > n for all v' > v
\\
\\ The code does this going down from biggest v and each new
\\ smallest best_n[v] is a term.
\\
\\ best_n[] might have holes where a given valuation never
\\ occurred.  Those are "oo" infinity and they have no effect
\\ on the records.  All n < 2^num_bits have been considered
\\ and the valuations occurring among those n are all that
\\ matter.


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