Self numbers
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.
Wikipedia Self numbers showing some tricks especially for the number above.
F#
<lang fsharp> // 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) </lang>
- 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
Go
Low memory
Simple-minded, uses very little memory (no sieve) but slow - over 2.5 minutes. <lang go>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))
}</lang>
- 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%. <lang go>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))
}</lang>
- 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. <lang go>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))
}</lang>
- 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
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. <lang julia>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()
</lang>
- 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. <lang julia>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()
</lang>
- 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)
Pascal
Just "sieving" with only one follower of every number
Extended to 10.23e9 <lang pascal>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.</lang>
- 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
Phix
Replacing the problematic sv[a+b+... line with a bit of dirty inline assembly improved performance by 90%
(Of course you lose bounds checking, type checking, negative subscripts, fraction handling, and all that jazz.)
We use a string of Y/N for the sieve to force one byte per element ('\0' and 1 would be equally valid).
<lang Phix>if machine_bits()=32 then crash("requires 64 bit") end if
function sieve()
string sv = repeat('N',2*1e9+9*9+1) 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>
- Output:
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 The 100,000,000th self number is 1,022,727,208 (11.0s)
generator dictionary
While this is dog-slow (see shocking estimate below), it is interesting to note that even by the time it generates the
10,000,000th number, it is only having to juggle a mere 27 generators.
Just a shame that we had to push over 10,000,000 generators onto that stack, and tried to push quite a few more.
Memory use is pretty low, around ~4MB.
[unlike the above, this is perfectly happy on both 32 and 64 bit]
Aside: the getd_index() check is often worth trying with phix dictionaries: if there is a high probability that
the key already exists, it will yield a win, but with a low probability it will just be unhelpful overhead.
<lang Phix>integer gd = new_dict(), g = 1, gdhead = 2, n = 0
function ng(integer n)
integer r = n while r do n += remainder(r,10) r = floor(r/10) end while return n
end function
function self(integer /*i*/) -- note: assumes sequentially invoked (arg i unused) integer nxt
n += 1 while n=gdhead do while g<=gdhead do nxt = ng(g) if getd_index(nxt, gd)=NULL then -- (~25% gain) setd(nxt,0,gd) end if g += 1 end while integer wasgdhead = gdhead while true do gdhead = pop_dict(gd)[1] if gdhead!=wasgdhead then exit end if
-- ?{"ding",wasgdhead} -- 2, once only...
end while nxt = ng(gdhead)
-- if getd_index(nxt, gd)=NULL then -- (~1% loss)
setd(nxt,0,gd)
-- end if
n += (n!=gdhead) end while return n
end function
atom t0 = time() printf(1,"The first 50 self numbers are:\n") pp(apply(tagset(50),self),{pp_IntFmt,"%3d",pp_IntCh,false})
constant limit = 10000000 integer chk = 100 printf(1,"\n i n size time\n") for i=51 to limit do
n = self(i) if i=chk then printf(1,"%,11d %,11d %6d %s\n",{i,n,dict_size(gd),elapsed(time()-t0)}) chk *= 10 end if
end for printf(1,"\nEstimated time for %,d :%s\n",{1e8,elapsed((time()-t0)*1e8/limit)})</lang>
- 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} i n size time 100 973 18 0.1s 1,000 10,188 13 0.2s 10,000 102,225 10 1.0s 100,000 1,022,675 20 9.3s 1,000,000 10,227,221 17 1 minute and 37s 10,000,000 102,272,662 27 16 minutes and 40s Estimated time for 100,000,000 :2 hours, 46 minutes and 37s
For the record, I would not be at all surprised should a translation of this beat 20 minutes (for 1e8)
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. <lang ecmascript> 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.")</lang>
- 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.