[seqfan] Re: A093893 Subsequence

[Hagen von EItzen]

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