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?

Improve this question
$\endgroup$
4
  • $\begingroup$ See also: groups.google.com/forum/#!topic/comp.soft-sys.math.mathematica/… $\endgroup$
    – Qianyi Guo
    Feb 5, 2015 at 5:52
  • 3
    $\begingroup$ The problem is that Subsets[Most[Divisors[720]]] eats up your RAM. $\endgroup$
    – Karsten7
    Feb 5, 2015 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$
    – Karsten7
    Feb 5, 2015 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, 2015 at 17:30

5 Answers 5

Reset to default
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).

Improve this answer
$\endgroup$
1
  • $\begingroup$ Beautiful. +1 .. $\endgroup$
    – Mr.Wizard
    Feb 6, 2015 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

Improve this answer
$\endgroup$
5
  • $\begingroup$ Hi, Brent ! Only one suggestion - please format the code according to the guidelines @ help centre ;__; $\endgroup$
    – Sektor
    Feb 6, 2015 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, 2015 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, 2015 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, 2015 at 5:52
  • $\begingroup$ Can someone supply me the email of Brent Baughn so that I may continue this line of programming. Bob Wilson [email protected] $\endgroup$
    – Robert G. Wilson v
    Jun 11, 2015 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} *)
Improve this answer
$\endgroup$
9
  • 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, 2015 at 14:47
  • 2
    $\begingroup$ @Mr.Wizard Been there, done that. Don't know why. $\endgroup$
    – Dr. belisarius
    Feb 5, 2015 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, 2015 at 14:55
  • 1
    $\begingroup$ @Mr.Wizard :D :D $\endgroup$
    – Dr. belisarius
    Feb 5, 2015 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, 2015 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.

Improve this answer
$\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[[email protected] == n, a]]; 
  n = n + 2];
2 Position[res, True]

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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