\\ Kevin Ryde, November 2022
\\
\\ Usage: gp <a358650.gp
\\
\\ This is some PARI/GP code to calculate terms of A358650, the
\\ Matula-Goebel number of the binomial tree of n vertices.

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


a(n) =
{
  my(max_k=logint(n,2),
     k=valuation(n,2),  \\ low 0-bits of n
     T=1);  \\ to become a(2^k)
  for(i=1,k, T*=prime(T));
  if(k==max_k, return(T));  \\ when n = 2^max_k

  my(pT=prime(T),
     pR=pT);  \\ pR = prime(a(n mod 2^k)), after k++
  while(k++ < max_k,
        T*=pT; pT=prime(T);
        if(bittest(n,k), pR=prime(T*pR)));
  T*pT*pR;   \\ most significant 1-bit of n
}

{
  print1("A358650 = ");
  for(n=1,12, print1(a(n),", "));
  print("...");
}

\\ This code uses the recurrences from the sequence
\\
\\     a(n) = a(2^k) * prime(a(r))
\\     where n = 2^k + r and 1 <= r <= 2^k
\\
\\     a(2^k) = a(2^(k-1)) * prime(a(2^(k-1)))
\\            = A033286(a(2^(k-1))
\\            = A076146(k+1)
\\     starting a(1) = 1
\\
\\ T = a(2^k) grows by the a(2^k) formula.
\\ If n is a power of 2 then T is the answer immediately.
\\ Otherwise T is applied as per the first formula to build up to n by
\\ 1-bits of n from least to most significant.
\\
\\ An implicit quantity R = a(n mod 2^k) is the bits accumulated so far.
\\ This is kept in the form pR = prime(R) since that and pT = prime(T) are
\\ ready for the final result which is
\\
\\     a(n) = a(2^(k-1)) * prime(a(2^(k-1))) * prime(a(r))
\\
\\ The reason for calculating pT = prime(T) ahead of its use is so only one
\\ prime(x) is called at the initial pR = pT at the least significant 1-bit.
\\
\\ The speed of PARI prime(x) is good within the "primelimit" but becomes
\\ slow much beyond that.  It's worth some care to call prime(x) only once
\\ for each necessary x.
\\
\\
\\ Limitations
\\ -----------
\\
\\ Speed is limited by prime(x) on ever larger x, and generally requires a a
\\ suitably large "primelimit" setting.  Circa PARI 2.13, a 32-bit build can
\\ only have primelimit <= about 4 billion.  This suffices (and in fact the
\\ default suffices too) up to n=59, but beyond that is unfeasably slow.


\\-----------------------------------------------------------------------------
\\ Vector of Terms
\\
\\ Like most recurrenceas, a vector of terms can be made by using the
\\ recurrence on various previously calculated terms.
\\
\\ In the following a_vector(), each prime(a(n)) is calculated if/when
\\ needed and then kept in pr[] for possible re-use.  The speed
\\ considerations noted above for large prime(x) still apply.

a_vector(len) =
{
  my(v=vector(len,i,1),   \\ v[1] = 1
     pr=vector(len),      \\ pr[n] = prime(v[n])
     k=0,r=1,
     n=1);  \\ v[n] was the last element stored
  while(n<len,
        if(!pr[r], pr[r] = prime(v[r]));
        v[n++] = v[2^k] * pr[r];
        if(r++ > 2^k, k++;r=1));
  v;
}

{
  \\ check a_vector() against a()
  my(v=a_vector(32));
  #v == 32                 || error();
  vector(#v,n, a(n)) == v  || error();
}

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