Calculating weird numbers

7

2

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?

Qianyi Guo

Posted 2015-02-05T05:49:02.057

Reputation: 173

3The problem is that Subsets[Most[Divisors[720]]] eats up your RAM. – Karsten 7. – 2015-02-05T06:46:50.047

1Does "how to optimize the code without prior knowledge of properties of weird numbers" mean, that you'll accept solely a brute force solution? – Karsten 7. – 2015-02-05T07:00:08.310

1You can avoid Subsets if you set up the check as a subset sum problem. FindInstance for example could be used for this. – Daniel Lichtblau – 2015-02-05T17:30:56.060

Answers

12

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).

Simon Woods

Posted 2015-02-05T05:49:02.057

Reputation: 81 905

Beautiful. +1 .. – Mr.Wizard – 2015-02-06T04:01:34.203

8

(*  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

Brent Baughn

Posted 2015-02-05T05:49:02.057

Reputation: 81

Hi, Brent ! Only one suggestion - please format the code according to the guidelines @ help centre ;__; – Sektor – 2015-02-06T21:36:07.667

7There 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. – Daniel Lichtblau – 2015-02-07T22:29:14.537

2PS 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. – Daniel Lichtblau – 2015-02-07T22:31:24.130

@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? – Mr.Wizard – 2015-02-09T05:52:43.520

Can someone supply me the email of Brent Baughn so that I may continue this line of programming. Bob Wilson rgwv@rgwv.com – Robert G. Wilson v – 2015-06-11T03:04:43.867

5

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} *)

Dr. belisarius

Posted 2015-02-05T05:49:02.057

Reputation: 112 848

2Why 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. – Mr.Wizard – 2015-02-05T14:47:43.870

2@Mr.Wizard Been there, done that. Don't know why. – Dr. belisarius – 2015-02-05T14:55:15.983

1I 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.) – Mr.Wizard – 2015-02-05T14:55:47.267

1@Mr.Wizard :D :D – Dr. belisarius – 2015-02-05T14:56:31.727

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? – Qianyi Guo – 2015-02-05T15:26:32.913

What is NthSubset in that code? – David G. Stork – 2015-02-05T17:57:40.367

@DavidG.Stork It returns the j-th subset of the set d in Lexicographical . It's a function defined in the Combinatorica package – Dr. belisarius – 2015-02-05T18:09:54.487

@Mr.Wizard that typically comes up in my code when I originally had a more complicated function application, or if I need to play quickly with multiple functions (some of which might take additional parameters while others don't). Why it would survive in a "finished" piece of code... well... I don't know. – evanb – 2015-02-06T00:57:31.087

@Mr.Wizard I remember I once wrote f[##]&@@ twice in the same Q&A. One of these times I was even trying to show somebody else how he could make his code shorter. Luckily J.M. was there to help, twice. :P.

– Jacob Akkerboom – 2015-02-08T22:11:21.737

5

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.

Daniel Lichtblau

Posted 2015-02-05T05:49:02.057

Reputation: 52 368

4

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}}

Dr. belisarius

Posted 2015-02-05T05:49:02.057

Reputation: 112 848