# [seqfan] Re: A093893 Subsequence

Hagen von EItzen math at von-eitzen.de
Sat May 30 00:14:00 CEST 2009

Bob wrote
> Dear Hagen,
>
> 	No, a(11) does not equal to 31^10. The divisors of 31^10 = 819628286980801
> are {1, 31, 961, 29791, 923521, 28629151, 887503681, 27512614111, 852891037441,
> 26439622160671, 819628286980801}. The 1486th subset {1, 31, 961, 29791, 923521,
> 28629151, 887503681} sums to 917087137, a prime. a(11)= 211^10.
>
> 	What I need is a value for a(10).
>
> Bob.
>
Oops, indeed my handcrafted attempts to optimize the search left out a
few subsets for n=11.
Meanwhile I tried the following definitions in PARI:

goodN(n) = goodv(divisors(n))
goodpp(p, n) = goodv(vector(n,k,p^(k-1)))  /* equivalent to
goodN(p^(n-1)) */
goodv(M) = local(i, ok=1, T, J=vector(#M,x,1));
for(i=3,2^#M-1,if(bitand(i,i-1),if(isprime(T=vecextract(M,i)*vecextract(J,i)~),ok=0;break)));ok
/* By the way: I'm still quite new to PARI. Is there no smoother way to
sum all elements of a vector? Something like plus@ in Mathematica? */

and found

//a(2) = 3, //////a(3) = //////49, //////a(4) = //////87, //////a(5) = //////130321, //////a(6) = //////4753, //////a(7) = //////139^6 = 7212549413161, //////a(8) = //////285541,
//////a(9) = //////211^2 * 421^2 = 7890946561,
//////a(10) = //////211^4 * 421 = 834472284661,
//////a(11) = //////211^10 = 174913992535407978606601,
//////a(12) <= 2311^11////////////
a(13) = //////2311^12 = 23205949656945057666311162427422570380321,
//////a(14) <= //////2311^13,
//////a(15) <= 92401^14 ~ 3.3*10^69; if of the form////// p^4*q^2, then p>10^4 or q>10^4.
//////a(16) = //////???, //////
a(17) //////> 10^80 (i.e. p^16 with p>10^5)//

My search for a(10) in PARI was as follows:
First look for numbers of the form p^9:
n=10;p=3;while(!goodpp(p,n),p=nextprime(p+1));lowest = p^(n-1)
Then find a lot of p such that no partial divisor sum of p^4 is prime:
g4=[];forprime(p=3,10000,if(goodpp(p,5),g4=concat(g4,p)))
Finally, try numbers of the form p^4*q with p taken from above list
(starting with big p to avoid need to test for large q beyond primelimit):

for(i=1,#g4,p=g4[#g4-i+1];p4=p^4;qlim=floor(lowest/p4);forprime(q=3,qlim,if(q!=p
&& goodN(p4*q), lowest=p4*q;print(p,"^4*",q," = ",lowest);break)))
This lowers the record in several steps down to
211^4*421////=834472284661//// and then gets stuck because I exceed the
primelimit.
However, it was fascinating that a whole lot of new records were of the
form p^4*211, i.e. the 211 from a(11) is quite ubiquitious.

For the 16 primes p<211, I tried a different approach and generated all
numbers p^4*q in ascending order, stopping when a good number was found.
In fact none was found.

The lack of findings esp. for a(17) may look like bad news for my simple
heuristical argument that a(n) is defined for (almost?) all n.
However, I only tested up to p=10^5. Even for p~10^100, the probability
for numbers n~p^16 being prime is still ~0.00027 and the probability of
no prime appearing
among 2^15 (naively assumed independent) partial sums of that size
(ignoring those trivially divisible by 2 or p) is still about
(1-0.00027)^(2^15) ~ 0.00014,
in other words we might have a(17) = p^16 for a *very* big p.

Hagen