Weird numbers: Difference between revisions
Ada! |
|||
Line 2,630:
::<syntaxhighlight lang="python"># anySum :: Int -> [Int] -> [Int]
from time import time▼
start = time()▼
from itertools import combinations
from math import prod
▲from time import time
primitivesp_nos = set(); primes = []
▲start = time()
def main():
weird_nos = []
n =
x = 1
while n > 0:
if isweird(x) == 1: weird_nos.append(x); n = n - 1
weird_nos.append(x)▼
x = x + 1
print("First", len(weird_nos), "weird nos:\n", weird_nos)
if x == [n]:
sum_fctrs = 1; fctrs
difference = sum_fctrs - 2 * n
if difference
if difference == 0: primitivesp_nos.add(n)▼
▲ while n % k == 0:
▲def isweird(n): # Checks if n is weird
▲ for i in fctrs:
▲ sum_fctrs = sum_fctrs * (i ** (x.count(i) + 1) - 1)//(i - 1)
▲ primitivesp_nos.add(n)
return 0
for i in range(2, len(x)):
for j in combinations(x, i):
product = prod(j)
if product not in primitivesp_nos:
else: return
x = 1; fctrs = {1}.union(set(i for i in fctrs if i <= difference))
ns = n - (difference + n - sum(fctrs))
else: primitivesp_nos.add(n); return 0
▲ primitivesp_nos.add(n)
main()
end = time()
print("Execution time: ", round(end - start, 2), "s")
-- -> [100,96]</syntaxhighlight>
|
Revision as of 14:37, 29 May 2024
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
In number theory, a weird number is a natural number that is abundant but not semiperfect (and therefore not perfect either).
In other words, the sum of the proper divisors of the number (divisors including 1 but not itself) is greater than the number itself (the number is abundant), but no subset of those divisors sums to the number itself (the number is not semiperfect).
For example:
- 12 is not a weird number.
- It is abundant; its proper divisors 1, 2, 3, 4, 6 sum to 16 (which is > 12),
- but it is semiperfect, e.g.: 6 + 4 + 2 == 12.
- 70 is a weird number.
- It is abundant; its proper divisors 1, 2, 5, 7, 10, 14, 35 sum to 74 (which is > 70),
- and there is no subset of proper divisors that sum to 70.
- Task
Find and display, here on this page, the first 25 weird numbers.
11l
F divisors(n)
V divs = [1]
[Int] divs2
V i = 2
L i * i <= n
I n % i == 0
V j = n I/ i
divs [+]= i
I i != j
divs2 [+]= j
i++
R divs2 [+] reversed(divs)
F abundant(n, divs)
R sum(divs) > n
F semiperfect(n, divs) -> Bool
I !divs.empty
V h = divs[0]
V t = divs[1..]
I n < h
R semiperfect(n, t)
E
R n == h | semiperfect(n - h, t) | semiperfect(n, t)
E
R 0B
F sieve(limit)
V w = [0B] * limit
L(i) (2 .< limit).step(2)
I w[i]
L.continue
V divs = divisors(i)
I !abundant(i, divs)
w[i] = 1B
E I semiperfect(i, divs)
L(j) (i .< limit).step(i)
w[j] = 1B
R w
V w = sieve(17'000)
V count = 0
print(‘The first 25 weird numbers:’)
L(n) (2..).step(2)
I !w[n]
print(n, end' ‘ ’)
count++
I count == 25
L.break
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Ada
This is not quite a translation of another version, in that:
- I first implemented it without looking at other implementations. This works, but is significantly slower than what others were reporting.
- I then compared to other implementations, and adapted what they had.
To be honest, I prefer the original implementation, even though it is much slower than many others, since:
- As far as I can tell, it's fairly consistent with timings, and doesn't much depend on the order of the factors.
- The "fast" approach is _highly_ dependent on the ordering of the factors. The implementers seem to have lucked into a good ordering that looks unnatural.
{{{out}}} When run on my machine, both versions produce the following output.
The first 25 weird numbers are 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Version 1 (slower; but fairly consistent)
pragma Ada_2022;
with Ada.Text_IO;
with Ada.Containers.Vectors;
with Ada.Containers.Ordered_Sets;
procedure Weird_Numbers is
package IO renames Ada.Text_IO;
-- allows us to type "IO.Put_Line" instead of "Ada.Text_IO.Put_Line"
package IntVecs is new Ada.Containers.Vectors
(Index_Type => Positive, Element_Type => Positive);
subtype IntVec is IntVecs.Vector;
use all type IntVec;
Factor_Cache : IntVec;
-- used to keep track of the factors on each iteration
package Intsets is new Ada.Containers.Ordered_Sets
(Element_Type => Positive);
subtype Intset is Intsets.Set;
Semiperfect_Cache : Intset;
Not_Abundant_Cache : Intset;
function Proper_Divisors (Value : Positive) return IntVec is
Cofactors : IntVec;
begin
Factor_Cache.Clear;
begin
Factor_Cache.Append (1);
for Factor in 2 .. Value / 2 loop
if Value rem Factor = 0 then
Factor_Cache.Append (Factor);
end if;
end loop;
return Factor_Cache;
end;
end Proper_Divisors;
function Is_Abundant (Value : Positive) return Boolean is
begin
if Not_Abundant_Cache.Contains (Value)
or else Semiperfect_Cache.Contains (Value)
then
return False;
end if;
if Proper_Divisors (Value)'Reduce ("+", 0) > Value then
return True;
end if;
Not_Abundant_Cache.Insert (Value);
return False;
end Is_Abundant;
function Subset_Of (Factors : IntVec; Ith : Positive) return IntVec is
-- returns the Ith subset of Factors
-- if Factors has N elements, then there are 2**N possibilities
-- this maps Ith to one of these possibilities using its binary representation
Result : IntVec;
Index : Natural := 1;
Ith_Remainder : Natural := Ith;
begin
while Ith_Remainder /= 0 loop
if Ith_Remainder rem 2 = 1 then
Result.Append (Factors (Index));
Ith_Remainder := @ - 1;
end if;
Ith_Remainder := @ / 2;
Index := @ + 1;
end loop;
return Result;
end Subset_Of;
function Is_Semiperfect (Value : Positive) return Boolean is
Factors : constant IntVec := Factor_Cache;
Subset : IntVec;
Sum : Natural;
begin
if (for some Previous of Semiperfect_Cache => Value rem Previous = 0)
then
return True;
end if;
for Ith in reverse 1 .. (2**Positive (Factors.Length) - 1) loop
Subset := Subset_Of (Factors, Positive (Ith));
Sum := Subset'Reduce ("+", 0);
if Sum = Value then
Semiperfect_Cache.Insert (Value);
return True;
end if;
end loop;
return False;
end Is_Semiperfect;
function Is_Weird (Value : Positive) return Boolean is
(Is_Abundant (Value) and then not Is_Semiperfect (Value));
Current : Positive := 2;
Number_Found : Natural := 0;
begin
IO.Put ("The first 25 weird numbers are");
while Number_Found < 25 loop
if Is_Weird (Current * 2) then
IO.Put (Integer'Image (Current * 2));
Number_Found := @ + 1;
end if;
Current := @ + 1;
end loop;
IO.New_Line;
end Weird_Numbers;
Version 2 (what most people are using; faster if you're lucky with the ordering of factors)
pragma Ada_2022;
with Ada.Text_IO;
with Ada.Containers.Vectors;
with Ada.Containers.Ordered_Sets;
procedure Weird_Numbers is
package IO renames Ada.Text_IO;
package IntVecs is new Ada.Containers.Vectors
(Index_Type => Positive, Element_Type => Positive);
subtype IntVec is IntVecs.Vector;
use all type IntVec;
Factor_Cache : IntVec;
package Intsets is new Ada.Containers.Ordered_Sets
(Element_Type => Positive);
subtype Intset is Intsets.Set;
Semiperfect_Cache : Intset;
Not_Abundant_Cache : Intset;
function Proper_Divisors (Value : Positive) return IntVec is
Cofactors : IntVec;
begin
Factor_Cache.Clear;
begin
Factor_Cache.Append (1);
for Factor in 2 .. Value / 2 loop
if Factor * Factor > Value then
exit;
end if;
if Value rem Factor = 0 then
Factor_Cache.Append (Factor);
if Factor * Factor /= Value then
Cofactors.Append (Value / Factor);
end if;
end if;
end loop;
Factor_Cache.Prepend_Vector (Cofactors);
-- this ordering is REALLY curious
-- for instance, the factors of 12 are ordered as 6, 4, 1, 2, 3
-- AND THIS MATTERS
-- both strictly ascending and strictly descending orders are much slower
return Factor_Cache;
end;
end Proper_Divisors;
function Is_Abundant (Value : Positive) return Boolean is
begin
if Not_Abundant_Cache.Contains (Value)
or else Semiperfect_Cache.Contains (Value)
then
return False;
end if;
if Proper_Divisors (Value)'Reduce ("+", 0) > Value then
return True;
end if;
Not_Abundant_Cache.Insert (Value);
return False;
end Is_Abundant;
function Alternate_Semiperfect (Value : Positive) return Boolean is
function Reroute
(Value : Positive; First : Positive) return Boolean
is
Head : Positive;
begin
if Factor_Cache.Last_Index >= First then
Head := Factor_Cache (First);
if Value < Head then
return Reroute (Value, First + 1);
else
return
Value = Head
or else Reroute (Value - Head, First + 1)
or else Reroute (Value, First + 1);
end if;
else
return False;
end if;
end Reroute;
begin
return
Reroute
(Value, Factor_Cache.First_Index);
end Alternate_Semiperfect;
function Is_Weird (Value : Positive) return Boolean is
(Is_Abundant (Value) and then not Alternate_Semiperfect (Value));
Current : Positive := 2;
Number_Found : Natural := 0;
begin
IO.Put ("The first 25 weird numbers are");
while Number_Found < 25 loop
if Is_Weird (Current * 2) then
IO.Put (Integer'Image (Current * 2));
Number_Found := @ + 1;
end if;
Current := @ + 1;
end loop;
IO.New_Line;
end Weird_Numbers;
ALGOL 68
Translation of the untweaked Go version 1 sample. Avoids creating separate array slices in the semiperfect routine, to save memory for Algol 68G version 2.8.3.
BEGIN # find wierd numbers - abundant but not semiperfect numbers - translation of Go #
# returns the divisors of n in descending order #
PROC divisors = ( INT n )[]INT:
BEGIN
INT max divs = 2 * ENTIER sqrt( n );
[ 1 : max divs ]INT divs;
[ 1 : max divs ]INT divs2;
INT d pos := 0, d2 pos := 0;
divs[ d pos +:= 1 ] := 1;
FOR i FROM 2 WHILE i * i <= n DO
IF n MOD i = 0 THEN
INT j = n OVER i;
divs[ d pos +:= 1 ] := i;
IF i /= j THEN divs2[ d2 pos +:= 1 ] := j FI
FI
OD;
FOR i FROM d pos BY -1 WHILE i > 0 DO
divs2[ d2 pos +:= 1 ] := divs[ i ]
OD;
divs2[ 1 : d2 pos ]
END # divisors # ;
# returns TRUE if n with divisors divs, is abundant, FALSE otherwise #
PROC abundant = ( INT n, []INT divs )BOOL:
BEGIN
INT sum := 0;
FOR i FROM LWB divs TO UPB divs DO sum +:= divs[ i ] OD;
sum > n
END # abundant # ;
# returns TRUE if n with divisors divs, is semiperfect, FALSE otherwise #
PROC semiperfect = ( INT n, []INT divs, INT lb, ub )BOOL:
IF ub < lb
THEN FALSE
ELIF INT h = divs[ lb ];
n < h
THEN semiperfect( n, divs, lb + 1, ub )
ELIF n = h
THEN TRUE
ELIF semiperfect( n - h, divs, lb + 1, ub )
THEN TRUE
ELSE semiperfect( n, divs, lb + 1, ub )
FI # semiperfect # ;
# returns a sieve where FALSE = abundant and not semiperfect #
PROC sieve = ( INT limit )[]BOOL:
BEGIN # Only interested in even numbers >= 2 #
[ 1 : limit ]BOOL w; FOR i FROM 1 TO limit DO w[ i ] := FALSE OD;
FOR i FROM 2 BY 2 TO limit DO
IF NOT w[ i ] THEN
[]INT divs = divisors( i );
IF NOT abundant( i, divs ) THEN
w[ i ] := TRUE
ELIF semiperfect( i, divs, LWB divs, UPB divs ) THEN
FOR j FROM i BY i TO limit DO w[ j ] := TRUE OD
FI
FI
OD;
w
END # sieve # ;
BEGIN # task #
[]BOOL w = sieve( 17 000 );
INT count := 0;
INT max = 25;
print( ( "The first 25 weird numbers are:", newline ) );
FOR n FROM 2 BY 2 WHILE count < max DO
IF NOT w[ n ] THEN
print( ( whole( n, 0 ), " " ) );
count +:= 1
FI
OD;
print( ( newline ) )
END
END
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
AppleScript
Functional
Applescript is not the recommended apparatus for this kind of experiment.
(Though after about 6 seconds (on this system) it does yield the first 25, and intermediates can be logged in the Messages channel of macOS Script Editor).
on run
take(25, weirds())
-- Gets there, but takes about 6 seconds on this system,
-- (logging intermediates through the Messages channel, for the impatient :-)
end run
-- weirds :: Gen [Int]
on weirds()
script
property x : 1
property v : 0
on |λ|()
repeat until isWeird(x)
set x to 1 + x
end repeat
set v to x
log v
set x to 1 + x
return v
end |λ|
end script
end weirds
-- isWeird :: Int -> Bool
on isWeird(n)
set ds to descProperDivisors(n)
set d to sum(ds) - n
0 < d and not hasSum(d, ds)
end isWeird
-- hasSum :: Int -> [Int] -> Bool
on hasSum(n, xs)
if {} ≠ xs then
set h to item 1 of xs
set t to rest of xs
if n < h then
hasSum(n, t)
else
n = h or hasSum(n - h, t) or hasSum(n, t)
end if
else
false
end if
end hasSum
-- GENERIC ------------------------------------------------
-- descProperDivisors :: Int -> [Int]
on descProperDivisors(n)
if n = 1 then
{1}
else
set realRoot to n ^ (1 / 2)
set intRoot to realRoot as integer
set blnPerfect to intRoot = realRoot
-- isFactor :: Int -> Bool
script isFactor
on |λ|(x)
n mod x = 0
end |λ|
end script
-- Factors up to square root of n,
set lows to filter(isFactor, enumFromTo(1, intRoot))
-- and cofactors of these beyond the square root,
-- integerQuotient :: Int -> Int
script integerQuotient
on |λ|(x)
(n / x) as integer
end |λ|
end script
set t to rest of lows
if blnPerfect then
set xs to t
else
set xs to lows
end if
map(integerQuotient, t) & (reverse of xs)
end if
end descProperDivisors
-- enumFromTo :: (Int, Int) -> [Int]
on enumFromTo(m, n)
if m ≤ n then
set lst to {}
repeat with i from m to n
set end of lst to i
end repeat
return lst
else
return {}
end if
end enumFromTo
-- filter :: (a -> Bool) -> [a] -> [a]
on filter(f, xs)
tell mReturn(f)
set lst to {}
set lng to length of xs
repeat with i from 1 to lng
set v to item i of xs
if |λ|(v, i, xs) then set end of lst to v
end repeat
return lst
end tell
end filter
-- foldl :: (a -> b -> a) -> a -> [b] -> a
on foldl(f, startValue, xs)
tell mReturn(f)
set v to startValue
set lng to length of xs
repeat with i from 1 to lng
set v to |λ|(v, item i of xs, i, xs)
end repeat
return v
end tell
end foldl
-- map :: (a -> b) -> [a] -> [b]
on map(f, xs)
tell mReturn(f)
set lng to length of xs
set lst to {}
repeat with i from 1 to lng
set end of lst to |λ|(item i of xs, i, xs)
end repeat
return lst
end tell
end map
-- sum :: [Num] -> Num
on sum(xs)
script add
on |λ|(a, b)
a + b
end |λ|
end script
foldl(add, 0, xs)
end sum
-- take :: Int -> Gen [a] -> [a]
on take(n, xs)
set ys to {}
repeat with i from 1 to n
set v to xs's |λ|()
if missing value is v then
return ys
else
set end of ys to v
end if
end repeat
return ys
end take
-- Lift 2nd class handler function into 1st class script wrapper
-- mReturn :: First-class m => (a -> b) -> m (a -> b)
on mReturn(f)
if script is class of f then
f
else
script
property |λ| : f
end script
end if
end mReturn
- Output:
{70, 836, 4030, 5830, 7192, 7912, 9272, 10430, 10570, 10792, 10990, 11410, 11690, 12110, 12530, 12670, 13370, 13510, 13790, 13930, 14770, 15610, 15890, 16030, 16310}
Idiomatic
0.69 seconds:
-- Sum n's proper divisors.
on aliquotSum(n)
if (n < 2) then return 0
set sum to 1
set sqrt to n ^ 0.5
set limit to sqrt div 1
if (limit = sqrt) then
set sum to sum + limit
set limit to limit - 1
end if
repeat with i from 2 to limit
if (n mod i is 0) then set sum to sum + i + n div i
end repeat
return sum
end aliquotSum
-- Return n's proper divisors.
on properDivisors(n)
set output to {}
if (n > 1) then
set sqrt to n ^ 0.5
set limit to sqrt div 1
if (limit = sqrt) then
set end of output to limit
set limit to limit - 1
end if
repeat with i from limit to 2 by -1
if (n mod i is 0) then
set beginning of output to i
set end of output to n div i
end if
end repeat
set beginning of output to 1
end if
return output
end properDivisors
-- Does a subset of the given list of numbers add up to the target value?
on subsetOf:numberList sumsTo:target
script o
property lst : numberList
property someNegatives : false
on ssp(target, i)
repeat while (i > 1)
set n to item i of my lst
set i to i - 1
if ((n = target) or (((n < target) or (someNegatives)) and (ssp(target - n, i)))) then return true
end repeat
return (target = beginning of my lst)
end ssp
end script
-- The search can be more efficient if it's known the list contains no negatives.
repeat with n in o's lst
if (n < 0) then
set o's someNegatives to true
exit repeat
end if
end repeat
return o's ssp(target, count o's lst)
end subsetOf:sumsTo:
-- Is n a weird number?
on isWeird(n)
-- Yes if its aliquot sum's greater than it and no subset of its proper divisors adds up to it.
-- Using aliquotSum() to get the divisor sum and then calling properDivisors() too if a list's actually
-- needed is generally faster than calling properDivisors() in the first place and summing the result.
set sum to aliquotSum(n)
if (sum > n) then
set divisors to properDivisors(n)
-- Check that no subset sums to the smaller (usually the latter) of n and sum - n.
tell (sum - n) to if (it < n) then set n to it
return (not (my subsetOf:divisors sumsTo:n))
else
return false
end if
end isWeird
-- Task code:
on weirdNumbers(target)
script o
property weirds : {}
end script
set n to 2
set counter to 0
repeat until (counter = target)
if (isWeird(n)) then
set end of o's weirds to n
set counter to counter + 1
end if
set n to n + 1
end repeat
return o's weirds
end weirdNumbers
weirdNumbers(25)
- Output:
{70, 836, 4030, 5830, 7192, 7912, 9272, 10430, 10570, 10792, 10990, 11410, 11690, 12110, 12530, 12670, 13370, 13510, 13790, 13930, 14770, 15610, 15890, 16030, 16310}
C
#include "stdio.h"
#include "stdlib.h"
#include "stdbool.h"
#include "string.h"
struct int_a {
int *ptr;
size_t size;
};
struct int_a divisors(int n) {
int *divs, *divs2, *out;
int i, j, c1 = 0, c2 = 0;
struct int_a array;
divs = malloc(n * sizeof(int) / 2);
divs2 = malloc(n * sizeof(int) / 2);
divs[c1++] = 1;
for (i = 2; i * i <= n; i++) {
if (n % i == 0) {
j = n / i;
divs[c1++] = i;
if (i != j) {
divs2[c2++] = j;
}
}
}
out = malloc((c1 + c2) * sizeof(int));
for (int i = 0; i < c2; i++) {
out[i] = divs2[i];
}
for (int i = 0; i < c1; i++) {
out[c2 + i] = divs[c1 - i - 1];
}
array.ptr = out;
array.size = c1 + c2;
free(divs);
free(divs2);
return array;
}
bool abundant(int n, struct int_a divs) {
int sum = 0;
int i;
for (i = 0; i < divs.size; i++) {
sum += divs.ptr[i];
}
return sum > n;
}
bool semiperfect(int n, struct int_a divs) {
if (divs.size > 0) {
int h = *divs.ptr;
int *t = divs.ptr + 1;
struct int_a ta;
ta.ptr = t;
ta.size = divs.size - 1;
if (n < h) {
return semiperfect(n, ta);
} else {
return n == h
|| semiperfect(n - h, ta)
|| semiperfect(n, ta);
}
} else {
return false;
}
}
bool *sieve(int limit) {
bool *w = calloc(limit, sizeof(bool));
struct int_a divs;
int i, j;
for (i = 2; i < limit; i += 2) {
if (w[i]) continue;
divs = divisors(i);
if (!abundant(i, divs)) {
w[i] = true;
} else if (semiperfect(i, divs)) {
for (j = i; j < limit; j += i) {
w[j] = true;
}
}
}
free(divs.ptr);
return w;
}
int main() {
bool *w = sieve(17000);
int count = 0;
int max = 25;
int n;
printf("The first 25 weird numbers:\n");
for (n = 2; count < max; n += 2) {
if (!w[n]) {
printf("%d ", n);
count++;
}
}
printf("\n");
free(w);
return 0;
}
- Output:
70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
C#
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace WeirdNumbers {
class Program {
static List<int> Divisors(int n) {
List<int> divs = new List<int> { 1 };
List<int> divs2 = new List<int>();
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
int j = n / i;
divs.Add(i);
if (i != j) {
divs2.Add(j);
}
}
}
divs.Reverse();
divs2.AddRange(divs);
return divs2;
}
static bool Abundant(int n, List<int> divs) {
return divs.Sum() > n;
}
static bool Semiperfect(int n, List<int> divs) {
if (divs.Count > 0) {
var h = divs[0];
var t = divs.Skip(1).ToList();
if (n < h) {
return Semiperfect(n, t);
} else {
return n == h
|| Semiperfect(n - h, t)
|| Semiperfect(n, t);
}
} else {
return false;
}
}
static List<bool> Sieve(int limit) {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
bool[] w = new bool[limit];
for (int i = 2; i < limit; i += 2) {
if (w[i]) continue;
var divs = Divisors(i);
if (!Abundant(i, divs)) {
w[i] = true;
} else if (Semiperfect(i, divs)) {
for (int j = i; j < limit; j += i) {
w[j] = true;
}
}
}
return w.ToList();
}
static void Main() {
var w = Sieve(17_000);
int count = 0;
int max = 25;
Console.WriteLine("The first 25 weird numbers:");
for (int n = 2; count < max; n += 2) {
if (!w[n]) {
Console.Write("{0} ", n);
count++;
}
}
Console.WriteLine();
}
}
}
- Output:
70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
C++
#include <algorithm>
#include <iostream>
#include <numeric>
#include <vector>
std::vector<int> divisors(int n) {
std::vector<int> divs = { 1 };
std::vector<int> divs2;
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
int j = n / i;
divs.push_back(i);
if (i != j) {
divs2.push_back(j);
}
}
}
std::copy(divs.cbegin(), divs.cend(), std::back_inserter(divs2));
return divs2;
}
bool abundant(int n, const std::vector<int> &divs) {
return std::accumulate(divs.cbegin(), divs.cend(), 0) > n;
}
template<typename IT>
bool semiperfect(int n, const IT &it, const IT &end) {
if (it != end) {
auto h = *it;
auto t = std::next(it);
if (n < h) {
return semiperfect(n, t, end);
} else {
return n == h
|| semiperfect(n - h, t, end)
|| semiperfect(n, t, end);
}
} else {
return false;
}
}
template<typename C>
bool semiperfect(int n, const C &c) {
return semiperfect(n, std::cbegin(c), std::cend(c));
}
std::vector<bool> sieve(int limit) {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
std::vector<bool> w(limit);
for (int i = 2; i < limit; i += 2) {
if (w[i]) continue;
auto divs = divisors(i);
if (!abundant(i, divs)) {
w[i] = true;
} else if (semiperfect(i, divs)) {
for (int j = i; j < limit; j += i) {
w[j] = true;
}
}
}
return w;
}
int main() {
auto w = sieve(17000);
int count = 0;
int max = 25;
std::cout << "The first 25 weird numbers:";
for (int n = 2; count < max; n += 2) {
if (!w[n]) {
std::cout << n << ' ';
count++;
}
}
std::cout << '\n';
return 0;
}
- Output:
The first 25 weird numbers:70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Crystal
def divisors(n : Int32) : Array(Int32)
divs = [1]
divs2 = [] of Int32
i = 2
while i * i < n
if n % i == 0
j = n // i
divs << i
divs2 << j if i != j
end
i += 1
end
i = divs.size - 1
# TODO: Use reverse
while i >= 0
divs2 << divs[i]
i -= 1
end
divs2
end
def abundant(n : Int32, divs : Array(Int32)) : Bool
divs.sum > n
end
def semiperfect(n : Int32, divs : Array(Int32)) : Bool
if divs.size > 0
h = divs[0]
t = divs[1..]
return n < h ? semiperfect(n, t) : n == h || semiperfect(n - h, t) || semiperfect(n, t)
end
return false
end
def sieve(limit : Int32) : Array(Bool)
# false denotes abundant and not semi-perfect.
# Only interested in even numbers >= 2
w = Array(Bool).new(limit, false) # An array filled with 'false'
i = 2
while i < limit
if !w[i]
divs = divisors i
if !abundant(i, divs)
w[i] = true
elsif semiperfect(i, divs)
j = i
while j < limit
w[j] = true
j += i
end
end
end
i += 2
end
w
end
def main
w = sieve 17000
count = 0
max = 25
print "The first 25 weird numbers are: "
n = 2
while count < max
if !w[n]
print "#{n} "
count += 1
end
n += 2
end
puts "\n"
end
require "benchmark"
puts Benchmark.measure { main }
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 # Benchmark with --release flag 0.046875 0.000000 0.046875 ( 0.040754)
D
Adding efficient "cut" condition in semiperfect recursive algorithm
import std.algorithm;
import std.array;
import std.stdio;
int[] divisors(int n) {
int[] divs = [1];
int[] divs2;
for (int i = 2; i * i <= n; i++) {
if (n % i == 0) {
int j = n / i;
divs ~= i;
if (i != j) {
divs2 ~= j;
}
}
}
divs2 ~= divs.reverse;
return divs2;
}
bool abundant(int n, int[] divs) {
return divs.sum() > n;
}
bool semiperfect(int n, int[] divs) {
// This algorithm is O(2^N) for N == divs.length when number is not semiperfect.
// Comparing with (divs.sum < n) instead (divs.length==0) removes unnecessary
// recursive binary tree branches.
auto s = divs.sum;
if(s == n)
return true;
else if ( s<n )
return false;
else {
auto h = divs[0];
auto t = divs[1..$];
if (n < h) {
return semiperfect(n, t);
} else {
return n == h
// Supossin h is part of the sum
|| semiperfect(n - h, t)
// Supossin h is not part of the sum
|| semiperfect(n, t);
}
}
}
bool[] sieve(int limit) {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
auto w = uninitializedArray!(bool[])(limit);
w[] = false;
for (int i = 2; i < limit; i += 2) {
if (w[i]) continue;
auto divs = divisors(i);
if (!abundant(i, divs)) {
w[i] = true;
} else if (semiperfect(i, divs)) {
for (int j = i; j < limit; j += i) {
w[j] = true;
}
}
}
return w;
}
void main() {
auto w = sieve(17_000);
int count = 0;
int max = 25;
writeln("The first 25 weird numbers:");
for (int n = 2; count < max; n += 2) {
if (!w[n]) {
write(n, ' ');
count++;
}
}
writeln;
}
- Output:
70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
F#
let divisors n = [1..n/2] |> List.filter (fun x->n % x = 0)
let abundant (n:int) divs = Seq.sum(divs) > n
let rec semiperfect (n:int) (divs:List<int>) =
if divs.Length > 0 then
let h = divs.Head
let t = divs.Tail
if n < h then
semiperfect n t
else
n = h || (semiperfect (n - h) t) || (semiperfect n t)
else false
let weird n =
let d = divisors n
if abundant n d then
not(semiperfect n d)
else
false
[<EntryPoint>]
let main _ =
let mutable i = 1
let mutable count = 0
while (count < 25) do
if (weird i) then
count <- count + 1
printf "%d -> %d\n" count i
i <- i + 1
0 // return an integer exit code
- Output:
1 -> 70 2 -> 836 3 -> 4030 4 -> 5830 5 -> 7192 6 -> 7912 7 -> 9272 8 -> 10430 9 -> 10570 10 -> 10792 11 -> 10990 12 -> 11410 13 -> 11690 14 -> 12110 15 -> 12530 16 -> 12670 17 -> 13370 18 -> 13510 19 -> 13790 20 -> 13930 21 -> 14770 22 -> 15610 23 -> 15890 24 -> 16030 25 -> 16310
Factor
The has-sum?
word is a translation of the Haskell function.
USING: combinators.short-circuit io kernel lists lists.lazy
locals math math.primes.factors prettyprint sequences ;
IN: rosetta-code.weird-numbers
:: has-sum? ( n seq -- ? )
seq [ f ] [
unclip-slice :> ( xs x )
n x < [ n xs has-sum? ] [
{
[ n x = ]
[ n x - xs has-sum? ]
[ n xs has-sum? ]
} 0||
] if
] if-empty ;
: weird? ( n -- ? )
dup divisors but-last reverse
{ [ sum < ] [ has-sum? not ] } 2&& ;
: weirds ( -- list ) 1 lfrom [ weird? ] lfilter ;
: weird-numbers-demo ( -- )
"First 25 weird numbers:" print
25 weirds ltake list>array . ;
MAIN: weird-numbers-demo
- Output:
First 25 weird numbers: { 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 }
FreeBASIC
Function GetFactors(n As Long,r() As Long) As Long
Redim r(0)
r(0)=1
Dim As Long count,acc
For z As Long=2 To n\2
If n Mod z=0 Then
count+=1:redim preserve r(0 to count)
r(count)=z
acc+=z
End If
Next z
Return 1+acc
End Function
sub sumcombinations(arr() As Long,n As Long,r As Long,index As Long,_data() As Long,i As Long,Byref ans As Long,ref As Long)
Dim As Long acc
If index=r Then
For j As Long=0 To r-1
acc+=_data(j)
If acc=ref Then ans=1:Return
If acc>ref then return
Next j
Return
End If
If i>=n Or ans<>0 Then Return
_data(index) = arr(i)
sumcombinations(arr(),n,r,index + 1,_data(),i+1,ans,ref)
sumcombinations(arr(),n,r,index,_data(),i+1,ans,ref)
End sub
Function IsWeird(u() As Long,num As Long) As Long
Redim As Long d()
Dim As Long ans
For r As Long=2 To Ubound(u)
Redim d(r)
ans=0
sumcombinations(u(),Ubound(u)+1,r,0,d(),0,ans,num)
If ans =1 Then Return 0
Next r
Return 1
End Function
Redim As Long u()
Dim As Long SumFactors,number=2,count
Do
number+=2
SumFactors=GetFactors(number,u())
If SumFactors>number Then
If IsWeird(u(),number) Then Print number;" ";:count+=1
End If
Loop Until count=25
Print
Print "first 25 done"
Sleep
- Output:
70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 first 25 done
Go
Version 1
This takes advantage of Hout's analysis (see talk page) when testing for primitive semi-perfect numbers.
It also uses a sieve so we can make use of the fact that all multiples of a semi-perfect number are themselves semi-perfect.
Runs in less than 10 ms on an Intel Core i7-8565U machine. The first fifty (with a sieve size of 27000) takes roughly double that.
When run on the same machine, the 'tweaked' version (linked to below), which was supplied by Enter your username, is almost 3 times faster than this.
package main
import "fmt"
func divisors(n int) []int {
divs := []int{1}
divs2 := []int{}
for i := 2; i*i <= n; i++ {
if n%i == 0 {
j := n / i
divs = append(divs, i)
if i != j {
divs2 = append(divs2, j)
}
}
}
for i := len(divs) - 1; i >= 0; i-- {
divs2 = append(divs2, divs[i])
}
return divs2
}
func abundant(n int, divs []int) bool {
sum := 0
for _, div := range divs {
sum += div
}
return sum > n
}
func semiperfect(n int, divs []int) bool {
le := len(divs)
if le > 0 {
h := divs[0]
t := divs[1:]
if n < h {
return semiperfect(n, t)
} else {
return n == h || semiperfect(n-h, t) || semiperfect(n, t)
}
} else {
return false
}
}
func sieve(limit int) []bool {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
w := make([]bool, limit)
for i := 2; i < limit; i += 2 {
if w[i] {
continue
}
divs := divisors(i)
if !abundant(i, divs) {
w[i] = true
} else if semiperfect(i, divs) {
for j := i; j < limit; j += i {
w[j] = true
}
}
}
return w
}
func main() {
w := sieve(17000)
count := 0
const max = 25
fmt.Println("The first 25 weird numbers are:")
for n := 2; count < max; n += 2 {
if !w[n] {
fmt.Printf("%d ", n)
count++
}
}
fmt.Println()
}
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Version 2 (Tweaked)
Link to a tweaked version at Try it Online!
Runs in under 6.0ms on the Tio server. The first fifty (with a sieve size of 26,533) takes under 12.0ms. Comments added where tweaks were applied
Haskell
weirds :: [Int]
weirds = filter abundantNotSemiperfect [1 ..]
abundantNotSemiperfect :: Int -> Bool
abundantNotSemiperfect n =
let ds = descProperDivisors n
d = sum ds - n
in 0 < d && not (hasSum d ds)
hasSum :: Int -> [Int] -> Bool
hasSum _ [] = False
hasSum n (x:xs)
| n < x = hasSum n xs
| otherwise = (n == x) || hasSum (n - x) xs || hasSum n xs
descProperDivisors
:: Integral a
=> a -> [a]
descProperDivisors n =
let root = (floor . sqrt) (fromIntegral n :: Double)
lows = filter ((0 ==) . rem n) [root,root - 1 .. 1]
factors
| n == root ^ 2 = tail lows
| otherwise = lows
in tail $ reverse (quot n <$> lows) ++ factors
main :: IO ()
main =
(putStrLn . unlines) $
zipWith (\i x -> show i ++ (" -> " ++ show x)) [1 ..] (take 25 weirds)
- Output:
1 -> 70 2 -> 836 3 -> 4030 4 -> 5830 5 -> 7192 6 -> 7912 7 -> 9272 8 -> 10430 9 -> 10570 10 -> 10792 11 -> 10990 12 -> 11410 13 -> 11690 14 -> 12110 15 -> 12530 16 -> 12670 17 -> 13370 18 -> 13510 19 -> 13790 20 -> 13930 21 -> 14770 22 -> 15610 23 -> 15890 24 -> 16030 25 -> 16310
J
This algorithm uses a sieve to eliminate multiples of semiperfect numbers from future testing.
factor=: [: }: [: , [: */&> [: { [: <@(^ i.@>:)/"1 [: |: __&q:
classify=: 3 : 0
weird =: perfect =: deficient =: abundant =: i. 0
a=: (i. -. 0 , deficient =: 1 , i.&.:(p:inv)) y NB. a are potential semi-perfect numbers
for_n. a do.
if. n e. a do.
factors=. factor n
sf =. +/ factors
if. sf < n do.
deficient =: deficient , n
else.
if. n < sf do.
abundant=: abundant , n
else.
perfect =: perfect , n
a =: a -. (2+i.)@<.&.(%&n) y NB. remove multiples of perfect numbers
continue.
end.
NB. compute sums of subsets to detect semiperfection
NB. the following algorithm correctly finds weird numbers less than 20000
NB. remove large terms necessary for the sum to reduce the Catalan tally of sets
factors =. /:~ factors NB. ascending sort
NB. if the sum of the length one outfixes is less n then the factor is required in the semiperfect set.
i_required =. n (1 i.~ (>(1+/\.]))) factors
target =. n - +/ i_required }. factors
t =. i_required {. factors
NB. work in chunks of 2^16 to reduce memory requirement
sp =. target e. ; (,:~2^16) <@([: +/"1 t #~ (_ ,(#t)) {. #:);.3 i. 2 ^ # t
if. sp do.
a =: a -. (2+i.)@<.&.(%&n) y NB. remove multiples of semi perfect numbers
else.
weird =: weird , n
a =: a -. n
end.
end.
end.
end.
a =: a -. deficient
weird
)
classify 20000 NB. the first 36 weird numbers 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 16730 16870 17272 17570 17990 18410 18830 18970 19390 19670 19810
Java
import java.util.ArrayList;
import java.util.List;
public class WeirdNumbers {
public static void main(String[] args) {
int n = 2;
// n += 2 : No odd weird numbers < 10^21
for ( int count = 1 ; count <= 25 ; n += 2 ) {
if ( isWeird(n) ) {
System.out.printf("w(%d) = %d%n", count, n);
count++;
}
}
}
private static boolean isWeird(int n) {
List<Integer> properDivisors = getProperDivisors(n);
return isAbundant(properDivisors, n) && ! isSemiPerfect(properDivisors, n);
}
private static boolean isAbundant(List<Integer> divisors, int n) {
int divisorSum = divisors.stream().mapToInt(i -> i.intValue()).sum();
return divisorSum > n;
}
// Use Dynamic Programming
private static boolean isSemiPerfect(List<Integer> divisors, int sum) {
int size = divisors.size();
// The value of subset[i][j] will be true if there is a subset of divisors[0..j-1] with sum equal to i
boolean subset[][] = new boolean[sum+1][size+1];
// If sum is 0, then answer is true
for (int i = 0; i <= size; i++) {
subset[0][i] = true;
}
// If sum is not 0 and set is empty, then answer is false
for (int i = 1; i <= sum; i++) {
subset[i][0] = false;
}
// Fill the subset table in bottom up manner
for ( int i = 1 ; i <= sum ; i++ ) {
for ( int j = 1 ; j <= size ; j++ ) {
subset[i][j] = subset[i][j-1];
int test = divisors.get(j-1);
if ( i >= test ) {
subset[i][j] = subset[i][j] || subset[i - test][j-1];
}
}
}
return subset[sum][size];
}
private static final List<Integer> getProperDivisors(int number) {
List<Integer> divisors = new ArrayList<Integer>();
long sqrt = (long) Math.sqrt(number);
for ( int i = 1 ; i <= sqrt ; i++ ) {
if ( number % i == 0 ) {
divisors.add(i);
int div = number / i;
if ( div != i && div != number ) {
divisors.add(div);
}
}
}
return divisors;
}
}
- Output:
w(1) = 70 w(2) = 836 w(3) = 4030 w(4) = 5830 w(5) = 7192 w(6) = 7912 w(7) = 9272 w(8) = 10430 w(9) = 10570 w(10) = 10792 w(11) = 10990 w(12) = 11410 w(13) = 11690 w(14) = 12110 w(15) = 12530 w(16) = 12670 w(17) = 13370 w(18) = 13510 w(19) = 13790 w(20) = 13930 w(21) = 14770 w(22) = 15610 w(23) = 15890 w(24) = 16030 w(25) = 16310
JavaScript
ES6
(() => {
'use strict';
// main :: IO ()
const main = () =>
take(25, weirds());
// weirds :: Gen [Int]
function* weirds() {
let
x = 1,
i = 1;
while (true) {
x = until(isWeird, succ, x)
console.log(i.toString() + ' -> ' + x)
yield x;
x = 1 + x;
i = 1 + i;
}
}
// isWeird :: Int -> Bool
const isWeird = n => {
const
ds = descProperDivisors(n),
d = sum(ds) - n;
return 0 < d && !hasSum(d, ds)
};
// hasSum :: Int -> [Int] -> Bool
const hasSum = (n, xs) => {
const go = (n, xs) =>
0 < xs.length ? (() => {
const
h = xs[0],
t = xs.slice(1);
return n < h ? (
go(n, t)
) : (
n == h || hasSum(n - h, t) || hasSum(n, t)
);
})() : false;
return go(n, xs);
};
// descProperDivisors :: Int -> [Int]
const descProperDivisors = n => {
const
rRoot = Math.sqrt(n),
intRoot = Math.floor(rRoot),
blnPerfect = rRoot === intRoot,
lows = enumFromThenTo(intRoot, intRoot - 1, 1)
.filter(x => (n % x) === 0);
return (
reverse(lows)
.slice(1)
.map(x => n / x)
).concat((blnPerfect ? tail : id)(lows))
};
// GENERIC FUNCTIONS ----------------------------
// enumFromThenTo :: Int -> Int -> Int -> [Int]
const enumFromThenTo = (x1, x2, y) => {
const d = x2 - x1;
return Array.from({
length: Math.floor(y - x2) / d + 2
}, (_, i) => x1 + (d * i));
};
// id :: a -> a
const id = x => x;
// reverse :: [a] -> [a]
const reverse = xs =>
'string' !== typeof xs ? (
xs.slice(0).reverse()
) : xs.split('').reverse().join('');
// succ :: Enum a => a -> a
const succ = x => 1 + x;
// sum :: [Num] -> Num
const sum = xs => xs.reduce((a, x) => a + x, 0);
// tail :: [a] -> [a]
const tail = xs => 0 < xs.length ? xs.slice(1) : [];
// take :: Int -> [a] -> [a]
// take :: Int -> String -> String
const take = (n, xs) =>
'GeneratorFunction' !== xs.constructor.constructor.name ? (
xs.slice(0, n)
) : [].concat.apply([], Array.from({
length: n
}, () => {
const x = xs.next();
return x.done ? [] : [x.value];
}));
// until :: (a -> Bool) -> (a -> a) -> a -> a
const until = (p, f, x) => {
let v = x;
while (!p(v)) v = f(v);
return v;
};
// MAIN ---
return main();
})();
- Output:
1 -> 70 2 -> 836 3 -> 4030 4 -> 5830 5 -> 7192 6 -> 7912 7 -> 9272 8 -> 10430 9 -> 10570 10 -> 10792 11 -> 10990 12 -> 11410 13 -> 11690 14 -> 12110 15 -> 12530 16 -> 12670 17 -> 13370 18 -> 13510 19 -> 13790 20 -> 13930 21 -> 14770 22 -> 15610 23 -> 15890 24 -> 16030 25 -> 16310
jq
Adapted from Wren
For an explanation, see the Go entry.
The following also works with gojq, the Go implementation of jq, though much more slowly.
# unordered
def proper_divisors:
. as $n
| if $n > 1 then 1,
( range(2; 1 + (sqrt|floor)) as $i
| if ($n % $i) == 0 then $i,
(($n / $i) | if . == $i then empty else . end)
else empty
end)
else empty
end;
# Is n semiperfect given that divs are the proper divisors
def semiperfect(n; divs):
(divs|length) as $le
| if $le == 0 then false
else divs[0] as $h
| if n == $h then true
elif $le == 1 then false
else divs[1:] as $t
| if n < $h then semiperfect(n; $t)
else semiperfect(n-$h; $t) or semiperfect(n; $t)
end
end
end ;
def sieve(limit):
# 'false' denotes abundant and not semi-perfect.
# Only interested in even numbers >= 2
(reduce range(6; limit; 6) as $j ([]; .[$j] = true)) # eliminates multiples of 3
| reduce range(2; limit; 2) as $i (.;
if (.[$i]|not)
then [$i|proper_divisors] as $divs
| ($divs | add) as $sum
| if $sum <= $i
then .[$i] = true
elif (semiperfect($sum-$i; $divs))
then reduce range($i; limit; $i) as $j (.; .[$j] = true)
else .
end
else .
end) ;
# Print up to $max weird numbers based on the given sieve size, $limit.
def task($limit; $max):
sieve($limit) as $w
| def weirds:
range(2; $w|length; 2) | select($w[.]|not);
# collect into an array for ease of counting
[limit($max; weirds)]
| "The first \(length) weird numbers are:", . ;
# The parameters should be set on the command line:
task($sieve; $limit)
Invocation:
jq -nrc --argjson sieve 16313 --argjson limit 25 -f weird.jq
- Output:
The first 25 weird numbers are: [70,836,4030,5830,7192,7912,9272,10430,10570,10792,10990,11410,11690,12110,12530,12670,13370,13510,13790,13930,14770,15610,15890,16030,16310]
Julia
using Primes
function nosuchsum(revsorted, num)
if sum(revsorted) < num
return true
end
for (i, n) in enumerate(revsorted)
if n > num
continue
elseif n == num
return false
elseif !nosuchsum(revsorted[i+1:end], num - n)
return false
end
end
true
end
function isweird(n)
if n < 70 || isodd(n)
return false
else
f = [one(n)]
for (p, x) in factor(n)
f = reduce(vcat, [f*p^i for i in 1:x], init=f)
end
pop!(f)
return sum(f) > n && nosuchsum(sort(f, rev=true), n)
end
end
function testweird(N)
println("The first $N weird numbers are: ")
count, n = 0, 69
while count < N
if isweird(n)
count += 1
print("$n ")
end
n += 1
end
end
testweird(25)
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Kotlin
// Version 1.3.21
fun divisors(n: Int): List<Int> {
val divs = mutableListOf(1)
val divs2 = mutableListOf<Int>()
var i = 2
while (i * i <= n) {
if (n % i == 0) {
val j = n / i
divs.add(i)
if (i != j) divs2.add(j)
}
i++
}
divs2.addAll(divs.asReversed())
return divs2
}
fun abundant(n: Int, divs: List<Int>) = divs.sum() > n
fun semiperfect(n: Int, divs: List<Int>): Boolean {
if (divs.size > 0) {
val h = divs[0]
val t = divs.subList(1, divs.size)
if (n < h) {
return semiperfect(n, t)
} else {
return n == h || semiperfect(n-h, t) || semiperfect(n, t)
}
} else {
return false
}
}
fun sieve(limit: Int): BooleanArray {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
val w = BooleanArray(limit)
for (i in 2 until limit step 2) {
if (w[i]) continue
val divs = divisors(i)
if (!abundant(i, divs)) {
w[i] = true
} else if (semiperfect(i, divs)) {
for (j in i until limit step i) w[j] = true
}
}
return w
}
fun main() {
val w = sieve(17000)
var count = 0
val max = 25
println("The first 25 weird numbers are:")
var n = 2
while (count < max) {
if (!w[n]) {
print("$n ")
count++
}
n += 2
}
println()
}
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Lua
function make(n, d)
local a = {}
for i=1,n do
table.insert(a, d)
end
return a
end
function reverse(t)
local n = #t
local i = 1
while i < n do
t[i],t[n] = t[n],t[i]
i = i + 1
n = n - 1
end
end
function tail(list)
return { select(2, unpack(list)) }
end
function divisors(n)
local divs = {}
table.insert(divs, 1)
local divs2 = {}
local i = 2
while i * i <= n do
if n % i == 0 then
local j = n / i
table.insert(divs, i)
if i ~= j then
table.insert(divs2, j)
end
end
i = i + 1
end
reverse(divs)
for i,v in pairs(divs) do
table.insert(divs2, v)
end
return divs2
end
function abundant(n, divs)
local sum = 0
for i,v in pairs(divs) do
sum = sum + v
end
return sum > n
end
function semiPerfect(n, divs)
if #divs > 0 then
local h = divs[1]
local t = tail(divs)
if n < h then
return semiPerfect(n, t)
else
return n == h
or semiPerfect(n - h, t)
or semiPerfect(n, t)
end
else
return false
end
end
function sieve(limit)
-- false denotes abundant and not semi-perfect.
-- Only interested in even numbers >= 2
local w = make(limit, false)
local i = 2
while i < limit do
if not w[i] then
local divs = divisors(i)
if not abundant(i, divs) then
w[i] = true
elseif semiPerfect(i, divs) then
local j = i
while j < limit do
w[j] = true
j = j + i
end
end
end
i = i + 1
end
return w
end
function main()
local w = sieve(17000)
local count = 0
local max = 25
print("The first 25 weird numbers:")
local n = 2
while count < max do
if not w[n] then
io.write(n, ' ')
count = count + 1
end
n = n + 2
end
print()
end
main()
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Mathematica / Wolfram Language
ClearAll[WeirdNumberQ, HasSumQ]
HasSumQ[n_Integer, xs_List] := HasSumHelperQ[n, ReverseSort[xs]]
HasSumHelperQ[n_Integer, xs_List] := Module[{h, t},
If[Length[xs] > 0,
h = First[xs];
t = Drop[xs, 1];
If[n < h,
HasSumHelperQ[n, t]
,
n == h \[Or] HasSumHelperQ[n - h, t] \[Or] HasSumHelperQ[n, t]
]
,
False
]
]
WeirdNumberQ[n_Integer] := Module[{divs},
divs = Most[Divisors[n]];
If[Total[divs] > n,
! HasSumQ[n, divs]
,
False
]
]
r = {};
n = 0;
While[
Length[r] < 25,
If[WeirdNumberQ[++n], AppendTo[r, n]]
]
Print[r]
- Output:
{70, 836, 4030, 5830, 7192, 7912, 9272, 10430, 10570, 10792, 10990, 11410, 11690, 12110, 12530, 12670, 13370, 13510, 13790, 13930, 14770, 15610, 15890, 16030, 16310}
Nim
import algorithm, math, strutils
func divisors(n: int): seq[int] =
var smallDivs = @[1]
for i in 2..sqrt(n.toFloat).int:
if n mod i == 0:
let j = n div i
smallDivs.add i
if i != j: result.add j
result.add reversed(smallDivs)
func abundant(n: int; divs: seq[int]): bool {.inline.}=
sum(divs) > n
func semiperfect(n: int; divs: seq[int]): bool =
if divs.len > 0:
let h = divs[0]
let t = divs[1..^1]
result = if n < h: semiperfect(n, t)
else: n == h or semiperfect(n - h, t) or semiperfect(n, t)
func sieve(limit: int): seq[bool] =
# False denotes abundant and not semi-perfect.
# Only interested in even numbers >= 2.
result.setLen(limit)
for i in countup(2, limit - 1, 2):
if result[i]: continue
let divs = divisors(i)
if not abundant(i, divs):
result[i] = true
elif semiperfect(i, divs):
for j in countup(i, limit - 1, i):
result[j] = true
const Max = 25
let w = sieve(17_000)
var list: seq[int]
echo "The first 25 weird numbers are:"
var n = 2
while list.len != Max:
if not w[n]: list.add n
inc n, 2
echo list.join(" ")
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Perl
use strict;
use feature 'say';
use List::Util 'sum';
use POSIX 'floor';
use Algorithm::Combinatorics 'subsets';
use ntheory <is_prime divisors>;
sub abundant {
my($x) = @_;
my $s = sum( my @l = is_prime($x) ? 1 : grep { $x != $_ } divisors($x) );
$s > $x ? ($s, sort { $b <=> $a } @l) : ();
}
my(@weird,$n);
while () {
$n++;
my ($sum, @div) = abundant($n);
next unless $sum; # Weird number must be abundant, skip it if it isn't.
next if $sum / $n > 1.1; # There aren't any weird numbers with a sum:number ratio greater than 1.08 or so.
if ($n >= 10430 and (! int $n%70) and is_prime(int $n/70)) {
# It's weird. All numbers of the form 70 * (a prime 149 or larger) are weird
} else {
my $next;
my $l = shift @div;
my $iter = subsets(\@div);
while (my $s = $iter->next) {
++$next and last if sum(@$s) == $n - $l;
}
next if $next;
}
push @weird, $n;
last if @weird == 25;
}
say "The first 25 weird numbers:\n" . join ' ', @weird;
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Simpler and faster solution:
use 5.010;
use strict;
use ntheory qw(vecsum divisors divisor_sum);
sub is_pseudoperfect {
my ($n, $d, $s, $m) = @_;
$d //= do { my @d = divisors($n); pop(@d); \@d };
$s //= vecsum(@$d);
$m //= $#$d;
return 0 if $m < 0;
while ($d->[$m] > $n) {
$s -= $d->[$m--];
}
return 1 if ($n == $s or $d->[$m] == $n);
is_pseudoperfect($n-$d->[$m], $d, $s-$d->[$m], $m - 1) ||
is_pseudoperfect($n, $d, $s-$d->[$m], $m - 1);
}
sub is_weird {
my ($n) = @_;
divisor_sum($n) > 2*$n and not is_pseudoperfect($n);
}
my @weird;
for (my $k = 1 ; @weird < 25 ; ++$k) {
push(@weird, $k) if is_weird($k);
}
say "The first 25 weird numbers:\n@weird";
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Phix
Sufficiently fast that I un-optimised it a bit to make it easier to follow.
with javascript_semantics function abundant(integer n, sequence divs) return sum(divs) > n end function function semiperfect(integer n, sequence divs) if length(divs)=0 then return false end if integer h = divs[1]; divs = divs[2..$] return n=h or (n>h and semiperfect(n-h, divs)) or semiperfect(n, divs) end function function sieve(integer limit) -- true denotes abundant and not semi-perfect. -- only interested in even numbers >= 2 sequence wierd := repeat(true,limit) for j=6 to limit by 6 do -- eliminate multiples of 3 wierd[j] = false end for for i=2 to limit by 2 do if wierd[i] then sequence divs := factors(i,-1) if not abundant(i,divs) then wierd[i] = false elsif semiperfect(i,divs) then for j=i to limit by i do wierd[j] = false end for end if end if end for return wierd end function --constant MAX = 25, sieve_limit = 16313 constant MAX = 50, sieve_limit = 26533 sequence wierd := sieve(sieve_limit), res = {} for i=2 to sieve_limit by 2 do if wierd[i] then res &= i if length(res)=MAX then exit end if end if end for printf(1,"%s\n",{join(shorten(res,"weird numbers",5,"%d"))})
- Output:
70 836 4030 5830 7192 ... 24710 25130 25690 26110 26530 (50 weird numbers)
Python
Functional
The first 50 seem to take c. 300 ms
'''Weird numbers'''
from itertools import chain, count, islice, repeat
from functools import reduce
from math import sqrt
from time import time
# weirds :: Gen [Int]
def weirds():
'''Non-finite stream of weird numbers.
(Abundant, but not semi-perfect)
OEIS: A006037
'''
def go(n):
ds = descPropDivs(n)
d = sum(ds) - n
return [n] if 0 < d and not hasSum(d, ds) else []
return concatMap(go)(count(1))
# hasSum :: Int -> [Int] -> Bool
def hasSum(n, xs):
'''Does any subset of xs sum to n ?
(Assuming xs to be sorted in descending
order of magnitude)'''
def go(n, xs):
if xs:
h, t = xs[0], xs[1:]
if n < h: # Head too big. Forget it. Tail ?
return go(n, t)
else:
# The head IS the target ?
# Or the tail contains a sum for the
# DIFFERENCE between the head and the target ?
# Or the tail contains some OTHER sum for the target ?
return n == h or go(n - h, t) or go(n, t)
else:
return False
return go(n, xs)
# descPropDivs :: Int -> [Int]
def descPropDivs(n):
'''Descending positive divisors of n,
excluding n itself.'''
root = sqrt(n)
intRoot = int(root)
blnSqr = root == intRoot
lows = [x for x in range(1, 1 + intRoot) if 0 == n % x]
return [
n // x for x in (
lows[1:-1] if blnSqr else lows[1:]
)
] + list(reversed(lows))
# --------------------------TEST---------------------------
# main :: IO ()
def main():
'''Test'''
start = time()
n = 50
xs = take(n)(weirds())
print(
(tabulated('First ' + str(n) + ' weird numbers:\n')(
lambda i: str(1 + i)
)(str)(5)(
index(xs)
)(range(0, n)))
)
print(
'\nApprox computation time: ' +
str(int(1000 * (time() - start))) + ' ms'
)
# -------------------------GENERIC-------------------------
# chunksOf :: Int -> [a] -> [[a]]
def chunksOf(n):
'''A series of lists of length n,
subdividing the contents of xs.
Where the length of xs is not evenly divible,
the final list will be shorter than n.'''
return lambda xs: reduce(
lambda a, i: a + [xs[i:n + i]],
range(0, len(xs), n), []
) if 0 < n else []
# compose (<<<) :: (b -> c) -> (a -> b) -> a -> c
def compose(g):
'''Right to left function composition.'''
return lambda f: lambda x: g(f(x))
# concatMap :: (a -> [b]) -> [a] -> [b]
def concatMap(f):
'''A concatenated list or string over which a function f
has been mapped.
The list monad can be derived by using an (a -> [b])
function which wraps its output in a list (using an
empty list to represent computational failure).
'''
return lambda xs: chain.from_iterable(map(f, xs))
# index (!!) :: [a] -> Int -> a
def index(xs):
'''Item at given (zero-based) index.'''
return lambda n: None if 0 > n else (
xs[n] if (
hasattr(xs, "__getitem__")
) else next(islice(xs, n, None))
)
# paddedMatrix :: a -> [[a]] -> [[a]]
def paddedMatrix(v):
''''A list of rows padded to equal length
(where needed) with instances of the value v.'''
def go(rows):
return paddedRows(
len(max(rows, key=len))
)(v)(rows)
return lambda rows: go(rows) if rows else []
# paddedRows :: Int -> a -> [[a]] -[[a]]
def paddedRows(n):
'''A list of rows padded (but never truncated)
to length n with copies of value v.'''
def go(v, xs):
def pad(x):
d = n - len(x)
return (x + list(repeat(v, d))) if 0 < d else x
return list(map(pad, xs))
return lambda v: lambda xs: go(v, xs) if xs else []
# showColumns :: Int -> [String] -> String
def showColumns(n):
'''A column-wrapped string
derived from a list of rows.'''
def go(xs):
def fit(col):
w = len(max(col, key=len))
def pad(x):
return x.ljust(4 + w, ' ')
return ''.join(map(pad, col))
q, r = divmod(len(xs), n)
return unlines(map(
fit,
transpose(paddedMatrix('')(
chunksOf(q + int(bool(r)))(
xs
)
))
))
return lambda xs: go(xs)
# succ :: Enum a => a -> a
def succ(x):
'''The successor of a value. For numeric types, (1 +).'''
return 1 + x if isinstance(x, int) else (
chr(1 + ord(x))
)
# tabulated :: String -> (a -> String) ->
# (b -> String) ->
# Int ->
# (a -> b) -> [a] -> String
def tabulated(s):
'''Heading -> x display function -> fx display function ->
number of columns -> f -> value list -> tabular string.'''
def go(xShow, fxShow, intCols, f, xs):
w = max(map(compose(len)(xShow), xs))
return s + '\n' + showColumns(intCols)([
xShow(x).rjust(w, ' ') + ' -> ' + fxShow(f(x)) for x in xs
])
return lambda xShow: lambda fxShow: lambda nCols: (
lambda f: lambda xs: go(
xShow, fxShow, nCols, f, xs
)
)
# take :: Int -> [a] -> [a]
# take :: Int -> String -> String
def take(n):
'''The prefix of xs of length n,
or xs itself if n > length xs.'''
return lambda xs: (
xs[0:n]
if isinstance(xs, list)
else list(islice(xs, n))
)
# transpose :: Matrix a -> Matrix a
def transpose(m):
'''The rows and columns of the argument transposed.
(The matrix containers and rows can be lists or tuples).'''
if m:
inner = type(m[0])
z = zip(*m)
return (type(m))(
map(inner, z) if tuple != inner else z
)
else:
return m
# unlines :: [String] -> String
def unlines(xs):
'''A single string derived by the intercalation
of a list of strings with the newline character.'''
return '\n'.join(xs)
# until :: (a -> Bool) -> (a -> a) -> a -> a
def until(p):
'''The result of repeatedly applying f until p holds.
The initial seed value is x.'''
def go(f, x):
v = x
while not p(v):
v = f(v)
return v
return lambda f: lambda x: go(f, x)
# MAIN ----------------------------------------------------
if __name__ == '__main__':
main()
- Output:
First 50 weird numbers: 1 -> 70 11 -> 10990 21 -> 14770 31 -> 18410 41 -> 22190 2 -> 836 12 -> 11410 22 -> 15610 32 -> 18830 42 -> 23170 3 -> 4030 13 -> 11690 23 -> 15890 33 -> 18970 43 -> 23590 4 -> 5830 14 -> 12110 24 -> 16030 34 -> 19390 44 -> 24290 5 -> 7192 15 -> 12530 25 -> 16310 35 -> 19670 45 -> 24430 6 -> 7912 16 -> 12670 26 -> 16730 36 -> 19810 46 -> 24710 7 -> 9272 17 -> 13370 27 -> 16870 37 -> 20510 47 -> 25130 8 -> 10430 18 -> 13510 28 -> 17272 38 -> 21490 48 -> 25690 9 -> 10570 19 -> 13790 29 -> 17570 39 -> 21770 49 -> 26110 10 -> 10792 20 -> 13930 30 -> 17990 40 -> 21910 50 -> 26530 Approx computation time: 284 ms
Alternate
Edit: I experimented a bit and created this second version. It's faster than the original, takes about 0.16 s for first 50.
# anySum :: Int -> [Int] -> [Int] from time import time start = time() from itertools import combinations from math import prod primitivesp_nos = set(); primes = [] def main(): weird_nos = [] n = 1766 x = 1 while n > 0: if isweird(x) == 1: weird_nos.append(x); n = n - 1 x = x + 1 print("First", len(weird_nos), "weird nos:\n", weird_nos) def isweird(n): global primes; global primitivesp_nos; x = []; a = n for i in primes: while a % i == 0: x.append(i); a = a//i if i * i > a: break if a > 1: x.append(a) if x == [n]: primes.append(n); return 0 sum_fctrs = 1; fctrs = set(x) for i in fctrs: sum_fctrs = sum_fctrs*(i**(x.count(i)+1)-1)//(i-1) difference = sum_fctrs - 2 * n if difference <= 0: if difference == 0: primitivesp_nos.add(n) return 0 for i in range(2, len(x)): for j in combinations(x, i): product = prod(j) if product not in primitivesp_nos: fctrs.add(product) else: return 0 x = 1; fctrs = {1}.union(set(i for i in fctrs if i <= difference)) ns = n - (difference + n - sum(fctrs)) if ns < 0: return 1 for d in fctrs: x |= x << d if not x >> ns & 1: return 1 else: primitivesp_nos.add(n); return 0 main() end = time() print("Execution time: ", round(end - start, 2), "s") -- -> [100,96]
Quackery
properdivisors
is defined at Proper divisors#Quackery.
[ stack ] is target ( --> s )
[ stack ] is success ( --> s )
[ stack ] is makeable ( --> s )
[ bit makeable take
2dup & 0 !=
dip [ | makeable put ] ] is made ( n --> b )
[ ' [ 0 ] swap
dup target put
properdivisors
0 over witheach +
target share > not iff
[ target release
2drop false ] done
true success put
0 makeable put
witheach
[ over witheach
[ over dip
[ +
dup target share = iff
[ false success replace
drop conclude ] done
dup target share < iff
[ dup made not iff
join else drop ]
else drop ] ]
success share not if conclude
drop ]
drop
target release
makeable release
success take ] is weird ( n --> b )
[] 0
[ 1+
dup weird if
[ tuck join swap ]
over size 25 = until ]
drop
echo
- Output:
[ 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 ]
Racket
#lang racket
(require math/number-theory)
(define (abundant? n proper-divisors)
(> (apply + proper-divisors) n))
(define (semi-perfect? n proper-divisors)
(let recur ((ds proper-divisors) (n n))
(or (zero? n)
(and (positive? n)
(pair? ds)
(or (recur (cdr ds) n)
(recur (cdr ds) (- n (car ds))))))))
(define (weird? n)
(let ((proper-divisors (drop-right (divisors n) 1))) ;; divisors includes n
(and (abundant? n proper-divisors) (not (semi-perfect? n proper-divisors)))))
(module+ main
(let recur ((i 0) (n 1) (acc null))
(cond [(= i 25) (reverse acc)]
[(weird? n) (recur (add1 i) (add1 n) (cons n acc))]
[else (recur i (add1 n) acc)])))
(module+ test
(require rackunit)
(check-true (weird? 70))
(check-false (weird? 12)))
- Output:
'(70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310)
Raku
(formerly Perl 6)
sub abundant (\x) {
my @l = x.is-prime ?? 1 !! flat
1, (2 .. x.sqrt.floor).map: -> \d {
my \y = x div d;
next if y * d !== x;
d !== y ?? (d, y) !! d
};
(my $s = @l.sum) > x ?? ($s, |@l.sort(-*)) !! ();
}
my @weird = (2, 4, {|($_ + 4, $_ + 6)} ... *).map: -> $n {
my ($sum, @div) = $n.&abundant;
next unless $sum; # Weird number must be abundant, skip it if it isn't.
next if $sum / $n > 1.1; # There aren't any weird numbers with a sum:number ratio greater than 1.08 or so.
if $n >= 10430 and ($n %% 70) and ($n div 70).is-prime {
# It's weird. All numbers of the form 70 * (a prime 149 or larger) are weird
} else {
my $next;
my $l = @div.shift;
++$next and last if $_.sum == $n - $l for @div.combinations;
next if $next;
}
$n
}
put "The first 25 weird numbers:\n", @weird[^25];
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
REXX
vanilla version
/*REXX program finds and displays N weird numbers in a vertical format (with index).*/
parse arg n cols . /*obtain optional arguments from the CL*/
if n=='' | n=="," then n= 25 /*Not specified? Then use the default.*/
if cols=='' | cols=="," then cols= 10 /* " " " " " " */
w= 10 /*width of a number in any column. */
if cols>0 then say ' index │'center(' weird numbers', 1 + cols*(w+1) )
if cols>0 then say '───────┼'center("" , 1 + cols*(w+1), '─')
idx= 1; $= /*index for the output list; $: 1 line*/
weirds= 0 /*the count of weird numbers (so far).*/
do j=2 by 2 until weirds==n /*examine even integers 'til have 'nuff*/
if \weird(j) then iterate /*Not a weird number? Then skip it. */
weirds= weirds + 1 /*bump the count of weird numbers. */
c= commas(j) /*maybe add commas to the number. */
$= $ right(c, max(w, length(c) ) ) /*add a nice prime ──► list, allow big#*/
if weirds//cols\==0 then iterate /*have we populated a line of output? */
say center(idx, 7)'│' substr($, 2); $= /*display what we have so far (cols). */
idx= idx + cols /*bump the index count for the output*/
end /*j*/
if $\=='' then say center(idx, 7)"│" substr($, 2) /*possible display residual output.*/
if cols>0 then say '───────┴'center("" , 1 + cols*(w+1), '─')
say
say 'Found ' commas(weirds) ' weird numbers'
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do ic=length(_)-3 to 1 by -3; _=insert(',', _, ic); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
DaS: procedure; parse arg x 1 z 1,b; a= 1 /*get X,Z,B (the 1st arg); init A list.*/
r= 0; q= 1 /* [↓] ══integer square root══ ___ */
do while q<=z; q=q*4; end /*R: an integer which will be √ X */
do while q>1; q=q%4; _= z-r-q; r=r%2; if _>=0 then do; z=_; r=r+q; end
end /*while q>1*/ /* [↑] compute the integer sqrt of X.*/
sig= a /*initialize the sigma so far. ___ */
do j=2 to r - (r*r==x) /*divide by some integers up to √ X */
if x//j==0 then do; a=a j; b= x%j b /*if ÷, add both divisors to α and ß. */
sig= sig +j +x%j /*bump the sigma (the sum of divisors).*/
end
end /*j*/ /* [↑] % is the REXX integer division*/
/* [↓] adjust for a square. ___*/
if j*j==x then return sig+j a j b /*Was X a square? If so, add √ X */
return sig a b /*return the divisors (both lists). */
/*──────────────────────────────────────────────────────────────────────────────────────*/
weird: procedure; parse arg x . /*obtain a # to be tested for weirdness*/
if x<70 | x//3==0 then return 0 /*test if X is too low or multiple of 3*/
parse value DaS(x) with sigma divs /*obtain sigma and the proper divisors.*/
if sigma<=x then return 0 /*X isn't abundant (sigma too small).*/
#= words(divs) /*count the number of divisors for X. */
if #<3 then return 0 /*Not enough divisors? " " */
if #>15 then return 0 /*number of divs > 15? It's not weird.*/
a.= /*initialize the A. stemmed array.*/
do i=1 for #; _= word(divs, i) /*obtain one of the divisors of X. */
@.i= _; a._= . /*assign proper divs──►@ array; also id*/
end /*i*/
df= sigma - x /*calculate difference between Σ and X.*/
if a.df==. then return 0 /*Any divisor is equal to DF? Not weird*/
c= 0 /*zero combo counter; calc. power of 2.*/
do p=1 for 2**#-2; c= c + 1 /*convert P──►binary with leading zeros*/
yy.c= strip( x2b( d2x(p) ), 'L', 0) /*store this particular combination. */
end /*p*/
/* [↓] decreasing partitions is faster*/
do part=c by -1 for c; s= 0 /*test of a partition add to the arg X.*/
_= yy.part; L= length(_) /*obtain one method of partitioning. */
do cp=L by -1 for L /*obtain a sum of a partition. */
if substr(_,cp,1) then do; s= s + @.cp /*1 bit? Then add ──►S*/
if s==x then return 0 /*Sum equal? Not weird*/
if s==df then return 0 /*Sum = DF? " " */
if s>x then iterate /*Sum too big? Try next*/
end
end /*cp*/
end /*part*/; return 1 /*no sum equal to X, so X is weird.*/
- output when using the default inputs:
index │ weird numbers ───────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ 70 836 4,030 5,830 7,192 7,912 9,272 10,430 10,570 10,792 11 │ 10,990 11,410 11,690 12,110 12,530 12,670 13,370 13,510 13,790 13,930 21 │ 14,770 15,610 15,890 16,030 16,310 ───────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────── Found 25 weird numbers
optimized version
This REXX program was optimized by finding primitive weird numbers (as in the 1st REXX version), and multiplying
them by prime numbers ≥ sigma(primitive weird numbers) to find higher weird numbers.
This version is about 300% faster than the 1st REXX version for larger amount of numbers.
/*REXX program finds and displays N weird numbers in a vertical format (with index).*/
parse arg n cols . /*obtain optional arguments from the CL*/
if n=='' | n=="," then n= 400 /*Not specified? Then use the default.*/
if cols=='' | cols=="," then cols= 10 /* " " " " " " */
w= 10 /*width of a number in any column. */
call genP /*generate primes just past Hp. */
if cols>0 then say ' index │'center(' weird numbers', 1 + cols*(w+1) )
if cols>0 then say '───────┼'center("" , 1 + cols*(w+1), '─')
weirds= 0; !!.= 0 /*the count of weird numbers (so far).*/
idx= 1; $= /*index for the output list; $: 1 line*/
do j=2 by 2 until weirds==n /*examine even integers 'til have 'nuff*/
if \weird(j) then iterate /*Not a weird number? Then skip it. */
weirds= weirds + 1 /*bump the count of weird numbers. */
do a=1 for # until _>hp; if @.a<sigma+j then iterate; _= j*@.a; !!._= 1
end /*a*/
c= commas(j) /*maybe add commas to the number. */
$= $ right(c, max(w, length(c) ) ) /*add a nice prime ──► list, allow big#*/
if weirds//cols\==0 then iterate /*have we populated a line of output? */
say center(idx, 7)'│' substr($, 2); $= /*display what we have so far (cols). */
idx= idx + cols /*bump the index count for the output*/
end /*j*/
if $\=='' then say center(idx, 7)"│" substr($, 2) /*possible display residual output.*/
if cols>0 then say '───────┴'center("" , 1 + cols*(w+1), '─')
say
say 'Found ' commas(weirds) ' weird numbers'
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do ic=length(_)-3 to 1 by -3; _=insert(',', _, ic); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
DaS: procedure; parse arg x 1 z 1,b; a= 1 /*get X,Z,B (the 1st arg); init A list.*/
r= 0; q= 1 /* [↓] ══integer square root══ ___ */
do while q<=z; q=q*4; end /*R: an integer which will be √ X */
do while q>1; q=q%4; _= z-r-q; r=r%2; if _>=0 then do; z=_; r=r+q; end
end /*while q>1*/ /* [↑] compute the integer sqrt of X.*/
sig = a /*initialize the sigma so far. ___ */
do j=2 to r - (r*r==x) /*divide by some integers up to √ X */
if x//j==0 then do; a=a j; b= x%j b /*if ÷, add both divisors to α & ß. */
sig= sig +j +x%j /*bump the sigma (the sum of Pdivisors)*/
end
end /*j*/ /* [↑] % is the REXX integer division*/
/* [↓] adjust for a square. ___*/
if j*j==x then return sig+j a j b /*Was X a square? If so, add √ X */
return sig a b /*return the divisors (both lists). */
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP: hp= 1000 * n /*high Prime limit; define 2 low primes*/
@.1=2; @.2=3; @.3=5; @.4=7; @.5=11 /*define some low primes. */
#=5; s.#= @.# **2 /*number of primes so far; prime². */
/* [↓] generate more primes ≤ high.*/
do j=@.#+2 by 2 for max(0, hp%2-@.#%2-1) /*find odd primes from here on. */
parse var j '' -1 _; if _==5 then iterate /*J divisible by 5? (right dig)*/
if j// 3==0 then iterate /*" " " 3? */
if j// 7==0 then iterate /*" " " 7? */
/* [↑] the above five lines saves time*/
do k=5 while s.k<=j /* [↓] divide by the known odd primes.*/
if j // @.k == 0 then iterate j /*Is J ÷ X? Then not prime. ___ */
end /*k*/ /* [↑] only process numbers ≤ √ J */
#= #+1; @.#= j; s.#= j*j /*bump # of Ps; assign next P; P²; P# */
end /*j*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
weird: procedure expose !!. sigma; parse arg x /*obtain a # to be tested for weirdness*/
if x<70 | x//3==0 then return 0 /*test if X is too low or multiple of 3*/
if !!.x then return 1 /*Is this a prime*previous #? Found one*/
parse value DaS(x) with sigma divs /*obtain sigma and the proper divisors.*/
if sigma<=x then return 0 /*X isn't abundant (sigma too small).*/
#= words(divs) /*count the number of divisors for X. */
if #<3 then return 0 /*Not enough divisors? " " */
if #>15 then return 0 /*number of divs > 15? It's not weird.*/
a.= /*initialize the A. stemmed array.*/
do i=1 for #; _= word(divs, i) /*obtain one of the divisors of X. */
@.i= _; a._= . /*assign proper divs──►@ array; also id*/
end /*i*/
df= sigma - x /*calculate difference between Σ and X.*/
if a.df==. then return 0 /*Any divisor is equal to DF? Not weird*/
c= 0; u= 2**# /*zero combo counter; calc. power of 2.*/
do p=1 for u-2; c= c + 1 /*convert P──►binary with leading zeros*/
yy.c= strip( x2b( d2x(p) ), 'L', 0) /*store this particular combination. */
end /*p*/
/* [↓] decreasing partitions is faster*/
do part=c by -1 for c; s= 0 /*test of a partition add to the arg X.*/
_= yy.part; L= length(_) /*obtain one method of partitioning. */
do cp=L by -1 for L /*obtain a sum of a partition. */
if substr(_,cp,1) then do; s= s + @.cp /*1 bit? Then add ──►S*/
if s==x then return 0 /*Sum equal? Not weird*/
if s==df then return 0 /*Sum = DF? " " */
if s>x then iterate /*Sum too big? Try next*/
end
end /*cp*/
end /*part*/
return 1 /*no sum equal to X, so X is weird.*/
- output when using the default inputs:
index │ weird numbers ───────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ 70 836 4,030 5,830 7,192 7,912 9,272 10,430 10,570 10,792 11 │ 10,990 11,410 11,690 12,110 12,530 12,670 13,370 13,510 13,790 13,930 21 │ 14,770 15,610 15,890 16,030 16,310 16,730 16,870 17,272 17,570 17,990 31 │ 18,410 18,830 18,970 19,390 19,670 19,810 20,510 21,490 21,770 21,910 41 │ 22,190 23,170 23,590 24,290 24,430 24,710 25,130 25,690 26,110 26,530 51 │ 26,810 27,230 27,790 28,070 28,630 29,330 29,470 30,170 30,310 30,730 61 │ 31,010 31,430 31,990 32,270 32,410 32,690 33,530 34,090 34,370 34,930 71 │ 35,210 35,630 36,470 36,610 37,870 38,290 38,990 39,410 39,830 39,970 81 │ 40,390 41,090 41,510 41,930 42,070 42,490 42,910 43,190 43,330 44,170 91 │ 44,870 45,010 45,290 45,710 46,130 46,270 47,110 47,390 47,810 48,370 101 │ 49,070 49,630 50,330 50,890 51,310 51,730 52,010 52,570 52,990 53,270 111 │ 53,830 54,110 55,090 55,790 56,630 56,770 57,470 57,610 57,890 58,030 121 │ 58,730 59,710 59,990 60,130 60,410 61,390 61,670 61,810 62,090 63,490 131 │ 63,770 64,330 65,030 65,590 65,870 66,290 66,710 67,690 67,970 68,390 141 │ 68,810 69,370 69,790 70,630 70,910 71,330 71,470 72,170 72,310 72,730 151 │ 73,430 73,570 74,270 74,410 74,830 76,090 76,370 76,510 76,790 77,210 161 │ 77,630 78,190 78,610 79,030 80,570 80,710 81,410 81,970 82,670 83,090 171 │ 83,510 84,070 84,910 85,190 85,610 86,030 86,170 86,590 87,430 88,130 181 │ 89,390 89,530 89,810 90,230 90,370 90,790 91,070 91,210 91,490 92,330 191 │ 92,470 92,890 95,270 95,690 96,110 96,670 97,930 98,630 99,610 99,890 201 │ 100,030 100,310 100,730 101,290 101,570 101,710 102,130 102,970 103,670 103,810 211 │ 104,090 104,230 104,510 104,930 105,770 106,610 107,170 108,010 108,430 108,710 221 │ 109,130 109,690 109,970 110,530 110,810 111,790 112,070 112,490 112,630 112,910 231 │ 113,330 113,470 113,890 114,590 115,990 116,410 116,690 116,830 118,510 118,790 241 │ 118,930 119,630 120,470 120,610 121,310 121,870 122,290 122,710 123,130 124,390 251 │ 124,810 125,090 125,230 126,070 126,770 127,610 128,170 129,290 130,270 130,690 261 │ 130,970 131,110 131,390 131,530 132,230 133,070 133,490 133,910 135,170 135,310 271 │ 136,430 136,570 138,110 138,530 139,090 139,510 139,790 139,930 140,210 140,770 281 │ 141,190 141,890 142,030 142,730 143,710 144,410 144,830 145,670 145,810 146,090 291 │ 146,230 146,930 147,770 147,910 149,030 149,170 149,590 149,870 150,010 150,710 301 │ 151,270 152,530 154,210 154,490 154,910 155,470 156,590 156,730 157,010 157,570 311 │ 158,690 158,830 159,110 159,670 160,090 160,510 160,790 161,630 161,770 163,310 321 │ 163,730 163,870 164,290 164,570 164,990 165,970 166,390 166,670 166,810 167,230 331 │ 167,510 167,930 168,770 169,190 169,610 170,590 170,870 171,290 172,130 172,690 341 │ 173,110 173,390 175,210 176,470 177,170 177,730 178,010 178,430 178,570 178,990 351 │ 180,530 181,370 181,510 182,630 183,190 183,470 184,310 185,290 185,990 186,130 361 │ 186,410 186,970 187,390 187,810 188,090 188,230 188,510 188,930 189,490 189,770 371 │ 189,910 190,330 191,030 191,170 191,870 192,430 192,710 193,690 194,390 195,230 381 │ 195,370 195,790 196,070 196,210 197,330 198,310 198,590 199,010 199,570 199,990 391 │ 200,270 201,530 202,090 202,790 203,210 203,630 204,190 204,890 205,730 206,710 ───────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────── Found 400 weird numbers
RPL
@ ------- Sieve manager --------- « 39 DUP STWS / CEIL « # 0b » 'x' 1 4 ROLL 1 SEQ » 'CSV' STO « SWAP 1 - 39 MOD LASTARG / IP 1 + ROT SWAP DUP2 GET 2 5 ROLL ^ R→B OR PUT » 'SSV' STO « SWAP 1 - 39 MOD LASTARG / IP 1 + ROT SWAP GET 2 ROT ^ R→B AND # 0b ≠ » 'SVS?' STO @ --- Weird numbers generator --- « → beg « DUP ←divs DOT ←n CASE DUP2 == THEN DROP2 1 SF END > THEN beg ←ndiv FOR j IF DUP j GET THEN DUP j 0 PUT j FINDSUM IF 1 FS? THEN ←ndiv 'j' STO END END NEXT END END DROP » » 'FINDSUM' STO @ ( [ bits ] 1st_bit → ) flag 1 is set if →n is semi-perfect « IF DUP 'Sieve' SVS? THEN DROP 0 ELSE DUP DIVIS SORT REVLIST TAIL AXL DUP SIZE EVAL → ←n ←divs ←ndiv « ←divs DUP 1 CON IF SWAP OVER DOT ←n < THEN DROP 0 ELSE 1 CF 1 FINDSUM 1 FC? IF DUP NOT ←n SieveMax 2 / ≤ AND THEN 2 SieveMax ←n / FLOOR FOR j j ←n * 'Sieve' SSV NEXT END END » END » 'WEIRD?' STO « 2 → n « { } 20000 'SieveMax' ▶ CSV 'Sieve' STO WHILE DUP SIZE 25 < REPEAT IF n WEIRD? THEN n + END 'n' 2 STO+ END » » 'TASK' STO
- Output:
1: { 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 }
It took 269 seconds for an emulator running on an iPhone XR => probably 3 hours for an HP-50 pocket calculator.
Ruby
def divisors(n)
divs = [1]
divs2 = []
i = 2
while i * i <= n
if n % i == 0 then
j = (n / i).to_i
divs.append(i)
if i != j then
divs2.append(j)
end
end
i = i + 1
end
divs2 += divs.reverse
return divs2
end
def abundant(n, divs)
return divs.sum > n
end
def semiperfect(n, divs)
if divs.length > 0 then
h = divs[0]
t = divs[1..-1]
if n < h then
return semiperfect(n, t)
else
return n == h || semiperfect(n - h, t) || semiperfect(n, t)
end
else
return false
end
end
def sieve(limit)
w = Array.new(limit, false)
i = 2
while i < limit
if not w[i] then
divs = divisors(i)
if not abundant(i, divs) then
w[i] = true
elsif semiperfect(i, divs) then
j = i
while j < limit
w[j] = true
j = j + i
end
end
end
i = i + 2
end
return w
end
def main
w = sieve(17000)
count = 0
max = 25
print "The first %d weird numbers:\n" % [max]
n = 2
while count < max
if not w[n] then
print n, " "
count = count + 1
end
n = n + 2
end
print "\n"
end
main()
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Sidef
func is_pseudoperfect(n, d = n.divisors.first(-1), s = d.sum, m = d.end) {
return false if (m < 0)
while (d[m] > n) {
s -= d[m--]
}
return true if (n == s)
return true if (d[m] == n)
__FUNC__(n-d[m], d, s-d[m], m-1) || __FUNC__(n, d, s-d[m], m-1)
}
func is_weird(n) {
(n.sigma > 2*n) && !is_pseudoperfect(n)
}
var w = (1..Inf -> lazy.grep(is_weird).first(25))
say "The first 25 weird numbers:\n#{w.join(' ')}"
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Visual Basic .NET
Performance is now on par with the python version, (but not quite up the Go version's performance), I applied what I could after reading the comments made by Hout on the discussion page.
This program is similar to the structure of the Go example. I found a couple of tweaks here and there to help with performance. For example, the divisors list is built on a single array instead of joining two, and it calculates the sum while creating the divisors list. The divisors list is headed by the difference between "n" and the sum of the divisors. The semiperfect() function checks for equality first (rather than chopping the head from the tail list first) to save a little more time. And of course, the parallel execution.
A new feature is that one can calculate weird numbers up to any reasonable number, just enter a command line parameter of more than zero. Another new feature is calculating weird numbers continuously until a key is pressed (like the spigot algorithm from the Pi task) - to do so, enter a command line parameter of less than 1.
This has no sieve cache, as one must "know" beforehand what number to cache up to, (for best results). Since there is no cache (runs slower), I added the parallel execution to make it run faster.
I haven't let it run long enough to see how high it can get before crashing, I suspect it should happen once the weird number being tested is around Int32.MaxValue (2,147,483,647). But long before that it will slow down quite a bit. It takes around 17 minutes to get to the 10,732nd weird number, which is the first over 7 million (7,000,210).
Module Module1
Dim resu As New List(Of Integer)
Function TestAbundant(n As Integer, ByRef divs As List(Of Integer)) As Boolean
divs = New List(Of Integer)
Dim sum As Integer = -n : For i As Integer = Math.Sqrt(n) To 1 Step -1
If n Mod i = 0 Then divs.Add(i) : Dim j As Integer = n / i : divs.Insert(0, j) : sum += i + j
Next : divs(0) = sum - divs(0) : Return divs(0) > 0
End Function
Function subList(src As List(Of Integer), Optional first As Integer = Integer.MinValue) As List(Of Integer)
subList = src.ToList : subList.RemoveAt(1)
End Function
Function semiperfect(divs As List(Of Integer)) As Boolean
If divs.Count < 2 Then Return False
Select Case divs.First.CompareTo(divs(1))
Case 0 : Return True
Case -1 : Return semiperfect(subList(divs))
Case 1 : Dim t As List(Of Integer) = subList(divs) : t(0) -= divs(1)
If semiperfect(t) Then Return True Else t(0) = divs.First : Return semiperfect(t)
End Select : Return False ' execution can't get here, just for compiler warning
End Function
Function Since(et As TimeSpan) As String ' big ugly routine to prettify the elasped time
If et > New TimeSpan(2000000) Then
Dim s As String = " " & et.ToString(), p As Integer = s.IndexOf(":"), q As Integer = s.IndexOf(".")
If q < p Then s = s.Insert(q, "Days") : s = s.Replace("Days.", "Days, ")
p = s.IndexOf(":") : s = s.Insert(p, "h") : s = s.Replace("h:", "h ")
p = s.IndexOf(":") : s = s.Insert(p, "m") : s = s.Replace("m:", "m ")
s = s.Replace(" 0", " ").Replace(" 0h", " ").Replace(" 0m", " ") & "s"
Return s.TrimStart()
Else
If et > New TimeSpan(1500) Then
Return et.TotalMilliseconds.ToString() & "ms"
Else
If et > New TimeSpan(15) Then
Return (et.TotalMilliseconds * 1000.0).ToString() & "µs"
Else
Return (et.TotalMilliseconds * 1000000.0).ToString() & "ns"
End If
End If
End If
End Function
Sub Main(args As String())
Dim sw As New Stopwatch, st As Integer = 2, stp As Integer = 1020, count As Integer = 0
Dim max As Integer = 25, halted As Boolean = False
If args.Length > 0 Then _
Dim t As Integer = Integer.MaxValue : If Integer.TryParse(args(0), t) Then max = If(t > 0, t, Integer.MaxValue)
If max = Integer.MaxValue Then
Console.WriteLine("Calculating weird numbers, press a key to halt.")
stp *= 10
Else
Console.WriteLine("The first {0} weird numbers:", max)
End If
If max < 25 Then stp = 140
sw.Start()
Do : Parallel.ForEach(Enumerable.Range(st, stp),
Sub(n)
Dim divs As List(Of Integer) = Nothing
If TestAbundant(n, divs) AndAlso Not semiperfect(divs) Then
SyncLock resu : resu.Add(n) : End SyncLock
End If
End Sub)
If resu.Count > 0 Then
resu.Sort()
If count + resu.Count > max Then
resu = resu.Take(max - count).ToList
End If
Console.Write(String.Join(" ", resu) & " ")
count += resu.Count : resu.Clear()
End If
If Console.KeyAvailable Then Console.ReadKey() : halted = True : Exit Do
st += stp
Loop Until count >= max
sw.Stop()
If max < Integer.MaxValue Then
Console.WriteLine(vbLf & "Computation time was {0}.", Since(sw.Elapsed))
If halted Then Console.WriteLine("Halted at number {0}.", count)
Else
Console.WriteLine(vbLf & "Computation time was {0} for the first {1} weird numbers.", Since(sw.Elapsed), count)
End If
End Sub
End Module
- Output:
Without any command line parameters:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 Computation time was 37.4931ms.
With command line arguments = 50
The first 50 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 16730 16870 17272 17570 17990 18410 18830 18970 19390 19670 19810 20510 21490 21770 21910 22190 23170 23590 24290 24430 24710 25130 25690 26110 26530 Computation time was 47.6589ms.
With command line arguments = 0
Calculating weird numbers, press a key to halt. 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 16730 16870 17272 17570 17990 18410 18830 18970 19390 19670 19810 20510 21490 21770 21910 22190 23170 23590 24290 24430 24710 25130 25690 26110 26530 26810 27230 27790 28070 28630 29330 29470 30170 30310 30730 31010 31430 31990 32270 32410 32690 33530 34090 34370 34930 35210 35630 36470 36610 37870 38290 38990 39410 39830 39970 40390 41090 41510 41930 42070 42490 42910 43190 43330 44170 44870 45010 45290 45356 45710 46130 46270 47110 47390 47810 48370 49070 49630 50330 50890 51310 51730 52010 52570 52990 53270 53830 54110 55090 55790 56630 56770 57470 57610 57890 58030 58730 59710 59990 60130 60410 61390 61670 61810 62090 63490 63770 64330 65030 65590 65870 66290 66710 67690 67970 68390 68810 69370 69790 70630 70910 71330 71470 72170 72310 72730 73430 73570 73616 74270 74410 74830 76090 76370 76510 76790 77210 77630 78190 78610 79030 80570 80710 81410 81970 82670 83090 83312 83510 84070 84910 85190 85610 86030 86170 86590 87430 88130 89390 89530 89810 90230 90370 90790 91070 91210 91388 91490 92330 92470 92890 95270 95690 96110 96670 97930 98630 99610 99890 100030 100310 100730 101290 101570 101710 102130 102970 103670 103810 104090 104230 104510 104930 105770 106610 107170 108010 108430 108710 109130 109690 109970 110530 110810 111790 112070 112490 112630 112910 113072 113330 113470 113890 114590 115990 116410 116690 116830 118510 118790 118930 119630 120470 120610 121310 121870 122290 122710 123130 124390 124810 125090 125230 126070 126770 127610 128170 129290 130270 130690 130970 131110 131390 131530 132230 133070 133490 133910 135170 135310 136430 136570 138110 138530 139090 139510 139790 139930 140210 140770 Computation time was 153.3649ms for the first 285 weird numbers.
Tail-end of a longer session:
6981310 6983108 6983270 6983690 6985090 6985510 6986630 6987190 6987610 6988030 6988310 6988730 6990130 6990970 6991390 6991468 6991670 6992930 6993070 6993490 6994610 6995030 6996484 6997270 6997970 6998110 6999230 6999370 7000210 7001330 7003010 7003172 7003430 7003990 7004830 7007210 7007630 7008890 7009030 Computation time was 17m 9.0062776s for the first 10742 weird numbers.
V (Vlang)
fn divisors(n int) []int {
mut divs := [1]
mut divs2 := []int{}
for i := 2; i*i <= n; i++ {
if n%i == 0 {
j := n / i
divs << i
if i != j {
divs2 << j
}
}
}
for i := divs.len - 1; i >= 0; i-- {
divs2 << divs[i]
}
return divs2
}
fn abundant(n int, divs []int) bool {
mut sum := 0
for div in divs {
sum += div
}
return sum > n
}
fn semiperfect(n int, divs []int) bool {
le := divs.len
if le > 0 {
h := divs[0]
t := divs[1..]
if n < h {
return semiperfect(n, t)
} else {
return n == h || semiperfect(n-h, t) || semiperfect(n, t)
}
} else {
return false
}
}
fn sieve(limit int) []bool {
// false denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
mut w := []bool{len: limit}
for i := 2; i < limit; i += 2 {
if w[i] {
continue
}
divs := divisors(i)
if !abundant(i, divs) {
w[i] = true
} else if semiperfect(i, divs) {
for j := i; j < limit; j += i {
w[j] = true
}
}
}
return w
}
fn main() {
w := sieve(17000)
mut count := 0
max := 25
println("The first 25 weird numbers are:")
for n := 2; count < max; n += 2 {
if !w[n] {
print("$n ")
count++
}
}
println('')
}
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
Wren
import "./math" for Int, Nums
import "./iterate" for Stepped
var semiperfect // recursive
semiperfect = Fn.new { |n, divs|
var le = divs.count
if (le == 0) return false
var h = divs[0]
if (n == h) return true
if (le == 1) return false
var t = divs[1..-1]
if (n < h) return semiperfect.call(n, t)
return semiperfect.call(n-h, t) || semiperfect.call(n, t)
}
var sieve = Fn.new { |limit|
// 'false' denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
var w = List.filled(limit, false)
for (j in Stepped.new(6...limit, 6)) w[j] = true // eliminate multiples of 3
for (i in Stepped.new(2...limit, 2)) {
if (!w[i]) {
var divs = Int.properDivisors(i)
var sum = Nums.sum(divs)
if (sum <= i) {
w[i] = true
} else if (semiperfect.call(sum-i, divs)) {
for (j in Stepped.new(i...limit, i)) w[j] = true
}
}
}
return w
}
var start = System.clock
var limit = 16313
var w = sieve.call(limit)
var count = 0
var max = 25
System.print("The first 25 weird numbers are:")
var n = 2
while (count < max) {
if (!w[n]) {
System.write("%(n) ")
count = count + 1
}
n = n + 2
}
System.print()
System.print("\nTook %(((System.clock-start)*1000).round) milliseconds")
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310 Took 144 milliseconds
XPL0
This runs on a Raspberry Pi. MAlloc in other versions of XPL0 work differently. Takes about 1.3 seconds.
def SizeOfInt = 4;
def \IntA\ Ptr, Size;
int Array(2);
func Divisors(N); \Returns a list of proper divisors for N
int N;
int Divs, Divs2, Out;
int I, J, C1, C2;
[C1:= 0; C2:= 0;
Divs:= MAlloc(N * SizeOfInt / 2);
Divs2:= MAlloc(N * SizeOfInt / 2);
Divs(C1):= 1; C1:= C1+1;
I:= 2;
while I*I <= N do
[if rem(N/I) = 0 then
[J:= N/I;
Divs(C1):= I; C1:= C1+1;
if I # J then
[Divs2(C2):= J; C2:= C2+1];
];
I:= I+1;
];
Out:= MAlloc((C1+C2) * SizeOfInt);
for I:= 0 to C2-1 do
Out(I):= Divs2(I);
for I:= 0 to C1-1 do
Out(C2+I):= Divs(C1-I-1);
Array(Ptr):= Out;
Array(Size):= C1 + C2;
Release(Divs);
Release(Divs2);
return Array;
];
func Abundant(N, Divs); \Returns 'true' if N is abundant
int N, Divs;
int Sum, I;
[Sum:= 0;
for I:= 0 to Divs(Size)-1 do
Sum:= Sum + Divs(Ptr,I);
return Sum > N;
];
func Semiperfect(N, Divs); \Returns 'true' if N is semiperfect
int N, Divs;
int H, T, TA(2);
[if Divs(Size) > 0 then
[H:= Divs(Ptr,0);
T:= Divs(Ptr)+SizeOfInt;
TA(Ptr):= T;
TA(Size):= Divs(Size)-1;
if N < H then
return Semiperfect(N, TA)
else return N = H or Semiperfect(N-H, TA) or Semiperfect(N, TA);
]
else return false;
];
func Sieve(Limit); \Return array of weird number indexes set 'false'
int Limit; \i.e. non-abundant and non-semiperfect
int W, Divs(2), I, J;
[W:= MAlloc(Limit * SizeOfInt);
for I:= 0 to Limit-1 do W(I):= 0; \for safety
I:= 2;
while I < Limit do
[if W(I) = 0 then
[Divs:= Divisors(I);
if not Abundant(I, Divs) then
W(I):= true
else if Semiperfect(I, Divs) then
[J:= I;
while J < Limit do
[W(J):= true;
J:= J+I;
];
];
];
I:= I+2;
];
Release(Divs(Ptr));
return W;
];
int W, Count, Max, N;
[W:= Sieve(17000);
Count:= 0;
Max:= 25;
Text(0, "The first 25 weird numbers:^m^j");
N:= 2;
while Count < Max do
[if not W(N) then
[IntOut(0, N); ChOut(0, ^ );
Count:= Count+1;
];
N:= N+2;
];
CrLf(0);
Release(W);
]
- Output:
The first 25 weird numbers: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310
zkl
fcn properDivs(n){
if(n==1) return(T);
( pd:=[1..(n).toFloat().sqrt()].filter('wrap(x){ n%x==0 }) )
.pump(pd,'wrap(pd){ if(pd!=1 and (y:=n/pd)!=pd ) y else Void.Skip })
}
fcn abundant(n,divs){ divs.sum(0) > n }
fcn semiperfect(n,divs){
if(divs){
h,t := divs[0], divs[1,*];
if(n<h) return(semiperfect(n,t));
return((n==h) or semiperfect(n - h, t) or semiperfect(n, t));
}
False
}
fcn sieve(limit){
// False denotes abundant and not semi-perfect.
// Only interested in even numbers >= 2
w:=List.createLong(limit,False);
foreach i in ([2..limit - 1, 2]){
if(w[i]) continue;
divs:=properDivs(i);
if(not abundant(i,divs)) w[i]=True;
else if(semiperfect(i,divs))
{ foreach j in ([i..limit - 1, i]){ w[j]=True; } }
}
w
}
w,count,max := sieve(17_000), 0, 25;
println("The first 25 weird numbers are:");
foreach n in ([2..* ,2]){
if(not w[n]){ print("%d ".fmt(n)); count+=1; }
if(count>=max) break;
}
println();
- Output:
The first 25 weird numbers are: 70 836 4030 5830 7192 7912 9272 10430 10570 10792 10990 11410 11690 12110 12530 12670 13370 13510 13790 13930 14770 15610 15890 16030 16310