[seqfan] Re: Who understands Granville numbers?
Maximilian Hasler
maximilian.hasler at gmail.com
Thu Oct 28 18:01:14 CEST 2010
> The comment on A118372 says:
> S={1}. Assume n>1 and that all numbers m<n have already been tested.
> Let s=sum{d: d|n, d<n and d in S}. If s<=n then n is now in S.
>
> I take that to mean:
> S_0 = {1}.
> s_n = sum{ d: d | n & d < n & d in S_{n-1} }
> S_n := { S_{n-1} s_n > n
> { S_{n-1} U {n} s_n <= n
>
>
I agree that the comment cannot be interpreted otherwise.
Anyway all divisors of n are <= n so the result of
this test (for n, and for d | n, d in S ?)
does not change subsequently if larger n's are added.
A118372(Nmax) = { S=Set(1); for(n=2,Nmax, my(s=sumdiv(n,d,if( d<n &
setsearch(S,d), d))); s > n & next; S=setunion(S,Set(n)); n==s &
print1(n","); )}
The sequence of numbers NOT in S after this procedure,
12,18,20,30,42,48,56,66,70,72,78,80,84,88,90,102,104,108,114,120,138,150,
162,174,180,186,192,196,200,210,220,222,246,252,258,260,264,270,272,280,
282,288,294,300,304,308,312,318,320,330,336,340,354,364,366,368,378,380,
390,396,402,408,420,426,432,438,448,450,456,460,462,464,468,474,476,480,
498,500,504,510,532,534,546,550,552,570,572,580,582,594,606,612,618,620,
630,642,644,648,650,654,660,672,678...
seems not in OEIS, either.
It is much more efficient to use & store this list rather than S itself,
(60ms vs 500ms, resp. 1.5s vs 40s for N=10^3 resp. 10^4 on my laptop)
i.e.
A118372(Nmax) = { C=[]; for(n=2,Nmax, my(s=sumdiv(n,d,if( d<n &
!setsearch(C,d), d))); s > n & C=setunion(C,Set(n)); n==s & print1(n","); )}
Maximilian
More information about the SeqFan
mailing list