Self numbers: Difference between revisions
m (→sliding sieve: removed blank line) |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
(58 intermediate revisions by 21 users not shown) | |||
Line 12: | Line 12: | ||
;*[[oeis:A003052|OEIS: A003052 - Self numbers or Colombian numbers]] |
;*[[oeis:A003052|OEIS: A003052 - Self numbers or Colombian numbers]] |
||
;*[[wp:Self_number|Wikipedia: Self numbers]] |
;*[[wp:Self_number|Wikipedia: Self numbers]] |
||
=={{header|ALGOL 68}}== |
|||
{{Trans|Go}} |
|||
<syntaxhighlight lang="algol68">BEGIN # find some self numbers numbers n such that there is no g such that g + sum of g's digits = n # |
|||
INT max number = 1 999 999 999 + 82; # maximum n plus g we will condifer # |
|||
# sieve the self numbers up to 1 999 999 999 # |
|||
[ 0 : max number ]BOOL self; FOR i TO UPB self DO self[ i ] := TRUE OD; |
|||
INT n := 0; |
|||
FOR s0 FROM 0 TO 1 DO |
|||
FOR d1 FROM 0 TO 9 DO |
|||
INT s1 = s0 + d1; |
|||
FOR d2 FROM 0 TO 9 DO |
|||
INT s2 = s1 + d2; |
|||
FOR d3 FROM 0 TO 9 DO |
|||
INT s3 = s2 + d3; |
|||
FOR d4 FROM 0 TO 9 DO |
|||
INT s4 = s3 + d4; |
|||
FOR d5 FROM 0 TO 9 DO |
|||
INT s5 = s4 + d5; |
|||
FOR d6 FROM 0 TO 9 DO |
|||
INT s6 = s5 + d6; |
|||
FOR d7 FROM 0 TO 9 DO |
|||
INT s7 = s6 + d7; |
|||
FOR d8 FROM 0 TO 9 DO |
|||
INT s8 = s7 + d8; |
|||
FOR d9 FROM 0 TO 9 DO |
|||
INT s9 = s8 + d9; |
|||
self[ s9 + n ] := FALSE; |
|||
n +:= 1 |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD |
|||
OD; |
|||
# show the first 50 self numbers # |
|||
INT s count := 0; |
|||
FOR i TO UPB self WHILE s count < 50 DO |
|||
IF self[ i ] THEN |
|||
print( ( " ", whole( i, -3 ) ) ); |
|||
IF ( s count +:= 1 ) MOD 18 = 0 THEN print( ( newline ) ) FI |
|||
FI |
|||
OD; |
|||
print( ( newline ) ); |
|||
# show the self numbers with power-of-10 indxes # |
|||
INT s show := 1; |
|||
s count := 0; |
|||
print( ( " nth self", newline ) ); |
|||
print( ( " n number", newline ) ); |
|||
FOR i TO UPB self DO |
|||
IF self[ i ] THEN |
|||
s count +:= 1; |
|||
IF s count = s show THEN |
|||
print( ( whole( s show, -9 ), " ", whole( i, -11 ), newline ) ); |
|||
s show *:= 10 |
|||
FI |
|||
FI |
|||
OD |
|||
END</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 |
|||
154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 |
|||
334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
nth self |
|||
n number |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1000 10188 |
|||
10000 102225 |
|||
100000 1022675 |
|||
1000000 10227221 |
|||
10000000 102272662 |
|||
100000000 1022727208 |
|||
</pre> |
|||
=={{header|AppleScript}}== |
=={{header|AppleScript}}== |
||
I couldn't follow the math in the Wikipedia entry, nor the discussion and code here so far. But an initial expedient of generating a list of all the integers from 1 to just over ten times the required number of results and then deleting those that ''could'' be derived |
I couldn't follow the math in the Wikipedia entry, nor the discussion and code here so far. But an initial expedient of generating a list of all the integers from 1 to just over ten times the required number of results and then deleting those that ''could'' be derived by the described method revealed the sequencing pattern on which the code below is based. On the test machine, it completes all three of the tests at the bottom in a total of around a millisecond. |
||
< |
<syntaxhighlight lang="applescript">(* |
||
Base-10 self numbers by index (single or range). |
Base-10 self numbers by index (single or range). |
||
Follows an observed sequence pattern whereby, after the initial single-digit odd numbers, self numbers are |
|||
grouped in runs whose members occur at numeric intervals of 11. Runs after the first one come in blocks of |
|||
ten: eight runs of ten numbers followed by two shorter runs. The numeric interval between runs is usually 2, |
|||
but that between shorter runs, and their length, depend on the highest-order digit change occurring in them. |
|||
This connection with significant digit change means every ten blocks form a higher-order block, every ten |
|||
digit to have changed in the number sequence since they were last set. |
|||
of these a higher-order-still block, and so on. |
|||
The code below appears to be good up to the last self number before 10^12 — ie. 999,999,999,997, which is |
|||
returned as the 97,777,777,792nd such number. After this, instead of zero-length shorter runs, the actual |
|||
pattern apparently starts again with a single run of 10, like the one at the beginning. |
|||
*) |
*) |
||
on selfNumbers(indexRange) |
on selfNumbers(indexRange) |
||
-- Subscript to monitor progress, store required results, and indicate when finished. |
|||
set indexRange to indexRange as list |
set indexRange to indexRange as list |
||
-- Script object with subhandlers and associated properties. |
|||
script storer |
|||
script |subscript| |
|||
property startIndex : beginning of indexRange |
property startIndex : beginning of indexRange |
||
property endIndex : end of indexRange |
property endIndex : end of indexRange |
||
property counter : 0 |
property counter : 0 |
||
property currentSelf : -1 |
|||
property output : {} |
property output : {} |
||
-- Advance to the next self number in the sequence, append it to the output if required, indicate if finished. |
|||
on doneAfter(thisSelf) |
|||
on doneAfterAdding(interval) |
|||
set currentSelf to currentSelf + interval |
|||
set counter to counter + 1 |
set counter to counter + 1 |
||
if (counter < startIndex) then return false |
if (counter < startIndex) then return false |
||
set end of my output to |
set end of my output to currentSelf |
||
return (counter = endIndex) |
return (counter = endIndex) |
||
end |
end doneAfterAdding |
||
-- If necessary, fast forward to the last self number before the lowest-order block containing the first number required. |
|||
on fastForward() |
|||
if (counter ≥ startIndex) then return |
|||
-- The highest-order blocks whose ends this script handles correctly contain 9,777,777,778 self numbers. |
|||
-- The difference between equivalently positioned numbers in these blocks is 100,000,000,001. |
|||
-- The figures for successively lower-order blocks have successively fewer 7s and 0s! |
|||
set indexDiff to 9.777777778E+9 |
|||
set numericDiff to 1.00000000001E+11 |
|||
repeat until ((indexDiff < 98) or (counter = startIndex)) |
|||
set test to counter + indexDiff |
|||
if (test < startIndex) then |
|||
set counter to test |
|||
set currentSelf to (currentSelf + numericDiff) |
|||
else |
|||
set indexDiff to (indexDiff + 2) div 10 |
|||
set numericDiff to numericDiff div 10 + 1 |
|||
end if |
|||
end repeat |
|||
end fastForward |
|||
-- Work out a shorter run length based on the most significant digit change about to happen. |
|||
on getShorterRunLength() |
|||
set shorterRunLength to 9 |
|||
set temp to (|subscript|'s currentSelf) div 1000 |
|||
repeat while (temp mod 10 is 9) |
|||
set shorterRunLength to shorterRunLength - 1 |
|||
set temp to temp div 10 |
|||
end repeat |
|||
return shorterRunLength |
|||
end getShorterRunLength |
|||
end script |
end script |
||
-- Start with the single-digit |
-- Main process. Start with the single-digit odd numbers and first run. |
||
repeat 5 times |
|||
if (|subscript|'s doneAfterAdding(2)) then return |subscript|'s output |
|||
repeat 4 times |
|||
end repeat |
|||
set thisSelf to thisSelf + 2 |
|||
repeat 9 times |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output |
|||
end repeat |
end repeat |
||
-- Fast forward if the start index hasn't yet been reached. |
|||
tell |subscript| to fastForward() |
|||
-- |
-- Sequencing loop, per lowest-order block. |
||
repeat |
|||
set shortGroupLength to 10 |
|||
-- Eight ten-number runs, each at a numeric interval of 2 from the end of the previous one. |
|||
repeat -- Per block. |
|||
-- First short group. Its length in the first block is 10, thereafter the same as that of the second |
|||
-- short group of the preceding block. The numeric interval preceding it is related to this length. |
|||
set thisSelf to thisSelf + 2 + (10 - shortGroupLength) * 13 |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
-- Consecutive numbers within groups are 11 apart. |
|||
repeat (shortGroupLength - 1) times |
|||
set thisSelf to thisSelf + 11 |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
end repeat |
|||
-- Eight ten-number groups, each preceded by an interval of 2. |
|||
repeat 8 times |
repeat 8 times |
||
if (|subscript|'s doneAfterAdding(2)) then return |subscript|'s output |
|||
set thisSelf to thisSelf + 2 |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
repeat 9 times |
repeat 9 times |
||
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
end repeat |
end repeat |
||
end repeat |
end repeat |
||
-- Two shorter runs, the second at an interval inversely related to their length. |
|||
set shorterRunLength to |subscript|'s getShorterRunLength() |
|||
-- Second short group. Its length depends on the most significant digit to tick over |
|||
repeat with interval in {2, 2 + (10 - shorterRunLength) * 13} |
|||
-- at the thousands boundary within it. The interval preceding it is 2. |
|||
if (|subscript|'s doneAfterAdding(interval)) then return |subscript|'s output |
|||
set thisSelf to thisSelf + 2 |
|||
repeat (shorterRunLength - 1) times |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output |
|||
set i to 2 |
|||
end repeat |
|||
repeat until (i > shortGroupLength) |
|||
set thisSelf to thisSelf + 11 |
|||
if (storer's doneAfter(thisSelf)) then return storer's output |
|||
if (thisSelf mod 1000 < 11) then |
|||
set temp to thisSelf div 1000 |
|||
repeat while (temp mod 10 is 0) |
|||
set shortGroupLength to shortGroupLength - 1 |
|||
set temp to temp div 10 |
|||
end repeat |
|||
end if |
|||
set i to i + 1 |
|||
end repeat |
end repeat |
||
end repeat |
end repeat |
||
return storer's output |
|||
end selfNumbers |
end selfNumbers |
||
Line 102: | Line 198: | ||
-- One hundred millionth: |
-- One hundred millionth: |
||
selfNumbers(100000000) |
selfNumbers(100000000) |
||
--> { |
--> {1.022727208E+9} |
||
-- 97,777,777,792nd: |
|||
selfNumbers(9.7777777792E+10) |
|||
--> {9.99999999997E+11}</syntaxhighlight> |
|||
=={{header|AWK}}== |
|||
<syntaxhighlight lang="awk"> |
|||
# syntax: GAWK -f SELF_NUMBERS.AWK |
|||
# converted from Go (low memory example) |
|||
BEGIN { |
|||
print("HH:MM:SS INDEX SELF") |
|||
print("-------- ---------- ----------") |
|||
count = 0 |
|||
digits = 1 |
|||
i = 1 |
|||
last_self = 0 |
|||
offset = 9 |
|||
pow = 10 |
|||
while (count < 1E8) { |
|||
is_self = 1 |
|||
start = max(i-offset,0) |
|||
sum = sum_digits(start) |
|||
for (j=start; j<i; j++) { |
|||
if (j + sum == i) { |
|||
is_self = 0 |
|||
break |
|||
} |
|||
sum = ((j+1) % 10 != 0) ? ++sum : sum_digits(j+1) |
|||
} |
|||
if (is_self) { |
|||
last_self = i |
|||
if (++count <= 50) { |
|||
selfs = selfs i " " |
|||
} |
|||
} |
|||
if (++i % pow == 0) { |
|||
pow *= 10 |
|||
digits++ |
|||
offset = digits * 9 |
|||
} |
|||
if (count ~ /^10*$/ && arr[count]++ == 0) { |
|||
printf("%8s %10s %10s\n",strftime("%H:%M:%S"),count,last_self) |
|||
} |
|||
} |
|||
printf("\nfirst 50 self numbers:\n%s\n",selfs) |
|||
exit(0) |
|||
} |
|||
function sum_digits(x, sum,y) { |
|||
while (x) { |
|||
y = x % 10 |
|||
sum += y |
|||
x = int(x/10) |
|||
} |
|||
return(sum) |
|||
} |
|||
function max(x,y) { return((x > y) ? x : y) } |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
HH:MM:SS INDEX SELF |
|||
-------- ---------- ---------- |
|||
00:36:35 1 1 |
|||
00:36:35 10 64 |
|||
00:36:35 100 973 |
|||
00:36:35 1000 10188 |
|||
00:36:36 10000 102225 |
|||
00:36:46 100000 1022675 |
|||
00:38:49 1000000 10227221 |
|||
01:03:01 10000000 102272662 |
|||
05:27:35 100000000 1022727208 |
|||
first 50 self numbers: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
Line 108: | Line 278: | ||
{{trans|Go}} |
{{trans|Go}} |
||
About 25% faster than Go (using GCC 7.5.0 -O3) mainly due to being able to iterate through the sieve using a pointer. |
About 25% faster than Go (using GCC 7.5.0 -O3) mainly due to being able to iterate through the sieve using a pointer. |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <time.h> |
#include <time.h> |
||
Line 171: | Line 341: | ||
printf("Took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC); |
printf("Took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 184: | Line 354: | ||
===Extended=== |
===Extended=== |
||
{{trans|Pascal}} |
{{trans|Pascal}} |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <time.h> |
#include <time.h> |
||
Line 260: | Line 430: | ||
printf("\nOverall took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC); |
printf("\nOverall took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 283: | Line 453: | ||
Overall took 11.574158 seconds. |
Overall took 11.574158 seconds. |
||
</pre> |
</pre> |
||
=={{header|C++}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="cpp">#include <array> |
|||
#include <iomanip> |
|||
#include <iostream> |
|||
const int MC = 103 * 1000 * 10000 + 11 * 9 + 1; |
|||
std::array<bool, MC + 1> SV; |
|||
void sieve() { |
|||
std::array<int, 10000> dS; |
|||
for (int a = 9, i = 9999; a >= 0; a--) { |
|||
for (int b = 9; b >= 0; b--) { |
|||
for (int c = 9, s = a + b; c >= 0; c--) { |
|||
for (int d = 9, t = s + c; d >= 0; d--) { |
|||
dS[i--] = t + d; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
for (int a = 0, n = 0; a < 103; a++) { |
|||
for (int b = 0, d = dS[a]; b < 1000; b++, n += 10000) { |
|||
for (int c = 0, s = d + dS[b] + n; c < 10000; c++) { |
|||
SV[dS[c] + s++] = true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
int main() { |
|||
sieve(); |
|||
std::cout << "The first 50 self numbers are:\n"; |
|||
for (int i = 0, count = 0; count <= 50; i++) { |
|||
if (!SV[i]) { |
|||
count++; |
|||
if (count <= 50) { |
|||
std::cout << i << ' '; |
|||
} else { |
|||
std::cout << "\n\n Index Self number\n"; |
|||
} |
|||
} |
|||
} |
|||
for (int i = 0, limit = 1, count = 0; i < MC; i++) { |
|||
if (!SV[i]) { |
|||
if (++count == limit) { |
|||
//System.out.printf("%,12d %,13d%n", count, i); |
|||
std::cout << std::setw(12) << count << " " << std::setw(13) << i << '\n'; |
|||
limit *= 10; |
|||
} |
|||
} |
|||
} |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The first 50 self numbers are: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
Index Self number |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1000 10188 |
|||
10000 102225 |
|||
100000 1022675 |
|||
1000000 10227221 |
|||
10000000 102272662 |
|||
100000000 1022727208</pre> |
|||
=={{header|C#|CSharp}}== |
=={{header|C#|CSharp}}== |
||
{{trans|Pascal}}via{{trans|Go}}(third version) Stripped down, as C# limits the size of an array to Int32.MaxValue, so the sieve isn't large enough to hit the 1,000,000,000th value. |
{{trans|Pascal}}via{{trans|Go}}(third version) Stripped down, as C# limits the size of an array to Int32.MaxValue, so the sieve isn't large enough to hit the 1,000,000,000th value. |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using static System.Console; |
using static System.Console; |
||
Line 318: | Line 558: | ||
WriteLine("\nOverall took {0}s", (DateTime.Now - st). TotalSeconds); |
WriteLine("\nOverall took {0}s", (DateTime.Now - st). TotalSeconds); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}}Timing from tio.run |
{{out}}Timing from tio.run |
||
<pre>Sieving took 3.4266187s |
<pre>Sieving took 3.4266187s |
||
Line 337: | Line 577: | ||
Overall took 7.0237244s</pre> |
Overall took 7.0237244s</pre> |
||
=={{header|Elixir}}== |
|||
<syntaxhighlight lang="elixir"> |
|||
defmodule SelfNums do |
|||
def digAndSum(number) when is_number(number) do |
|||
Integer.digits(number) |> |
|||
Enum.reduce( 0, fn(num, acc) -> num + acc end ) |> |
|||
(fn(x) -> x + number end).() |
|||
end |
|||
def selfFilter(list, firstNth) do |
|||
numRange = Enum.to_list 1..firstNth |
|||
numRange -- list |
|||
end |
|||
end |
|||
defmodule SelfTest do |
|||
import SelfNums |
|||
stop = 1000 |
|||
Enum.to_list 1..stop |> |
|||
Enum.map(&digAndSum/1) |> |
|||
SelfNums.selfFilter(stop) |> |
|||
Enum.take(50) |> |
|||
IO.inspect |
|||
end |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 110, 121, 132, 143, 154, |
|||
165, 176, 187, 198, 209, 211, 222, 233, 244, 255, 266, 277, 288, 299, 310, 312, |
|||
323, 334, 345, 356, 367, 378, 389, 400, 411, 413, 424, 435, 446, 457, 468] |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |
||
< |
<syntaxhighlight lang="fsharp"> |
||
// Self numbers. Nigel Galloway: October 6th., 2020 |
// Self numbers. Nigel Galloway: October 6th., 2020 |
||
let fN g=let rec fG n g=match n/10 with 0->n+g |i->fG i (g+(n%10)) in fG g g |
let fN g=let rec fG n g=match n/10 with 0->n+g |i->fG i (g+(n%10)) in fG g g |
||
Line 346: | Line 623: | ||
Self |> Seq.take 50 |> Seq.iter(printf "%d "); printfn "" |
Self |> Seq.take 50 |> Seq.iter(printf "%d "); printfn "" |
||
printfn "\n%d" (Seq.item 99999999 Self) |
printfn "\n%d" (Seq.item 99999999 Self) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 353: | Line 630: | ||
1022727208 |
1022727208 |
||
</pre> |
</pre> |
||
=={{header|FreeBASIC}}== |
|||
{{trans|Ring}} |
|||
<syntaxhighlight lang="freebasic">Print "The first 50 self numbers are:" |
|||
Dim As Boolean flag |
|||
Dim As Ulong m, p, sum, number(), sum2 |
|||
Dim As Ulong n = 0 |
|||
Dim As Ulong num = 0 |
|||
Dim As Ulong limit = 51 |
|||
Dim As Ulong limit2 = 100000000 |
|||
Dim As String strnum |
|||
Do |
|||
n += 1 |
|||
For m = 1 To n |
|||
flag = True |
|||
sum = 0 |
|||
strnum = Str(m) |
|||
For p = 1 To Len(strnum) |
|||
sum += Val(Mid(strnum,p,1)) |
|||
Next p |
|||
sum2 = m + sum |
|||
If sum2 = n Then |
|||
flag = False |
|||
Exit For |
|||
Else |
|||
flag = True |
|||
End If |
|||
Next m |
|||
If flag = True Then |
|||
num += 1 |
|||
If num < limit Then Print ""; num; ". "; n |
|||
If num >= limit2 Then |
|||
Print "The "; limit2; "th self number is: "; n |
|||
Exit Do |
|||
End If |
|||
End If |
|||
Loop |
|||
Sleep</syntaxhighlight> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
===Low memory=== |
===Low memory=== |
||
Simple-minded, uses very little memory (no sieve) but slow - over 2.5 minutes. |
Simple-minded, uses very little memory (no sieve) but slow - over 2.5 minutes. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 424: | Line 742: | ||
fmt.Println("\nThe 100 millionth self number is", lastSelf) |
fmt.Println("\nThe 100 millionth self number is", lastSelf) |
||
fmt.Println("Took", time.Since(st)) |
fmt.Println("Took", time.Since(st)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 441: | Line 759: | ||
Have also incorporated Enter your username's suggestion (see Talk page) of using partial sums for each loop which improves performance by about 25%. |
Have also incorporated Enter your username's suggestion (see Talk page) of using partial sums for each loop which improves performance by about 25%. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 503: | Line 821: | ||
} |
} |
||
fmt.Println("Took", time.Since(st)) |
fmt.Println("Took", time.Since(st)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 517: | Line 835: | ||
{{trans|Pascal}} |
{{trans|Pascal}} |
||
This uses horst.h's ideas (see Talk page) to find up to the 1 billionth self number in a reasonable time and using less memory than the simple 'sieve based' approach above would have needed. |
This uses horst.h's ideas (see Talk page) to find up to the 1 billionth self number in a reasonable time and using less memory than the simple 'sieve based' approach above would have needed. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 592: | Line 910: | ||
} |
} |
||
fmt.Println("\nOverall took", time.Since(st)) |
fmt.Println("\nOverall took", time.Since(st)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 615: | Line 933: | ||
Overall took 14.647314803s |
Overall took 14.647314803s |
||
</pre> |
</pre> |
||
=={{header|Haskell}}== |
|||
The solution is quite straightforward. The length of the foreseeing window in filtering procedure (81) is chosen empirically and doesn't have any theoretical background. |
|||
<syntaxhighlight lang="haskell">import Control.Monad (forM_) |
|||
import Text.Printf |
|||
selfs :: [Integer] |
|||
selfs = sieve (sumFs [0..]) [0..] |
|||
where |
|||
sumFs = zipWith (+) [ a+b+c+d+e+f+g+h+i+j |
|||
| a <- [0..9] , b <- [0..9] |
|||
, c <- [0..9] , d <- [0..9] |
|||
, e <- [0..9] , f <- [0..9] |
|||
, g <- [0..9] , h <- [0..9] |
|||
, i <- [0..9] , j <- [0..9] ] |
|||
-- More idiomatic list generator is about three times slower |
|||
-- sumFs = zipWith (+) $ sum <$> replicateM 10 [0..9] |
|||
sieve (f:fs) (n:ns) |
|||
| n > f = sieve fs (n:ns) |
|||
| n `notElem` take 81 (f:fs) = n : sieve (f:fs) ns |
|||
| otherwise = sieve (f:fs) ns |
|||
main = do |
|||
print $ take 50 selfs |
|||
forM_ [1..8] $ \i -> |
|||
printf "1e%v\t%v\n" (i :: Int) (selfs !! (10^i-1))</syntaxhighlight> |
|||
<pre>$ ghc -O2 SelfNum.hs && time ./SelfNum |
|||
[1,3,5,7,9,20,31,42,53,64,75,86,97,108,110,121,132,143,154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323,334,345,356,367,378,389,400,411,413,424,435,446,457,468] |
|||
1e1 64 |
|||
1e2 973 |
|||
1e3 10188 |
|||
1e4 102225 |
|||
1e5 1022675 |
|||
1e6 10227221 |
|||
1e7 102272662 |
|||
1e8 1022727208 |
|||
275.98 user 3.11 system 4:41.02 elapsed</pre> |
|||
=={{header|Java}}== |
|||
{{trans|C#}} |
|||
<syntaxhighlight lang="java">public class SelfNumbers { |
|||
private static final int MC = 103 * 1000 * 10000 + 11 * 9 + 1; |
|||
private static final boolean[] SV = new boolean[MC + 1]; |
|||
private static void sieve() { |
|||
int[] dS = new int[10_000]; |
|||
for (int a = 9, i = 9999; a >= 0; a--) { |
|||
for (int b = 9; b >= 0; b--) { |
|||
for (int c = 9, s = a + b; c >= 0; c--) { |
|||
for (int d = 9, t = s + c; d >= 0; d--) { |
|||
dS[i--] = t + d; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
for (int a = 0, n = 0; a < 103; a++) { |
|||
for (int b = 0, d = dS[a]; b < 1000; b++, n += 10000) { |
|||
for (int c = 0, s = d + dS[b] + n; c < 10000; c++) { |
|||
SV[dS[c] + s++] = true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
public static void main(String[] args) { |
|||
sieve(); |
|||
System.out.println("The first 50 self numbers are:"); |
|||
for (int i = 0, count = 0; count <= 50; i++) { |
|||
if (!SV[i]) { |
|||
count++; |
|||
if (count <= 50) { |
|||
System.out.printf("%d ", i); |
|||
} else { |
|||
System.out.printf("%n%n Index Self number%n"); |
|||
} |
|||
} |
|||
} |
|||
for (int i = 0, limit = 1, count = 0; i < MC; i++) { |
|||
if (!SV[i]) { |
|||
if (++count == limit) { |
|||
System.out.printf("%,12d %,13d%n", count, i); |
|||
limit *= 10; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The first 50 self numbers are: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
Index Self number |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1,000 10,188 |
|||
10,000 102,225 |
|||
100,000 1,022,675 |
|||
1,000,000 10,227,221 |
|||
10,000,000 102,272,662 |
|||
100,000,000 1,022,727,208</pre> |
|||
=={{header|jq}}== |
|||
'''Adapted from [[#Julia|Julia]]''' |
|||
{{works with|jq}} |
|||
<syntaxhighlight lang="jq"> |
|||
def sumdigits: tostring | explode | map([.]|implode|tonumber) | add; |
|||
def gsum: . + sumdigits; |
|||
def isnonself: |
|||
def ndigits: tostring|length; |
|||
. as $i |
|||
| ($i - ($i|ndigits)*9) as $n |
|||
| any( range($i-1; [0,$n]|max; -1); |
|||
gsum == $i); |
|||
# an array |
|||
def last81: |
|||
[limit(81; range(1; infinite) | select(isnonself))]; |
|||
# output an unbounded stream |
|||
def selfnumbers: |
|||
foreach range(1; infinite) as $i ( [0, last81]; |
|||
.[0] = false |
|||
| .[1] as $last81 |
|||
| if $last81 | bsearch($i) < 0 |
|||
then .[0] = $i |
|||
| if $i > $last81[-1] then .[1] = $last81[1:] + [$i | gsum ] else . end |
|||
else . |
|||
end; |
|||
.[0] | select(.) ); |
|||
"The first 50 self numbers are:", last81[:50], |
|||
"", |
|||
(nth(100000000 - 1; selfnumbers) |
|||
| if . == 1022727208 |
|||
then "Yes, \(.) is the 100,000,000th self number." |
|||
else "No, \(.i) is actually the 100,000,000th self number." |
|||
end)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
The first 50 self numbers are: |
|||
[2,4,6,8,10,11,12,13,14,15,16,17,18,19,21,22,23,24,25,26,27,28,29,30,32,33,34,35,36,37,38,39,40,41,43,44,45,46,47,48,49,50,51,52,54,55,56,57,58,59] |
|||
Yes, 1022727208 is the 100,000,000th self number. |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 620: | Line 1,093: | ||
Note that 81 is the window size because the sum of digits of 999,999,999 |
Note that 81 is the window size because the sum of digits of 999,999,999 |
||
(the largest digit sum of a counting number less than 1022727208) is 81. |
(the largest digit sum of a counting number less than 1022727208) is 81. |
||
< |
<syntaxhighlight lang="julia">gsum(i) = sum(digits(i)) + i |
||
isnonself(i) = any(x -> gsum(x) == i, i-1:-1:i-max(1, ndigits(i)*9)) |
isnonself(i) = any(x -> gsum(x) == i, i-1:-1:i-max(1, ndigits(i)*9)) |
||
const last81 = filter(isnonself, 1:5000)[1:81] |
const last81 = filter(isnonself, 1:5000)[1:81] |
||
Line 646: | Line 1,119: | ||
checkselfnumbers() |
checkselfnumbers() |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
||
Line 655: | Line 1,128: | ||
{{trans|Pascal}} |
{{trans|Pascal}} |
||
Contains tweaks peculiar to the "10 to the nth" self number. Timings include compilation times. |
Contains tweaks peculiar to the "10 to the nth" self number. Timings include compilation times. |
||
< |
<syntaxhighlight lang="julia">const MAXCOUNT = 103 * 10000 * 10000 + 11 * 9 + 1 |
||
function dosieve!(sieve, digitsum9999) |
function dosieve!(sieve, digitsum9999) |
||
Line 698: | Line 1,171: | ||
@time findselves() |
@time findselves() |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
Sieve time: |
Sieve time: |
||
Line 715: | Line 1,188: | ||
16.999383 seconds (42.92 k allocations: 9.595 GiB, 0.01% gc time) |
16.999383 seconds (42.92 k allocations: 9.595 GiB, 0.01% gc time) |
||
</pre> |
</pre> |
||
=={{header|Kotlin}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="scala">private const val MC = 103 * 1000 * 10000 + 11 * 9 + 1 |
|||
private val SV = BooleanArray(MC + 1) |
|||
private fun sieve() { |
|||
val dS = IntArray(10000) |
|||
run { |
|||
var a = 9 |
|||
var i = 9999 |
|||
while (a >= 0) { |
|||
for (b in 9 downTo 0) { |
|||
var c = 9 |
|||
val s = a + b |
|||
while (c >= 0) { |
|||
var d = 9 |
|||
val t = s + c |
|||
while (d >= 0) { |
|||
dS[i--] = t + d |
|||
d-- |
|||
} |
|||
c-- |
|||
} |
|||
} |
|||
a-- |
|||
} |
|||
} |
|||
var a = 0 |
|||
var n = 0 |
|||
while (a < 103) { |
|||
var b = 0 |
|||
val d = dS[a] |
|||
while (b < 1000) { |
|||
var c = 0 |
|||
var s = d + dS[b] + n |
|||
while (c < 10000) { |
|||
SV[dS[c] + s++] = true |
|||
c++ |
|||
} |
|||
b++ |
|||
n += 10000 |
|||
} |
|||
a++ |
|||
} |
|||
} |
|||
fun main() { |
|||
sieve() |
|||
println("The first 50 self numbers are:") |
|||
run { |
|||
var i = 0 |
|||
var count = 0 |
|||
while (count <= 50) { |
|||
if (!SV[i]) { |
|||
count++ |
|||
if (count <= 50) { |
|||
print("$i ") |
|||
} else { |
|||
println() |
|||
println() |
|||
println(" Index Self number") |
|||
} |
|||
} |
|||
i++ |
|||
} |
|||
} |
|||
var i = 0 |
|||
var limit = 1 |
|||
var count = 0 |
|||
while (i < MC) { |
|||
if (!SV[i]) { |
|||
if (++count == limit) { |
|||
println("%,12d %,13d".format(count, i)) |
|||
limit *= 10 |
|||
} |
|||
} |
|||
i++ |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The first 50 self numbers are: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
Index Self number |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1,000 10,188 |
|||
10,000 102,225 |
|||
100,000 1,022,675 |
|||
1,000,000 10,227,221 |
|||
10,000,000 102,272,662 |
|||
100,000,000 1,022,727,208</pre> |
|||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica"> |
|||
sum[g_] := g + Total@IntegerDigits@g |
|||
ming[n_] := n - IntegerLength[n]*9 |
|||
self[n_] := NoneTrue [Range[ming[n], n - 1], sum[#] == n &] |
|||
Module[{t = 1, x = 1}, |
|||
Reap[ |
|||
While[t <= 50, |
|||
If[self[x], Sow[x]; t++]; x++] |
|||
][[2, 1]]] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{1,3,5,7,9,20,31,42,53,64,75,86,97,108,110,121,132,143,154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323,334,345,356,367,378,389,400,411,413,424,435,446,457,468}</pre> |
|||
=={{header|Nim}}== |
|||
In order to use less memory, we have chosen to use indexing at bit level. So, our sieve is a custom object defined by its length in bits and its value which is a sequence of bytes. With bit indexing, the sieve uses eight times less memory than with byte indexing. |
|||
Of course, there is a trade off to this strategy: reading values from and writing values to the sieve are significantly slower. |
|||
===Simple sieve=== |
|||
{{trans|Go}} |
|||
We use the Go algorithm with bit indexing. As a consequence the sieve uses about 250MB instead of 1 GB. And the program is about five times slower. |
|||
Note that we used a sequence of ten nested loops as in the Go solution but we have not memorized the intermediate sums as the C compiler does a good job to detect the loop invariants (remember, Nim produces C code and this code has proved to be quite optimizable by the C compiler). The ten loops looks a lot better this way, too 🙂. |
|||
<syntaxhighlight lang="nim">import bitops, strutils, std/monotimes, times |
|||
const MaxCount = 2 * 1_000_000_000 + 9 * 9 |
|||
# Bit string used to represent an array of booleans. |
|||
type BitString = object |
|||
len: Natural # length in bits. |
|||
values: seq[byte] # Sequence containing the bits. |
|||
proc newBitString(n: Natural): BitString = |
|||
## Return a new bit string of length "n" bits. |
|||
result.len = n |
|||
result.values.setLen((n + 7) shr 3) |
|||
template checkIndex(i, length: Natural) {.used.} = |
|||
## Check if index "i" is less than the array length. |
|||
if i >= length: |
|||
raise newException(IndexDefect, "index $1 not in 0 .. $2".format(i, length)) |
|||
proc `[]`(bs: BitString; i: Natural): bool = |
|||
## Return the value of bit at index "i" as a boolean. |
|||
when compileOption("boundchecks"): |
|||
checkIndex(i, bs.len) |
|||
result = bs.values[i shr 3].testbit(i and 0x07) |
|||
proc `[]=`(bs: var BitString; i: Natural; value: bool) = |
|||
## Set the bit at index "i" to the given value. |
|||
when compileOption("boundchecks"): |
|||
checkIndex(i, bs.len) |
|||
if value: bs.values[i shr 3].setBit(i and 0x07) |
|||
else: bs.values[i shr 3].clearBit(i and 0x07) |
|||
proc fill(sieve: var BitString) = |
|||
## Fill a sieve. |
|||
var n = 0 |
|||
for a in 0..1: |
|||
for b in 0..9: |
|||
for c in 0..9: |
|||
for d in 0..9: |
|||
for e in 0..9: |
|||
for f in 0..9: |
|||
for g in 0..9: |
|||
for h in 0..9: |
|||
for i in 0..9: |
|||
for j in 0..9: |
|||
sieve[a + b + c + d + e + f + g + h + i + j + n] = true |
|||
inc n |
|||
let t0 = getMonoTime() |
|||
var sieve = newBitString(MaxCount + 1) |
|||
sieve.fill() |
|||
echo "Sieve time: ", getMonoTime() - t0 |
|||
# Find first 50. |
|||
echo "\nFirst 50 self numbers:" |
|||
var count = 0 |
|||
var line = "" |
|||
for n in 0..MaxCount: |
|||
if not sieve[n]: |
|||
inc count |
|||
line.addSep(" ") |
|||
line.add $n |
|||
if count == 50: break |
|||
echo line |
|||
# Find 1st, 10th, 100th, ..., 100_000_000th. |
|||
echo "\n Rank Value" |
|||
var limit = 1 |
|||
count = 0 |
|||
for n in 0..MaxCount: |
|||
if not sieve[n]: inc count |
|||
if count == limit: |
|||
echo ($count).align(10), ($n).align(12) |
|||
limit *= 10 |
|||
echo "Total time: ", getMonoTime() - t0</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Sieve time: 2 seconds, 59 milliseconds, 67 microseconds, and 152 nanoseconds |
|||
First 50 self numbers: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
Rank Value |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1000 10188 |
|||
10000 102225 |
|||
100000 1022675 |
|||
1000000 10227221 |
|||
10000000 102272662 |
|||
100000000 1022727208 |
|||
Total time: 7 seconds, 903 milliseconds, 752 microseconds, and 944 nanoseconds</pre> |
|||
===Improved sieve=== |
|||
{{trans|Pascal}} |
|||
Using bit indexing is very useful here as with byte indexing, the sieve needs 10GB. On a computer with only 8 GB , as this is the case of the laptop I use to run these programs, it fails to execute (I have a very small swap and don’t want to use the swap anyway). With bit indexing, the sieve needs only 1,25GB which is more reasonable. |
|||
Of course, the program is slower but not in the same proportion that in the previous program: it is about twice slower than the Pascal version. Note that the execution time varies significantly according to the way statements are written. For instance, writing <code>if not sieve[n]: inc count</code> has proved to be more efficient than writing <code>inc count, ord(not sieve[n])</code> or <code>inc count, 1 - ord(sieve[n])</code> which is surprising as the latter forms avoid a jump. Maybe changing some other statements could give better results, but current time is already satisfying. |
|||
<syntaxhighlight lang="nim">import bitops, strutils, std/monotimes, times |
|||
const MaxCount = 103 * 10_000 * 10_000 + 11 * 9 + 1 |
|||
# Bit string used to represent an array of booleans. |
|||
type BitString = object |
|||
len: Natural |
|||
values: seq[byte] |
|||
proc newBitString(n: Natural): BitString = |
|||
## Return a new bit string of length "n" bits. |
|||
result.len = n |
|||
result.values.setLen((n + 7) shr 3) |
|||
template checkIndex(i, length: Natural) {.used.} = |
|||
## Check if index "i" is less than the array length. |
|||
if i >= length: |
|||
raise newException(IndexDefect, "index $1 not in 0 .. $2".format(i, length)) |
|||
proc `[]`(bs: BitString; i: Natural): bool = |
|||
## Return the value of bit at index "i" as a boolean. |
|||
when compileOption("boundchecks"): |
|||
checkIndex(i, bs.len) |
|||
result = bs.values[i shr 3].testbit(i and 0x07) |
|||
proc `[]=`(bs: var BitString; i: Natural; value: bool) = |
|||
## Set the bit at index "i" to the given value. |
|||
when compileOption("boundchecks"): |
|||
checkIndex(i, bs.len) |
|||
if value: bs.values[i shr 3].setBit(i and 0x07) |
|||
else: bs.values[i shr 3].clearBit(i and 0x07) |
|||
proc initDigitSum9999(): array[10000, byte] {.compileTime.} = |
|||
## Return the array of the digit sums for numbers 0 to 9999. |
|||
var i = 0 |
|||
for a in 0..9: |
|||
for b in 0..9: |
|||
for c in 0..9: |
|||
for d in 0..9: |
|||
result[i] = byte(a + b + c + d) |
|||
inc i |
|||
const DigitSum9999 = initDigitSum9999() |
|||
proc fill(sieve: var BitString) = |
|||
## Fill a sieve. |
|||
var n = 0 |
|||
for a in 0..102: |
|||
for b in 0..9999: |
|||
var s = DigitSum9999[a].int + DigitSum9999[b].int + n |
|||
for c in 0..9999: |
|||
sieve[DigitSum9999[c].int + s] = true |
|||
inc s |
|||
inc n, 10_000 |
|||
let t0 = getMonoTime() |
|||
var sieve = newBitString(MaxCount + 1) |
|||
sieve.fill() |
|||
echo "Sieve time: ", getMonoTime() - t0 |
|||
# Find first 50. |
|||
echo "\nFirst 50 self numbers:" |
|||
var count = 0 |
|||
var line = "" |
|||
for n in 0..MaxCount: |
|||
if not sieve[n]: |
|||
inc count |
|||
line.addSep(" ") |
|||
line.add $n |
|||
if count == 50: break |
|||
echo line |
|||
# Find 1st, 10th, 100th, ..., 1_000_000_000th. |
|||
echo "\n Rank Value" |
|||
var limit = 1 |
|||
count = 0 |
|||
for n in 0..MaxCount: |
|||
if not sieve[n]: inc count |
|||
if count == limit: |
|||
echo ($count).align(10), ($n).align(12) |
|||
limit *= 10 |
|||
echo "Total time: ", getMonoTime() - t0</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Sieve time: 13 seconds, 340 milliseconds, 45 microseconds, and 528 nanoseconds |
|||
First 50 self numbers: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
Rank Value |
|||
1 1 |
|||
10 64 |
|||
100 973 |
|||
1000 10188 |
|||
10000 102225 |
|||
100000 1022675 |
|||
1000000 10227221 |
|||
10000000 102272662 |
|||
100000000 1022727208 |
|||
1000000000 10227272649 |
|||
Total time: 28 seconds, 135 milliseconds, 481 microseconds, and 697 nanoseconds</pre> |
|||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
Line 720: | Line 1,533: | ||
Just "sieving" with only one follower of every number {{trans|Go}} |
Just "sieving" with only one follower of every number {{trans|Go}} |
||
Extended to 10.23e9 |
Extended to 10.23e9 |
||
< |
<syntaxhighlight lang="pascal">program selfnumb; |
||
{$IFDEF FPC} |
{$IFDEF FPC} |
||
{$MODE Delphi} |
{$MODE Delphi} |
||
Line 816: | Line 1,629: | ||
end; |
end; |
||
end; |
end; |
||
END.</ |
END.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> time ./selfnumb |
<pre> time ./selfnumb |
||
Line 834: | Line 1,647: | ||
real 0m13,252s</pre> |
real 0m13,252s</pre> |
||
=={{header| |
=={{header|Perl}}== |
||
{{trans| |
{{trans|Raku}} |
||
<syntaxhighlight lang="perl">use strict; |
|||
Replacing the problematic sv[a+b+... line with a bit of dirty inline assembly improved performance by 90%<br> |
|||
use warnings; |
|||
(Of course you lose bounds checking, type checking, negative subscripts, fraction handling, and all that jazz.)<br> |
|||
use feature 'say'; |
|||
We use a string of Y/N for the sieve to force one byte per element ('\0' and 1 would be equally valid). |
|||
use List::Util qw(max sum); |
|||
<lang Phix>if machine_bits()=32 then crash("requires 64 bit") end if |
|||
my ($i, $pow, $digits, $offset, $lastSelf, @selfs) |
|||
function sieve() |
|||
= ( 1, 10, 1, 9, 0, ); |
|||
string sv = repeat('N',2*1e9+9*9+1) -- (1.86GB) |
|||
integer n = 0 |
|||
for a=0 to 1 do |
|||
for b=0 to 9 do |
|||
for c=0 to 9 do |
|||
for d=0 to 9 do |
|||
for e=0 to 9 do |
|||
for f=0 to 9 do |
|||
for g=0 to 9 do |
|||
for h=0 to 9 do |
|||
for i=0 to 9 do |
|||
for j=0 to 9 do |
|||
-- n += 1 |
|||
-- sv[a+b+c+d+e+f+g+h+i+j+n] = 'Y' |
|||
#ilASM{ |
|||
[32] -- (allows clean compilation on 32 bit, before crash as above) |
|||
[64] |
|||
mov rax,[a] |
|||
mov r12,[b] |
|||
mov r13,[c] |
|||
mov r14,[d] |
|||
mov r15,[e] |
|||
add r12,r13 |
|||
add r14,r15 |
|||
add rax,r12 |
|||
mov rdi,[sv] |
|||
add rax,r14 |
|||
mov r12,[f] |
|||
mov r13,[g] |
|||
mov r14,[h] |
|||
mov r15,[i] |
|||
add r12,r13 |
|||
shl rdi,2 |
|||
mov rcx,[n] |
|||
mov r13,[j] |
|||
add r14,r15 |
|||
add rax,r12 |
|||
add rax,r14 |
|||
add r13,rcx |
|||
add rax,r13 |
|||
add rcx,1 |
|||
mov byte[rdi+rax],'Y' |
|||
mov [n],rcx } |
|||
end for |
|||
end for |
|||
end for |
|||
end for |
|||
end for |
|||
end for |
|||
end for |
|||
end for |
|||
printf(1,"%d,%d\r",{a,b}) -- (show progress) |
|||
end for |
|||
end for |
|||
return sv |
|||
end function |
|||
atom t0 = time() |
|||
string sv = sieve() |
|||
printf(1,"sieve build took %s\n",{elapsed(time()-t0)}) |
|||
integer count = 0 |
|||
printf(1,"The first 50 self numbers are:\n") |
|||
for i=1 to length(sv) do |
|||
if sv[i]='N' then |
|||
count += 1 |
|||
if count <= 50 then |
|||
printf(1,"%3d ",i-1) |
|||
if remainder(count,10)=0 then |
|||
printf(1,"\n") |
|||
end if |
|||
end if |
|||
if count == 1e8 then |
|||
printf(1,"\nThe %,dth self number is %,d (%s)\n", |
|||
{count,i-1,elapsed(time()-t0)}) |
|||
exit |
|||
end if |
|||
end if |
|||
end for</lang> |
|||
{{out}} |
|||
<pre> |
|||
sieve build took 6.6s |
|||
The first 50 self numbers are: |
|||
1 3 5 7 9 20 31 42 53 64 |
|||
75 86 97 108 110 121 132 143 154 165 |
|||
176 187 198 209 211 222 233 244 255 266 |
|||
277 288 299 310 312 323 334 345 356 367 |
|||
378 389 400 411 413 424 435 446 457 468 |
|||
my $final = 50; |
|||
The 100,000,000th self number is 1,022,727,208 (11.0s) |
|||
</pre> |
|||
while () { |
|||
=== generator dictionary === |
|||
my $isSelf = 1; |
|||
While this is dog-slow (see shocking estimate below), it is interesting to note that even by the time it generates the |
|||
my $sum = my $start = sum split //, max(($i-$offset), 0); |
|||
10,000,000th number, it is only having to juggle a mere 27 generators. |
|||
for ( my $j = $start; $j < $i; $j++ ) { |
|||
Just a shame that we had to push over 10,000,000 generators onto that stack, and tried to push quite a few more. |
|||
if ($j+$sum == $i) { $isSelf = 0 ; last } |
|||
Memory use is pretty low, around ~4MB.<br> |
|||
($j+1)%10 != 0 ? $sum++ : ( $sum = sum split '', ($j+1) ); |
|||
[unlike the above, this is perfectly happy on both 32 and 64 bit]<br> |
|||
} |
|||
Long story short: this works much the same as a prime sieve, in which you only need to eliminate multiples |
|||
of previous primes. Here, you only need to eliminate digital additions of previous safe numbers. |
|||
Only after writing this did I understand how to write a sliding sieve, which it turns out is a much better idea (see below). |
|||
Still, this is pretty interesting and quite neat. |
|||
<lang Phix>integer gd = new_dict(), gdhead = 2, n = 0 |
|||
if ($isSelf) { |
|||
function ng(integer n) |
|||
push @selfs, $lastSelf = $i; |
|||
integer r = n |
|||
last if @selfs == $final; |
|||
while r do |
|||
} |
|||
n += remainder(r,10) |
|||
r = floor(r/10) |
|||
end while |
|||
return n |
|||
end function |
|||
next unless ++$i % $pow == 0; |
|||
function self(integer /*i*/) |
|||
$pow *= 10; |
|||
-- note: assumes sequentially invoked (arg i unused) |
|||
$offset = 9 * $digits++ |
|||
} |
|||
while n=gdhead do |
|||
gdhead = pop_dict(gd)[1] |
|||
setd(ng(gdhead),0,gd) |
|||
n += (n!=gdhead) |
|||
end while |
|||
setd(ng(n),0,gd) |
|||
return n |
|||
end function |
|||
say "The first 50 self numbers are:\n" . join ' ', @selfs;</syntaxhighlight> |
|||
atom t0 = time() |
|||
{{out}} |
|||
printf(1,"The first 50 self numbers are:\n") |
|||
<pre>The first 50 self numbers are: |
|||
pp(apply(tagset(50),self),{pp_IntFmt,"%3d",pp_IntCh,false}) |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468</pre> |
|||
=={{header|Phix}}== |
|||
constant limit = 10000000 |
|||
{{trans|AppleScript}} |
|||
integer chk = 100 |
|||
Certainly puts my previous rubbish attempts ([[Self_numbers\Phix|archived here]]) to shame.<br> |
|||
printf(1,"\n i n size time\n") |
|||
The precise nature of the difference-pattern eludes me, I will admit. |
|||
for i=51 to limit do |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
n = self(i) |
|||
<span style="color: #000080;font-style:italic;">-- |
|||
if i=chk then |
|||
-- Base-10 self numbers by index (single or range). |
|||
printf(1,"%,11d %,11d %6d %s\n",{i,n,dict_size(gd),elapsed(time()-t0)}) |
|||
-- Follows an observed sequence pattern whereby, after the initial single-digit odd numbers, self numbers are |
|||
chk *= 10 |
|||
-- grouped in runs whose members occur at numeric intervals of 11. Runs after the first one come in blocks of |
|||
end if |
|||
-- ten: eight runs of ten numbers followed by two shorter runs. The numeric interval between runs is usually 2, |
|||
end for |
|||
-- but that between shorter runs, and their length, depend on the highest-order digit change occurring in them. |
|||
printf(1,"\nEstimated time for %,d :%s\n",{1e8,elapsed((time()-t0)*1e8/limit)})</lang> |
|||
-- This connection with significant digit change means every ten blocks form a higher-order block, every ten |
|||
-- of these a higher-order-still block, and so on. |
|||
-- |
|||
-- The code below appears to be good up to the last self number before 10^12 — ie. 999,999,999,997, which is |
|||
-- returned as the 97,777,777,792nd such number. After this, instead of zero-length shorter runs, the actual |
|||
-- pattern apparently starts again with a single run of 10, like the one at the beginning. |
|||
--</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">startIndex</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">endIndex</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">counter</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">currentSelf</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">output</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">interval</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">-- Advance to the next self number in the sequence, append it to the output if required, indicate if finished.</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">currentSelf</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">interval</span> |
|||
<span style="color: #000000;">counter</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">counter</span> <span style="color: #0000FF;">>=</span> <span style="color: #000000;">startIndex</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">output</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">currentSelf</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">counter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">endIndex</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #004600;">true</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #004600;">false</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">selfNumbers</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">indexRange</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">startIndex</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">indexRange</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">endIndex</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">indexRange</span><span style="color: #0000FF;">[$]</span> |
|||
<span style="color: #000000;">counter</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
<span style="color: #000000;">currentSelf</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
|||
<span style="color: #000000;">output</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> |
|||
<span style="color: #000080;font-style:italic;">-- Main process. Start with the single-digit odd numbers and first run.</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #000080;font-style:italic;">-- If necessary, fast forward to last self number before the lowest-order block containing first number rqd.</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">counter</span><span style="color: #0000FF;"><</span><span style="color: #000000;">startIndex</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000080;font-style:italic;">-- The highest-order blocks whose ends this handles correctly contain 9,777,777,778 self numbers. |
|||
-- The difference between equivalently positioned numbers in these blocks is 100,000,000,001. |
|||
-- The figures for successively lower-order blocks have successively fewer 7s and 0s!</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">indexDiff</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">9777777778</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">numericDiff</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">100000000001</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">indexDiff</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">98</span> <span style="color: #008080;">and</span> <span style="color: #000000;">counter</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">startIndex</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">counter</span><span style="color: #0000FF;">+</span><span style="color: #000000;">indexDiff</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">startIndex</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">counter</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">indexDiff</span> |
|||
<span style="color: #000000;">currentSelf</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">numericDiff</span> |
|||
<span style="color: #008080;">else</span> |
|||
<span style="color: #000000;">indexDiff</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">indexDiff</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">10</span> <span style="color: #000080;font-style:italic;">-- (..78->80->8)</span> |
|||
<span style="color: #000000;">numericDiff</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">numericDiff</span><span style="color: #0000FF;">+</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">10</span> <span style="color: #000080;font-style:italic;">-- (..01->10->1)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #000080;font-style:italic;">-- Sequencing loop, per lowest-order block.</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000080;font-style:italic;">-- Eight ten-number runs, each at a numeric interval of 2 from the end of the previous one.</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">8</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #000080;font-style:italic;">-- Two shorter runs, the second at an interval inversely related to their length.</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">shorterRunLength</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> |
|||
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">currentSelf</span><span style="color: #0000FF;">/</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">-- Work out a shorter run length based on the most significant digit change about to happen.</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">temp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">9</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">shorterRunLength</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #000000;">temp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">temp</span><span style="color: #0000FF;">/</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">interval</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">2</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">interval</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">doneAfterAdding</span><span style="color: #0000FF;">(</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span><span style="color: #000000;">shorterRunLength</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">output</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #000000;">interval</span> <span style="color: #0000FF;">+=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">9</span><span style="color: #0000FF;">-</span><span style="color: #000000;">shorterRunLength</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">13</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The first 50 self numbers are:\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">selfNumbers</span><span style="color: #0000FF;">({</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">50</span><span style="color: #0000FF;">}),{</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%3d"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntCh</span><span style="color: #0000FF;">,</span><span style="color: #004600;">false</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">=</span><span style="color: #000000;">8</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The %,dth safe number is %,d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">selfNumbers</span><span style="color: #0000FF;">({</span><span style="color: #000000;">n</span><span style="color: #0000FF;">})[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"completed in %s\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">))</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 987: | Line 1,786: | ||
154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323, |
154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323, |
||
334,345,356,367,378,389,400,411,413,424,435,446,457,468} |
334,345,356,367,378,389,400,411,413,424,435,446,457,468} |
||
The 100,000,000th safe number is 1,022,727,208 |
|||
The 1,000,000,000th safe number is 10,227,272,649 |
|||
completed in 0.1s |
|||
</pre> |
|||
=={{header|Python}}== |
|||
i n size time |
|||
{{works with|Python|2.7}} |
|||
100 973 18 0.1s |
|||
<syntaxhighlight lang="python">class DigitSumer : |
|||
1,000 10,188 13 0.2s |
|||
def __init__(self): |
|||
10,000 102,225 10 1.0s |
|||
sumdigit = lambda n : sum( map( int,str( n ))) |
|||
100,000 1,022,675 20 9.3s |
|||
self.t = [sumdigit( i ) for i in xrange( 10000 )] |
|||
1,000,000 10,227,221 17 1 minute and 37s |
|||
def __call__ ( self,n ): |
|||
10,000,000 102,272,662 27 16 minutes and 40s |
|||
r = 0 |
|||
while n >= 10000 : |
|||
n,q = divmod( n,10000 ) |
|||
r += self.t[q] |
|||
return r + self.t[n] |
|||
Estimated time for 100,000,000 :2 hours, 46 minutes and 37s |
|||
</pre> |
|||
For the record, I would not be at all surprised should a translation of this beat 20 minutes (for 1e8) |
|||
def self_numbers (): |
|||
=== sliding sieve === |
|||
d = DigitSumer() |
|||
Mid-speed, perhaps helped by a slightly smarter way of calculating/updating the digit sums<br> |
|||
s = set([]) |
|||
Similar to some other entries, we (only) need a sieve of 9*9, +1 here as I test an entry after the slide. |
|||
i = 1 |
|||
<lang Phix>--sequence sieve = repeat(0,82), -- (~25% slower) |
|||
while 1 : |
|||
sequence sieve = repeat(0,8192), |
|||
n = i + d( i ) |
|||
if i in s : |
|||
integer offset = 0, |
|||
s.discard( i ) |
|||
else: |
|||
yield i |
|||
s.add( n ) |
|||
procedure next_self() |
|||
i += 1 |
|||
n += 1 |
|||
for i=length(digits) to 1 by -1 do |
|||
integer d = digits[i] |
|||
if d!=9 then |
|||
digits[i] = d+1 |
|||
digit_sum += 1 |
|||
exit |
|||
end if |
|||
digits[i] = 0 |
|||
digit_sum -= 9 |
|||
end for |
|||
integer k = n+digit_sum-offset |
|||
if k>length(sieve) then |
|||
integer j = 1 |
|||
for i=n-offset to length(sieve) do |
|||
sieve[j] = sieve[i] |
|||
j += 1 |
|||
end for |
|||
sieve[j..$] = 0 |
|||
offset = n-1 |
|||
k = digit_sum+1 |
|||
end if |
|||
sieve[k] = 1 |
|||
if sieve[n-offset]=0 then exit end if |
|||
end while |
|||
-- (result in n) |
|||
end procedure |
|||
import time |
|||
constant limit = 100000000 |
|||
p = 100 |
|||
atom t0 = time() |
|||
t = time.time() |
|||
printf(1,"The first 50 self numbers are:\n") |
|||
for i,s in enumerate( self_numbers(),1 ): |
|||
integer chk = 100 |
|||
if i <= 50 : |
|||
print s, |
|||
if i |
if i == 50 : print |
||
if i == p : |
|||
printf(1," %3d%s",{n,iff(mod(i,25)=0?"\n":"")}) |
|||
print '%7.1f sec %9dth = %d'%( time.time()-t,i,s ) |
|||
elsif i=chk then |
|||
p *= 10</syntaxhighlight> |
|||
printf(1,"\n i n time\n") |
|||
end if |
|||
printf(1,"%,12d %,13d %s\n",{i,n,elapsed(time()-t0)}) |
|||
chk *= 10 |
|||
end if |
|||
end for |
|||
printf(1,"\nEstimated time for %,d :%s\n",{1e9,elapsed((time()-t0)*1e9/limit)})</lang> |
|||
{{out}} |
{{out}} |
||
<pre>1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
<pre> |
|||
0.0 sec 100th = 973 |
|||
The first 50 self numbers are: |
|||
0.0 sec 1000th = 10188 |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 |
|||
0.1 sec 10000th = 102225 |
|||
222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
1.0 sec 100000th = 1022675 |
|||
11.4 sec 1000000th = 10227221 |
|||
143.4 sec 10000000th = 102272662 |
|||
1812.0 sec 100000000th = 1022727208</pre> |
|||
=={{header|Raku}}== |
|||
i n time |
|||
Translated the low memory version of the Go entry but showed only the first 50 self numbers. The machine for running this task (a Xeon E3110+8GB memory) is showing its age as, 1) took over 6 minutes to complete the Go entry, 2) not even able to run the other two Go alternative entries and 3) needed over 47 minutes to reach 1e6 iterations here. Anyway I will try this on an i5 box later to see how it goes. |
|||
100 973 0.1s |
|||
{{trans|Go}} |
|||
1,000 10,188 0.1s |
|||
<syntaxhighlight lang="raku" line># 20201127 Raku programming solution |
|||
10,000 102,225 0.1s |
|||
100,000 1,022,675 0.1s |
|||
1,000,000 10,227,221 0.5s |
|||
10,000,000 102,272,662 4.5s |
|||
100,000,000 1,022,727,208 44.5s |
|||
my ( $st, $count, $i, $pow, $digits, $offset, $lastSelf, $done, @selfs) = |
|||
Estimated time for 1,000,000,000 :7 minutes and 25s |
|||
now, 0, 1, 10, 1, 9, 0, False; |
|||
</pre> |
|||
# while ( $count < 1e8 ) { |
|||
until $done { |
|||
my $isSelf = True; |
|||
my $sum = (my $start = max ($i-$offset), 0).comb.sum; |
|||
loop ( my $j = $start; $j < $i; $j+=1 ) { |
|||
if $j+$sum == $i { $isSelf = False and last } |
|||
($j+1)%10 != 0 ?? ( $sum+=1 ) !! ( $sum = ($j+1).comb.sum ) ; |
|||
} |
|||
if $isSelf { |
|||
$count+=1; |
|||
$lastSelf = $i; |
|||
if $count ≤ 50 { |
|||
@selfs.append: $i; |
|||
if $count == 50 { |
|||
say "The first 50 self numbers are:\n", @selfs; |
|||
$done = True; |
|||
} |
|||
} |
|||
} |
|||
$i+=1; |
|||
if $i % $pow == 0 { |
|||
$pow *= 10; |
|||
$digits+=1 ; |
|||
$offset = $digits * 9 |
|||
} |
|||
} |
|||
# say "The 100 millionth self number is ", $lastSelf; |
|||
# say "Took ", now - $st, " seconds."</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The first 50 self numbers are: |
|||
[1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468]</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
=== first 50 self numbers === |
=== first 50 self numbers === |
||
< |
<syntaxhighlight lang="rexx">/*REXX program displays N self numbers (aka Colombian or Devlali numbers). OEIS A3052.*/ |
||
parse arg n . /*obtain optional argument from the CL.*/ |
parse arg n . /*obtain optional argument from the CL.*/ |
||
if n=='' | n=="," then n= 50 /*Not specified? Then use the default.*/ |
if n=='' | n=="," then n= 50 /*Not specified? Then use the default.*/ |
||
tell = n>0; n= abs(n) /*TELL: show the self numbers if N>0 */ |
|||
@.= . /*initialize the array of self numbers.*/ |
|||
do j=1 for n*10 /*scan through ten times the #s wanted.*/ |
do j=1 for n*10 /*scan through ten times the #s wanted.*/ |
||
$= j /*1st part of sum is the number itself.*/ |
$= j /*1st part of sum is the number itself.*/ |
||
do k=1 for length(j) |
do k=1 for length(j) /*sum the decimal digits in the number.*/ |
||
$= $ + substr(j, k, 1) |
$= $ + substr(j, k, 1) /*add a particular digit to the sum. */ |
||
end /*k*/ |
end /*k*/ |
||
@.$= /*mark J as not being a self number. */ |
@.$= /*mark J as not being a self number. */ |
||
end /*j*/ |
end /*j*/ /* ─── */ |
||
list= 1 /*initialize list to the |
list= 1 /*initialize the list to the 1st number*/ |
||
#= 1 /*the count of self numbers (so far). */ |
|||
do i=3 until #==n; if @.i=='' then iterate /*Not a self number? Then skip it. */ |
|||
#= # + 1; list= list i /*bump counter of self #'s; add to list*/ |
|||
end /*i*/ |
|||
#= # + 1 /*bump the counter of self numbers. */ |
|||
/*stick a fork in it, we're all done. */ |
|||
say n " self numbers were found." /*display the title for the output list*/ |
|||
end /*i*/ |
|||
if tell then say list /*display list of self numbers ──►term.*/</syntaxhighlight> |
|||
say 'the first ' n " self numbers are:" /*display the title for the output list*/ |
|||
say list /*display list of self numbers ──►term.*/ |
|||
exit 0 /*stick a fork in it, we're all done. */</lang> |
|||
{{out|output|text= when using the default input:}} |
{{out|output|text= when using the default input:}} |
||
<pre> |
<pre> |
||
50 self numbers were found. |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
||
</pre> |
</pre> |
||
Line 1,107: | Line 1,909: | ||
=== ten millionth self number === |
=== ten millionth self number === |
||
{{trans|Go (low memory)}} |
{{trans|Go (low memory)}} |
||
< |
<syntaxhighlight lang="rexx">/*REXX pgm displays the Nth self number, aka Colombian or Devlali numbers. OEIS A3052.*/ |
||
numeric digits 20 /*ensure enough decimal digits for #. */ |
numeric digits 20 /*ensure enough decimal digits for #. */ |
||
parse arg high . /*obtain optional argument from the CL.*/ |
parse arg high . /*obtain optional argument from the CL.*/ |
||
if high=='' | high=="," then high= |
if high=='' | high=="," then high= 10000000 /*Not specified? Then use 10M default.*/ |
||
i= 1; pow= 10; digs= 1; offset= 9; $= 0 /*$: the last self number found. */ |
i= 1; pow= 10; digs= 1; offset= 9; $= 0 /*$: the last self number found. */ |
||
#= 0 /*count of self numbers (so far). */ |
#= 0 /*count of self numbers (so far). */ |
||
do while #<high; isSelf= 1 /*assume a self number (so far). */ |
do while #<high; isSelf= 1 /*assume a self number (so far). */ |
||
start= max(i-offset, 0) |
start= max(i-offset, 0) /*find start #; must be non-negative. */ |
||
sum= sumDigs(start) |
sum= sumDigs(start) /*obtain the sum of the decimal digits.*/ |
||
do j=start to i-1 |
do j=start to i-1 |
||
Line 1,121: | Line 1,923: | ||
iterate /*keep looking for more self numbers. */ |
iterate /*keep looking for more self numbers. */ |
||
end |
end |
||
if (j+1)//10==0 then sum= sumDigs(j+1) /*obtain the sum of the decimal digits.*/ |
|||
jp= j + 1 /*shortcut variable for next statement.*/ |
|||
else sum= sum + 1 /*bump " " " " " " */ |
|||
if jp//10==0 then sum= sumDigs(jp) |
|||
else sum= sum + 1 |
|||
end /*j*/ |
end /*j*/ |
||
Line 1,129: | Line 1,930: | ||
$= i /*save the last self number found. */ |
$= i /*save the last self number found. */ |
||
end |
end |
||
i= i + 1 /*bump the self number by unity. */ |
|||
i= i + 1 |
|||
if i//pow==0 then do; |
if i//pow==0 then do; digs= digs + 1 /* " " number of decimal digits. */ |
||
pow= pow * 10 /* " " power " a factor of ten. */ |
|||
offset= digs * 9 /* |
offset= digs * 9 /* " " offset " " " " nine. */ |
||
end |
end |
||
end /*while*/ |
end /*while*/ |
||
Line 1,144: | Line 1,945: | ||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
commas: parse arg _; do c=length(_)-3 to 1 by -3; _=insert(',', _, c); end; return _ |
commas: parse arg _; do c=length(_)-3 to 1 by -3; _=insert(',', _, c); end; return _ |
||
th: parse arg th; |
th: parse arg th; return word('th st nd rd', 1 +(th//10)*(th//100%10\==1)*(th//10<4))</syntaxhighlight> |
||
{{out|output|text= when using the default input:}} |
{{out|output|text= when using the default input:}} |
||
<pre> |
<pre> |
||
the 100,000,000th self number is: 1,022,727,208 |
the 100,000,000th self number is: 1,022,727,208 |
||
</pre> |
|||
=={{header|Ring}}== |
|||
<syntaxhighlight lang="ring"> |
|||
load "stdlib.ring" |
|||
see "working..." + nl |
|||
see "The first 50 self numbers are:" + nl |
|||
n = 0 |
|||
num = 0 |
|||
limit = 51 |
|||
limit2 = 10000000 |
|||
while true |
|||
n = n + 1 |
|||
for m = 1 to n |
|||
flag = 1 |
|||
sum = 0 |
|||
strnum = string(m) |
|||
for p = 1 to len(strnum) |
|||
sum = sum + number(strnum[p]) |
|||
next |
|||
sum2 = m + sum |
|||
if sum2 = n |
|||
flag = 0 |
|||
exit |
|||
else |
|||
flag = 1 |
|||
ok |
|||
next |
|||
if flag = 1 |
|||
num = num + 1 |
|||
if num < limit |
|||
see "" + num + ". " + n + nl |
|||
ok |
|||
if num = limit2 |
|||
see "The " + limit2 + "th self number is: " + n + nl |
|||
ok |
|||
if num > limit2 |
|||
exit |
|||
ok |
|||
ok |
|||
end |
|||
see "done..." + nl |
|||
</syntaxhighlight> |
|||
Output: |
|||
<pre> |
|||
working... |
|||
The first 50 self numbers are: |
|||
1. 1 |
|||
2. 3 |
|||
3. 5 |
|||
4. 7 |
|||
5. 9 |
|||
6. 20 |
|||
7. 31 |
|||
8. 42 |
|||
9. 53 |
|||
10. 64 |
|||
11. 75 |
|||
12. 86 |
|||
13. 97 |
|||
14. 108 |
|||
15. 110 |
|||
16. 121 |
|||
17. 132 |
|||
18. 143 |
|||
19. 154 |
|||
20. 165 |
|||
21. 176 |
|||
22. 187 |
|||
23. 198 |
|||
24. 209 |
|||
25. 211 |
|||
26. 222 |
|||
27. 233 |
|||
28. 244 |
|||
29. 255 |
|||
30. 266 |
|||
31. 277 |
|||
32. 288 |
|||
33. 299 |
|||
34. 310 |
|||
35. 312 |
|||
36. 323 |
|||
37. 334 |
|||
38. 345 |
|||
39. 356 |
|||
40. 367 |
|||
41. 378 |
|||
42. 389 |
|||
43. 400 |
|||
44. 411 |
|||
45. 413 |
|||
46. 424 |
|||
47. 435 |
|||
48. 446 |
|||
49. 457 |
|||
50. 468 |
|||
The 10000000th self number is: 1022727208 |
|||
done... |
|||
</pre> |
|||
=={{header|RPL}}== |
|||
====Brute force==== |
|||
Using a sieve: |
|||
« 0 |
|||
'''WHILE''' OVER '''REPEAT''' |
|||
SWAP 10 MOD LASTARG / IP |
|||
ROT ROT + |
|||
'''END''' + |
|||
» '<span style="color:blue">DIGSUM</span>' STO |
|||
« 500 0 → max n |
|||
« { } DUP max + 0 CON |
|||
1 CF |
|||
'''DO''' 'n' INCR |
|||
DUP <span style="color:blue">DIGSUM</span> + |
|||
'''IFERR''' 1 PUT '''THEN''' DROP2 1 SF '''END''' |
|||
'''UNTIL''' 1 FS? '''END''' |
|||
1 |
|||
'''WHILE''' 3 PICK SIZE 50 < '''REPEAT''' |
|||
'''IF''' DUP2 GET NOT '''THEN''' ROT OVER + ROT ROT '''END''' |
|||
1 + |
|||
'''END''' DROP2 |
|||
» » '<span style="color:blue">TASK</span>' STO |
|||
{{out}} |
|||
<pre> |
|||
1: { 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 } |
|||
</pre> |
|||
Runs in 2 minutes 8 seconds on a HP-48SX. |
|||
====Maximilian F. Hasler's algorithm==== |
|||
Translation of the PARI/GP formula on the OEIS page: |
|||
« → n |
|||
« { } 1 SF |
|||
1 n 2 / IP n XPON 1 + 9 * MIN '''FOR''' j |
|||
'''IF''' n j - <span style="color:blue">DIGSUM</span> j == '''THEN''' 1 CF n 'j' STO '''END''' |
|||
'''NEXT''' |
|||
1 FS? |
|||
» '<span style="color:blue">SELF?</span>' STO |
|||
« { } |
|||
'''WHILE''' OVER SIZE 50 < REPEAT |
|||
'''IF''' DUP <span style="color:blue">SELF?</span> '''THEN''' SWAP OVER + SWAP '''END''' |
|||
1 + |
|||
'''END''' DROP |
|||
» '<span style="color:blue">TASK</span>' STO |
|||
Same result as above. No need for sieve, but much slower: 10 minutes 52 seconds on a HP-48SX. |
|||
=={{header|Sidef}}== |
|||
Algorithm by David A. Corneth (OEIS: A003052). |
|||
<syntaxhighlight lang="ruby">func is_self_number(n) { |
|||
if (n < 30) { |
|||
return (((n < 10) && (n.is_odd)) || (n == 20)) |
|||
} |
|||
var qd = (1 + n.ilog10) |
|||
var r = (1 + (n-1)%9) |
|||
var h = (r + 9*(r%2))/2 |
|||
var ld = 10 |
|||
while (h + 9*qd >= n%ld) { |
|||
ld *= 10 |
|||
} |
|||
var vs = idiv(n, ld).sumdigits |
|||
n %= ld |
|||
0..qd -> none { |i| |
|||
vs + sumdigits(n - h - 9*i) == (h + 9*i) |
|||
} |
|||
} |
|||
say is_self_number.first(50).join(' ')</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
</pre> |
|||
Simpler algorithm (by M. F. Hasler): |
|||
<syntaxhighlight lang="ruby">func is_self_number(n) { |
|||
1..min(n>>1, 9*n.len) -> none {|i| sumdigits(n-i) == i } && (n > 0) |
|||
}</syntaxhighlight> |
|||
=={{header|Standard ML}}== |
|||
<syntaxhighlight lang="ocaml"> |
|||
open List; |
|||
val rec selfNumberNr = fn NR => |
|||
let |
|||
val rec sumdgt = fn 0 => 0 | n => Int.rem (n, 10) + sumdgt (Int.quot(n ,10)); |
|||
val rec isSelf = fn ([],l1,l2) => [] |
|||
| (x::tt,l1,l2) => if exists (fn i=>i=x) l1 orelse exists (fn i=>i=x) l2 |
|||
then ( isSelf (tt,l1,l2)) else x::isSelf (tt,l1,l2) ; |
|||
val rec partcount = fn (n, listIn , count , selfs) => |
|||
if count >= NR then nth (selfs, length selfs + NR - count -1) |
|||
else |
|||
let |
|||
val range = tabulate (81 , fn i => 81*n +i+1) ; |
|||
val listOut = map (fn i => i + sumdgt i ) range ; |
|||
val selfs = isSelf (range,listIn,listOut) |
|||
in |
|||
partcount ( n+1 , listOut , count+length (selfs) , selfs ) |
|||
end; |
|||
in |
|||
partcount (0,[],0,[]) |
|||
end; |
|||
app ((fn s => print (s ^ " ")) o Int.toString o selfNumberNr) (tabulate (50,fn i=>i+1)) ; |
|||
selfNumberNr 100000000 ; |
|||
</syntaxhighlight> |
|||
output |
|||
<pre> |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
1022727208 |
|||
</pre> |
</pre> |
||
Line 1,157: | Line 2,178: | ||
Unsurprisingly, very slow compared to the Go version as Wren is interpreted and uses floating point arithmetic for all numerical work. |
Unsurprisingly, very slow compared to the Go version as Wren is interpreted and uses floating point arithmetic for all numerical work. |
||
< |
<syntaxhighlight lang="wren">var sieve = Fn.new { |
||
var sv = List.filled(2*1e9+9*9+1, false) |
var sv = List.filled(2*1e9+9*9+1, false) |
||
var n = 0 |
var n = 0 |
||
Line 1,208: | Line 2,229: | ||
} |
} |
||
} |
} |
||
System.print("Took %(System.clock-st) seconds.")</ |
System.print("Took %(System.clock-st) seconds.")</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,217: | Line 2,238: | ||
The 100 millionth self number is 1022727208 |
The 100 millionth self number is 1022727208 |
||
Took 222.789713 seconds. |
Took 222.789713 seconds. |
||
</pre> |
|||
=={{header|XPL0}}== |
|||
{{trans|Go}} |
|||
As run on Raspberry Pi4. |
|||
<syntaxhighlight lang "XPL0">def LenSV = 2_000_000_000 + 9*9 + 1; |
|||
func Sieve; |
|||
char SV; |
|||
int N, S(8), A, B, C, D, E, F, G, H, I, J; |
|||
[SV:= MAlloc(LenSV); |
|||
N:= 0; |
|||
for A:= 0 to 1 do |
|||
[for B:= 0 to 9 do |
|||
[S(0):= A + B; |
|||
for C:= 0 to 9 do |
|||
[S(1):= S(0) + C; |
|||
for D:= 0 to 9 do |
|||
[S(2):= S(1) + D; |
|||
for E:= 0 to 9 do |
|||
[S(3):= S(2) + E; |
|||
for F:= 0 to 9 do |
|||
[S(4):= S(3) + F; |
|||
for G:= 0 to 9 do |
|||
[S(5):= S(4) + G; |
|||
for H:= 0 to 9 do |
|||
[S(6):= S(5) + H; |
|||
for I:= 0 to 9 do |
|||
[S(7):= S(6) + I; |
|||
for J:= 0 to 9 do |
|||
[SV(S(7)+J+N):= true; |
|||
N:= N+1; |
|||
] |
|||
] |
|||
] |
|||
] |
|||
] |
|||
] |
|||
] |
|||
] |
|||
] |
|||
]; |
|||
return SV; |
|||
]; |
|||
char SV; |
|||
int ST, Count, I; |
|||
[ST:= GetTime; |
|||
SV:= Sieve; |
|||
Count:= 0; |
|||
Text(0, "The first 50 self numbers are:^m^j"); |
|||
for I:= 0 to LenSV-1 do |
|||
[if SV(I) = false then |
|||
[Count:= Count+1; |
|||
if Count <= 50 then |
|||
[IntOut(0, I); ChOut(0, ^ )]; |
|||
if Count = 100_000_000 then |
|||
[Text(0, "^m^j^m^jThe 100 millionth self number is "); |
|||
IntOut(0, I); |
|||
CrLf(0); |
|||
I:= LenSV; |
|||
]; |
|||
]; |
|||
]; |
|||
Text(0, "Took "); RlOut(0, float(GetTime-ST) / 1e6); CrLf(0); |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
The first 50 self numbers are: |
|||
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 |
|||
The 100 millionth self number is 1022727208 |
|||
Took 29.46575 |
|||
</pre> |
</pre> |
Revision as of 10:55, 4 February 2024
You are encouraged to solve this task according to the task description, using any language you may know.
A number n is a self number if there is no number g such that g + the sum of g's digits = n. So 18 is not a self number because 9+9=18, 43 is not a self number because 35+5+3=43.
The task is:
Display the first 50 self numbers; I believe that the 100000000th self number is 1022727208. You should either confirm or dispute my conjecture.
224036583-1 is a Mersenne prime, claimed to also be a self number. Extra credit to anyone proving it.
- See also
ALGOL 68
BEGIN # find some self numbers numbers n such that there is no g such that g + sum of g's digits = n #
INT max number = 1 999 999 999 + 82; # maximum n plus g we will condifer #
# sieve the self numbers up to 1 999 999 999 #
[ 0 : max number ]BOOL self; FOR i TO UPB self DO self[ i ] := TRUE OD;
INT n := 0;
FOR s0 FROM 0 TO 1 DO
FOR d1 FROM 0 TO 9 DO
INT s1 = s0 + d1;
FOR d2 FROM 0 TO 9 DO
INT s2 = s1 + d2;
FOR d3 FROM 0 TO 9 DO
INT s3 = s2 + d3;
FOR d4 FROM 0 TO 9 DO
INT s4 = s3 + d4;
FOR d5 FROM 0 TO 9 DO
INT s5 = s4 + d5;
FOR d6 FROM 0 TO 9 DO
INT s6 = s5 + d6;
FOR d7 FROM 0 TO 9 DO
INT s7 = s6 + d7;
FOR d8 FROM 0 TO 9 DO
INT s8 = s7 + d8;
FOR d9 FROM 0 TO 9 DO
INT s9 = s8 + d9;
self[ s9 + n ] := FALSE;
n +:= 1
OD
OD
OD
OD
OD
OD
OD
OD
OD
OD;
# show the first 50 self numbers #
INT s count := 0;
FOR i TO UPB self WHILE s count < 50 DO
IF self[ i ] THEN
print( ( " ", whole( i, -3 ) ) );
IF ( s count +:= 1 ) MOD 18 = 0 THEN print( ( newline ) ) FI
FI
OD;
print( ( newline ) );
# show the self numbers with power-of-10 indxes #
INT s show := 1;
s count := 0;
print( ( " nth self", newline ) );
print( ( " n number", newline ) );
FOR i TO UPB self DO
IF self[ i ] THEN
s count +:= 1;
IF s count = s show THEN
print( ( whole( s show, -9 ), " ", whole( i, -11 ), newline ) );
s show *:= 10
FI
FI
OD
END
- Output:
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 nth self n number 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208
AppleScript
I couldn't follow the math in the Wikipedia entry, nor the discussion and code here so far. But an initial expedient of generating a list of all the integers from 1 to just over ten times the required number of results and then deleting those that could be derived by the described method revealed the sequencing pattern on which the code below is based. On the test machine, it completes all three of the tests at the bottom in a total of around a millisecond.
(*
Base-10 self numbers by index (single or range).
Follows an observed sequence pattern whereby, after the initial single-digit odd numbers, self numbers are
grouped in runs whose members occur at numeric intervals of 11. Runs after the first one come in blocks of
ten: eight runs of ten numbers followed by two shorter runs. The numeric interval between runs is usually 2,
but that between shorter runs, and their length, depend on the highest-order digit change occurring in them.
This connection with significant digit change means every ten blocks form a higher-order block, every ten
of these a higher-order-still block, and so on.
The code below appears to be good up to the last self number before 10^12 — ie. 999,999,999,997, which is
returned as the 97,777,777,792nd such number. After this, instead of zero-length shorter runs, the actual
pattern apparently starts again with a single run of 10, like the one at the beginning.
*)
on selfNumbers(indexRange)
set indexRange to indexRange as list
-- Script object with subhandlers and associated properties.
script |subscript|
property startIndex : beginning of indexRange
property endIndex : end of indexRange
property counter : 0
property currentSelf : -1
property output : {}
-- Advance to the next self number in the sequence, append it to the output if required, indicate if finished.
on doneAfterAdding(interval)
set currentSelf to currentSelf + interval
set counter to counter + 1
if (counter < startIndex) then return false
set end of my output to currentSelf
return (counter = endIndex)
end doneAfterAdding
-- If necessary, fast forward to the last self number before the lowest-order block containing the first number required.
on fastForward()
if (counter ≥ startIndex) then return
-- The highest-order blocks whose ends this script handles correctly contain 9,777,777,778 self numbers.
-- The difference between equivalently positioned numbers in these blocks is 100,000,000,001.
-- The figures for successively lower-order blocks have successively fewer 7s and 0s!
set indexDiff to 9.777777778E+9
set numericDiff to 1.00000000001E+11
repeat until ((indexDiff < 98) or (counter = startIndex))
set test to counter + indexDiff
if (test < startIndex) then
set counter to test
set currentSelf to (currentSelf + numericDiff)
else
set indexDiff to (indexDiff + 2) div 10
set numericDiff to numericDiff div 10 + 1
end if
end repeat
end fastForward
-- Work out a shorter run length based on the most significant digit change about to happen.
on getShorterRunLength()
set shorterRunLength to 9
set temp to (|subscript|'s currentSelf) div 1000
repeat while (temp mod 10 is 9)
set shorterRunLength to shorterRunLength - 1
set temp to temp div 10
end repeat
return shorterRunLength
end getShorterRunLength
end script
-- Main process. Start with the single-digit odd numbers and first run.
repeat 5 times
if (|subscript|'s doneAfterAdding(2)) then return |subscript|'s output
end repeat
repeat 9 times
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output
end repeat
-- Fast forward if the start index hasn't yet been reached.
tell |subscript| to fastForward()
-- Sequencing loop, per lowest-order block.
repeat
-- Eight ten-number runs, each at a numeric interval of 2 from the end of the previous one.
repeat 8 times
if (|subscript|'s doneAfterAdding(2)) then return |subscript|'s output
repeat 9 times
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output
end repeat
end repeat
-- Two shorter runs, the second at an interval inversely related to their length.
set shorterRunLength to |subscript|'s getShorterRunLength()
repeat with interval in {2, 2 + (10 - shorterRunLength) * 13}
if (|subscript|'s doneAfterAdding(interval)) then return |subscript|'s output
repeat (shorterRunLength - 1) times
if (|subscript|'s doneAfterAdding(11)) then return |subscript|'s output
end repeat
end repeat
end repeat
end selfNumbers
-- Demo calls:
-- First to fiftieth self numbers.
selfNumbers({1, 50})
--> {1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 110, 121, 132, 143, 154, 165, 176, 187, 198, 209, 211, 222, 233, 244, 255, 266, 277, 288, 299, 310, 312, 323, 334, 345, 356, 367, 378, 389, 400, 411, 413, 424, 435, 446, 457, 468}
-- One hundred millionth:
selfNumbers(100000000)
--> {1.022727208E+9}
-- 97,777,777,792nd:
selfNumbers(9.7777777792E+10)
--> {9.99999999997E+11}
AWK
# syntax: GAWK -f SELF_NUMBERS.AWK
# converted from Go (low memory example)
BEGIN {
print("HH:MM:SS INDEX SELF")
print("-------- ---------- ----------")
count = 0
digits = 1
i = 1
last_self = 0
offset = 9
pow = 10
while (count < 1E8) {
is_self = 1
start = max(i-offset,0)
sum = sum_digits(start)
for (j=start; j<i; j++) {
if (j + sum == i) {
is_self = 0
break
}
sum = ((j+1) % 10 != 0) ? ++sum : sum_digits(j+1)
}
if (is_self) {
last_self = i
if (++count <= 50) {
selfs = selfs i " "
}
}
if (++i % pow == 0) {
pow *= 10
digits++
offset = digits * 9
}
if (count ~ /^10*$/ && arr[count]++ == 0) {
printf("%8s %10s %10s\n",strftime("%H:%M:%S"),count,last_self)
}
}
printf("\nfirst 50 self numbers:\n%s\n",selfs)
exit(0)
}
function sum_digits(x, sum,y) {
while (x) {
y = x % 10
sum += y
x = int(x/10)
}
return(sum)
}
function max(x,y) { return((x > y) ? x : y) }
- Output:
HH:MM:SS INDEX SELF -------- ---------- ---------- 00:36:35 1 1 00:36:35 10 64 00:36:35 100 973 00:36:35 1000 10188 00:36:36 10000 102225 00:36:46 100000 1022675 00:38:49 1000000 10227221 01:03:01 10000000 102272662 05:27:35 100000000 1022727208 first 50 self numbers: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468
C
Sieve based
About 25% faster than Go (using GCC 7.5.0 -O3) mainly due to being able to iterate through the sieve using a pointer.
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef unsigned char bool;
#define TRUE 1
#define FALSE 0
#define MILLION 1000000
#define BILLION 1000 * MILLION
#define MAX_COUNT 2*BILLION + 9*9 + 1
void sieve(bool *sv) {
int n = 0, s[8], a, b, c, d, e, f, g, h, i, j;
for (a = 0; a < 2; ++a) {
for (b = 0; b < 10; ++b) {
s[0] = a + b;
for (c = 0; c < 10; ++c) {
s[1] = s[0] + c;
for (d = 0; d < 10; ++d) {
s[2] = s[1] + d;
for (e = 0; e < 10; ++e) {
s[3] = s[2] + e;
for (f = 0; f < 10; ++f) {
s[4] = s[3] + f;
for (g = 0; g < 10; ++g) {
s[5] = s[4] + g;
for (h = 0; h < 10; ++h) {
s[6] = s[5] + h;
for (i = 0; i < 10; ++i) {
s[7] = s[6] + i;
for (j = 0; j < 10; ++j) {
sv[s[7] + j+ n++] = TRUE;
}
}
}
}
}
}
}
}
}
}
}
int main() {
int count = 0;
clock_t begin = clock();
bool *p, *sv = (bool*) calloc(MAX_COUNT, sizeof(bool));
sieve(sv);
printf("The first 50 self numbers are:\n");
for (p = sv; p < sv + MAX_COUNT; ++p) {
if (!*p) {
if (++count <= 50) printf("%ld ", p-sv);
if (count == 100 * MILLION) {
printf("\n\nThe 100 millionth self number is %ld\n", p-sv);
break;
}
}
}
free(sv);
printf("Took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC);
return 0;
}
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 The 100 millionth self number is 1022727208 Took 1.521486 seconds.
Extended
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef unsigned char bool;
#define TRUE 1
#define FALSE 0
#define MILLION 1000000LL
#define BILLION 1000 * MILLION
#define MAX_COUNT 103LL*10000*10000 + 11*9 + 1
int digitSum[10000];
void init() {
int i = 9999, s, t, a, b, c, d;
for (a = 9; a >= 0; --a) {
for (b = 9; b >= 0; --b) {
s = a + b;
for (c = 9; c >= 0; --c) {
t = s + c;
for (d = 9; d >= 0; --d) {
digitSum[i] = t + d;
--i;
}
}
}
}
}
void sieve(bool *sv) {
int a, b, c;
long long s, n = 0;
for (a = 0; a < 103; ++a) {
for (b = 0; b < 10000; ++b) {
s = digitSum[a] + digitSum[b] + n;
for (c = 0; c < 10000; ++c) {
sv[digitSum[c]+s] = TRUE;
++s;
}
n += 10000;
}
}
}
int main() {
long long count = 0, limit = 1;
clock_t begin = clock(), end;
bool *p, *sv = (bool*) calloc(MAX_COUNT, sizeof(bool));
init();
sieve(sv);
printf("Sieving took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC);
printf("\nThe first 50 self numbers are:\n");
for (p = sv; p < sv + MAX_COUNT; ++p) {
if (!*p) {
if (++count <= 50) {
printf("%ld ", p-sv);
} else {
printf("\n\n Index Self number\n");
break;
}
}
}
count = 0;
for (p = sv; p < sv + MAX_COUNT; ++p) {
if (!*p) {
if (++count == limit) {
printf("%10lld %11ld\n", count, p-sv);
limit *= 10;
if (limit == 10 * BILLION) break;
}
}
}
free(sv);
printf("\nOverall took %lf seconds.\n", (double)(clock() - begin) / CLOCKS_PER_SEC);
return 0;
}
- Output:
Sieving took 7.429969 seconds. The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 1000000000 10227272649 Overall took 11.574158 seconds.
C++
#include <array>
#include <iomanip>
#include <iostream>
const int MC = 103 * 1000 * 10000 + 11 * 9 + 1;
std::array<bool, MC + 1> SV;
void sieve() {
std::array<int, 10000> dS;
for (int a = 9, i = 9999; a >= 0; a--) {
for (int b = 9; b >= 0; b--) {
for (int c = 9, s = a + b; c >= 0; c--) {
for (int d = 9, t = s + c; d >= 0; d--) {
dS[i--] = t + d;
}
}
}
}
for (int a = 0, n = 0; a < 103; a++) {
for (int b = 0, d = dS[a]; b < 1000; b++, n += 10000) {
for (int c = 0, s = d + dS[b] + n; c < 10000; c++) {
SV[dS[c] + s++] = true;
}
}
}
}
int main() {
sieve();
std::cout << "The first 50 self numbers are:\n";
for (int i = 0, count = 0; count <= 50; i++) {
if (!SV[i]) {
count++;
if (count <= 50) {
std::cout << i << ' ';
} else {
std::cout << "\n\n Index Self number\n";
}
}
}
for (int i = 0, limit = 1, count = 0; i < MC; i++) {
if (!SV[i]) {
if (++count == limit) {
//System.out.printf("%,12d %,13d%n", count, i);
std::cout << std::setw(12) << count << " " << std::setw(13) << i << '\n';
limit *= 10;
}
}
}
return 0;
}
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208
C#
via
(third version) Stripped down, as C# limits the size of an array to Int32.MaxValue, so the sieve isn't large enough to hit the 1,000,000,000th value.
using System;
using static System.Console;
class Program {
const int mc = 103 * 1000 * 10000 + 11 * 9 + 1;
static bool[] sv = new bool[mc + 1];
static void sieve() { int[] dS = new int[10000];
for (int a = 9, i = 9999; a >= 0; a--)
for (int b = 9; b >= 0; b--)
for (int c = 9, s = a + b; c >= 0; c--)
for (int d = 9, t = s + c; d >= 0; d--)
dS[i--] = t + d;
for (int a = 0, n = 0; a < 103; a++)
for (int b = 0, d = dS[a]; b < 1000; b++, n += 10000)
for (int c = 0, s = d + dS[b] + n; c < 10000; c++)
sv[dS[c] + s++] = true; }
static void Main() { DateTime st = DateTime.Now; sieve();
WriteLine("Sieving took {0}s", (DateTime.Now - st).TotalSeconds);
WriteLine("\nThe first 50 self numbers are:");
for (int i = 0, count = 0; count <= 50; i++) if (!sv[i]) {
count++; if (count <= 50) Write("{0} ", i);
else WriteLine("\n\n Index Self number"); }
for (int i = 0, limit = 1, count = 0; i < mc; i++)
if (!sv[i]) if (++count == limit) {
WriteLine("{0,12:n0} {1,13:n0}", count, i);
if (limit == 1e9) break; limit *= 10; }
WriteLine("\nOverall took {0}s", (DateTime.Now - st). TotalSeconds);
}
}
- Output:
Timing from tio.run
Sieving took 3.4266187s The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1,000 10,188 10,000 102,225 100,000 1,022,675 1,000,000 10,227,221 10,000,000 102,272,662 100,000,000 1,022,727,208 Overall took 7.0237244s
Elixir
defmodule SelfNums do
def digAndSum(number) when is_number(number) do
Integer.digits(number) |>
Enum.reduce( 0, fn(num, acc) -> num + acc end ) |>
(fn(x) -> x + number end).()
end
def selfFilter(list, firstNth) do
numRange = Enum.to_list 1..firstNth
numRange -- list
end
end
defmodule SelfTest do
import SelfNums
stop = 1000
Enum.to_list 1..stop |>
Enum.map(&digAndSum/1) |>
SelfNums.selfFilter(stop) |>
Enum.take(50) |>
IO.inspect
end
- Output:
[1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 110, 121, 132, 143, 154, 165, 176, 187, 198, 209, 211, 222, 233, 244, 255, 266, 277, 288, 299, 310, 312, 323, 334, 345, 356, 367, 378, 389, 400, 411, 413, 424, 435, 446, 457, 468]
F#
// Self numbers. Nigel Galloway: October 6th., 2020
let fN g=let rec fG n g=match n/10 with 0->n+g |i->fG i (g+(n%10)) in fG g g
let Self=let rec Self n i g=seq{let g=g@([n..i]|>List.map fN) in yield! List.except g [n..i]; yield! Self (n+100) (i+100) (List.filter(fun n->n>i) g)} in Self 0 99 []
Self |> Seq.take 50 |> Seq.iter(printf "%d "); printfn ""
printfn "\n%d" (Seq.item 99999999 Self)
- Output:
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 1022727208
FreeBASIC
Print "The first 50 self numbers are:"
Dim As Boolean flag
Dim As Ulong m, p, sum, number(), sum2
Dim As Ulong n = 0
Dim As Ulong num = 0
Dim As Ulong limit = 51
Dim As Ulong limit2 = 100000000
Dim As String strnum
Do
n += 1
For m = 1 To n
flag = True
sum = 0
strnum = Str(m)
For p = 1 To Len(strnum)
sum += Val(Mid(strnum,p,1))
Next p
sum2 = m + sum
If sum2 = n Then
flag = False
Exit For
Else
flag = True
End If
Next m
If flag = True Then
num += 1
If num < limit Then Print ""; num; ". "; n
If num >= limit2 Then
Print "The "; limit2; "th self number is: "; n
Exit Do
End If
End If
Loop
Sleep
Go
Low memory
Simple-minded, uses very little memory (no sieve) but slow - over 2.5 minutes.
package main
import (
"fmt"
"time"
)
func sumDigits(n int) int {
sum := 0
for n > 0 {
sum += n % 10
n /= 10
}
return sum
}
func max(x, y int) int {
if x > y {
return x
}
return y
}
func main() {
st := time.Now()
count := 0
var selfs []int
i := 1
pow := 10
digits := 1
offset := 9
lastSelf := 0
for count < 1e8 {
isSelf := true
start := max(i-offset, 0)
sum := sumDigits(start)
for j := start; j < i; j++ {
if j+sum == i {
isSelf = false
break
}
if (j+1)%10 != 0 {
sum++
} else {
sum = sumDigits(j + 1)
}
}
if isSelf {
count++
lastSelf = i
if count <= 50 {
selfs = append(selfs, i)
if count == 50 {
fmt.Println("The first 50 self numbers are:")
fmt.Println(selfs)
}
}
}
i++
if i%pow == 0 {
pow *= 10
digits++
offset = digits * 9
}
}
fmt.Println("\nThe 100 millionth self number is", lastSelf)
fmt.Println("Took", time.Since(st))
}
- Output:
The first 50 self numbers are: [1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468] The 100 millionth self number is 1022727208 Took 2m35.531949399s
Sieve based
Simple sieve, requires a lot of memory but quick - around 2 seconds.
Nested 'for's used rather than a recursive function for extra speed.
Have also incorporated Enter your username's suggestion (see Talk page) of using partial sums for each loop which improves performance by about 25%.
package main
import (
"fmt"
"time"
)
func sieve() []bool {
sv := make([]bool, 2*1e9+9*9 + 1)
n := 0
var s [8]int
for a := 0; a < 2; a++ {
for b := 0; b < 10; b++ {
s[0] = a + b
for c := 0; c < 10; c++ {
s[1] = s[0] + c
for d := 0; d < 10; d++ {
s[2] = s[1] + d
for e := 0; e < 10; e++ {
s[3] = s[2] + e
for f := 0; f < 10; f++ {
s[4] = s[3] + f
for g := 0; g < 10; g++ {
s[5] = s[4] + g
for h := 0; h < 10; h++ {
s[6] = s[5] + h
for i := 0; i < 10; i++ {
s[7] = s[6] + i
for j := 0; j < 10; j++ {
sv[s[7]+j+n] = true
n++
}
}
}
}
}
}
}
}
}
}
return sv
}
func main() {
st := time.Now()
sv := sieve()
count := 0
fmt.Println("The first 50 self numbers are:")
for i := 0; i < len(sv); i++ {
if !sv[i] {
count++
if count <= 50 {
fmt.Printf("%d ", i)
}
if count == 1e8 {
fmt.Println("\n\nThe 100 millionth self number is", i)
break
}
}
}
fmt.Println("Took", time.Since(st))
}
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 The 100 millionth self number is 1022727208 Took 1.984969034s
Extended
This uses horst.h's ideas (see Talk page) to find up to the 1 billionth self number in a reasonable time and using less memory than the simple 'sieve based' approach above would have needed.
package main
import (
"fmt"
"time"
)
const MAX_COUNT = 103*1e4*1e4 + 11*9 + 1
var sv = make([]bool, MAX_COUNT+1)
var digitSum = make([]int, 1e4)
func init() {
i := 9999
var s, t int
for a := 9; a >= 0; a-- {
for b := 9; b >= 0; b-- {
s = a + b
for c := 9; c >= 0; c-- {
t = s + c
for d := 9; d >= 0; d-- {
digitSum[i] = t + d
i--
}
}
}
}
}
func sieve() {
n := 0
for a := 0; a < 103; a++ {
for b := 0; b < 1e4; b++ {
s := digitSum[a] + digitSum[b] + n
for c := 0; c < 1e4; c++ {
sv[digitSum[c]+s] = true
s++
}
n += 1e4
}
}
}
func main() {
st := time.Now()
sieve()
fmt.Println("Sieving took", time.Since(st))
count := 0
fmt.Println("\nThe first 50 self numbers are:")
for i := 0; i < len(sv); i++ {
if !sv[i] {
count++
if count <= 50 {
fmt.Printf("%d ", i)
} else {
fmt.Println("\n\n Index Self number")
break
}
}
}
count = 0
limit := 1
for i := 0; i < len(sv); i++ {
if !sv[i] {
count++
if count == limit {
fmt.Printf("%10d %11d\n", count, i)
limit *= 10
if limit == 1e10 {
break
}
}
}
}
fmt.Println("\nOverall took", time.Since(st))
}
- Output:
Sieving took 8.286841692s The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 1000000000 10227272649 Overall took 14.647314803s
Haskell
The solution is quite straightforward. The length of the foreseeing window in filtering procedure (81) is chosen empirically and doesn't have any theoretical background.
import Control.Monad (forM_)
import Text.Printf
selfs :: [Integer]
selfs = sieve (sumFs [0..]) [0..]
where
sumFs = zipWith (+) [ a+b+c+d+e+f+g+h+i+j
| a <- [0..9] , b <- [0..9]
, c <- [0..9] , d <- [0..9]
, e <- [0..9] , f <- [0..9]
, g <- [0..9] , h <- [0..9]
, i <- [0..9] , j <- [0..9] ]
-- More idiomatic list generator is about three times slower
-- sumFs = zipWith (+) $ sum <$> replicateM 10 [0..9]
sieve (f:fs) (n:ns)
| n > f = sieve fs (n:ns)
| n `notElem` take 81 (f:fs) = n : sieve (f:fs) ns
| otherwise = sieve (f:fs) ns
main = do
print $ take 50 selfs
forM_ [1..8] $ \i ->
printf "1e%v\t%v\n" (i :: Int) (selfs !! (10^i-1))
$ ghc -O2 SelfNum.hs && time ./SelfNum [1,3,5,7,9,20,31,42,53,64,75,86,97,108,110,121,132,143,154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323,334,345,356,367,378,389,400,411,413,424,435,446,457,468] 1e1 64 1e2 973 1e3 10188 1e4 102225 1e5 1022675 1e6 10227221 1e7 102272662 1e8 1022727208 275.98 user 3.11 system 4:41.02 elapsed
Java
public class SelfNumbers {
private static final int MC = 103 * 1000 * 10000 + 11 * 9 + 1;
private static final boolean[] SV = new boolean[MC + 1];
private static void sieve() {
int[] dS = new int[10_000];
for (int a = 9, i = 9999; a >= 0; a--) {
for (int b = 9; b >= 0; b--) {
for (int c = 9, s = a + b; c >= 0; c--) {
for (int d = 9, t = s + c; d >= 0; d--) {
dS[i--] = t + d;
}
}
}
}
for (int a = 0, n = 0; a < 103; a++) {
for (int b = 0, d = dS[a]; b < 1000; b++, n += 10000) {
for (int c = 0, s = d + dS[b] + n; c < 10000; c++) {
SV[dS[c] + s++] = true;
}
}
}
}
public static void main(String[] args) {
sieve();
System.out.println("The first 50 self numbers are:");
for (int i = 0, count = 0; count <= 50; i++) {
if (!SV[i]) {
count++;
if (count <= 50) {
System.out.printf("%d ", i);
} else {
System.out.printf("%n%n Index Self number%n");
}
}
}
for (int i = 0, limit = 1, count = 0; i < MC; i++) {
if (!SV[i]) {
if (++count == limit) {
System.out.printf("%,12d %,13d%n", count, i);
limit *= 10;
}
}
}
}
}
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1,000 10,188 10,000 102,225 100,000 1,022,675 1,000,000 10,227,221 10,000,000 102,272,662 100,000,000 1,022,727,208
jq
Adapted from Julia
def sumdigits: tostring | explode | map([.]|implode|tonumber) | add;
def gsum: . + sumdigits;
def isnonself:
def ndigits: tostring|length;
. as $i
| ($i - ($i|ndigits)*9) as $n
| any( range($i-1; [0,$n]|max; -1);
gsum == $i);
# an array
def last81:
[limit(81; range(1; infinite) | select(isnonself))];
# output an unbounded stream
def selfnumbers:
foreach range(1; infinite) as $i ( [0, last81];
.[0] = false
| .[1] as $last81
| if $last81 | bsearch($i) < 0
then .[0] = $i
| if $i > $last81[-1] then .[1] = $last81[1:] + [$i | gsum ] else . end
else .
end;
.[0] | select(.) );
"The first 50 self numbers are:", last81[:50],
"",
(nth(100000000 - 1; selfnumbers)
| if . == 1022727208
then "Yes, \(.) is the 100,000,000th self number."
else "No, \(.i) is actually the 100,000,000th self number."
end)
- Output:
The first 50 self numbers are: [2,4,6,8,10,11,12,13,14,15,16,17,18,19,21,22,23,24,25,26,27,28,29,30,32,33,34,35,36,37,38,39,40,41,43,44,45,46,47,48,49,50,51,52,54,55,56,57,58,59] Yes, 1022727208 is the 100,000,000th self number.
Julia
The code first bootstraps a sliding window of size 81 and then uses this as a sieve. Note that 81 is the window size because the sum of digits of 999,999,999 (the largest digit sum of a counting number less than 1022727208) is 81.
gsum(i) = sum(digits(i)) + i
isnonself(i) = any(x -> gsum(x) == i, i-1:-1:i-max(1, ndigits(i)*9))
const last81 = filter(isnonself, 1:5000)[1:81]
function checkselfnumbers()
i, selfcount = 1, 0
while selfcount <= 100_000_000 && i <= 1022727208
if !(i in last81)
selfcount += 1
if selfcount < 51
print(i, " ")
elseif selfcount == 51
println()
elseif selfcount == 100_000_000
println(i == 1022727208 ?
"Yes, $i is the 100,000,000th self number." :
"No, instead $i is the 100,000,000th self number.")
end
end
popfirst!(last81)
push!(last81, gsum(i))
i += 1
end
end
checkselfnumbers()
- Output:
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Yes, 1022727208 is the 100,000,000th self number.
Faster version
Contains tweaks peculiar to the "10 to the nth" self number. Timings include compilation times.
const MAXCOUNT = 103 * 10000 * 10000 + 11 * 9 + 1
function dosieve!(sieve, digitsum9999)
n = 1
for a in 1:103, b in 1:10000
s = digitsum9999[a] + digitsum9999[b] + n
for c in 1:10000
sieve[digitsum9999[c] + s] = true
s += 1
end
n += 10000
end
end
initdigitsum() = reverse!(vec([sum(k) for k in Iterators.product(9:-1:0, 9:-1:0, 9:-1:0, 9:-1:0)]))
function findselves()
sieve = zeros(Bool, MAXCOUNT+1)
println("Sieve time:")
@time begin
digitsum = initdigitsum()
dosieve!(sieve, digitsum)
end
cnt = 1
for i in 1:MAXCOUNT+1
if !sieve[i]
cnt > 50 && break
print(i, " ")
cnt += 1
end
end
println()
limit, cnt = 1, 0
for i in 0:MAXCOUNT
cnt += 1 - sieve[i + 1]
if cnt == limit
println(lpad(cnt, 10), lpad(i, 12))
limit *= 10
end
end
end
@time findselves()
- Output:
Sieve time: 7.187635 seconds (2 allocations: 78.203 KiB) 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 1000000000 10227272649 16.999383 seconds (42.92 k allocations: 9.595 GiB, 0.01% gc time)
Kotlin
private const val MC = 103 * 1000 * 10000 + 11 * 9 + 1
private val SV = BooleanArray(MC + 1)
private fun sieve() {
val dS = IntArray(10000)
run {
var a = 9
var i = 9999
while (a >= 0) {
for (b in 9 downTo 0) {
var c = 9
val s = a + b
while (c >= 0) {
var d = 9
val t = s + c
while (d >= 0) {
dS[i--] = t + d
d--
}
c--
}
}
a--
}
}
var a = 0
var n = 0
while (a < 103) {
var b = 0
val d = dS[a]
while (b < 1000) {
var c = 0
var s = d + dS[b] + n
while (c < 10000) {
SV[dS[c] + s++] = true
c++
}
b++
n += 10000
}
a++
}
}
fun main() {
sieve()
println("The first 50 self numbers are:")
run {
var i = 0
var count = 0
while (count <= 50) {
if (!SV[i]) {
count++
if (count <= 50) {
print("$i ")
} else {
println()
println()
println(" Index Self number")
}
}
i++
}
}
var i = 0
var limit = 1
var count = 0
while (i < MC) {
if (!SV[i]) {
if (++count == limit) {
println("%,12d %,13d".format(count, i))
limit *= 10
}
}
i++
}
}
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Index Self number 1 1 10 64 100 973 1,000 10,188 10,000 102,225 100,000 1,022,675 1,000,000 10,227,221 10,000,000 102,272,662 100,000,000 1,022,727,208
Mathematica/Wolfram Language
sum[g_] := g + Total@IntegerDigits@g
ming[n_] := n - IntegerLength[n]*9
self[n_] := NoneTrue [Range[ming[n], n - 1], sum[#] == n &]
Module[{t = 1, x = 1},
Reap[
While[t <= 50,
If[self[x], Sow[x]; t++]; x++]
][[2, 1]]]
- Output:
{1,3,5,7,9,20,31,42,53,64,75,86,97,108,110,121,132,143,154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323,334,345,356,367,378,389,400,411,413,424,435,446,457,468}
Nim
In order to use less memory, we have chosen to use indexing at bit level. So, our sieve is a custom object defined by its length in bits and its value which is a sequence of bytes. With bit indexing, the sieve uses eight times less memory than with byte indexing.
Of course, there is a trade off to this strategy: reading values from and writing values to the sieve are significantly slower.
Simple sieve
We use the Go algorithm with bit indexing. As a consequence the sieve uses about 250MB instead of 1 GB. And the program is about five times slower.
Note that we used a sequence of ten nested loops as in the Go solution but we have not memorized the intermediate sums as the C compiler does a good job to detect the loop invariants (remember, Nim produces C code and this code has proved to be quite optimizable by the C compiler). The ten loops looks a lot better this way, too 🙂.
import bitops, strutils, std/monotimes, times
const MaxCount = 2 * 1_000_000_000 + 9 * 9
# Bit string used to represent an array of booleans.
type BitString = object
len: Natural # length in bits.
values: seq[byte] # Sequence containing the bits.
proc newBitString(n: Natural): BitString =
## Return a new bit string of length "n" bits.
result.len = n
result.values.setLen((n + 7) shr 3)
template checkIndex(i, length: Natural) {.used.} =
## Check if index "i" is less than the array length.
if i >= length:
raise newException(IndexDefect, "index $1 not in 0 .. $2".format(i, length))
proc `[]`(bs: BitString; i: Natural): bool =
## Return the value of bit at index "i" as a boolean.
when compileOption("boundchecks"):
checkIndex(i, bs.len)
result = bs.values[i shr 3].testbit(i and 0x07)
proc `[]=`(bs: var BitString; i: Natural; value: bool) =
## Set the bit at index "i" to the given value.
when compileOption("boundchecks"):
checkIndex(i, bs.len)
if value: bs.values[i shr 3].setBit(i and 0x07)
else: bs.values[i shr 3].clearBit(i and 0x07)
proc fill(sieve: var BitString) =
## Fill a sieve.
var n = 0
for a in 0..1:
for b in 0..9:
for c in 0..9:
for d in 0..9:
for e in 0..9:
for f in 0..9:
for g in 0..9:
for h in 0..9:
for i in 0..9:
for j in 0..9:
sieve[a + b + c + d + e + f + g + h + i + j + n] = true
inc n
let t0 = getMonoTime()
var sieve = newBitString(MaxCount + 1)
sieve.fill()
echo "Sieve time: ", getMonoTime() - t0
# Find first 50.
echo "\nFirst 50 self numbers:"
var count = 0
var line = ""
for n in 0..MaxCount:
if not sieve[n]:
inc count
line.addSep(" ")
line.add $n
if count == 50: break
echo line
# Find 1st, 10th, 100th, ..., 100_000_000th.
echo "\n Rank Value"
var limit = 1
count = 0
for n in 0..MaxCount:
if not sieve[n]: inc count
if count == limit:
echo ($count).align(10), ($n).align(12)
limit *= 10
echo "Total time: ", getMonoTime() - t0
- Output:
Sieve time: 2 seconds, 59 milliseconds, 67 microseconds, and 152 nanoseconds First 50 self numbers: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Rank Value 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 Total time: 7 seconds, 903 milliseconds, 752 microseconds, and 944 nanoseconds
Improved sieve
Using bit indexing is very useful here as with byte indexing, the sieve needs 10GB. On a computer with only 8 GB , as this is the case of the laptop I use to run these programs, it fails to execute (I have a very small swap and don’t want to use the swap anyway). With bit indexing, the sieve needs only 1,25GB which is more reasonable.
Of course, the program is slower but not in the same proportion that in the previous program: it is about twice slower than the Pascal version. Note that the execution time varies significantly according to the way statements are written. For instance, writing if not sieve[n]: inc count
has proved to be more efficient than writing inc count, ord(not sieve[n])
or inc count, 1 - ord(sieve[n])
which is surprising as the latter forms avoid a jump. Maybe changing some other statements could give better results, but current time is already satisfying.
import bitops, strutils, std/monotimes, times
const MaxCount = 103 * 10_000 * 10_000 + 11 * 9 + 1
# Bit string used to represent an array of booleans.
type BitString = object
len: Natural
values: seq[byte]
proc newBitString(n: Natural): BitString =
## Return a new bit string of length "n" bits.
result.len = n
result.values.setLen((n + 7) shr 3)
template checkIndex(i, length: Natural) {.used.} =
## Check if index "i" is less than the array length.
if i >= length:
raise newException(IndexDefect, "index $1 not in 0 .. $2".format(i, length))
proc `[]`(bs: BitString; i: Natural): bool =
## Return the value of bit at index "i" as a boolean.
when compileOption("boundchecks"):
checkIndex(i, bs.len)
result = bs.values[i shr 3].testbit(i and 0x07)
proc `[]=`(bs: var BitString; i: Natural; value: bool) =
## Set the bit at index "i" to the given value.
when compileOption("boundchecks"):
checkIndex(i, bs.len)
if value: bs.values[i shr 3].setBit(i and 0x07)
else: bs.values[i shr 3].clearBit(i and 0x07)
proc initDigitSum9999(): array[10000, byte] {.compileTime.} =
## Return the array of the digit sums for numbers 0 to 9999.
var i = 0
for a in 0..9:
for b in 0..9:
for c in 0..9:
for d in 0..9:
result[i] = byte(a + b + c + d)
inc i
const DigitSum9999 = initDigitSum9999()
proc fill(sieve: var BitString) =
## Fill a sieve.
var n = 0
for a in 0..102:
for b in 0..9999:
var s = DigitSum9999[a].int + DigitSum9999[b].int + n
for c in 0..9999:
sieve[DigitSum9999[c].int + s] = true
inc s
inc n, 10_000
let t0 = getMonoTime()
var sieve = newBitString(MaxCount + 1)
sieve.fill()
echo "Sieve time: ", getMonoTime() - t0
# Find first 50.
echo "\nFirst 50 self numbers:"
var count = 0
var line = ""
for n in 0..MaxCount:
if not sieve[n]:
inc count
line.addSep(" ")
line.add $n
if count == 50: break
echo line
# Find 1st, 10th, 100th, ..., 1_000_000_000th.
echo "\n Rank Value"
var limit = 1
count = 0
for n in 0..MaxCount:
if not sieve[n]: inc count
if count == limit:
echo ($count).align(10), ($n).align(12)
limit *= 10
echo "Total time: ", getMonoTime() - t0
- Output:
Sieve time: 13 seconds, 340 milliseconds, 45 microseconds, and 528 nanoseconds First 50 self numbers: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 Rank Value 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 1000000000 10227272649 Total time: 28 seconds, 135 milliseconds, 481 microseconds, and 697 nanoseconds
Pascal
Just "sieving" with only one follower of every number
Extended to 10.23e9
program selfnumb;
{$IFDEF FPC}
{$MODE Delphi}
{$Optimization ON,ALL}
{$IFEND}
{$IFDEF DELPHI} {$APPTYPE CONSOLE} {$IFEND}
uses
sysutils;
const
MAXCOUNT =103*10000*10000+11*9+ 1;
type
tDigitSum9999 = array[0..9999] of Uint8;
tpDigitSum9999 = ^tDigitSum9999;
var
DigitSum9999 : tDigitSum9999;
sieve : array of boolean;
procedure dosieve;
var
pSieve : pBoolean;
pDigitSum :tpDigitSum9999;
n,c,b,a,s : NativeInt;
Begin
pSieve := @sieve[0];
pDigitSum := @DigitSum9999[0];
n := 0;
for a := 0 to 102 do
for b := 0 to 9999 do
Begin
s := pDigitSum^[a]+pDigitSum^[b]+n;
for c := 0 to 9999 do
Begin
pSieve[pDigitSum^[c]+s] := true;
s+=1;
end;
inc(n,10000);
end;
end;
procedure InitDigitSum;
var
i,d,c,b,a : NativeInt;
begin
i := 9999;
for a := 9 downto 0 do
for b := 9 downto 0 do
for c := 9 downto 0 do
for d := 9 downto 0 do
Begin
DigitSum9999[i] := a+b+c+d;
dec(i);
end;
end;
procedure OutPut(cnt,i:NativeUint);
Begin
writeln(cnt:10,i:12);
end;
var
pSieve : pboolean;
T0 : Uint64;
i,cnt,limit,One: NativeUInt;
BEGIN
setlength(sieve,MAXCOUNT);
pSieve := @sieve[0];
T0 := GetTickCount64;
InitDigitSum;
dosieve;
writeln('Sievetime : ',(GetTickCount64-T0 )/1000:8:3,' sec');
//find first 50
cnt := 0;
for i := 0 to MAXCOUNT do
Begin
if NOT(pSieve[i]) then
Begin
inc(cnt);
if cnt <= 50 then
write(i:4)
else
BREAK;
end;
end;
writeln;
One := 1;
limit := One;
cnt := 0;
for i := 0 to MAXCOUNT do
Begin
inc(cnt,One-Ord(pSieve[i]));
if cnt = limit then
Begin
OutPut(cnt,i);
limit := limit*10;
end;
end;
END.
- Output:
time ./selfnumb Sievetime : 6.579 sec 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 1 1 10 64 100 973 1000 10188 10000 102225 100000 1022675 1000000 10227221 10000000 102272662 100000000 1022727208 1000000000 10227272649 real 0m13,252s
Perl
use strict;
use warnings;
use feature 'say';
use List::Util qw(max sum);
my ($i, $pow, $digits, $offset, $lastSelf, @selfs)
= ( 1, 10, 1, 9, 0, );
my $final = 50;
while () {
my $isSelf = 1;
my $sum = my $start = sum split //, max(($i-$offset), 0);
for ( my $j = $start; $j < $i; $j++ ) {
if ($j+$sum == $i) { $isSelf = 0 ; last }
($j+1)%10 != 0 ? $sum++ : ( $sum = sum split '', ($j+1) );
}
if ($isSelf) {
push @selfs, $lastSelf = $i;
last if @selfs == $final;
}
next unless ++$i % $pow == 0;
$pow *= 10;
$offset = 9 * $digits++
}
say "The first 50 self numbers are:\n" . join ' ', @selfs;
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468
Phix
Certainly puts my previous rubbish attempts (archived here) to shame.
The precise nature of the difference-pattern eludes me, I will admit.
-- -- Base-10 self numbers by index (single or range). -- Follows an observed sequence pattern whereby, after the initial single-digit odd numbers, self numbers are -- grouped in runs whose members occur at numeric intervals of 11. Runs after the first one come in blocks of -- ten: eight runs of ten numbers followed by two shorter runs. The numeric interval between runs is usually 2, -- but that between shorter runs, and their length, depend on the highest-order digit change occurring in them. -- This connection with significant digit change means every ten blocks form a higher-order block, every ten -- of these a higher-order-still block, and so on. -- -- The code below appears to be good up to the last self number before 10^12 — ie. 999,999,999,997, which is -- returned as the 97,777,777,792nd such number. After this, instead of zero-length shorter runs, the actual -- pattern apparently starts again with a single run of 10, like the one at the beginning. -- integer startIndex, endIndex, counter atom currentSelf sequence output function doneAfterAdding(integer interval, n) -- Advance to the next self number in the sequence, append it to the output if required, indicate if finished. for i=1 to n do currentSelf += interval counter += 1 if counter >= startIndex then output &= currentSelf end if if counter = endIndex then return true end if end for return false end function function selfNumbers(sequence indexRange) startIndex = indexRange[1] endIndex = indexRange[$] counter = 0 currentSelf = -1 output = {} -- Main process. Start with the single-digit odd numbers and first run. if doneAfterAdding(2,5) then return output end if if doneAfterAdding(11,9) then return output end if -- If necessary, fast forward to last self number before the lowest-order block containing first number rqd. if counter<startIndex then -- The highest-order blocks whose ends this handles correctly contain 9,777,777,778 self numbers. -- The difference between equivalently positioned numbers in these blocks is 100,000,000,001. -- The figures for successively lower-order blocks have successively fewer 7s and 0s! atom indexDiff = 9777777778, numericDiff = 100000000001 while indexDiff>=98 and counter!=startIndex do if counter+indexDiff < startIndex then counter += indexDiff currentSelf += numericDiff else indexDiff = (indexDiff+2)/10 -- (..78->80->8) numericDiff = (numericDiff+9)/10 -- (..01->10->1) end if end while end if -- Sequencing loop, per lowest-order block. while true do -- Eight ten-number runs, each at a numeric interval of 2 from the end of the previous one. for i=1 to 8 do if doneAfterAdding(2,1) then return output end if if doneAfterAdding(11,9) then return output end if end for -- Two shorter runs, the second at an interval inversely related to their length. integer shorterRunLength = 8, temp = floor(currentSelf/1000) -- Work out a shorter run length based on the most significant digit change about to happen. while remainder(temp,10)=9 do shorterRunLength -= 1 temp = floor(temp/10) end while integer interval = 2 for i=1 to 2 do if doneAfterAdding(interval,1) then return output end if if doneAfterAdding(11,shorterRunLength) then return output end if interval += (9-shorterRunLength)*13 end for end while end function atom t0 = time() printf(1,"The first 50 self numbers are:\n") pp(selfNumbers({1, 50}),{pp_IntFmt,"%3d",pp_IntCh,false}) for p=8 to 9 do integer n = power(10,p) printf(1,"The %,dth safe number is %,d\n",{n,selfNumbers({n})[1]}) end for printf(1,"completed in %s\n",elapsed(time()-t0))
- Output:
The first 50 self numbers are: { 1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97,108,110,121,132,143, 154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323, 334,345,356,367,378,389,400,411,413,424,435,446,457,468} The 100,000,000th safe number is 1,022,727,208 The 1,000,000,000th safe number is 10,227,272,649 completed in 0.1s
Python
class DigitSumer :
def __init__(self):
sumdigit = lambda n : sum( map( int,str( n )))
self.t = [sumdigit( i ) for i in xrange( 10000 )]
def __call__ ( self,n ):
r = 0
while n >= 10000 :
n,q = divmod( n,10000 )
r += self.t[q]
return r + self.t[n]
def self_numbers ():
d = DigitSumer()
s = set([])
i = 1
while 1 :
n = i + d( i )
if i in s :
s.discard( i )
else:
yield i
s.add( n )
i += 1
import time
p = 100
t = time.time()
for i,s in enumerate( self_numbers(),1 ):
if i <= 50 :
print s,
if i == 50 : print
if i == p :
print '%7.1f sec %9dth = %d'%( time.time()-t,i,s )
p *= 10
- Output:
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 0.0 sec 100th = 973 0.0 sec 1000th = 10188 0.1 sec 10000th = 102225 1.0 sec 100000th = 1022675 11.4 sec 1000000th = 10227221 143.4 sec 10000000th = 102272662 1812.0 sec 100000000th = 1022727208
Raku
Translated the low memory version of the Go entry but showed only the first 50 self numbers. The machine for running this task (a Xeon E3110+8GB memory) is showing its age as, 1) took over 6 minutes to complete the Go entry, 2) not even able to run the other two Go alternative entries and 3) needed over 47 minutes to reach 1e6 iterations here. Anyway I will try this on an i5 box later to see how it goes.
# 20201127 Raku programming solution
my ( $st, $count, $i, $pow, $digits, $offset, $lastSelf, $done, @selfs) =
now, 0, 1, 10, 1, 9, 0, False;
# while ( $count < 1e8 ) {
until $done {
my $isSelf = True;
my $sum = (my $start = max ($i-$offset), 0).comb.sum;
loop ( my $j = $start; $j < $i; $j+=1 ) {
if $j+$sum == $i { $isSelf = False and last }
($j+1)%10 != 0 ?? ( $sum+=1 ) !! ( $sum = ($j+1).comb.sum ) ;
}
if $isSelf {
$count+=1;
$lastSelf = $i;
if $count ≤ 50 {
@selfs.append: $i;
if $count == 50 {
say "The first 50 self numbers are:\n", @selfs;
$done = True;
}
}
}
$i+=1;
if $i % $pow == 0 {
$pow *= 10;
$digits+=1 ;
$offset = $digits * 9
}
}
# say "The 100 millionth self number is ", $lastSelf;
# say "Took ", now - $st, " seconds."
- Output:
The first 50 self numbers are: [1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468]
REXX
first 50 self numbers
/*REXX program displays N self numbers (aka Colombian or Devlali numbers). OEIS A3052.*/
parse arg n . /*obtain optional argument from the CL.*/
if n=='' | n=="," then n= 50 /*Not specified? Then use the default.*/
tell = n>0; n= abs(n) /*TELL: show the self numbers if N>0 */
@.= . /*initialize the array of self numbers.*/
do j=1 for n*10 /*scan through ten times the #s wanted.*/
$= j /*1st part of sum is the number itself.*/
do k=1 for length(j) /*sum the decimal digits in the number.*/
$= $ + substr(j, k, 1) /*add a particular digit to the sum. */
end /*k*/
@.$= /*mark J as not being a self number. */
end /*j*/ /* ─── */
list= 1 /*initialize the list to the 1st number*/
#= 1 /*the count of self numbers (so far). */
do i=3 until #==n; if @.i=='' then iterate /*Not a self number? Then skip it. */
#= # + 1; list= list i /*bump counter of self #'s; add to list*/
end /*i*/
/*stick a fork in it, we're all done. */
say n " self numbers were found." /*display the title for the output list*/
if tell then say list /*display list of self numbers ──►term.*/
- output when using the default input:
50 self numbers were found. 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468
ten millionth self number
/*REXX pgm displays the Nth self number, aka Colombian or Devlali numbers. OEIS A3052.*/
numeric digits 20 /*ensure enough decimal digits for #. */
parse arg high . /*obtain optional argument from the CL.*/
if high=='' | high=="," then high= 10000000 /*Not specified? Then use 10M default.*/
i= 1; pow= 10; digs= 1; offset= 9; $= 0 /*$: the last self number found. */
#= 0 /*count of self numbers (so far). */
do while #<high; isSelf= 1 /*assume a self number (so far). */
start= max(i-offset, 0) /*find start #; must be non-negative. */
sum= sumDigs(start) /*obtain the sum of the decimal digits.*/
do j=start to i-1
if j+sum==i then do; isSelf= 0 /*found a non self number. */
iterate /*keep looking for more self numbers. */
end
if (j+1)//10==0 then sum= sumDigs(j+1) /*obtain the sum of the decimal digits.*/
else sum= sum + 1 /*bump " " " " " " */
end /*j*/
if isSelf then do; #= # + 1 /*bump the count of self numbers. */
$= i /*save the last self number found. */
end
i= i + 1 /*bump the self number by unity. */
if i//pow==0 then do; digs= digs + 1 /* " " number of decimal digits. */
pow= pow * 10 /* " " power " a factor of ten. */
offset= digs * 9 /* " " offset " " " " nine. */
end
end /*while*/
say
say 'the ' commas(high)th(high) " self number is: " commas($)
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
sumDigs: parse arg s 2 x; do k=1 for length(x) /*get 1st dig, & also get the rest.*/
s= s + substr(x, k, 1) /*add a particular digit to the sum.*/
end /*k*/; return s
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do c=length(_)-3 to 1 by -3; _=insert(',', _, c); end; return _
th: parse arg th; return word('th st nd rd', 1 +(th//10)*(th//100%10\==1)*(th//10<4))
- output when using the default input:
the 100,000,000th self number is: 1,022,727,208
Ring
load "stdlib.ring"
see "working..." + nl
see "The first 50 self numbers are:" + nl
n = 0
num = 0
limit = 51
limit2 = 10000000
while true
n = n + 1
for m = 1 to n
flag = 1
sum = 0
strnum = string(m)
for p = 1 to len(strnum)
sum = sum + number(strnum[p])
next
sum2 = m + sum
if sum2 = n
flag = 0
exit
else
flag = 1
ok
next
if flag = 1
num = num + 1
if num < limit
see "" + num + ". " + n + nl
ok
if num = limit2
see "The " + limit2 + "th self number is: " + n + nl
ok
if num > limit2
exit
ok
ok
end
see "done..." + nl
Output:
working... The first 50 self numbers are: 1. 1 2. 3 3. 5 4. 7 5. 9 6. 20 7. 31 8. 42 9. 53 10. 64 11. 75 12. 86 13. 97 14. 108 15. 110 16. 121 17. 132 18. 143 19. 154 20. 165 21. 176 22. 187 23. 198 24. 209 25. 211 26. 222 27. 233 28. 244 29. 255 30. 266 31. 277 32. 288 33. 299 34. 310 35. 312 36. 323 37. 334 38. 345 39. 356 40. 367 41. 378 42. 389 43. 400 44. 411 45. 413 46. 424 47. 435 48. 446 49. 457 50. 468 The 10000000th self number is: 1022727208 done...
RPL
Brute force
Using a sieve:
« 0 WHILE OVER REPEAT SWAP 10 MOD LASTARG / IP ROT ROT + END + » 'DIGSUM' STO « 500 0 → max n « { } DUP max + 0 CON 1 CF DO 'n' INCR DUP DIGSUM + IFERR 1 PUT THEN DROP2 1 SF END UNTIL 1 FS? END 1 WHILE 3 PICK SIZE 50 < REPEAT IF DUP2 GET NOT THEN ROT OVER + ROT ROT END 1 + END DROP2 » » 'TASK' STO
- Output:
1: { 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 }
Runs in 2 minutes 8 seconds on a HP-48SX.
Maximilian F. Hasler's algorithm
Translation of the PARI/GP formula on the OEIS page:
« → n « { } 1 SF 1 n 2 / IP n XPON 1 + 9 * MIN FOR j IF n j - DIGSUM j == THEN 1 CF n 'j' STO END NEXT 1 FS? » 'SELF?' STO « { } WHILE OVER SIZE 50 < REPEAT IF DUP SELF? THEN SWAP OVER + SWAP END 1 + END DROP » 'TASK' STO
Same result as above. No need for sieve, but much slower: 10 minutes 52 seconds on a HP-48SX.
Sidef
Algorithm by David A. Corneth (OEIS: A003052).
func is_self_number(n) {
if (n < 30) {
return (((n < 10) && (n.is_odd)) || (n == 20))
}
var qd = (1 + n.ilog10)
var r = (1 + (n-1)%9)
var h = (r + 9*(r%2))/2
var ld = 10
while (h + 9*qd >= n%ld) {
ld *= 10
}
var vs = idiv(n, ld).sumdigits
n %= ld
0..qd -> none { |i|
vs + sumdigits(n - h - 9*i) == (h + 9*i)
}
}
say is_self_number.first(50).join(' ')
- Output:
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468
Simpler algorithm (by M. F. Hasler):
func is_self_number(n) {
1..min(n>>1, 9*n.len) -> none {|i| sumdigits(n-i) == i } && (n > 0)
}
Standard ML
open List;
val rec selfNumberNr = fn NR =>
let
val rec sumdgt = fn 0 => 0 | n => Int.rem (n, 10) + sumdgt (Int.quot(n ,10));
val rec isSelf = fn ([],l1,l2) => []
| (x::tt,l1,l2) => if exists (fn i=>i=x) l1 orelse exists (fn i=>i=x) l2
then ( isSelf (tt,l1,l2)) else x::isSelf (tt,l1,l2) ;
val rec partcount = fn (n, listIn , count , selfs) =>
if count >= NR then nth (selfs, length selfs + NR - count -1)
else
let
val range = tabulate (81 , fn i => 81*n +i+1) ;
val listOut = map (fn i => i + sumdgt i ) range ;
val selfs = isSelf (range,listIn,listOut)
in
partcount ( n+1 , listOut , count+length (selfs) , selfs )
end;
in
partcount (0,[],0,[])
end;
app ((fn s => print (s ^ " ")) o Int.toString o selfNumberNr) (tabulate (50,fn i=>i+1)) ;
selfNumberNr 100000000 ;
output
1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 1022727208
Wren
Just the sieve based version as the low memory version would take too long to run in Wren.
Note that you need a lot of memory to run this as Bools in Wren require 8 bytes of storage compared to 1 byte in Go.
Unsurprisingly, very slow compared to the Go version as Wren is interpreted and uses floating point arithmetic for all numerical work.
var sieve = Fn.new {
var sv = List.filled(2*1e9+9*9+1, false)
var n = 0
var s = [0] * 8
for (a in 0..1) {
for (b in 0..9) {
s[0] = a + b
for (c in 0..9) {
s[1] = s[0] + c
for (d in 0..9) {
s[2] = s[1] + d
for (e in 0..9) {
s[3] = s[2] + e
for (f in 0..9) {
s[4] = s[3] + f
for (g in 0..9) {
s[5] = s[4] + g
for (h in 0..9) {
s[6] = s[5] + h
for (i in 0..9) {
s[7] = s[6] + i
for (j in 0..9) {
sv[s[7] + j + n] = true
n = n + 1
}
}
}
}
}
}
}
}
}
}
return sv
}
var st = System.clock
var sv = sieve.call()
var count = 0
System.print("The first 50 self numbers are:")
for (i in 0...sv.count) {
if (!sv[i]) {
count = count + 1
if (count <= 50) System.write("%(i) ")
if (count == 1e8) {
System.print("\n\nThe 100 millionth self number is %(i)")
break
}
}
}
System.print("Took %(System.clock-st) seconds.")
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 The 100 millionth self number is 1022727208 Took 222.789713 seconds.
XPL0
As run on Raspberry Pi4.
def LenSV = 2_000_000_000 + 9*9 + 1;
func Sieve;
char SV;
int N, S(8), A, B, C, D, E, F, G, H, I, J;
[SV:= MAlloc(LenSV);
N:= 0;
for A:= 0 to 1 do
[for B:= 0 to 9 do
[S(0):= A + B;
for C:= 0 to 9 do
[S(1):= S(0) + C;
for D:= 0 to 9 do
[S(2):= S(1) + D;
for E:= 0 to 9 do
[S(3):= S(2) + E;
for F:= 0 to 9 do
[S(4):= S(3) + F;
for G:= 0 to 9 do
[S(5):= S(4) + G;
for H:= 0 to 9 do
[S(6):= S(5) + H;
for I:= 0 to 9 do
[S(7):= S(6) + I;
for J:= 0 to 9 do
[SV(S(7)+J+N):= true;
N:= N+1;
]
]
]
]
]
]
]
]
]
];
return SV;
];
char SV;
int ST, Count, I;
[ST:= GetTime;
SV:= Sieve;
Count:= 0;
Text(0, "The first 50 self numbers are:^m^j");
for I:= 0 to LenSV-1 do
[if SV(I) = false then
[Count:= Count+1;
if Count <= 50 then
[IntOut(0, I); ChOut(0, ^ )];
if Count = 100_000_000 then
[Text(0, "^m^j^m^jThe 100 millionth self number is ");
IntOut(0, I);
CrLf(0);
I:= LenSV;
];
];
];
Text(0, "Took "); RlOut(0, float(GetTime-ST) / 1e6); CrLf(0);
]
- Output:
The first 50 self numbers are: 1 3 5 7 9 20 31 42 53 64 75 86 97 108 110 121 132 143 154 165 176 187 198 209 211 222 233 244 255 266 277 288 299 310 312 323 334 345 356 367 378 389 400 411 413 424 435 446 457 468 The 100 millionth self number is 1022727208 Took 29.46575