\\ A short PARI-program to compute OEIS-sequence A351029 "Number of integers whose arithmetic derivative is equal to the n-th primorial."
\\ and A366890 "Numbers k whose arithmetic derivative (A003415) is a primorial number (A002110) > 1 and A001222(k) > 2."
\\ First draft, 2024-01-09 by Antti Karttunen
\\ Second editing, 2024-01-16 by Antti Karttunen
\\ Load into PARI with command: read("A369000.gp.txt");
\\ This algorithm is based on the following facts:
\\ A327978 (numbers whose arithmetic derivative is a primorial number > 1)
\\ is a subsequence of A327862 (numbers whose arithmetic derivative is of the form 4k+2),
\\ which in turn is a subsequence of A235992 (numbers with an even arithmetic derivative),
\\ which is an union of multiples of 4 and odd numbers with an even number of prime factors with multiplicity.
\\ Because arithmetic derivative of any multiple of 4 is still multiple of 4, it leaves only the latter choice.
\\ Furthermore, A327978 is a subsequence of A004709 (cubefree numbers: numbers that are not divisible by any cube > 1),
\\ therefore we need to consider only weakly increasing sequences of odd primes (3, 5, 7, 11, ...) of even length,
\\ where any prime is duplicated at most once, but no more.
\\ For increasing prime p, A003415(p*n) = n + p*A003415(n) is a strictly monotonic sequence for any n > 1,
\\ so we can optimize the last prime off from the prime-list and iterate over one element shorter sequences.
\\ As A003415(p*n) = n + p*A003415(n), so p = (A003415(p*n)-n) / A003415(n).
vminsols = vector(101); \\ A368703 o=0: a(n) is the least integer k for which A003415(k) = A002110(n), and 0 if no such k exists.
vmaxsols = vector(101); \\ A368704 o=1: a(n) is the greatest integer k for which A003415(k) = A002110(n), and 0 if no such k exists.
A002110(n) = prod(i=1,n,prime(i));
A003415vec(tv) = { my(n=factorback(tv), s=0, m=1, spf); for(i=1,#tv,spf = tv[i]; n /= spf; s += m*n; m *= spf); (s); }; \\ Compute Arithmetic derivative from the vector of primes.
\\ The following function returns:
\\ 0 if max(pv) > (lim-factorback(pv))/A003415(factorback(pv)) [indicating that it's useless to continue and we should backtrack, i.e., call nextvector, or maybe even stop]
\\ p when lim = p*A003415(pv) + 1*pv, with some p >= max(pv), i.e., when lim-factorback(pv) is a divisible by A003415(factorback(pv)) and the divisor is a prime p.
\\ otherwise returns 1 if max(pv) < (lim-factorback(pv))/A003415(factorback(pv)) [indicating that we should just bump up the rightmost prime]
A003415vrl(pv,lim) = { my(n=factorback(pv), x=lim-n, s=0, m=1, spf, u=n); for(i=1,#pv,spf = pv[i]; u /= spf; s += m*u; m *= spf); if(((x/s)>1)));
nextvector(tv,ip) = { my(nv=tv); nv[ip]=nextprime(1+nv[ip]); for(i=1+ip,#nv,if(i>2 && nv[i-1]==nv[i-2], nv[i]=nextprime(1+nv[i-1]), nv[i]=nv[i-1])); (nv); };
\\ nextvector(tv,ip) = { my(nv=tv); nv[ip]=nextprime(1+nv[ip]); for(i=1+ip,#nv,if(i>2 && nv[i-1]==nv[i-2], nv[i]=nextprime(1+nv[i-1]), nv[i]=nv[i-1])); print(nv); (nv); };
\\ For vl=3 and size=5 (A002210(5) = 2310) we iterate over the following 24 prime-vectors:
\\ tv=[3, 3, 5] ip=3
\\ tv=[3, 3, 7] ip=3
\\ tv=[3, 3, 11] ip=3
\\ tv=[3, 3, 13] ip=3
\\ tv=[3, 3, 17] ip=3
\\ tv=[3, 3, 19] ip=3
\\ tv=[3, 5, 5] ip=2
\\ tv=[3, 5, 7] ip=3
\\ tv=[3, 5, 11] ip=3
\\ tv=[3, 5, 13] ip=3
\\ tv=[3, 5, 17] ip=3
\\ tv=[3, 7, 7] ip=2
\\ tv=[3, 7, 11] ip=3
\\ tv=[3, 7, 13] ip=3
\\ tv=[3, 7, 17] ip=3
\\ tv=[3, 11, 11] ip=2
\\ tv=[5, 5, 7] ip=1
\\ tv=[5, 5, 11] ip=3
\\ tv=[5, 5, 13] ip=3
\\ tv=[5, 7, 7] ip=2
\\ tv=[5, 7, 11] ip=3
\\ tv=[5, 7, 13] ip=3
\\ tv=[5, 11, 11] ip=2
\\ tv=[7, 7, 11] ip=1
loop_over_prime_vectors_of_length(vl,size,deb) = {
my(f, u = A002110(size), ip=vl, c=0, totits=0, previts=0, tv = firstvector(vl), sol);
if(2==deb,print("tv=",tv, " ip=",ip));
while(1,totits++;
f = A003415vrl(tv,u);
if(f, \\ Still makes sense to keep on incrementing the last prime in tv?
ip = vl; \\ So set or keep the ip pointing at the rightmost prime.
if(f>1, c++; \\ Found a solution (f is then our missing prime).
sol = f*factorback(tv);
if(!vminsols[size] || sol < vminsols[size], vminsols[size] = sol);
if(sol > vmaxsols[size], vmaxsols[size] = sol);
if(vl>2,
write("b366890_terms_found_by_the_search_order.txt",sol);
if(deb, print("Found a ",if(vl>2,"non",""),"semiprime solution: ",concat(tv,[f])," = ",sol))
);
),
/* Otherwise (when A003415vrl returned 0) we should backtrack towards left: */
ip = ip-1;
if(!ip,
if(deb, print("Size=",size,": Loop over vectors of length ",vl," finished. Number of iterations: ",totits,". Number of solutions: ", c));
return(c);
);
);
if(ip==vl, tv[vl] = nextprime(1+tv[vl]), \\ Keep on bumping up the rightmost (largest) prime in the vector (just an optimization)
tv = nextvector(tv,ip) \\ As this would suffice in both cases.
);
if(2==deb,print("tv=",tv, " ip=",ip));
if((3==deb && tv[1]!=tv[2])
||
(4==deb && 1==ip),
print("tv=",tv, " its=", totits-previts); \\ Show the test vector always when we have backtracked to its start.
previts = totits;
)
);
};
A351029_or_A369000(n,vls) = {
my(u=A002110(n), s=0);
forstep(vl=vls,oo,2,
if(A003415vec(firstvector(vl)) > u,
break,
s += loop_over_prime_vectors_of_length(vl,n,4+(vl==1));
)
);
write("b368703.txt", n, " ", vminsols[n]);
write("b368704.txt", n, " ", vmaxsols[n]);
return(s);
};
A351029(n) = A351029_or_A369000(n,1);
A369000(n) = A351029_or_A369000(n,3); \\ Exclude the Goldbachian solutions, start straight from the test-prime vectors of length 3.
for(n=1,11,write("b351029.txt", n, " ", A351029(n)));
print("A368703 = ",concat([2], vminsols));
print("A368704 = ",vmaxsols);
for(n=1,13,write("b369000.txt", n, " ", A369000(n))); \\ Takes days.
return(0);
\\ For vl=3 and size=5 (A002110(6) = 30030) we iterate over the following 217 prime-vectors:
\\ tv=[3, 3, 5] ip=3
\\ tv=[3, 3, 7] ip=3
\\ tv=[3, 3, 11] ip=3
\\ tv=[3, 3, 13] ip=3
\\ tv=[3, 3, 17] ip=3
\\ tv=[3, 3, 19] ip=3
\\ tv=[3, 3, 23] ip=3
\\ tv=[3, 3, 29] ip=3
\\ tv=[3, 3, 31] ip=3
\\ tv=[3, 3, 37] ip=3
\\ tv=[3, 3, 41] ip=3
\\ tv=[3, 3, 43] ip=3
\\ tv=[3, 3, 47] ip=3
\\ tv=[3, 3, 53] ip=3
\\ tv=[3, 3, 59] ip=3
\\ tv=[3, 3, 61] ip=3
\\ tv=[3, 3, 67] ip=3
\\ tv=[3, 3, 71] ip=3
\\ tv=[3, 5, 5] ip=2
\\ tv=[3, 5, 7] ip=3
\\ tv=[3, 5, 11] ip=3
\\ tv=[3, 5, 13] ip=3
\\ tv=[3, 5, 17] ip=3
\\ tv=[3, 5, 19] ip=3
\\ tv=[3, 5, 23] ip=3
\\ tv=[3, 5, 29] ip=3
\\ tv=[3, 5, 31] ip=3
\\ tv=[3, 5, 37] ip=3
\\ tv=[3, 5, 41] ip=3
\\ tv=[3, 5, 43] ip=3
\\ tv=[3, 5, 47] ip=3
\\ tv=[3, 5, 53] ip=3
\\ tv=[3, 5, 59] ip=3
\\ tv=[3, 5, 61] ip=3
\\ tv=[3, 7, 7] ip=2
\\ tv=[3, 7, 11] ip=3
\\ tv=[3, 7, 13] ip=3
\\ tv=[3, 7, 17] ip=3
\\ tv=[3, 7, 19] ip=3
\\ tv=[3, 7, 23] ip=3
\\ tv=[3, 7, 29] ip=3
\\ tv=[3, 7, 31] ip=3
\\ tv=[3, 7, 37] ip=3
\\ tv=[3, 7, 41] ip=3
\\ tv=[3, 7, 43] ip=3
\\ tv=[3, 7, 47] ip=3
\\ tv=[3, 7, 53] ip=3
\\ tv=[3, 11, 11] ip=2
\\ tv=[3, 11, 13] ip=3
\\ tv=[3, 11, 17] ip=3
\\ tv=[3, 11, 19] ip=3
\\ tv=[3, 11, 23] ip=3
\\ tv=[3, 11, 29] ip=3
\\ tv=[3, 11, 31] ip=3
\\ tv=[3, 11, 37] ip=3
\\ tv=[3, 11, 41] ip=3
\\ tv=[3, 11, 43] ip=3
\\ tv=[3, 11, 47] ip=3
\\ tv=[3, 13, 13] ip=2
\\ tv=[3, 13, 17] ip=3
\\ tv=[3, 13, 19] ip=3
\\ tv=[3, 13, 23] ip=3
\\ tv=[3, 13, 29] ip=3
\\ tv=[3, 13, 31] ip=3
\\ tv=[3, 13, 37] ip=3
\\ tv=[3, 13, 41] ip=3
\\ tv=[3, 17, 17] ip=2
\\ tv=[3, 17, 19] ip=3
\\ tv=[3, 17, 23] ip=3
\\ tv=[3, 17, 29] ip=3
\\ tv=[3, 17, 31] ip=3
\\ tv=[3, 17, 37] ip=3
\\ tv=[3, 19, 19] ip=2
\\ tv=[3, 19, 23] ip=3
\\ tv=[3, 19, 29] ip=3
\\ tv=[3, 19, 31] ip=3
\\ tv=[3, 19, 37] ip=3
\\ tv=[3, 23, 23] ip=2
\\ tv=[3, 23, 29] ip=3
\\ tv=[3, 23, 31] ip=3
\\ tv=[3, 23, 37] ip=3
\\ tv=[3, 29, 29] ip=2
\\ tv=[5, 5, 7] ip=1
\\ tv=[5, 5, 11] ip=3
\\ tv=[5, 5, 13] ip=3
\\ tv=[5, 5, 17] ip=3
\\ tv=[5, 5, 19] ip=3
\\ tv=[5, 5, 23] ip=3
\\ tv=[5, 5, 29] ip=3
\\ tv=[5, 5, 31] ip=3
\\ tv=[5, 5, 37] ip=3
\\ tv=[5, 5, 41] ip=3
\\ tv=[5, 5, 43] ip=3
\\ tv=[5, 5, 47] ip=3
\\ tv=[5, 5, 53] ip=3
\\ tv=[5, 7, 7] ip=2
\\ tv=[5, 7, 11] ip=3
\\ tv=[5, 7, 13] ip=3
\\ tv=[5, 7, 17] ip=3
\\ tv=[5, 7, 19] ip=3
\\ tv=[5, 7, 23] ip=3
\\ tv=[5, 7, 29] ip=3
\\ tv=[5, 7, 31] ip=3
\\ tv=[5, 7, 37] ip=3
\\ tv=[5, 7, 41] ip=3
\\ tv=[5, 7, 43] ip=3
\\ tv=[5, 7, 47] ip=3
\\ tv=[5, 7, 53] ip=3
\\ tv=[5, 11, 11] ip=2
\\ tv=[5, 11, 13] ip=3
\\ tv=[5, 11, 17] ip=3
\\ tv=[5, 11, 19] ip=3
\\ tv=[5, 11, 23] ip=3
\\ tv=[5, 11, 29] ip=3
\\ tv=[5, 11, 31] ip=3
\\ tv=[5, 11, 37] ip=3
\\ tv=[5, 11, 41] ip=3
\\ tv=[5, 13, 13] ip=2
\\ tv=[5, 13, 17] ip=3
\\ tv=[5, 13, 19] ip=3
\\ tv=[5, 13, 23] ip=3
\\ tv=[5, 13, 29] ip=3
\\ tv=[5, 13, 31] ip=3
\\ tv=[5, 13, 37] ip=3
\\ tv=[5, 13, 41] ip=3
\\ tv=[5, 17, 17] ip=2
\\ tv=[5, 17, 19] ip=3
\\ tv=[5, 17, 23] ip=3
\\ tv=[5, 17, 29] ip=3
\\ tv=[5, 17, 31] ip=3
\\ tv=[5, 17, 37] ip=3
\\ tv=[5, 19, 19] ip=2
\\ tv=[5, 19, 23] ip=3
\\ tv=[5, 19, 29] ip=3
\\ tv=[5, 19, 31] ip=3
\\ tv=[5, 19, 37] ip=3
\\ tv=[5, 23, 23] ip=2
\\ tv=[5, 23, 29] ip=3
\\ tv=[5, 29, 29] ip=2
\\ tv=[7, 7, 11] ip=1
\\ tv=[7, 7, 13] ip=3
\\ tv=[7, 7, 17] ip=3
\\ tv=[7, 7, 19] ip=3
\\ tv=[7, 7, 23] ip=3
\\ tv=[7, 7, 29] ip=3
\\ tv=[7, 7, 31] ip=3
\\ tv=[7, 7, 37] ip=3
\\ tv=[7, 7, 41] ip=3
\\ tv=[7, 7, 43] ip=3
\\ tv=[7, 11, 11] ip=2
\\ tv=[7, 11, 13] ip=3
\\ tv=[7, 11, 17] ip=3
\\ tv=[7, 11, 19] ip=3
\\ tv=[7, 11, 23] ip=3
\\ tv=[7, 11, 29] ip=3
\\ tv=[7, 11, 31] ip=3
\\ tv=[7, 11, 37] ip=3
\\ tv=[7, 13, 13] ip=2
\\ tv=[7, 13, 17] ip=3
\\ tv=[7, 13, 19] ip=3
\\ tv=[7, 13, 23] ip=3
\\ tv=[7, 13, 29] ip=3
\\ tv=[7, 13, 31] ip=3
\\ tv=[7, 13, 37] ip=3
\\ tv=[7, 17, 17] ip=2
\\ tv=[7, 17, 19] ip=3
\\ tv=[7, 17, 23] ip=3
\\ tv=[7, 17, 29] ip=3
\\ tv=[7, 17, 31] ip=3
\\ tv=[7, 19, 19] ip=2
\\ tv=[7, 19, 23] ip=3
\\ tv=[7, 19, 29] ip=3
\\ tv=[7, 19, 31] ip=3
\\ tv=[7, 23, 23] ip=2
\\ tv=[7, 23, 29] ip=3
\\ tv=[7, 29, 29] ip=2
\\ tv=[11, 11, 13] ip=1
\\ tv=[11, 11, 17] ip=3
\\ tv=[11, 11, 19] ip=3
\\ tv=[11, 11, 23] ip=3
\\ tv=[11, 11, 29] ip=3
\\ tv=[11, 11, 31] ip=3
\\ tv=[11, 11, 37] ip=3
\\ tv=[11, 13, 13] ip=2
\\ tv=[11, 13, 17] ip=3
\\ tv=[11, 13, 19] ip=3
\\ tv=[11, 13, 23] ip=3
\\ tv=[11, 13, 29] ip=3
\\ tv=[11, 13, 31] ip=3
\\ tv=[11, 17, 17] ip=2
\\ tv=[11, 17, 19] ip=3
\\ tv=[11, 17, 23] ip=3
\\ tv=[11, 17, 29] ip=3
\\ tv=[11, 19, 19] ip=2
\\ tv=[11, 19, 23] ip=3
\\ tv=[11, 19, 29] ip=3
\\ tv=[11, 23, 23] ip=2
\\ tv=[11, 23, 29] ip=3
\\ tv=[11, 29, 29] ip=2
\\ tv=[13, 13, 17] ip=1
\\ tv=[13, 13, 19] ip=3
\\ tv=[13, 13, 23] ip=3
\\ tv=[13, 13, 29] ip=3
\\ tv=[13, 17, 17] ip=2
\\ tv=[13, 17, 19] ip=3
\\ tv=[13, 17, 23] ip=3
\\ tv=[13, 17, 29] ip=3
\\ tv=[13, 19, 19] ip=2
\\ tv=[13, 19, 23] ip=3
\\ tv=[13, 19, 29] ip=3
\\ tv=[13, 23, 23] ip=2
\\ tv=[17, 17, 19] ip=1
\\ tv=[17, 17, 23] ip=3
\\ tv=[17, 19, 19] ip=2
\\ tv=[17, 19, 23] ip=3
\\ tv=[17, 23, 23] ip=2
\\ tv=[19, 19, 23] ip=1
\\ For vl=5 and size=5 (A002110(6) = 30030) we iterate over the following 12 prime-vectors:
\\ tv=[3, 3, 5, 5, 7] ip=5
\\ tv=[3, 3, 5, 5, 11] ip=5
\\ tv=[3, 3, 5, 7, 7] ip=4
\\ tv=[3, 3, 5, 7, 11] ip=5
\\ tv=[3, 3, 5, 11, 11] ip=4
\\ tv=[3, 3, 7, 7, 11] ip=3
\\ tv=[3, 5, 5, 7, 7] ip=2
\\ tv=[3, 5, 5, 7, 11] ip=5
\\ tv=[3, 5, 5, 11, 11] ip=4
\\ tv=[3, 5, 7, 7, 11] ip=3
\\ tv=[3, 7, 7, 11, 11] ip=2
\\ tv=[5, 5, 7, 7, 11] ip=1
\\ For vl=7 and size=8 (A002110(8) = 9699690) we iterate over the following 43 prime-vectors:
\\ tv=[3, 3, 5, 5, 7, 7, 11] ip=7
\\ tv=[3, 3, 5, 5, 7, 7, 13] ip=7
\\ tv=[3, 3, 5, 5, 7, 7, 17] ip=7
\\ tv=[3, 3, 5, 5, 7, 7, 19] ip=7
\\ tv=[3, 3, 5, 5, 7, 7, 23] ip=7
\\ tv=[3, 3, 5, 5, 7, 7, 29] ip=7
\\ tv=[3, 3, 5, 5, 7, 11, 11] ip=6
\\ tv=[3, 3, 5, 5, 7, 11, 13] ip=7
\\ tv=[3, 3, 5, 5, 7, 11, 17] ip=7
\\ tv=[3, 3, 5, 5, 7, 11, 19] ip=7
\\ tv=[3, 3, 5, 5, 7, 11, 23] ip=7
\\ tv=[3, 3, 5, 5, 7, 13, 13] ip=6
\\ tv=[3, 3, 5, 5, 7, 13, 17] ip=7
\\ tv=[3, 3, 5, 5, 7, 13, 19] ip=7
\\ tv=[3, 3, 5, 5, 7, 17, 17] ip=6
\\ tv=[3, 3, 5, 5, 11, 11, 13] ip=5
\\ tv=[3, 3, 5, 5, 11, 11, 17] ip=7
\\ tv=[3, 3, 5, 5, 11, 13, 13] ip=6
\\ tv=[3, 3, 5, 5, 11, 13, 17] ip=7
\\ tv=[3, 3, 5, 5, 11, 17, 17] ip=6
\\ tv=[3, 3, 5, 5, 13, 13, 17] ip=5
\\ tv=[3, 3, 5, 7, 7, 11, 11] ip=4
\\ tv=[3, 3, 5, 7, 7, 11, 13] ip=7
\\ tv=[3, 3, 5, 7, 7, 11, 17] ip=7
\\ tv=[3, 3, 5, 7, 7, 11, 19] ip=7
\\ tv=[3, 3, 5, 7, 7, 13, 13] ip=6
\\ tv=[3, 3, 5, 7, 7, 13, 17] ip=7
\\ tv=[3, 3, 5, 7, 7, 17, 17] ip=6
\\ tv=[3, 3, 5, 7, 11, 11, 13] ip=5
\\ tv=[3, 3, 5, 7, 11, 11, 17] ip=7
\\ tv=[3, 3, 5, 7, 11, 13, 13] ip=6
\\ tv=[3, 3, 5, 7, 13, 13, 17] ip=5
\\ tv=[3, 3, 5, 11, 11, 13, 13] ip=4
\\ tv=[3, 3, 7, 7, 11, 11, 13] ip=3
\\ tv=[3, 5, 5, 7, 7, 11, 11] ip=2 A003415(3*5*5*7*7*11*11*13) = 7386610 < A002110(8) = 9699690
\\ tv=[3, 5, 5, 7, 7, 11, 13] ip=7 A003415(3*5*5*7*7*11*13*13) = 8634080 < A002110(8) = 9699690
\\ tv=[3, 5, 5, 7, 7, 11, 17] ip=7
\\ tv=[3, 5, 5, 7, 7, 13, 13] ip=6
\\ tv=[3, 5, 5, 7, 11, 11, 13] ip=5
\\ tv=[3, 5, 5, 11, 11, 13, 13] ip=4
\\ tv=[3, 5, 7, 7, 11, 11, 13] ip=3
\\ tv=[3, 7, 7, 11, 11, 13, 13] ip=2
\\ tv=[5, 5, 7, 7, 11, 11, 13] ip=1 A003415(5*5*7*7*11*11*13*13) = 25585560 > A002110(8) = 9699690