7
$\begingroup$

A weird number is a number such that the sum of the proper divisors (divisors including 1 but not itself) of the number is greater than the number, but no subset of these divisors sums to to number itself. Find all the weird numbers up to 10000.

I wrote the following code

Select[Range[100],( ((Plus @@ Most[Divisors[#]])>#) && 
  (Times @@ (DeleteDuplicates[Plus @@@ Subsets[Most[Divisors[#]]]]-#))!=0)&]

Although it works fine in the range (0,100), the program crashed when it turns to (0,10000), after several minutes of computation. I know this direct method isn't efficient, but how to optimize the code without prior knowledge of properties of weird numbers?

$\endgroup$
  • $\begingroup$ See also: groups.google.com/forum/#!topic/comp.soft-sys.math.mathematica/… $\endgroup$ – Qianyi Guo Feb 5 '15 at 5:52
  • 3
    $\begingroup$ The problem is that Subsets[Most[Divisors[720]]] eats up your RAM. $\endgroup$ – Karsten 7. Feb 5 '15 at 6:46
  • 1
    $\begingroup$ Does "how to optimize the code without prior knowledge of properties of weird numbers" mean, that you'll accept solely a brute force solution? $\endgroup$ – Karsten 7. Feb 5 '15 at 7:00
  • 1
    $\begingroup$ You can avoid Subsets if you set up the check as a subset sum problem. FindInstance for example could be used for this. $\endgroup$ – Daniel Lichtblau Feb 5 '15 at 17:30
12
$\begingroup$

This is quite good:

f[d_, n_] := With[{x = Pick[d, Thread[Accumulate[d] - n < 0]]},
  Scan[f[x, n - Total[#]] &, Subsets[Complement[d, x], {1, Infinity}]]]

f[_, 0] := Throw @ True;

weird[n_] := (DivisorSigma[1, n] > 2 n) && ! TrueQ @ Catch @ f[Most @ Divisors[n], n]

Select[Range[10000], weird] // AbsoluteTiming
(* {1.450002, {70, 836, 4030, 5830, 7192, 7912, 9272}} *)

The basic idea is that, for example, to make a total of $70$ from its divisors $\{1, 2, 5, 7, 10, 14, 35\}$, the $35$ must be included because the total of all the others is insufficient. The problem is then reduced to making a total of $70 - 35 = 35$ from $\{1, 2, 5, 7, 10, 14\}$. The logic can be applied recursively until you reach a target of $0$ (in which case the number is non-weird) or run out of divisors (in which case it is weird).

$\endgroup$
  • $\begingroup$ Beautiful. +1 .. $\endgroup$ – Mr.Wizard Feb 6 '15 at 4:01
8
$\begingroup$
(*  Input: Range of even numbers --- Output: Primitive weird numbers *)
Block[{$RecursionLimit=Infinity},
       subOfSum[ss_, kk_, rr_]:= Module[{s=ss, k=kk, r=rr},
        If[ s+w[[k]] >=mm && s +w[[k]] <=m,  t=False;  Goto[ done](*Found*),
            If[s +w[[k]]+w[[k +1]] <=m, subOfSum[s +w[[k]], k+1, r-w[[k]]]];
        If[s+r -w[[k]] >= m && s+w[[k+1]] <= m, subOfSum[s,k+1, r-w[[k]]]]]; t](* end subOfSum *);
        greedyQ[ab_]:= Module[{ abn=ab, v, sum, s, j, jj, k}, tt=False;
        jj = Length[w]; (*start search*)
        Do[s = r; sum = 0; Do[ v=w[[j]]; sum = sum + v;
        If[sum > abn, sum = sum -v; Goto[nxt]];
        If[sum == abn, tt = True; Goto[ doneG]]; s=s-v; Label[nxt], {j, jj, 1, -1}];
           jj= jj-1, {k, 1, jj-1}] ; Label[doneG];
        (* True means found, False not found *) tt] (*end  greedyQ*);
        cnt=0; Do[ If[Mod[n,3]==0, Goto[agn]]; r=DivisorSigma[1, n]; m=r- 2*n;
        If[m > 0, fi = FactorInteger[n]; largestP = fi[[Length[fi]]][[1]];
        nn = n/largestP; If[m > 2*nn || Length[fi] < 3,  Goto[agn]];
        If[DivisorSigma[1, nn] > 2*nn, Goto[agn]]; t = True; r = r-n;
        ww = Divisors[n]; lenW = Length[ww];
        Do[If[ww[[i]] <= m, w = Drop[ww, i- lenW]; Break[], r = r-ww[[i]]],{i, lenW -1, 1, -1}];
        If[r >= m, If[greedyQ[ m], t = False, (*Powers of 2 dropped*)
            exp2 = fi[[1]][[2]]; sig2= 2^(exp2+1) -1; mm = m - sig2;
            lenW= Length[w]; ww={};
            If[exp2 > 1, Do[Do[If[w[[i]] == 2^ii, ww = AppendTo[ww, w[[i]]]], {i, 1, lenW}],
                { ii, 0, exp2}]; w= Complement[w, ww]
          (*end T if*),
              w= Drop[w, 2]];
          (*end Pwr2*)
        t =subOfSum[0, 1, r]]]; Label[done];
        If[t, Print[++cnt, "   ", n, "  ", t]]];
        Label[agn], {n, 2, 1000000, 2}]]

The above is an algorithm I wrote for finding primitive even weird numbers. The code is long and I'm sure, the coding is not very well written, but it is fast. The range of numbers to check for is at the bottom where I have {n, 2, 1000000, 2}. This is for finding all even weird numbers in range of 2 to 10^6. Any range can be put in and is only limited by Mathematica to factor the numbers.

Brent

How the algorithm works. It calculates primitive weird numbers which are always primitive abundant numbers. These are weird numbers that are not multiples of other weird numbers.

The algorithm eliminates numbers which can't be primitive weird until a finial recursive search is done to check if a sum of divisors can make up the abundance of the number. The abundance of the number is m = Sigma[n]-2*n. After several simple checks, the number is checked to see if it's not a primitive abundant number. For all primitive abundant numbers, Sigma[n]/n <= 2*(1+1/p), where p is the largest prime of the number n. If m > (2*n)/p, it can't be primitive abundant. Also, n/p needs to be a deficient number to be primitive weird.

All the divisors of n are eliminated that are greater than m. Let r equal to the sum of all the remaining divisors. If r < m, then n is weird and done. If not, do a quick search using a greedy type algorithm trying to form a sum of the remaining divisors to equal the abundance m. If this fails, the number may be weird.

Since n is even, remove all the divisors of Sigma[2^x] and put the reminding divisors into a list w . Note, the divisors of a number 2^x can be summed to any number from 1 to Sigma[2^x]. Let mm equal to the abundance minus Sigma[2^x]. If a sum can be made of the list w that fall in the range between m to mm, then the number is not weird.

Let r equal to the sum of divisors in w before the powers of 2 are removed. Call the recursive algorithm subOfSum[ 0, 1, r] to try to find a subset of divisors w to be between mm and m. If none can be found, then n is primitive weird.

Sorry my code is so bad. I hope the description can be followed. If you have any questions on anything, let me know.

Brent

$\endgroup$
  • $\begingroup$ Hi, Brent ! Only one suggestion - please format the code according to the guidelines @ help centre ;__; $\endgroup$ – Sektor Feb 6 '15 at 21:36
  • 7
    $\begingroup$ There should be an award for this answer. It is the most deplorably bad Mathematica code I have seen in a long time. And the most blindingly fast for the problem at hand. That latter in my mind makes it worthy of being the accepted answer (not that I object to the one that was actually chosen, just giving an opinion). But the ratio of speed to code awfulness is something almost surreal. Now I just need to spend a bit of time trying to understand how it works so well. $\endgroup$ – Daniel Lichtblau Feb 7 '15 at 22:29
  • 2
    $\begingroup$ PS I mean no offense. The speed to my mind trumps the style (although the latter would have a considerable effect on maintainability, were this production code). I was simply amazed by that disparity. $\endgroup$ – Daniel Lichtblau Feb 7 '15 at 22:31
  • $\begingroup$ @Daniel Those are some of the funniest comment's I've ever read. A brilliant algorithm hacked out any way possible perhaps? Brent, would you please consider providing a written description of the algorithm you use here? $\endgroup$ – Mr.Wizard Feb 9 '15 at 5:52
  • $\begingroup$ Can someone supply me the email of Brent Baughn so that I may continue this line of programming. Bob Wilson rgwv@rgwv.com $\endgroup$ – Robert G. Wilson v Jun 11 '15 at 3:04
5
$\begingroup$

Copied from OEIS without thinking about it:

Needs["Combinatorica`"] (*then*) 
fQ[n_] := Block[{d, l, t, i},
   If[DivisorSigma[1, n] > 2 n && Mod[n, 6] != 0,
    d = Take[Divisors[n], {1, -2}]; l = 2^Length[d]; 
    t = Table[NthSubset[j, d], {j, l - 1}];
    i = 2;
    While[i < l && Plus @@ t[[i]] != n, i++]
    ];
   If[i == l, True, False]];
Select[Range[100], fQ[#] &] 
(*from Robert G.Wilson v,May 20 2005*)

(* {70} *)
$\endgroup$
  • 2
    $\begingroup$ Why do people (not you) keep using e.g. fQ[#] &? I've created this redundancy myself a few times without thinking but it is disturbingly prevalent. $\endgroup$ – Mr.Wizard Feb 5 '15 at 14:47
  • 2
    $\begingroup$ @Mr.Wizard Been there, done that. Don't know why. $\endgroup$ – Dr. belisarius Feb 5 '15 at 14:55
  • 1
    $\begingroup$ I don't know if you intended that to be funny but once again you make me laugh. (I hope you take this as kindly as I mean it.) $\endgroup$ – Mr.Wizard Feb 5 '15 at 14:55
  • 1
    $\begingroup$ @Mr.Wizard :D :D $\endgroup$ – Dr. belisarius Feb 5 '15 at 14:56
  • $\begingroup$ The program still uses a huge amount of memory for large numbers. My concern is that this problem comes from an introductory book on Mathematica, and almost all problems can be evaluated in several seconds. Are there more efficient ways to do this? $\endgroup$ – Qianyi Guo Feb 5 '15 at 15:26
5
$\begingroup$

Faster however:

Timing[n = 1;
 Reap[While[n < 10000, d = Most@Divisors@n;
    a = Array[q, Length@d];
    If[Mod[n, 6] != 0 && Tr@d > n && 
      Length[FindInstance[a.d == n && (And @@ Thread[0 <= a <= 1]), a,
          Integers]] == 0,
     Sow[n]];
    n = n + 1];]
 ]

(* {14.188934, {Null, {{70, 836, 4030, 5830, 7192, 7912, 9272}}}} *)

If I use the fact that the smallest odd weird has to be fairly large (larger than 791,000 at least) it can be made a bit faster.

$\endgroup$
4
$\begingroup$

The following is not too slow and uses moderate resources:

n = 2; res = {};
Dynamic@n
While[n < 10000,
  d = Most@Divisors@n;
  a = Array[q, Length@d];
  AppendTo[res, Tr@d > n && Mod[n, 6] != 0 && ! SatisfiableQ[Boole@a.d == n, a]]; 
  n = n + 2];
2 Position[res, True]

{ {70}, {836}, {4030}, {5830}, {7192}, {7912}, {9272}}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.