Humble numbers: Difference between revisions
→{{header|Nim}}: add much faster versions...
(add RPL) |
GordonBGood (talk | contribs) (→{{header|Nim}}: add much faster versions...) |
||
Line 4,660:
17 10815
18 12759</pre>
===Faster Implementation Iterating Over Logarithms===
While the above submission is adequate for the tasts as defined, it is incredibly slow and doesn't scale linearily with increasing range as Hamming/Humble/N-smooth sequences should. The following code adapts [[Hamming_numbers#Much_faster_iterating_version_using_logarithmic_calculations|one of the Nim Hamming submissions (the fastest sequence using a queue)]], adding the extra prime factor of seven and simplifying the code, as well only using the logarithmic approximation which is adequate to the task of only showing the first 50 humble numbers and the digit counts up to a very high number of digits, as follows:
<syntaxhighlight lang="nim">import std/math, std/strutils
from std/times import inMilliseconds
from std/monotimes import getMonoTime, `-`
let lgshft = 50; let lg10 = 1'u64 shl lgshft; let fac = 2'f64.pow lgshft.float64
let lg7 = (7'f64.log10 * fac).round.uint64
let lg5 = (5'f64.log10 * fac).round.uint64
let lg3 = (3'f64.log10 * fac).round.uint64; let lg2 = lg10 - lg5
proc lr2UInt64(lr: uint64): uint64 = pow(10, (lr.float64 / fac)).round.uint64
iterator humblesLogQ(): uint64 = # {.closure.} =
var hd2 = lg2; var hd3 = lg3; var hd5 = lg5; var hd7 = lg7
var s2msk, s2hdi, s2tli, s3msk, s3hdi, s3tli, s5msk, s5hdi, s5tli = 0
var s2 = newSeq[uint64] 0
var s3 = newSeq[uint64] 0
var s5 = newSeq[uint64] 0
yield 0'u64
while true:
yield hd2
if s2tli == s2hdi:
let osz = if s2msk == 0: 512 else: s2msk + 1
s2.setLen (osz + osz); s2msk = s2.len - 1
if osz > 512:
if s2hdi == 0: s2tli = osz
else: # put extra space between tail and head...
copyMem(addr(s2[s2hdi + osz]), addr(s2[s2hdi]),
sizeof(uint64) * (osz - s2hdi)); s2hdi += osz
s2[s2tli] = hd2 + lg2; s2tli = (s2tli + 1) and s2msk
hd2 = s2[s2hdi]
if hd2 < hd3: s2hdi = (s2hdi + 1) and s2msk
else:
hd2 = hd3
if s3tli == s3hdi:
let osz = if s3msk == 0: 512 else: s3msk + 1
s3.setLen (osz + osz); s3msk = s3.len - 1
if osz > 512:
if s3hdi == 0: s3tli = osz
else: # put extra space between tail and head...
copyMem(addr(s3[s3hdi + osz]), addr(s3[s3hdi]),
sizeof(uint64) * (osz - s3hdi)); s3hdi += osz
s3[s3tli] = hd3 + lg3; s3tli = (s3tli + 1) and s3msk
hd3 = s3[s3hdi]
if hd3 < hd5: s3hdi = (s3hdi + 1) and s3msk
else:
hd3 = hd5
if s5tli == s5hdi:
let osz = if s5msk == 0: 512 else: s5msk + 1
s5.setLen (osz + osz); s5msk = s5.len - 1
if osz > 512:
if s5hdi == 0: s5tli = osz
else: # put extra space between tail and head...
copyMem(addr(s5[s5hdi + osz]), addr(s5[s5hdi]),
sizeof(uint64) * (osz - s5hdi)); s5hdi += osz
s5[s5tli] = hd5 + lg5; s5tli = (s5tli + 1) and s5msk
hd5 = s5[s5hdi]
if hd5 < hd7: s5hdi = (s5hdi + 1) and s5msk
else: hd5 = hd7; hd7 += lg7
proc commaString(s: string): string =
let sz = s.len; let sqlen = (sz + 2) div 3
var sq = newSeq[string](sqlen); var ndx = sqlen - 1
for i in countdown(sz - 1, 0, 3): sq[ndx] = s[(max(i-2, 0) .. i)]; ndx -= 1
sq.join ","
# testing it...
let numdigits = 255.uint64
echo "First 50 Humble numbers:"
var cnt = 0
for h in humblesLogQ():
stdout.write $h.lr2UInt64, " "; cnt += 1
if cnt >= 50: break
echo "\r\nCount of humble numbers for each digit length 1-", $numdigits, ":"
echo "Digits Count Accum"
let lmt = lg10 * numdigits
let strt = getMonoTime()
var cdigits = 0'u64; cnt = 0; var acnt = 0.int
for h in humblesLogQ():
if (h shr lgshft) <= cdigits: cnt += 1
else:
cdigits += 1; acnt += cnt
echo ($cdigits).align(4) & ($cnt).commaString.align(14) &
($acnt).commaString.align(19)
cnt = 1
if cdigits >= numdigits: break
let elpsd = (getMonoTime() - strt).inMilliseconds
echo "Counting took ", elpsd, " milliseconds."</syntaxhighlight>
{{out}}
<pre>First 50 Humble numbers:
1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 28 30 32 35 36 40 42 45 48 49 50 54 56 60 63 64 70 72 75 80 81 84 90 96 98 100 105 108 112 120
Count of humble numbers for each digit length 1-255:
Digits Count Accum
1 9 9
2 36 45
3 95 140
4 197 337
5 356 693
6 579 1,272
7 882 2,154
8 1,272 3,426
9 1,767 5,193
10 2,381 7,574
11 3,113 10,687
12 3,984 14,671
13 5,002 19,673
14 6,187 25,860
15 7,545 33,405
16 9,081 42,486
17 10,815 53,301
18 12,759 66,060
19 14,927 80,987
20 17,323 98,310
21 19,960 118,270
22 22,853 141,123
23 26,015 167,138
24 29,458 196,596
25 33,188 229,784
.
.
. results as for the C++ or Pascal versions...
.
.
.
250 30,938,881 1,954,289,627
251 31,310,645 1,985,600,272
252 31,685,379 2,017,285,651
253 32,063,093 2,049,348,744
254 32,443,792 2,081,792,536
255 32,827,496 2,114,620,032
Counting took 6096 milliseconds.</pre>
The above code as run on an Intel i5-6500 at 3.6 GHz (boost single-threaded) is faster than the Pascal version doing about the same as to algorithm, and would be able to show all the digit counts up to 877 digits in about 15 minutes (about six times faster than the Pascal version) except that this machine doesn't have enough memory (about eight Gigabytes required plus as required by the operating system and run-time, and this machine only has eight Gigabytes total), showing that this code is many times faster than the Pascal version while simpler to read.
===Direct Calculation Without Sorting the Sequence===
As noticed in the C++ Direct Generation - Variant submission, there is no need for the complexity and execution time cost of processing the humble numbers in order in order to show the digit counts; the following Nim code implements that algorithm, with the first fifty humble numbers able to be produced by any simple sequence generator, here using the generator from the submission above this:
<syntaxhighlight lang="nim">import std/math, std/strutils
from std/times import inMilliseconds
from std/monotimes import getMonoTime, `-`
let lgshft = 50; let lg10 = 1'u64 shl lgshft; let fac = 2'f64.pow lgshft.float64
let lg7 = (7'f64.log10 * fac).round.uint64
let lg5 = (5'f64.log10 * fac).round.uint64
let lg3 = (3'f64.log10 * fac).round.uint64; let lg2 = lg10 - lg5
proc lr2UInt64(lr: uint64): uint64 = pow(10, (lr.float64 / fac)).round.uint64
iterator humblesLogQ(): uint64 = # {.closure.} =
var hd2 = lg2; var hd3 = lg3; var hd5 = lg5; var hd7 = lg7
var s2msk, s2hdi, s2tli, s3msk, s3hdi, s3tli, s5msk, s5hdi, s5tli = 0
var s2 = newSeq[uint64] 0
var s3 = newSeq[uint64] 0
var s5 = newSeq[uint64] 0
yield 0'u64
while true:
yield hd2
if s2tli == s2hdi:
let osz = if s2msk == 0: 512 else: s2msk + 1
s2.setLen (osz + osz); s2msk = s2.len - 1
if osz > 512:
if s2hdi == 0: s2tli = osz
else: # put extra space between tail and head...
copyMem(addr(s2[s2hdi + osz]), addr(s2[s2hdi]),
sizeof(uint64) * (osz - s2hdi)); s2hdi += osz
s2[s2tli] = hd2 + lg2; s2tli = (s2tli + 1) and s2msk
hd2 = s2[s2hdi]
if hd2 < hd3: s2hdi = (s2hdi + 1) and s2msk
else:
hd2 = hd3
if s3tli == s3hdi:
let osz = if s3msk == 0: 512 else: s3msk + 1
s3.setLen (osz + osz); s3msk = s3.len - 1
if osz > 512:
if s3hdi == 0: s3tli = osz
else: # put extra space between tail and head...
copyMem(addr(s3[s3hdi + osz]), addr(s3[s3hdi]),
sizeof(uint64) * (osz - s3hdi)); s3hdi += osz
s3[s3tli] = hd3 + lg3; s3tli = (s3tli + 1) and s3msk
hd3 = s3[s3hdi]
if hd3 < hd5: s3hdi = (s3hdi + 1) and s3msk
else:
hd3 = hd5
if s5tli == s5hdi:
let osz = if s5msk == 0: 512 else: s5msk + 1
s5.setLen (osz + osz); s5msk = s5.len - 1
if osz > 512:
if s5hdi == 0: s5tli = osz
else: # put extra space between tail and head...
copyMem(addr(s5[s5hdi + osz]), addr(s5[s5hdi]),
sizeof(uint64) * (osz - s5hdi)); s5hdi += osz
s5[s5tli] = hd5 + lg5; s5tli = (s5tli + 1) and s5msk
hd5 = s5[s5hdi]
if hd5 < hd7: s5hdi = (s5hdi + 1) and s5msk
else: hd5 = hd7; hd7 += lg7
proc commaString(s: string): string =
let sz = s.len; let sqlen = (sz + 2) div 3
var sq = newSeq[string](sqlen); var ndx = sqlen - 1
for i in countdown(sz - 1, 0, 3): sq[ndx] = s[(max(i-2, 0) .. i)]; ndx -= 1
sq.join ","
# testing it...
let numdigits = 877.uint64
echo "First 50 Humble numbers:"
var cnt = 0
for h in humblesLogQ():
stdout.write $h.lr2UInt64, " "; cnt += 1
if cnt >= 50: break
echo "\r\nCount of humble numbers for each digit length 1-", $numdigits, ":"
echo "Digits Count Accum"
let lmt = lg10 * numdigits
let strt = getMonoTime()
var bins = newSeq[int](numdigits)
for w in countup(0'u64, lmt, lg7):
for x in countup(w, lmt, lg5):
for y in countup(x, lmt, lg3):
for z in countup(y, lmt, lg2):
bins[z shr lgshft] += 1
var a = 0
for i, c in bins:
a += c
echo ($(i + 1)).align(4) & ($c).commaString.align(14) &
($a).commaString.align(19)
let elpsd = (getMonoTime() - strt).inMilliseconds
echo "Counting took ", elpsd, " milliseconds."</syntaxhighlight>
{{out}}
<pre>First 50 Humble numbers:
1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 28 30 32 35 36 40 42 45 48 49 50 54 56 60 63 64 70 72 75 80 81 84 90 96 98 100 105 108 112 120
Count of humble numbers for each digit length 1-877:
Digits Count Accum
1 9 9
2 36 45
3 95 140
4 197 337
5 356 693
6 579 1,272
7 882 2,154
8 1,272 3,426
9 1,767 5,193
10 2,381 7,574
11 3,113 10,687
12 3,984 14,671
13 5,002 19,673
14 6,187 25,860
15 7,545 33,405
16 9,081 42,486
17 10,815 53,301
18 12,759 66,060
19 14,927 80,987
20 17,323 98,310
21 19,960 118,270
22 22,853 141,123
23 26,015 167,138
24 29,458 196,596
25 33,188 229,784
.
.
. results as for the C++ or Pascal versions...
.
.
.
860 1,252,394,180 270,098,254,942
861 1,256,764,708 271,355,019,650
862 1,261,145,413 272,616,165,063
863 1,265,536,277 273,881,701,340
864 1,269,937,307 275,151,638,647
865 1,274,348,541 276,425,987,188
866 1,278,769,968 277,704,757,156
867 1,283,201,615 278,987,958,771
868 1,287,643,503 280,275,602,274
869 1,292,095,618 281,567,697,892
870 1,296,557,975 282,864,255,867
871 1,301,030,613 284,165,286,480
872 1,305,513,506 285,470,799,986
873 1,310,006,698 286,780,806,684
874 1,314,510,190 288,095,316,874
875 1,319,023,979 289,414,340,853
876 1,323,548,095 290,737,888,948
877 1,328,082,553 292,065,971,501
Counting took 171076 milliseconds.</pre>
The above time as run on an Intel i5-6500 at 3.6 GHz with single-threaded boost is about the same speed as any of the fast language implementations of the same algorithm, including Haskell.
=={{header|Pascal}}==
| |||