possible error in A005835

Maximilian Hasler Maximilian.Hasler at martinique.univ-ag.fr
Sun Apr 6 17:31:43 CEST 2008


> > Should the defn of A005835 be
> > Pseudoperfect (or semiperfect) numbers: some set of _distinct_ proper
> > divisors of n sums to n (or: not all divisors are allowed to be equal?)
>  The fact that it's a >set< of divisors implies distinctness.
>  If duplicates were permissible, multiset would be used in place of set.

In fact, using "subset" instead of "set" would create less ambiguities IMHO.

Based on my code for A006037 (of which I'm no more sure why it never
loops infinitely... IMO a vecextract(,"^-1") is missing in the last
instruction, but visibly this is never reached for abundant
numbers....)
I suggest the following PARI code (which in turn should be OK - I put
lots of comments which could of course be deleted if considered too
verbose for a %o, reducing it to about 2-3 lines)

It prints the values up to 1000 in < 50 msec on my laptop, while
checking blindly all subsets takes quite long.

%o A005835 (PARI)
isA005835(n, d=0)={ local(t);
 /* Return nonzero iff n is the sum of a subset of d which defaults to
the set of proper divisors of n */
 if( !d, /* Initialize d */ d=vecextract(divisors(n), "^-1")
 ,/*else check if n equals to one element of d */ setsearch( Set(d), n
) & return(1));
 /* Remove terms > n */ while( #d>1 & d[ #d]>n, d=vecextract(d, "^-1"));
 /* If n is not smaller than the sum of all terms, we're done */ n >=
(t = sum(i=1, #d, d[i])) & return( n==t );
 /* If n is larger than M=max(d), then try to write n-M in terms of d \ { M } */
 n > d[#d] & isA005835( n - d[#d], vecextract( d, "^-1") ) & return(1);
 /* else only d \ {M} is needed */ isA005835( n, vecextract( d, "^-1" ))}

%o A005835 (PARI) for(n=1,1000,isA005835(n)&print1(n","))

I get the same a(1)..a(1000) as T.D.Noe in 0.3 sec.
Maximilian





More information about the SeqFan mailing list