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