greedy computation of pseudoperfect numbers?

franktaw at netscape.net franktaw at netscape.net
Wed Apr 9 00:06:23 CEST 2008


For PARI, you get an advantage by taking the product of 1 + x^d[i] 
+x*O(x^n);
this keeps PARI from wasting time computing sums larger than n.  I 
don't know
whether you can do the same thing in Mathematica.

Incidently, when you look at what is going on behind the scenes, this 
is really
a brute force approach to computing all sums of numbers in d.  A greedy
approach with backtracking is theoretically faster.  (You want to 
actually use
an approach called "brute force and ignorance" - keep track of the sum 
of
the remaining smaller values, and regard the current partial sum as a 
failure
when this is too small.)  Given the built-in support for polynomial 
arithmetic in
any of these products, you may need to have quite highly divisible 
numbers
before it's actually faster.

Franklin T. Adams-Watters

P.S. Regarding comprehensibility of code: though I use PARI, I 
generally find
Mathematica to be quite readable.  Maple, on the other hand, tends to 
be
almost completely opaque.

-----Original Message-----
From: T. D. Noe <noe at sspectra.com>

At 8:45 AM -0400 4/8/08, Maximilian Hasler wrote:
>> For both sorts of pseudoperfect numbers, allowing divisor one 
A005835,
>>  and forbidding divisor one A136446, a greedy method _seems_ to 
always
>>  work for deciding whther a given n in in the sequence:

Let d be the list of divisors of n (we can include or exclude 1).  
Compute
terms of the polynomial product_i (1+x^d[i]) up to the x^n term.  
(Perhaps
your PARI code does this. I don't read PARI.) If the x^n term is 
nonzero,
then n is the sum of a subset of the divisors in d.  This is fairly fast
and deterministic.  Using this method, the first 10 odd terms of 
A136446 are

945, 1575, 2205, 2835, 3465, 4095, 4725, 5355, 5775, 5985

which starts the same as A005231, odd abundant numbers.  Here is the
Mathematica code:

t = {}; n = 1; While[Length[t] < 10, n = n + 2;
 d = Rest[Most[Divisors[n]]];
 c = SeriesCoefficient[
   Series[Product[1 + x^d[[i]], {i, Length[d]}], {x, 0, n}], n];
 If[c > 0, AppendTo[t, n]]]; t

Best regards,

Tony










More information about the SeqFan mailing list