\\ Kevin Ryde, June 2023 \\ This is some PARI/GP code to calculate terms of A363369 \\ which is the number of steps x -> x+1 or x/prime required \\ to go from n to 1. \\ \\ A363369(n) individual term \\ A363369_vector(len) vector of initial terms \\ \\ A vector of terms is fast when all those terms are wanted. \\ An individual term can go to much larger n than a vector, \\ up to sizes factor(n) and various factor(x+1) can handle. default(recover,0); default(strictargs,1); A363369(n) = { if(n<=2,return(n-1)); my(xlist = List([ n ]), \\ x values flist = List([ factor(n) ]), \\ corresponding factor(x) ret=2); \\ result if loop finds somewhere 2 steps to 1 if(vecsum(flist[1][,2])==1, return(1)); \\ n prime while(1, \\ 2 prime factors is 2 steps to reach 1 for(i=1,#flist, if(vecsum(flist[i][,2])==2, return(ret))); my(xnew=List([]), fnew=List([]), pos); \\ divisions /prime for(i=1,#xlist, for(j=1,#flist[i][,2], \\ each prime factor f[j,1] my(f=flist[i]); f[j,2] || next; \\ must exponent remaining my(x = xlist[i]/f[j,1]); \\ new element my(pos=setsearch(xnew,x,1)); if(pos, \\ if not already in xlist[] f[j,2]--; listinsert(xnew,x,pos); listinsert(fnew,f,pos)))); \\ additions +1, \\ and if any is prime then that's 2 steps to reach 1 for(i=1,#xlist, my(x=xlist[i]+1); \\ new term my(pos=setsearch(xnew,x,1)); if(pos, \\ if not already in xlist[] my(f=factor(x)); if(vecsum(f[,2])==1, return(ret)); \\ prime listinsert(xnew,x,pos); listinsert(fnew,f,pos))); xlist=xnew; flist=fnew; ret++); } A363369_vector(len) = { my(a=vector(if(len,nextprime(len)), n, n>1), lo=3,hi=5); while(lo < len, a[hi-1] = 2; \\ prime-1 forstep(n=hi-2, lo+1, -1, my(f=factor(n), m=a[n+1]); \\ by n+1 for(i=1,#f[,1], \\ by each n/prime m = min(m, a[n/f[i,1]])); a[n] = m+1); lo=hi; hi=nextprime(hi+1)); Vec(a,len); } { print1("A363369 = "); for(n=1,20, print1(A363369(n),", ")); print("..."); } \\ consistency of A363369() and A363369_vector() A363369_vector(1000) == vector(1000,n, A363369(n)) || error(); \\----------------------------------------------------------------------------- \\ A363369() Implementation Notes \\ ------------------------------ \\ \\ A363369() starts from a vector [n] and applies all descents \\ to all elements of that vector, \\ \\ x -> x+1 and all x/p \\ \\ xlist[] is the current list of x values (no duplicates). \\ flist[] is their factorisations. \\ \\ flist[i] = factor(xlist[i]) \\ but may have some 0 exponents left in flist[i][j,2] \\ \\ xlist[] has no primes, so everything is at least 2 steps \\ to reach 1. \\ \\ An xlist[] element with exactly 2 prime factors (possibly \\ 2 copies of the same) is 2 steps to reach 1, and is the \\ minimum possible. \\ \\ "ret" has these 2 steps counted already, so is the return \\ value when a 2 steps to reach 1 is found. \\ \\ Descents x -> x/prime are as easy as dividing each prime \\ factor using the factor() matrix in flist[], and \\ decrementing the exponent in that matrix. This might \\ create many new terms (in xnew[]) but is little work. \\ \\ Descent x -> x+1 requires a factor(x+1) to make the new \\ fnew[] entry and to identify possible x+1 prime. If x+1 is \\ prime then that's 2 steps from x to 1 and so return "ret". \\ \\ x+1 steps are considered after x/prime since maybe one of \\ them is an xnew[] entry already, so no need to factor(x+1). \\ \\ x+1 steps are considered in ascending order since smaller \\ x+1 is likely faster to factor() and if it reveals x+1 is \\ prime then never need to factor() larger x+1. \\ \\ In x/prime steps, no attention is paid to the order primes \\ are divided out. This means steps may divide p then q and \\ also q then p. The xnew[] duplicate detection discards the \\ latter. Don't think any sophisticated enumeration is \\ needed, since time for various factor() will dominate any \\ sizeable n. \\ A363369_vector() Implementation Notes \\ ------------------------------------- \\ \\ Vector terms a[n] = a(n) are found firstly by \\ \\ a[1] = 0 \\ a[prime] = 1 \\ a[prime-1] = 2 for prime >= 5 \\ \\ and other terms by \\ \\ a[n] = 1 + min a[x] \\ x = each n+1 and n/p, p prime factor of n \\ \\ n is run in downward ranges between primes so that a[n+1] \\ is available, \\ \\ prime prime \\ "lo" lo+1 ... hi-2 hi-1 "hi" \\ | <--- n --- | \\ \\ Bertrand's postulate shows there is at least one prime \\ between n and n/2. This means n/2 <= lo and so all terms \\ a[n/prime] are already set. \\ forfactored() Not Used \\ ---------------------- \\ \\ A363369_vector() could use forfactored() instead of \\ individual factor(n). But that measures a little slower, \\ presumably mostly due to interpreter overheads. \\ \\ forfactored() goes by ascending n so a[n] can only be set \\ first using its a[n/prime] descents. a[n+1] is brought in \\ by having the store of a[n+1] go back to revise a[n]. \\ If a[n] > a[n+1] + 1 then set a[n]=a[n+1]+1 is the smallest. \\ And if a[n] changes then a[n-1] might be improved by a[n]+1, \\ and so on back. \\ \\ The divide steps a[n/prime] used by a given n are never \\ revised this way, again due to Bertrand's postulate \\ described above. There's always some prime s between n/2 \\ and n and it has a[s]=1 which is never improved so revisions \\ stop there. \\----------------------------------------------------------------------------- print("end");