\\ Kevin Ryde, July 2022 \\ \\ Usage: gp <a064002.gp \\ \\ This is some PARI/GP code to calculate terms of A064002. \\ a(n) is an individual term by descents, and suits large or small n. \\ a_vector(len) is initial terms using just additions. default(recover,0); default(strictargs,1); \\----------------------------------------------------------------------------- \\ Individual Term { \\ return a 2-element vector [i,j], for n>=2 my(pair(n) = my(r, s=sqrtint((n-1)<<1,&r)); if(r>s, r-=s;s++, r+=s); [r>>1, s]); a(n) = if(n<=1, n, my(i,j); [i,j]=pair(n); self()(i) + self()(j) + 1); } { print1("A064002 = "); for(n=1,12, print1(a(n),", ")); print("..."); } \\ a(n) is a recursive implementation of the defining formula \\ \\ a(n) = a(i) + a(j) + 1 \\ \\ pair(n) returns the relevant pair i,j for use at n \\ \\ i = A002260(n-1) \\ j = A002024(n-1) \\ \\ calculated together by a sqrtint() and its "remainder". \\ \\ Recursion is easy, and it measures faster than an explicit List() of \\ pending descents. \\ \\ Peak recursion depth is modest since the smallest n requiring depth k is \\ A108225(k) and that becomes very large very quickly. Eg. depth 30 is a \\ roughly 10 megabyte number. \\----------------------------------------------------------------------------- \\ Vector of Terms \\ A vector of initial terms from the sequence can follow the defining \\ formula above, but iterating through i,j pairs \\ \\ all 1 <= i <= j \\ in order ascending j and within j ascending i \\ \\ The first pair i=1,j=1 is for n=2. a_vector(len) = { my(v=vector(len,i,1), i=1, j=1); \\ 1 <= j <= i for(n=2,len, v[n] = v[i] + v[j] + 1; if(j++ > i, i++;j=1)); v; } { \\ check a_vector() against a() my(v=a_vector(500)); #v == 500 || error(); vector(#v,n, a(n)) == v || error(); } \\----------------------------------------------------------------------------- print("end");