Humble numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Humble numbers are positive integers which have no prime factors > 7.
Humble numbers are also called 7-smooth numbers, and sometimes called highly composite,
although this conflicts with another meaning of highly composite numbers.
Another way to express the above is:
humble = 2i × 3j × 5k × 7m
where i, j, k, m ≥ 0
- Task
-
- show the first 50 humble numbers (in a horizontal list)
- show the number of humble numbers that have x decimal digits for all x's up to n (inclusive).
- show (as many as feasible or reasonable for above) on separate lines
- show all output here on this page
- Related tasks
- References
11l
F is_humble(i)
I i <= 1
R 1B
I i % 2 == 0 {R is_humble(i I/ 2)}
I i % 3 == 0 {R is_humble(i I/ 3)}
I i % 5 == 0 {R is_humble(i I/ 5)}
I i % 7 == 0 {R is_humble(i I/ 7)}
R 0B
DefaultDict[Int, Int] humble
V limit = 7F'FF
V count = 0
V num = 1
L count < limit
I is_humble(num)
humble[String(num).len]++
I count < 50
print(num, end' ‘ ’)
count++
num++
print()
print()
print(‘Of the first ’count‘ humble numbers:’)
L(num) 1 .< humble.len - 1
I num !C humble
L.break
print(‘#5 have #. digits’.format(humble[num], num))- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Action!
...which is based on the Algol 68 sample.
As with the PL/M sample, we only consider up to 4 digit Humble Numbers, as Action! integers are limited to 16 bit signed or unsiged.
;;; Find some humble numbers - numbers with no prime factors above 7
PROC humbleStat( CARD s, d ) ;;; displays a statistic about humble numbers
Print( "There are " )
IF s < 10 THEN Put(' ) FI
IF s < 100 THEN Put(' ) FI
PrintC( s )Print( " humble numbers with " )PrintC( d )Print( " digit" )
IF d > 1 THEN Put('s) FI
PutE()
RETURN
PROC Main() ;;; find and print humble numbers
DEFINE MAX_HUMBLE = "400"
CARD ARRAY H( MAX_HUMBLE )
CARD h1, h2, h3, h4, h5, h6, hPos, m
CARD p2, p3, p5, p7
CARD last2, last3, last5, last7
; 1 is the first humble number
H( 0 ) = 1
h1 = 0 h2 = 0 h3 = 0 h4 = 0 h5 = 0 h6 = 0
last2 = 0 last3 = 0 last5 = 0 last7 = 0
p2 = 2 p3 = 3 p5 = 5 p7 = 7
FOR hPos = 1 TO MAX_HUMBLE - 1 DO
; the next humble number is the lowest of the
; next multiples of 2, 3, 5, 7
IF p2 < p3 THEN m = p2 ELSE m = p3 FI
IF p5 < m THEN m = p5 FI
IF p7 < m THEN m = p7 FI
H( hPos ) = m
IF m = p2 THEN last2 = last2 + 1 p2 = 2 * H( last2 ) FI
IF m = p3 THEN last3 = last3 + 1 p3 = 3 * H( last3 ) FI
IF m = p5 THEN last5 = last5 + 1 p5 = 5 * H( last5 ) FI
IF m = p7 THEN last7 = last7 + 1 p7 = 7 * H( last7 ) FI
OD
FOR hPos = 0 TO 49 DO
Put(' )PrintC( H( hPos ) )
OD
PutE()
FOR hPos = 0 TO MAX_HUMBLE - 1 DO
m = H( hPos )
IF m < 10 THEN h1 = h1 + 1
ELSEIF m < 100 THEN h2 = h2 + 1
ELSEIF m < 1000 THEN h3 = h3 + 1
ELSEIF m < 10000 THEN h4 = h4 + 1
FI
OD
humbleStat( h1, 1 )
humbleStat( h2, 2 )
humbleStat( h3, 3 )
humbleStat( h4, 4 )
RETURN- Output:
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 There are 9 humble numbers with 1 digit There are 36 humble numbers with 2 digits There are 95 humble numbers with 3 digits There are 197 humble numbers with 4 digits
Ada
with Ada.Text_IO;
procedure Show_Humble is
type Positive is range 1 .. 2**63 - 1;
First : constant Positive := Positive'First;
Last : constant Positive := 999_999_999;
function Is_Humble (I : in Positive) return Boolean is
begin
if I <= 1 then return True;
elsif I mod 2 = 0 then return Is_Humble (I / 2);
elsif I mod 3 = 0 then return Is_Humble (I / 3);
elsif I mod 5 = 0 then return Is_Humble (I / 5);
elsif I mod 7 = 0 then return Is_Humble (I / 7);
else return False;
end if;
end Is_Humble;
subtype Digit_Range is Natural range First'Image'Length - 1 .. Last'Image'Length - 1;
Digit_Count : array (Digit_Range) of Natural := (others => 0);
procedure Count_Humble_Digits is
use Ada.Text_IO;
Humble_Count : Natural := 0;
Len : Natural;
begin
Put_Line ("The first 50 humble numbers:");
for N in First .. Last loop
if Is_Humble (N) then
Len := N'Image'Length - 1;
Digit_Count (Len) := Digit_Count (Len) + 1;
if Humble_Count < 50 then
Put (N'Image);
Put (" ");
end if;
Humble_Count := Humble_Count + 1;
end if;
end loop;
New_Line (2);
end Count_Humble_Digits;
procedure Show_Digit_Counts is
package Natural_IO is
new Ada.Text_IO.Integer_IO (Natural);
use Ada.Text_IO;
use Natural_IO;
Placeholder : String := "Digits Count";
Image_Digit : String renames Placeholder (1 .. 6);
Image_Count : String renames Placeholder (8 .. Placeholder'Last);
begin
Put_Line ("The digit counts of humble numbers:");
Put_Line (Placeholder);
for Digit in Digit_Count'Range loop
Put (Image_Digit, Digit);
Put (Image_Count, Digit_Count (Digit));
Put_Line (Placeholder);
end loop;
end Show_Digit_Counts;
begin
Count_Humble_Digits;
Show_Digit_Counts;
end Show_Humble;
- Output:
The 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
The digit counts of humble numbers:
Digits Count
1 9
2 36
3 95
4 197
5 356
6 579
7 882
8 1272
9 1767
ALGOL 68
BEGIN # find some Humble numbers - numbers with no prime factors above 7 #
INT max humble = 2048;
INT max shown humble = 50;
PROC min = ( INT a, b )INT: IF a < b THEN a ELSE b FI;
[ 1 : max humble ]INT h;
[ 0 : 6 ]INT h count;
FOR i FROM LWB h count TO UPB h count DO h count[ i ] := 0 OD;
INT p2 := 2, p3 := 3, p5 := 5, p7 := 7;
INT last2 := 1, last3 := 1, last5 := 1, last7 := 1;
# 1 is the first humble number ( 2^0 * 3^0 * 5^0 * 7^0 ) and has 1 digit #
h[ 1 ] := 1;
h count[ 1 ] := 1;
print( ( "1" ) );
FOR n FROM 2 TO max humble DO
# the next humble number is the lowest of the next multiples of #
# 2, 3, 5, 7 #
INT m = min( min( min( p2, p3 ), p5 ), p7 );
h[ n ] := m;
IF n <= max shown humble THEN print( ( " ", whole( m, 0 ) ) ) FI;
IF m = p2 THEN p2 := 2 * h[ last2 := last2 + 1 ] FI;
IF m = p3 THEN p3 := 3 * h[ last3 := last3 + 1 ] FI;
IF m = p5 THEN p5 := 5 * h[ last5 := last5 + 1 ] FI;
IF m = p7 THEN p7 := 7 * h[ last7 := last7 + 1 ] FI;
h count[ IF m < 10 THEN 1
ELIF m < 100 THEN 2
ELIF m < 1 000 THEN 3
ELIF m < 10 000 THEN 4
ELIF m < 100 000 THEN 5
ELIF m < 1 000 000 THEN 6
ELSE 0
FI
] +:= 1
OD;
print( ( newline ) );
FOR i TO 6 DO
print( ( "There are "
, whole( h count[ i ], -4 )
, " Humble numbers with "
, whole( i, 0 )
, " digits"
, newline
)
)
OD
END- Output:
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 There are 9 Humble numbers with 1 digits There are 36 Humble numbers with 2 digits There are 95 Humble numbers with 3 digits There are 197 Humble numbers with 4 digits There are 356 Humble numbers with 5 digits There are 579 Humble numbers with 6 digits
ALGOL W
As noted by other samples, this is similar to the Hamming Numbers task. This is a modified version of the Algol W solution for Hamming Numbers. The numbers are generated in sequence.
begin % find some Humble numbers - numbers with no prime factors above 7 %
% returns the minimum of a and b %
integer procedure min ( integer value a, b ) ; if a < b then a else b;
% find and print Humble Numbers %
integer MAX_HUMBLE;
MAX_HUMBLE := 2048;
begin
integer array H( 1 :: MAX_HUMBLE );
integer p2, p3, p5, p7, last2, last3, last5, last7, h1, h2, h3, h4, h5, h6;
i_w := 1; s_w := 1; % output formatting %
% 1 is the first Humble number %
H( 1 ) := 1;
h1 := h2 := h3 := h4 := h5 := h6 := 0;
last2 := last3 := last5 := last7 := 1;
p2 := 2;
p3 := 3;
p5 := 5;
p7 := 7;
for hPos := 2 until MAX_HUMBLE do begin
integer m;
% the next Humble number is the lowest of the next multiple of 2, 3, 5, 7 %
m := min( min( min( p2, p3 ), p5 ), p7 );
H( hPos ) := m;
if m = p2 then begin
% the Humble number was the next multiple of 2 %
% the next multiple of 2 will now be twice the Humble number following %
% the previous multple of 2 %
last2 := last2 + 1;
p2 := 2 * H( last2 )
end if_used_power_of_2 ;
if m = p3 then begin
last3 := last3 + 1;
p3 := 3 * H( last3 )
end if_used_power_of_3 ;
if m = p5 then begin
last5 := last5 + 1;
p5 := 5 * H( last5 )
end if_used_power_of_5 ;
if m = p7 then begin
last7 := last7 + 1;
p7 := 7 * H( last7 )
end if_used_power_of_5 ;
end for_hPos ;
i_w := 1; s_w := 1; % output formatting %
write( H( 1 ) );
for hPos := 2 until 50 do writeon( H( hPos ) );
for hPos := 1 until MAX_HUMBLE do begin
integer m;
m := H( hPos );
if m < 10 then h1 := h1 + 1
else if m < 100 then h2 := h2 + 1
else if m < 1000 then h3 := h3 + 1
else if m < 10000 then h4 := h4 + 1
else if m < 100000 then h5 := h5 + 1
else if m < 1000000 then h6 := h6 + 1
end for_hPos ;
i_w := 5; s_w := 0;
write( "there are", h1, " Humble numbers with 1 digit" );
write( "there are", h2, " Humble numbers with 2 digits" );
write( "there are", h3, " Humble numbers with 3 digits" );
write( "there are", h4, " Humble numbers with 4 digits" );
write( "there are", h5, " Humble numbers with 5 digits" );
write( "there are", h6, " Humble numbers with 6 digits" )
end
end.- Output:
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 there are 9 Humble numbers with 1 digit there are 36 Humble numbers with 2 digits there are 95 Humble numbers with 3 digits there are 197 Humble numbers with 4 digits there are 356 Humble numbers with 5 digits there are 579 Humble numbers with 6 digits
AutoHotkey
n := 1, c := 0
while (c < 50)
{
if isHumbleNumbers(prime_numbers(n))
c++, result .= n " "
n++
}
n := 1, l := 0, c:=[]
loop
{
if (l:=StrLen(n)) > 5
break
if isHumbleNumbers(prime_numbers(n))
c[l] := c[l] ? c[l] + 1 : 1
n++
}
for i, v in c
result .= "`n" i ":`t" v
MsgBox, 262144, ,% result
return
isHumbleNumbers(x){
for i, v in x
if v > 7
return false
return true
}
prime_numbers(n) {
if (n <= 3)
return [n]
ans := [], done := false
while !done {
if !Mod(n,2)
ans.push(2), n/=2
else if !Mod(n,3)
ans.push(3), n/=3
else if (n = 1)
return ans
else {
sr:=sqrt(n), done:=true, i:=6
; try to divide the checked number by all numbers till its square root.
while (i <= sr+6) {
if !Mod(n, i-1) { ; is n divisible by i-1?
ans.push(i-1), n/=i-1, done:=false
break
}
if !Mod(n, i+1) { ; is n divisible by i+1?
ans.push(i+1), n/=i+1, done:=false
break
}
i += 6
}
}
}
ans.push(n)
return ans
}
- Output:
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 1: 9 2: 36 3: 95 4: 197 5: 356
AWK
# syntax: GAWK -f HUMBLE_NUMBERS.AWK
#
# sorting:
# PROCINFO["sorted_in"] is used by GAWK
# SORTTYPE is used by Thompson Automation's TAWK
#
BEGIN {
PROCINFO["sorted_in"] = "@ind_num_asc" ; SORTTYPE = 1
n = 1
for (; count<5193; n++) {
if (is_humble(n)) {
arr[length(n)]++
if (count++ < 50) {
printf("%d ",n)
}
}
}
printf("\nCount Digits of the first %d humble numbers:\n",count)
for (i in arr) {
printf("%5d %6d\n",arr[i],i)
}
exit(0)
}
function is_humble(i) {
if (i <= 1) { return(1) }
if (i % 2 == 0) { return(is_humble(i/2)) }
if (i % 3 == 0) { return(is_humble(i/3)) }
if (i % 5 == 0) { return(is_humble(i/5)) }
if (i % 7 == 0) { return(is_humble(i/7)) }
return(0)
}
- Output:
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 Digits of the first 5193 humble numbers:
9 1
36 2
95 3
197 4
356 5
579 6
882 7
1272 8
1767 9
C
#include <limits.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>
bool isHumble(int i) {
if (i <= 1) return true;
if (i % 2 == 0) return isHumble(i / 2);
if (i % 3 == 0) return isHumble(i / 3);
if (i % 5 == 0) return isHumble(i / 5);
if (i % 7 == 0) return isHumble(i / 7);
return false;
}
int main() {
int limit = SHRT_MAX;
int humble[16];
int count = 0;
int num = 1;
char buffer[16];
memset(humble, 0, sizeof(humble));
for (; count < limit; num++) {
if (isHumble(num)) {
size_t len;
sprintf_s(buffer, sizeof(buffer), "%d", num);
len = strlen(buffer);
if (len >= 16) {
break;
}
humble[len]++;
if (count < 50) {
printf("%d ", num);
}
count++;
}
}
printf("\n\n");
printf("Of the first %d humble numbers:\n", count);
for (num = 1; num < 10; num++) {
printf("%5d have %d digits\n", humble[num], num);
}
return 0;
}
- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
C#
using System;
using System.Collections.Generic;
namespace HumbleNumbers {
class Program {
static bool IsHumble(int i) {
if (i <= 1) return true;
if (i % 2 == 0) return IsHumble(i / 2);
if (i % 3 == 0) return IsHumble(i / 3);
if (i % 5 == 0) return IsHumble(i / 5);
if (i % 7 == 0) return IsHumble(i / 7);
return false;
}
static void Main() {
var limit = short.MaxValue;
Dictionary<int, int> humble = new Dictionary<int, int>();
var count = 0;
var num = 1;
while (count < limit) {
if (IsHumble(num)) {
var str = num.ToString();
var len = str.Length;
if (humble.ContainsKey(len)) {
humble[len]++;
} else {
humble[len] = 1;
}
if (count < 50) Console.Write("{0} ", num);
count++;
}
num++;
}
Console.WriteLine("\n");
Console.WriteLine("Of the first {0} humble numbers:", count);
num = 1;
while (num < humble.Count - 1) {
if (humble.ContainsKey(num)) {
var c = humble[num];
Console.WriteLine("{0,5} have {1,2} digits", c, num);
num++;
} else {
break;
}
}
}
}
}
- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Direct Generation
#define BI
using System;
using System.Linq;
using System.Collections.Generic;
#if BI
using UI = System.Numerics.BigInteger;
#else
using UI = System.UInt64;
#endif
class Program {
static void Main(string[] args) {
#if BI
const int max = 100;
#else
const int max = 19;
#endif
List<UI> h = new List<UI> { 1 };
UI x2 = 2, x3 = 3, x5 = 5, x7 = 7, hm = 2, lim = 10;
int i = 0, j = 0, k = 0, l = 0, lc = 0, d = 1;
Console.WriteLine("Digits Count Time Mb used");
var elpsd = -DateTime.Now.Ticks;
do {
h.Add(hm);
if (hm == x2) x2 = h[++i] << 1;
if (hm == x3) x3 = (h[++j] << 1) + h[j];
if (hm == x5) x5 = (h[++k] << 2) + h[k];
if (hm == x7) x7 = (h[++l] << 3) - h[l];
hm = x2; if (x3 < hm) hm = x3; if (x5 < hm) hm = x5; if (x7 < hm) hm = x7;
if (hm >= lim) {
Console.WriteLine("{0,3} {1,9:n0} {2,9:n0} ms {3,9:n0}", d, h.Count - lc,
(elpsd + DateTime.Now.Ticks) / 10000, GC.GetTotalMemory(false) / 1000000);
lc = h.Count; if (++d > max) break; lim *= 10;
}
} while (true);
Console.WriteLine("{0,13:n0} Total", lc);
int firstAmt = 50;
Console.WriteLine("The first {0} humble numbers are: {1}", firstAmt, string.Join(" ",h.Take(firstAmt)));
}
}
- Output:
Results from a core i7-7700 @ 3.6Ghz.
BigIntegers: (tabulates up to 100 digits in about 3/4 of a minute, but a lot of memory is consumed - 4.2 GB)
Digits Count Time Mb used 1 9 5 ms 0 2 36 7 ms 0 3 95 7 ms 0 4 197 7 ms 0 5 356 7 ms 0 6 579 8 ms 0 7 882 8 ms 0 8 1,272 9 ms 0 9 1,767 10 ms 1 10 2,381 11 ms 2 11 3,113 14 ms 3 12 3,984 23 ms 1 13 5,002 27 ms 4 14 6,187 34 ms 2 15 7,545 39 ms 6 16 9,081 60 ms 4 17 10,815 75 ms 3 18 12,759 88 ms 11 19 14,927 105 ms 4 20 17,323 116 ms 7 21 19,960 157 ms 6 22 22,853 183 ms 16 23 26,015 217 ms 10 24 29,458 241 ms 14 25 33,188 279 ms 14 26 37,222 327 ms 24 27 41,568 365 ms 28 28 46,245 408 ms 30 29 51,254 472 ms 23 30 56,618 526 ms 26 31 62,338 607 ms 49 32 68,437 697 ms 39 33 74,917 762 ms 47 34 81,793 819 ms 48 35 89,083 894 ms 53 36 96,786 1,017 ms 56 37 104,926 1,114 ms 58 38 113,511 1,245 ms 99 39 122,546 1,350 ms 104 40 132,054 1,473 ms 107 41 142,038 1,640 ms 101 42 152,515 1,772 ms 106 43 163,497 1,902 ms 113 44 174,986 2,040 ms 121 45 187,004 2,240 ms 165 46 199,565 2,371 ms 178 47 212,675 2,501 ms 187 48 226,346 2,640 ms 194 49 240,590 2,792 ms 209 50 255,415 2,977 ms 223 51 270,843 3,246 ms 236 52 286,880 3,463 ms 248 53 303,533 3,745 ms 395 54 320,821 3,988 ms 414 55 338,750 4,203 ms 428 56 357,343 4,447 ms 443 57 376,599 4,734 ms 460 58 396,533 5,127 ms 418 59 417,160 5,442 ms 438 60 438,492 5,782 ms 464 61 460,533 6,139 ms 489 62 483,307 6,519 ms 514 63 506,820 6,918 ms 545 64 531,076 7,560 ms 706 65 556,104 7,986 ms 740 66 581,902 8,591 ms 771 67 608,483 9,056 ms 805 68 635,864 9,470 ms 843 69 664,053 9,988 ms 876 70 693,065 10,597 ms 918 71 722,911 11,102 ms 961 72 753,593 11,880 ms 1,000 73 785,141 12,471 ms 1,047 74 817,554 13,056 ms 1,092 75 850,847 13,767 ms 1,140 76 885,037 14,551 ms 1,724 77 920,120 15,362 ms 1,776 78 956,120 16,131 ms 1,834 79 993,058 16,986 ms 1,901 80 1,030,928 17,776 ms 1,967 81 1,069,748 18,658 ms 2,037 82 1,109,528 19,866 ms 1,839 83 1,150,287 20,780 ms 1,911 84 1,192,035 21,748 ms 1,985 85 1,234,774 22,715 ms 2,067 86 1,278,527 23,799 ms 2,147 87 1,323,301 24,826 ms 2,235 88 1,369,106 25,953 ms 2,322 89 1,415,956 27,119 ms 2,411 90 1,463,862 29,195 ms 3,041 91 1,512,840 30,487 ms 3,138 92 1,562,897 31,732 ms 3,241 93 1,614,050 32,995 ms 3,339 94 1,666,302 34,338 ms 3,451 95 1,719,669 35,809 ms 3,560 96 1,774,166 37,386 ms 3,673 97 1,829,805 38,912 ms 3,800 98 1,886,590 40,474 ms 3,933 99 1,944,540 42,073 ms 4,077 100 2,003,661 43,808 ms 4,222 51,428,827 Total The first 50 humble numbers are: 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
UInt64s: (comment out "#define BI" at the top of the code)
Digits Count Time Mb used
1 9 4 ms 0
2 36 5 ms 0
3 95 6 ms 0
4 197 6 ms 0
5 356 6 ms 0
6 579 6 ms 0
7 882 6 ms 0
8 1,272 6 ms 0
9 1,767 6 ms 0
10 2,381 6 ms 0
11 3,113 6 ms 0
12 3,984 6 ms 0
13 5,002 6 ms 0
14 6,187 6 ms 0
15 7,545 7 ms 1
16 9,081 7 ms 1
17 10,815 7 ms 1
18 12,759 7 ms 2
19 14,927 7 ms 2
80,987 Total
The first 50 humble numbers are: 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
Direct Generation via Logarithms
Similar to one of the design elements of the Pascal entry (and a few others), add logarithms together instead of multiplying big numbers. Surprisingly, only about 10-11 digits of precision is needed, so the fixed-point logs fit in an UInt64. It's a very bad memory hog though (17GB!), so the Pascal entry is much better in that respect. It's quick, doing 255 digits in 17 seconds (about 60x faster than the Direct Generation BigInteger version above), comparable to the speed of the Pascal vesion. It does double duty (range 1..nth and digits tabulation), which slows performance a little. When the code that generates the range (1..nth) is removed, it can execute in about 15 seconds. (on the core i7-7700 @ 3.6Ghz)
It does have an issue with reporting humble numbers greater than ~1e19, as the code that returns the fixed-point log cannot fit the result into a UInt64. This is a non-issue for this task because the largest humble number asked for is 120. (Heh, on the Hamming Number task, that would be an issue.) However, the count of humble numbers in each decade is correct. If it is necessary to report large humble numbers, the System.Numerics library could be used and a function written to provide an arbitrary precision BigInteger.Exp() result.
Why use fixed-point logarithms of UIint64 instead of double? Because the rounding of the doubles when added together throws the sums off a bit so they don't match properly when incrementing the i, j, k, & l variables. If one were to change the 'fac" variable to a larger number, such as 1e15, there is too much "noise" on the least significant bits and the ijkl variables advance unevenly enough to throw off some of the counts. Some of the solutions presented here implement an "error banding" check to defeat this issue, but it seems a bit over complicated.
using System;
using UI = System.UInt64;
class Program {
// write a range (1..num) to the console when num < 0, just write the nth when num > 0, otherwise write the digits tabulation
// note: when doing range or nth, if num > ~1e19 the results will appear incorrect as UInt64 can't express numbers that large
static void humLog(int digs, int num = 0) {
bool range = num < 0, nth = num > 0, small = range | nth; num = Math.Abs(num);
int maxdim = num;
if (range | nth) digs = num.ToString().Length; // calculate number of digits when range or nth is specified
//const int maxdim = 2_147_483_647; // 2GB limit (Int32.MaxValue), causes out of memory error
//const int maxdim = 2_146_435_071; // max practical amount
//const int maxdim = 2_114_620_032; // amount needed for 255 digits
else maxdim = 2_114_620_032;
const double fac = 1e11;
UI lb2 = (UI)Math.Round(fac * Math.Log(2)), lb3 = (UI)Math.Round(fac * Math.Log(3)), lb5 = (UI)Math.Round(fac * Math.Log(5)),
lb7 = (UI)Math.Round(fac * Math.Log(7)), lb0 = (UI)Math.Round(fac * Math.Log(10)), hm,
x2 = lb2, x3 = lb3, x5 = lb5, x7 = lb7, lim = lb0;
int i = 0, j = 0, k = 0, l = 0, lc = 0, d = 1, hi = 1;
UI[] h = new UI[maxdim]; h[0] = 1;
var st = DateTime.Now.Ticks;
if (range) Console.Write("The first {0} humble numbers are: 1 ", num);
else if (nth) Console.Write("The {0}{1} humble number is ", num, (num % 10) switch { 1 => "st", 2 => "nd", 3 => "rd", _ => "th", });
else Console.WriteLine("\nDigits Dig Count Tot Count Time Mb used");
do { hm = x2; if (x3 < hm) hm = x3; if (x5 < hm) hm = x5; if (x7 < hm) hm = x7; // select the minimum
if (hm >= lim && !small) { // passed another decade, so output results
Console.WriteLine("{0,3} {1,13:n0} {4,16:n0} {2,9:n3}s {3,9:n0}", d, hi - lc,
((DateTime.Now.Ticks - st) / 10000)/1000.0, GC.GetTotalMemory(false) / 1000000, hi);
lc = hi; if (++d > digs) break; lim += lb0; }
h[hi++] = (hm); if (small) { if (nth && hi == num) { Console.WriteLine(Math.Round(Math.Exp(hm / fac))); break; }
if (range) { Console.Write("{0} ", Math.Round(Math.Exp(hm / fac))); if (hi == num) { Console.WriteLine(); break; } } }
if (hm == x2) x2 = h[++i] + lb2; if (hm == x3) x3 = h[++j] + lb3;
if (hm == x5) x5 = h[++k] + lb5; if (hm == x7) x7 = h[++l] + lb7;
} while (true);
if (!(range | nth)) Console.WriteLine("{0,17:n0} Total", lc);
}
static void Main(string[] args) {
humLog(0, -50); // see the range 1..50
humLog(255); // see tabulation for digits 1 to 255
}
}
- Output:
verified results against the Pascal entry:
The first 50 humble numbers are: 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
Digits Dig Count Tot Count Time Mb used
1 9 9 0.000s 16,917
2 36 45 0.038s 16,917
3 95 140 0.039s 16,917
4 197 337 0.039s 16,917
5 356 693 0.039s 16,917
6 579 1,272 0.039s 16,917
7 882 2,154 0.039s 16,917
8 1,272 3,426 0.039s 16,917
9 1,767 5,193 0.040s 16,917
10 2,381 7,574 0.040s 16,917
11 3,113 10,687 0.040s 16,917
12 3,984 14,671 0.040s 16,917
13 5,002 19,673 0.040s 16,917
14 6,187 25,860 0.041s 16,917
15 7,545 33,405 0.041s 16,917
16 9,081 42,486 0.041s 16,917
17 10,815 53,301 0.041s 16,917
18 12,759 66,060 0.041s 16,917
19 14,927 80,987 0.042s 16,917
20 17,323 98,310 0.042s 16,917
21 19,960 118,270 0.043s 16,917
22 22,853 141,123 0.044s 16,917
23 26,015 167,138 0.045s 16,917
24 29,458 196,596 0.045s 16,917
25 33,188 229,784 0.047s 16,917
26 37,222 267,006 0.047s 16,917
27 41,568 308,574 0.048s 16,917
28 46,245 354,819 0.049s 16,917
29 51,254 406,073 0.050s 16,917
30 56,618 462,691 0.051s 16,917
31 62,338 525,029 0.052s 16,917
32 68,437 593,466 0.088s 16,917
33 74,917 668,383 0.094s 16,917
34 81,793 750,176 0.095s 16,917
35 89,083 839,259 0.096s 16,917
36 96,786 936,045 0.097s 16,917
37 104,926 1,040,971 0.098s 16,917
38 113,511 1,154,482 0.099s 16,917
39 122,546 1,277,028 0.100s 16,917
40 132,054 1,409,082 0.101s 16,917
41 142,038 1,551,120 0.102s 16,917
42 152,515 1,703,635 0.104s 16,917
43 163,497 1,867,132 0.106s 16,917
44 174,986 2,042,118 0.109s 16,917
45 187,004 2,229,122 0.111s 16,917
46 199,565 2,428,687 0.113s 16,917
47 212,675 2,641,362 0.116s 16,917
48 226,346 2,867,708 0.118s 16,917
49 240,590 3,108,298 0.120s 16,917
50 255,415 3,363,713 0.123s 16,917
51 270,843 3,634,556 0.145s 16,917
52 286,880 3,921,436 0.169s 16,917
53 303,533 4,224,969 0.172s 16,917
54 320,821 4,545,790 0.176s 16,917
55 338,750 4,884,540 0.181s 16,917
56 357,343 5,241,883 0.225s 16,917
57 376,599 5,618,482 0.229s 16,917
58 396,533 6,015,015 0.233s 16,917
59 417,160 6,432,175 0.238s 16,917
60 438,492 6,870,667 0.275s 16,917
61 460,533 7,331,200 0.279s 16,917
62 483,307 7,814,507 0.285s 16,917
63 506,820 8,321,327 0.290s 16,917
64 531,076 8,852,403 0.299s 16,917
65 556,104 9,408,507 0.319s 16,917
66 581,902 9,990,409 0.341s 16,917
67 608,483 10,598,892 0.348s 16,917
68 635,864 11,234,756 0.370s 16,917
69 664,053 11,898,809 0.393s 16,917
70 693,065 12,591,874 0.413s 16,917
71 722,911 13,314,785 0.436s 16,917
72 753,593 14,068,378 0.457s 16,917
73 785,141 14,853,519 0.480s 16,917
74 817,554 15,671,073 0.488s 16,917
75 850,847 16,521,920 0.497s 16,917
76 885,037 17,406,957 0.506s 16,917
77 920,120 18,327,077 0.514s 16,917
78 956,120 19,283,197 0.522s 16,917
79 993,058 20,276,255 0.531s 16,917
80 1,030,928 21,307,183 0.540s 16,917
81 1,069,748 22,376,931 0.548s 16,917
82 1,109,528 23,486,459 0.559s 16,917
83 1,150,287 24,636,746 0.570s 16,917
84 1,192,035 25,828,781 0.580s 16,917
85 1,234,774 27,063,555 0.597s 16,917
86 1,278,527 28,342,082 0.608s 16,917
87 1,323,301 29,665,383 0.622s 16,917
88 1,369,106 31,034,489 0.637s 16,917
89 1,415,956 32,450,445 0.652s 16,917
90 1,463,862 33,914,307 0.665s 16,917
91 1,512,840 35,427,147 0.678s 16,917
92 1,562,897 36,990,044 0.692s 16,917
93 1,614,050 38,604,094 0.705s 16,917
94 1,666,302 40,270,396 0.720s 16,917
95 1,719,669 41,990,065 0.734s 16,917
96 1,774,166 43,764,231 0.749s 16,917
97 1,829,805 45,594,036 0.767s 16,917
98 1,886,590 47,480,626 0.783s 16,917
99 1,944,540 49,425,166 0.799s 16,917
100 2,003,661 51,428,827 0.816s 16,917
101 2,063,972 53,492,799 0.834s 16,917
102 2,125,486 55,618,285 0.855s 16,917
103 2,188,204 57,806,489 0.873s 16,917
104 2,252,146 60,058,635 0.897s 16,917
105 2,317,319 62,375,954 0.922s 16,917
106 2,383,733 64,759,687 0.944s 16,917
107 2,451,413 67,211,100 0.965s 16,917
108 2,520,360 69,731,460 0.990s 16,917
109 2,590,584 72,322,044 1.010s 16,917
110 2,662,102 74,984,146 1.031s 16,917
111 2,734,927 77,719,073 1.055s 16,917
112 2,809,069 80,528,142 1.082s 16,917
113 2,884,536 83,412,678 1.107s 16,917
114 2,961,346 86,374,024 1.134s 16,917
115 3,039,502 89,413,526 1.158s 16,917
116 3,119,022 92,532,548 1.188s 16,917
117 3,199,918 95,732,466 1.216s 16,917
118 3,282,203 99,014,669 1.243s 16,917
119 3,365,883 102,380,552 1.273s 16,917
120 3,450,981 105,831,533 1.303s 16,917
121 3,537,499 109,369,032 1.332s 16,917
122 3,625,444 112,994,476 1.371s 16,917
123 3,714,838 116,709,314 1.406s 16,917
124 3,805,692 120,515,006 1.435s 16,917
125 3,898,015 124,413,021 1.466s 16,917
126 3,991,818 128,404,839 1.503s 16,917
127 4,087,110 132,491,949 1.536s 16,917
128 4,183,914 136,675,863 1.573s 16,917
129 4,282,228 140,958,091 1.614s 16,917
130 4,382,079 145,340,170 1.654s 16,917
131 4,483,467 149,823,637 1.693s 16,917
132 4,586,405 154,410,042 1.732s 16,917
133 4,690,902 159,100,944 1.770s 16,917
134 4,796,979 163,897,923 1.808s 16,917
135 4,904,646 168,802,569 1.848s 16,917
136 5,013,909 173,816,478 1.887s 16,917
137 5,124,783 178,941,261 1.928s 16,917
138 5,237,275 184,178,536 1.969s 16,917
139 5,351,407 189,529,943 2.012s 16,917
140 5,467,187 194,997,130 2.055s 16,917
141 5,584,624 200,581,754 2.098s 16,917
142 5,703,728 206,285,482 2.143s 16,917
143 5,824,512 212,109,994 2.189s 16,917
144 5,946,992 218,056,986 2.236s 16,917
145 6,071,177 224,128,163 2.284s 16,917
146 6,197,080 230,325,243 2.331s 16,917
147 6,324,708 236,649,951 2.377s 16,917
148 6,454,082 243,104,033 2.426s 16,917
149 6,585,205 249,689,238 2.475s 16,917
150 6,718,091 256,407,329 2.524s 16,917
151 6,852,749 263,260,078 2.577s 16,917
152 6,989,204 270,249,282 2.629s 16,917
153 7,127,454 277,376,736 2.681s 16,917
154 7,267,511 284,644,247 2.736s 16,917
155 7,409,395 292,053,642 2.799s 16,917
156 7,553,112 299,606,754 2.863s 16,917
157 7,698,677 307,305,431 2.926s 16,917
158 7,846,103 315,151,534 2.989s 16,917
159 7,995,394 323,146,928 3.054s 16,917
160 8,146,567 331,293,495 3.125s 16,917
161 8,299,638 339,593,133 3.190s 16,917
162 8,454,607 348,047,740 3.257s 16,917
163 8,611,505 356,659,245 3.324s 16,917
164 8,770,324 365,429,569 3.394s 16,917
165 8,931,081 374,360,650 3.476s 16,917
166 9,093,797 383,454,447 3.546s 16,917
167 9,258,476 392,712,923 3.619s 16,917
168 9,425,127 402,138,050 3.693s 16,917
169 9,593,778 411,731,828 3.766s 16,917
170 9,764,417 421,496,245 3.843s 16,917
171 9,937,068 431,433,313 3.927s 16,917
172 10,111,745 441,545,058 4.003s 16,917
173 10,288,458 451,833,516 4.087s 16,917
174 10,467,215 462,300,731 4.167s 16,917
175 10,648,032 472,948,763 4.245s 16,917
176 10,830,920 483,779,683 4.329s 16,917
177 11,015,896 494,795,579 4.409s 16,917
178 11,202,959 505,998,538 4.493s 16,917
179 11,392,128 517,390,666 4.580s 16,917
180 11,583,420 528,974,086 4.672s 16,917
181 11,776,838 540,750,924 4.760s 16,917
182 11,972,395 552,723,319 4.855s 16,917
183 12,170,108 564,893,427 4.954s 16,917
184 12,369,985 577,263,412 5.047s 16,917
185 12,572,037 589,835,449 5.143s 16,917
186 12,776,285 602,611,734 5.249s 16,917
187 12,982,725 615,594,459 5.348s 16,917
188 13,191,377 628,785,836 5.454s 16,917
189 13,402,256 642,188,092 5.558s 16,917
190 13,615,367 655,803,459 5.667s 16,917
191 13,830,730 669,634,189 5.772s 16,917
192 14,048,347 683,682,536 5.883s 16,917
193 14,268,236 697,950,772 6.002s 16,917
194 14,490,415 712,441,187 6.109s 16,917
195 14,714,880 727,156,067 6.215s 16,917
196 14,941,651 742,097,718 6.324s 16,917
197 15,170,748 757,268,466 6.436s 16,917
198 15,402,165 772,670,631 6.559s 16,917
199 15,635,928 788,306,559 6.684s 16,917
200 15,872,045 804,178,604 6.809s 16,917
201 16,110,527 820,289,131 6.941s 16,917
202 16,351,384 836,640,515 7.068s 16,917
203 16,594,632 853,235,147 7.199s 16,917
204 16,840,283 870,075,430 7.332s 16,917
205 17,088,342 887,163,772 7.462s 16,917
206 17,338,826 904,502,598 7.596s 16,917
207 17,591,739 922,094,337 7.737s 16,917
208 17,847,107 939,941,444 7.868s 16,917
209 18,104,934 958,046,378 8.002s 16,917
210 18,365,234 976,411,612 8.138s 16,917
211 18,628,013 995,039,625 8.275s 16,917
212 18,893,289 1,013,932,914 8.426s 16,917
213 19,161,068 1,033,093,982 8.578s 16,917
214 19,431,375 1,052,525,357 8.727s 16,917
215 19,704,205 1,072,229,562 8.878s 16,917
216 19,979,576 1,092,209,138 9.042s 16,917
217 20,257,500 1,112,466,638 9.199s 16,917
218 20,537,988 1,133,004,626 9.357s 16,917
219 20,821,062 1,153,825,688 9.522s 16,917
220 21,106,720 1,174,932,408 9.680s 16,917
221 21,394,982 1,196,327,390 9.836s 16,917
222 21,685,859 1,218,013,249 9.997s 16,917
223 21,979,347 1,239,992,596 10.162s 16,917
224 22,275,484 1,262,268,080 10.331s 16,917
225 22,574,265 1,284,842,345 10.510s 16,917
226 22,875,700 1,307,718,045 10.700s 16,917
227 23,179,816 1,330,897,861 10.888s 16,917
228 23,486,609 1,354,384,470 11.079s 16,917
229 23,796,098 1,378,180,568 11.262s 16,917
230 24,108,300 1,402,288,868 11.453s 16,917
231 24,423,216 1,426,712,084 11.629s 16,917
232 24,740,870 1,451,452,954 11.809s 16,917
233 25,061,260 1,476,514,214 11.994s 16,917
234 25,384,397 1,501,898,611 12.202s 16,917
235 25,710,307 1,527,608,918 12.406s 16,917
236 26,038,994 1,553,647,912 12.616s 16,917
237 26,370,474 1,580,018,386 12.831s 16,917
238 26,704,760 1,606,723,146 13.049s 16,917
239 27,041,843 1,633,764,989 13.256s 16,917
240 27,381,757 1,661,146,746 13.453s 16,917
241 27,724,512 1,688,871,258 13.655s 16,917
242 28,070,118 1,716,941,376 13.871s 16,917
243 28,418,579 1,745,359,955 14.094s 16,917
244 28,769,910 1,774,129,865 14.315s 16,917
245 29,124,123 1,803,253,988 14.540s 16,917
246 29,481,235 1,832,735,223 14.768s 16,917
247 29,841,260 1,862,576,483 15.005s 16,917
248 30,204,196 1,892,780,679 15.231s 16,917
249 30,570,067 1,923,350,746 15.453s 16,917
250 30,938,881 1,954,289,627 15.694s 16,917
251 31,310,645 1,985,600,272 15.941s 16,917
252 31,685,379 2,017,285,651 16.208s 16,917
253 32,063,093 2,049,348,744 16.456s 16,917
254 32,443,792 2,081,792,536 16.702s 16,917
255 32,827,496 2,114,620,032 16.952s 16,917
2,114,620,032 Total
C++
#include <iomanip>
#include <iostream>
#include <map>
#include <sstream>
bool isHumble(int i) {
if (i <= 1) return true;
if (i % 2 == 0) return isHumble(i / 2);
if (i % 3 == 0) return isHumble(i / 3);
if (i % 5 == 0) return isHumble(i / 5);
if (i % 7 == 0) return isHumble(i / 7);
return false;
}
auto toString(int n) {
std::stringstream ss;
ss << n;
return ss.str();
}
int main() {
auto limit = SHRT_MAX;
std::map<int, int> humble;
auto count = 0;
auto num = 1;
while (count < limit) {
if (isHumble(num)) {
auto str = toString(num);
auto len = str.length();
auto it = humble.find(len);
if (it != humble.end()) {
it->second++;
} else {
humble[len] = 1;
}
if (count < 50) std::cout << num << ' ';
count++;
}
num++;
}
std::cout << "\n\n";
std::cout << "Of the first " << count << " humble numbers:\n";
num = 1;
while (num < humble.size() - 1) {
auto it = humble.find(num);
if (it != humble.end()) {
auto c = *it;
std::cout << std::setw(5) << c.second << " have " << std::setw(2) << num << " digits\n";
num++;
} else {
break;
}
}
return 0;
}
- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Direct Generation - Variant
A direct generation variant. Rather quick, as the humble numbers are not generated in order. And the digits are not counted individually, the log representation of each humble number is just binned into the decade tally with a simple division by log(10). Note: g++ compiler options: -O3 -std=c++17
#include <chrono>
#include <cmath>
#include <locale>
using UI = unsigned long long;
using namespace std;
using namespace chrono;
int limc = 877;
const double l2 = log(2), l3 = log(3), l5 = log(5), l7 = log(7), l0 = log(10), fac = 1e12;
static bool IsHum(int num) { if (num <= 1) return true; for (int j : { 2, 3, 5, 7 })
if (num % j == 0) return IsHum(num / j); return false; }
// slow way to determine whether numbers are humble numbers
static void First_Slow(int firstAmt) { printf("The first %d humble numbers are: ", firstAmt);
for (int gg = 0, g = 1; gg < firstAmt; g++) if (IsHum(g)) { printf("%d ", g); gg++; }
printf("\n\n"); }
int main(int argc, char **argv) {
First_Slow(50); setlocale(LC_ALL, ""); auto st = steady_clock::now();
if (argc > 1) limc = stoi(argv[1]);
UI *bins = new UI[limc], lb0 = (UI)round(fac * l0),
lb2 = (UI)round(fac * l2), lb3 = (UI)round(fac * l3),
lb5 = (UI)round(fac * l5), lb7 = (UI)round(fac * l7),
tot = 0, lmt = limc * lb0, lm2 = lb5 * 3;
printf("Digits Count Accum\n");
for (int g = 0; g < limc; g++) bins[g] = 0;
for (UI i = 0; i < lmt; i += lb2) for (UI j = i; j < lmt; j += lb3)
for (UI k = j; k < lmt; k += lb5) for (UI l = k; l < lmt; l += lb7)
bins[l / lb0]++;
for (int f = 0, g = 1; f < limc; f = g++) { tot += bins[f];
//if (g < 110 || g % 100 == 0 || (g < 200 && g % 10 == 0)) // uncomment to emulate pascal output
printf ("%4d %'13llu %'18llu\n", g, bins[f], tot); }
delete [] bins;
printf("Counting took %8f seconds\n", duration<double>(steady_clock::now() - st).count());
}
- Output:
Seems to give correct values as compared to the pascal (modification of hamming numbers fast alternative) version. And goes noticeably faster, up to 877 digits in about 3 1/4 minutes, where as pascal takes 1 1/3 hours to get to 877 digits.
The first 50 humble numbers are: 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 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 26 37,222 267,006 27 41,568 308,574 28 46,245 354,819 29 51,254 406,073 30 56,618 462,691 31 62,338 525,029 32 68,437 593,466 33 74,917 668,383 34 81,793 750,176 35 89,083 839,259 36 96,786 936,045 37 104,926 1,040,971 38 113,511 1,154,482 39 122,546 1,277,028 40 132,054 1,409,082 41 142,038 1,551,120 42 152,515 1,703,635 43 163,497 1,867,132 44 174,986 2,042,118 45 187,004 2,229,122 46 199,565 2,428,687 47 212,675 2,641,362 48 226,346 2,867,708 49 240,590 3,108,298 50 255,415 3,363,713 51 270,843 3,634,556 52 286,880 3,921,436 53 303,533 4,224,969 54 320,821 4,545,790 55 338,750 4,884,540 56 357,343 5,241,883 57 376,599 5,618,482 58 396,533 6,015,015 59 417,160 6,432,175 60 438,492 6,870,667 61 460,533 7,331,200 62 483,307 7,814,507 63 506,820 8,321,327 64 531,076 8,852,403 65 556,104 9,408,507 66 581,902 9,990,409 67 608,483 10,598,892 68 635,864 11,234,756 69 664,053 11,898,809 70 693,065 12,591,874 71 722,911 13,314,785 72 753,593 14,068,378 73 785,141 14,853,519 74 817,554 15,671,073 75 850,847 16,521,920 76 885,037 17,406,957 77 920,120 18,327,077 78 956,120 19,283,197 79 993,058 20,276,255 80 1,030,928 21,307,183 81 1,069,748 22,376,931 82 1,109,528 23,486,459 83 1,150,287 24,636,746 84 1,192,035 25,828,781 85 1,234,774 27,063,555 86 1,278,527 28,342,082 87 1,323,301 29,665,383 88 1,369,106 31,034,489 89 1,415,956 32,450,445 90 1,463,862 33,914,307 91 1,512,840 35,427,147 92 1,562,897 36,990,044 93 1,614,050 38,604,094 94 1,666,302 40,270,396 95 1,719,669 41,990,065 96 1,774,166 43,764,231 97 1,829,805 45,594,036 98 1,886,590 47,480,626 99 1,944,540 49,425,166 100 2,003,661 51,428,827 101 2,063,972 53,492,799 102 2,125,486 55,618,285 103 2,188,204 57,806,489 104 2,252,146 60,058,635 105 2,317,319 62,375,954 106 2,383,733 64,759,687 107 2,451,413 67,211,100 108 2,520,360 69,731,460 109 2,590,584 72,322,044 110 2,662,102 74,984,146 111 2,734,927 77,719,073 112 2,809,069 80,528,142 113 2,884,536 83,412,678 114 2,961,346 86,374,024 115 3,039,502 89,413,526 116 3,119,022 92,532,548 117 3,199,918 95,732,466 118 3,282,203 99,014,669 119 3,365,883 102,380,552 120 3,450,981 105,831,533 121 3,537,499 109,369,032 122 3,625,444 112,994,476 123 3,714,838 116,709,314 124 3,805,692 120,515,006 125 3,898,015 124,413,021 126 3,991,818 128,404,839 127 4,087,110 132,491,949 128 4,183,914 136,675,863 129 4,282,228 140,958,091 130 4,382,079 145,340,170 131 4,483,467 149,823,637 132 4,586,405 154,410,042 133 4,690,902 159,100,944 134 4,796,979 163,897,923 135 4,904,646 168,802,569 136 5,013,909 173,816,478 137 5,124,783 178,941,261 138 5,237,275 184,178,536 139 5,351,407 189,529,943 140 5,467,187 194,997,130 141 5,584,624 200,581,754 142 5,703,728 206,285,482 143 5,824,512 212,109,994 144 5,946,992 218,056,986 145 6,071,177 224,128,163 146 6,197,080 230,325,243 147 6,324,708 236,649,951 148 6,454,082 243,104,033 149 6,585,205 249,689,238 150 6,718,091 256,407,329 151 6,852,749 263,260,078 152 6,989,204 270,249,282 153 7,127,454 277,376,736 154 7,267,511 284,644,247 155 7,409,395 292,053,642 156 7,553,112 299,606,754 157 7,698,677 307,305,431 158 7,846,103 315,151,534 159 7,995,394 323,146,928 160 8,146,567 331,293,495 161 8,299,638 339,593,133 162 8,454,607 348,047,740 163 8,611,505 356,659,245 164 8,770,324 365,429,569 165 8,931,081 374,360,650 166 9,093,797 383,454,447 167 9,258,476 392,712,923 168 9,425,127 402,138,050 169 9,593,778 411,731,828 170 9,764,417 421,496,245 171 9,937,068 431,433,313 172 10,111,745 441,545,058 173 10,288,458 451,833,516 174 10,467,215 462,300,731 175 10,648,032 472,948,763 176 10,830,920 483,779,683 177 11,015,896 494,795,579 178 11,202,959 505,998,538 179 11,392,128 517,390,666 180 11,583,420 528,974,086 181 11,776,838 540,750,924 182 11,972,395 552,723,319 183 12,170,108 564,893,427 184 12,369,985 577,263,412 185 12,572,037 589,835,449 186 12,776,285 602,611,734 187 12,982,725 615,594,459 188 13,191,377 628,785,836 189 13,402,256 642,188,092 190 13,615,367 655,803,459 191 13,830,730 669,634,189 192 14,048,347 683,682,536 193 14,268,236 697,950,772 194 14,490,415 712,441,187 195 14,714,880 727,156,067 196 14,941,651 742,097,718 197 15,170,748 757,268,466 198 15,402,165 772,670,631 199 15,635,928 788,306,559 200 15,872,045 804,178,604 201 16,110,527 820,289,131 202 16,351,384 836,640,515 203 16,594,632 853,235,147 204 16,840,283 870,075,430 205 17,088,342 887,163,772 206 17,338,826 904,502,598 207 17,591,739 922,094,337 208 17,847,107 939,941,444 209 18,104,934 958,046,378 210 18,365,234 976,411,612 211 18,628,013 995,039,625 212 18,893,289 1,013,932,914 213 19,161,068 1,033,093,982 214 19,431,375 1,052,525,357 215 19,704,205 1,072,229,562 216 19,979,576 1,092,209,138 217 20,257,500 1,112,466,638 218 20,537,988 1,133,004,626 219 20,821,062 1,153,825,688 220 21,106,720 1,174,932,408 221 21,394,982 1,196,327,390 222 21,685,859 1,218,013,249 223 21,979,347 1,239,992,596 224 22,275,484 1,262,268,080 225 22,574,265 1,284,842,345 226 22,875,700 1,307,718,045 227 23,179,816 1,330,897,861 228 23,486,609 1,354,384,470 229 23,796,098 1,378,180,568 230 24,108,300 1,402,288,868 231 24,423,216 1,426,712,084 232 24,740,870 1,451,452,954 233 25,061,260 1,476,514,214 234 25,384,397 1,501,898,611 235 25,710,307 1,527,608,918 236 26,038,994 1,553,647,912 237 26,370,474 1,580,018,386 238 26,704,760 1,606,723,146 239 27,041,843 1,633,764,989 240 27,381,757 1,661,146,746 241 27,724,512 1,688,871,258 242 28,070,118 1,716,941,376 243 28,418,579 1,745,359,955 244 28,769,910 1,774,129,865 245 29,124,123 1,803,253,988 246 29,481,235 1,832,735,223 247 29,841,260 1,862,576,483 248 30,204,196 1,892,780,679 249 30,570,067 1,923,350,746 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 256 33,214,214 2,147,834,246 257 33,603,951 2,181,438,197 258 33,996,731 2,215,434,928 259 34,392,563 2,249,827,491 260 34,791,447 2,284,618,938 261 35,193,409 2,319,812,347 262 35,598,451 2,355,410,798 263 36,006,588 2,391,417,386 264 36,417,835 2,427,835,221 265 36,832,207 2,464,667,428 266 37,249,701 2,501,917,129 267 37,670,346 2,539,587,475 268 38,094,140 2,577,681,615 269 38,521,101 2,616,202,716 270 38,951,249 2,655,153,965 271 39,384,582 2,694,538,547 272 39,821,110 2,734,359,657 273 40,260,860 2,774,620,517 274 40,703,832 2,815,324,349 275 41,150,044 2,856,474,393 276 41,599,505 2,898,073,898 277 42,052,225 2,940,126,123 278 42,508,220 2,982,634,343 279 42,967,499 3,025,601,842 280 43,430,073 3,069,031,915 281 43,895,955 3,112,927,870 282 44,365,160 3,157,293,030 283 44,837,693 3,202,130,723 284 45,313,573 3,247,444,296 285 45,792,800 3,293,237,096 286 46,275,398 3,339,512,494 287 46,761,379 3,386,273,873 288 47,250,752 3,433,524,625 289 47,743,528 3,481,268,153 290 48,239,713 3,529,507,866 291 48,739,321 3,578,247,187 292 49,242,372 3,627,489,559 293 49,748,878 3,677,238,437 294 50,258,835 3,727,497,272 295 50,772,270 3,778,269,542 296 51,289,181 3,829,558,723 297 51,809,596 3,881,368,319 298 52,333,528 3,933,701,847 299 52,860,969 3,986,562,816 300 53,391,941 4,039,954,757 301 53,926,465 4,093,881,222 302 54,464,535 4,148,345,757 303 55,006,178 4,203,351,935 304 55,551,405 4,258,903,340 305 56,100,219 4,315,003,559 306 56,652,634 4,371,656,193 307 57,208,666 4,428,864,859 308 57,768,319 4,486,633,178 309 58,331,610 4,544,964,788 310 58,898,562 4,603,863,350 311 59,469,169 4,663,332,519 312 60,043,451 4,723,375,970 313 60,621,404 4,783,997,374 314 61,203,070 4,845,200,444 315 61,788,445 4,906,988,889 316 62,377,542 4,969,366,431 317 62,970,363 5,032,336,794 318 63,566,920 5,095,903,714 319 64,167,245 5,160,070,959 320 64,771,341 5,224,842,300 321 65,379,217 5,290,221,517 322 65,990,881 5,356,212,398 323 66,606,345 5,422,818,743 324 67,225,617 5,490,044,360 325 67,848,732 5,557,893,092 326 68,475,682 5,626,368,774 327 69,106,479 5,695,475,253 328 69,741,139 5,765,216,392 329 70,379,662 5,835,596,054 330 71,022,084 5,906,618,138 331 71,668,401 5,978,286,539 332 72,318,631 6,050,605,170 333 72,972,775 6,123,577,945 334 73,630,849 6,197,208,794 335 74,292,868 6,271,501,662 336 74,958,848 6,346,460,510 337 75,628,800 6,422,089,310 338 76,302,729 6,498,392,039 339 76,980,651 6,575,372,690 340 77,662,572 6,653,035,262 341 78,348,501 6,731,383,763 342 79,038,476 6,810,422,239 343 79,732,485 6,890,154,724 344 80,430,535 6,970,585,259 345 81,132,652 7,051,717,911 346 81,838,839 7,133,556,750 347 82,549,112 7,216,105,862 348 83,263,504 7,299,369,366 349 83,981,986 7,383,351,352 350 84,704,585 7,468,055,937 351 85,431,326 7,553,487,263 352 86,162,202 7,639,649,465 353 86,897,253 7,726,546,718 354 87,636,466 7,814,183,184 355 88,379,852 7,902,563,036 356 89,127,430 7,991,690,466 357 89,879,218 8,081,569,684 358 90,635,220 8,172,204,904 359 91,395,454 8,263,600,358 360 92,159,923 8,355,760,281 361 92,928,637 8,448,688,918 362 93,701,621 8,542,390,539 363 94,478,883 8,636,869,422 364 95,260,421 8,732,129,843 365 96,046,268 8,828,176,111 366 96,836,421 8,925,012,532 367 97,630,889 9,022,643,421 368 98,429,696 9,121,073,117 369 99,232,843 9,220,305,960 370 100,040,353 9,320,346,313 371 100,852,238 9,421,198,551 372 101,668,499 9,522,867,050 373 102,489,144 9,625,356,194 374 103,314,199 9,728,670,393 375 104,143,674 9,832,814,067 376 104,977,576 9,937,791,643 377 105,815,912 10,043,607,555 378 106,658,705 10,150,266,260 379 107,505,954 10,257,772,214 380 108,357,683 10,366,129,897 381 109,213,906 10,475,343,803 382 110,074,614 10,585,418,417 383 110,939,841 10,696,358,258 384 111,809,588 10,808,167,846 385 112,683,864 10,920,851,710 386 113,562,690 11,034,414,400 387 114,446,071 11,148,860,471 388 115,334,026 11,264,194,497 389 116,226,561 11,380,421,058 390 117,123,693 11,497,544,751 391 118,025,424 11,615,570,175 392 118,931,765 11,734,501,940 393 119,842,745 11,854,344,685 394 120,758,360 11,975,103,045 395 121,678,631 12,096,781,676 396 122,603,562 12,219,385,238 397 123,533,164 12,342,918,402 398 124,467,455 12,467,385,857 399 125,406,460 12,592,792,317 400 126,350,163 12,719,142,480 401 127,298,591 12,846,441,071 402 128,251,745 12,974,692,816 403 129,209,646 13,103,902,462 404 130,172,320 13,234,074,782 405 131,139,753 13,365,214,535 406 132,111,971 13,497,326,506 407 133,088,973 13,630,415,479 408 134,070,789 13,764,486,268 409 135,057,421 13,899,543,689 410 136,048,875 14,035,592,564 411 137,045,182 14,172,637,746 412 138,046,338 14,310,684,084 413 139,052,346 14,449,736,430 414 140,063,238 14,589,799,668 415 141,079,018 14,730,878,686 416 142,099,700 14,872,978,386 417 143,125,297 15,016,103,683 418 144,155,805 15,160,259,488 419 145,191,242 15,305,450,730 420 146,231,641 15,451,682,371 421 147,276,992 15,598,959,363 422 148,327,320 15,747,286,683 423 149,382,623 15,896,669,306 424 150,442,914 16,047,112,220 425 151,508,210 16,198,620,430 426 152,578,540 16,351,198,970 427 153,653,894 16,504,852,864 428 154,734,281 16,659,587,145 429 155,819,722 16,815,406,867 430 156,910,218 16,972,317,085 431 158,005,799 17,130,322,884 432 159,106,484 17,289,429,368 433 160,212,251 17,449,641,619 434 161,323,131 17,610,964,750 435 162,439,131 17,773,403,881 436 163,560,269 17,936,964,150 437 164,686,558 18,101,650,708 438 165,818,012 18,267,468,720 439 166,954,623 18,434,423,343 440 168,096,417 18,602,519,760 441 169,243,407 18,771,763,167 442 170,395,599 18,942,158,766 443 171,553,023 19,113,711,789 444 172,715,671 19,286,427,460 445 173,883,550 19,460,311,010 446 175,056,690 19,635,367,700 447 176,235,089 19,811,602,789 448 177,418,768 19,989,021,557 449 178,607,751 20,167,629,308 450 179,802,013 20,347,431,321 451 181,001,587 20,528,432,908 452 182,206,494 20,710,639,402 453 183,416,728 20,894,056,130 454 184,632,321 21,078,688,451 455 185,853,274 21,264,541,725 456 187,079,580 21,451,621,305 457 188,311,276 21,639,932,581 458 189,548,370 21,829,480,951 459 190,790,870 22,020,271,821 460 192,038,791 22,212,310,612 461 193,292,138 22,405,602,750 462 194,550,919 22,600,153,669 463 195,815,161 22,795,968,830 464 197,084,876 22,993,053,706 465 198,360,060 23,191,413,766 466 199,640,733 23,391,054,499 467 200,926,905 23,591,981,404 468 202,218,584 23,794,199,988 469 203,515,800 23,997,715,788 470 204,818,544 24,202,534,332 471 206,126,832 24,408,661,164 472 207,440,685 24,616,101,849 473 208,760,119 24,824,861,968 474 210,085,119 25,034,947,087 475 211,415,716 25,246,362,803 476 212,751,926 25,459,114,729 477 214,093,747 25,673,208,476 478 215,441,203 25,888,649,679 479 216,794,302 26,105,443,981 480 218,153,048 26,323,597,029 481 219,517,464 26,543,114,493 482 220,887,563 26,764,002,056 483 222,263,344 26,986,265,400 484 223,644,828 27,209,910,228 485 225,032,028 27,434,942,256 486 226,424,942 27,661,367,198 487 227,823,599 27,889,190,797 488 229,228,006 28,118,418,803 489 230,638,169 28,349,056,972 490 232,054,102 28,581,111,074 491 233,475,830 28,814,586,904 492 234,903,336 29,049,490,240 493 236,336,654 29,285,826,894 494 237,775,798 29,523,602,692 495 239,220,764 29,762,823,456 496 240,671,580 30,003,495,036 497 242,128,250 30,245,623,286 498 243,590,771 30,489,214,057 499 245,059,183 30,734,273,240 500 246,533,493 30,980,806,733 501 248,013,699 31,228,820,432 502 249,499,809 31,478,320,241 503 250,991,840 31,729,312,081 504 252,489,811 31,981,801,892 505 253,993,740 32,235,795,632 506 255,503,630 32,491,299,262 507 257,019,487 32,748,318,749 508 258,541,318 33,006,860,067 509 260,069,150 33,266,929,217 510 261,602,998 33,528,532,215 511 263,142,858 33,791,675,073 512 264,688,761 34,056,363,834 513 266,240,686 34,322,604,520 514 267,798,667 34,590,403,187 515 269,362,724 34,859,765,911 516 270,932,863 35,130,698,774 517 272,509,090 35,403,207,864 518 274,091,421 35,677,299,285 519 275,679,853 35,952,979,138 520 277,274,402 36,230,253,540 521 278,875,111 36,509,128,651 522 280,481,965 36,789,610,616 523 282,094,975 37,071,705,591 524 283,714,160 37,355,419,751 525 285,339,516 37,640,759,267 526 286,971,071 37,927,730,338 527 288,608,851 38,216,339,189 528 290,252,842 38,506,592,031 529 291,903,058 38,798,495,089 530 293,559,520 39,092,054,609 531 295,222,219 39,387,276,828 532 296,891,211 39,684,168,039 533 298,566,487 39,982,734,526 534 300,248,040 40,282,982,566 535 301,935,890 40,584,918,456 536 303,630,050 40,888,548,506 537 305,330,541 41,193,879,047 538 307,037,382 41,500,916,429 539 308,750,566 41,809,666,995 540 310,470,102 42,120,137,097 541 312,196,005 42,432,333,102 542 313,928,307 42,746,261,409 543 315,667,000 43,061,928,409 544 317,412,108 43,379,340,517 545 319,163,637 43,698,504,154 546 320,921,580 44,019,425,734 547 322,685,975 44,342,111,709 548 324,456,823 44,666,568,532 549 326,234,148 44,992,802,680 550 328,017,955 45,320,820,635 551 329,808,237 45,650,628,872 552 331,605,020 45,982,233,892 553 333,408,329 46,315,642,221 554 335,218,158 46,650,860,379 555 337,034,532 46,987,894,911 556 338,857,458 47,326,752,369 557 340,686,929 47,667,439,298 558 342,522,977 48,009,962,275 559 344,365,622 48,354,327,897 560 346,214,855 48,700,542,752 561 348,070,697 49,048,613,449 562 349,933,166 49,398,546,615 563 351,802,250 49,750,348,865 564 353,677,982 50,104,026,847 565 355,560,389 50,459,587,236 566 357,449,446 50,817,036,682 567 359,345,186 51,176,381,868 568 361,247,622 51,537,629,490 569 363,156,743 51,900,786,233 570 365,072,599 52,265,858,832 571 366,995,178 52,632,854,010 572 368,924,482 53,001,778,492 573 370,860,550 53,372,639,042 574 372,803,377 53,745,442,419 575 374,752,970 54,120,195,389 576 376,709,351 54,496,904,740 577 378,672,538 54,875,577,278 578 380,642,517 55,256,219,795 579 382,619,329 55,638,839,124 580 384,602,972 56,023,442,096 581 386,593,454 56,410,035,550 582 388,590,789 56,798,626,339 583 390,595,010 57,189,221,349 584 392,606,094 57,581,827,443 585 394,624,074 57,976,451,517 586 396,648,961 58,373,100,478 587 398,680,746 58,771,781,224 588 400,719,473 59,172,500,697 589 402,765,138 59,575,265,835 590 404,817,751 59,980,083,586 591 406,877,325 60,386,960,911 592 408,943,877 60,795,904,788 593 411,017,405 61,206,922,193 594 413,097,937 61,620,020,130 595 415,185,486 62,035,205,616 596 417,280,043 62,452,485,659 597 419,381,640 62,871,867,299 598 421,490,275 63,293,357,574 599 423,605,964 63,716,963,538 600 425,728,730 64,142,692,268 601 427,858,584 64,570,550,852 602 429,995,514 65,000,546,366 603 432,139,549 65,432,685,915 604 434,290,700 65,866,976,615 605 436,448,976 66,303,425,591 606 438,614,403 66,742,039,994 607 440,786,981 67,182,826,975 608 442,966,711 67,625,793,686 609 445,153,610 68,070,947,296 610 447,347,703 68,518,294,999 611 449,549,001 68,967,844,000 612 451,757,501 69,419,601,501 613 453,973,228 69,873,574,729 614 456,196,175 70,329,770,904 615 458,426,365 70,788,197,269 616 460,663,826 71,248,861,095 617 462,908,554 71,711,769,649 618 465,160,561 72,176,930,210 619 467,419,860 72,644,350,070 620 469,686,448 73,114,036,518 621 471,960,360 73,585,996,878 622 474,241,618 74,060,238,496 623 476,530,204 74,536,768,700 624 478,826,132 75,015,594,832 625 481,129,428 75,496,724,260 626 483,440,089 75,980,164,349 627 485,758,145 76,465,922,494 628 488,083,618 76,954,006,112 629 490,416,477 77,444,422,589 630 492,756,760 77,937,179,349 631 495,104,481 78,432,283,830 632 497,459,641 78,929,743,471 633 499,822,266 79,429,565,737 634 502,192,373 79,931,758,110 635 504,569,937 80,436,328,047 636 506,954,994 80,943,283,041 637 509,347,557 81,452,630,598 638 511,747,640 81,964,378,238 639 514,155,264 82,478,533,502 640 516,570,420 82,995,103,922 641 518,993,113 83,514,097,035 642 521,423,376 84,035,520,411 643 523,861,218 84,559,381,629 644 526,306,648 85,085,688,277 645 528,759,684 85,614,447,961 646 531,220,323 86,145,668,284 647 533,688,567 86,679,356,851 648 536,164,458 87,215,521,309 649 538,648,000 87,754,169,309 650 541,139,199 88,295,308,508 651 543,638,069 88,838,946,577 652 546,144,610 89,385,091,187 653 548,658,840 89,933,750,027 654 551,180,786 90,484,930,813 655 553,710,454 91,038,641,267 656 556,247,840 91,594,889,107 657 558,792,969 92,153,682,076 658 561,345,844 92,715,027,920 659 563,906,480 93,278,934,400 660 566,474,914 93,845,409,314 661 569,051,115 94,414,460,429 662 571,635,117 94,986,095,546 663 574,226,937 95,560,322,483 664 576,826,567 96,137,149,050 665 579,434,030 96,716,583,080 666 582,049,366 97,298,632,446 667 584,672,536 97,883,304,982 668 587,303,581 98,470,608,563 669 589,942,514 99,060,551,077 670 592,589,322 99,653,140,399 671 595,244,049 100,248,384,448 672 597,906,695 100,846,291,143 673 600,577,260 101,446,868,403 674 603,255,766 102,050,124,169 675 605,942,236 102,656,066,405 676 608,636,654 103,264,703,059 677 611,339,057 103,876,042,116 678 614,049,454 104,490,091,570 679 616,767,836 105,106,859,406 680 619,494,237 105,726,353,643 681 622,228,669 106,348,582,312 682 624,971,111 106,973,553,423 683 627,721,619 107,601,275,042 684 630,480,188 108,231,755,230 685 633,246,807 108,865,002,037 686 636,021,523 109,501,023,560 687 638,804,337 110,139,827,897 688 641,595,238 110,781,423,135 689 644,394,268 111,425,817,403 690 647,201,430 112,073,018,833 691 650,016,725 112,723,035,558 692 652,840,181 113,375,875,739 693 655,671,798 114,031,547,537 694 658,511,580 114,690,059,117 695 661,359,560 115,351,418,677 696 664,215,752 116,015,634,429 697 667,080,136 116,682,714,565 698 669,952,749 117,352,667,314 699 672,833,595 118,025,500,909 700 675,722,681 118,701,223,590 701 678,620,041 119,379,843,631 702 681,525,676 120,061,369,307 703 684,439,577 120,745,808,884 704 687,361,771 121,433,170,655 705 690,292,278 122,123,462,933 706 693,231,099 122,816,694,032 707 696,178,256 123,512,872,288 708 699,133,759 124,212,006,047 709 702,097,596 124,914,103,643 710 705,069,795 125,619,173,438 711 708,050,393 126,327,223,831 712 711,039,368 127,038,263,199 713 714,036,747 127,752,299,946 714 717,042,543 128,469,342,489 715 720,056,742 129,189,399,231 716 723,079,380 129,912,478,611 717 726,110,482 130,638,589,093 718 729,150,037 131,367,739,130 719 732,198,064 132,099,937,194 720 735,254,569 132,835,191,763 721 738,319,559 133,573,511,322 722 741,393,060 134,314,904,382 723 744,475,094 135,059,379,476 724 747,565,659 135,806,945,135 725 750,664,744 136,557,609,879 726 753,772,388 137,311,382,267 727 756,888,595 138,068,270,862 728 760,013,385 138,828,284,247 729 763,146,783 139,591,431,030 730 766,288,759 140,357,719,789 731 769,439,338 141,127,159,127 732 772,598,558 141,899,757,685 733 775,766,400 142,675,524,085 734 778,942,907 143,454,466,992 735 782,128,077 144,236,595,069 736 785,321,894 145,021,916,963 737 788,524,396 145,810,441,359 738 791,735,604 146,602,176,963 739 794,955,517 147,397,132,480 740 798,184,157 148,195,316,637 741 801,421,525 148,996,738,162 742 804,667,618 149,801,405,780 743 807,922,471 150,609,328,251 744 811,186,096 151,420,514,347 745 814,458,498 152,234,972,845 746 817,739,685 153,052,712,530 747 821,029,676 153,873,742,206 748 824,328,464 154,698,070,670 749 827,636,092 155,525,706,762 750 830,952,560 156,356,659,322 751 834,277,871 157,190,937,193 752 837,612,047 158,028,549,240 753 840,955,082 158,869,504,322 754 844,306,994 159,713,811,316 755 847,667,815 160,561,479,131 756 851,037,549 161,412,516,680 757 854,416,192 162,266,932,872 758 857,803,768 163,124,736,640 759 861,200,279 163,985,936,919 760 864,605,743 164,850,542,662 761 868,020,193 165,718,562,855 762 871,443,606 166,590,006,461 763 874,876,004 167,464,882,465 764 878,317,410 168,343,199,875 765 881,767,824 169,224,967,699 766 885,227,258 170,110,194,957 767 888,695,753 170,998,890,710 768 892,173,277 171,891,063,987 769 895,659,858 172,786,723,845 770 899,155,529 173,685,879,374 771 902,660,259 174,588,539,633 772 906,174,094 175,494,713,727 773 909,697,050 176,404,410,777 774 913,229,102 177,317,639,879 775 916,770,298 178,234,410,177 776 920,320,644 179,154,730,821 777 923,880,129 180,078,610,950 778 927,448,783 181,006,059,733 779 931,026,635 181,937,086,368 780 934,613,655 182,871,700,023 781 938,209,881 183,809,909,904 782 941,815,324 184,751,725,228 783 945,429,974 185,697,155,202 784 949,053,881 186,646,209,083 785 952,687,045 187,598,896,128 786 956,329,448 188,555,225,576 787 959,981,129 189,515,206,705 788 963,642,105 190,478,848,810 789 967,312,357 191,446,161,167 790 970,991,935 192,417,153,102 791 974,680,833 193,391,833,935 792 978,379,049 194,370,212,984 793 982,086,613 195,352,299,597 794 985,803,536 196,338,103,133 795 989,529,810 197,327,632,943 796 993,265,478 198,320,898,421 797 997,010,545 199,317,908,966 798 1,000,764,994 200,318,673,960 799 1,004,528,858 201,323,202,818 800 1,008,302,151 202,331,504,969 801 1,012,084,878 203,343,589,847 802 1,015,877,069 204,359,466,916 803 1,019,678,727 205,379,145,643 804 1,023,489,826 206,402,635,469 805 1,027,310,417 207,429,945,886 806 1,031,140,520 208,461,086,406 807 1,034,980,118 209,496,066,524 808 1,038,829,250 210,534,895,774 809 1,042,687,909 211,577,583,683 810 1,046,556,085 212,624,139,768 811 1,050,433,829 213,674,573,597 812 1,054,321,168 214,728,894,765 813 1,058,218,065 215,787,112,830 814 1,062,124,559 216,849,237,389 815 1,066,040,649 217,915,278,038 816 1,069,966,339 218,985,244,377 817 1,073,901,670 220,059,146,047 818 1,077,846,651 221,136,992,698 819 1,081,801,268 222,218,793,966 820 1,085,765,541 223,304,559,507 821 1,089,739,494 224,394,299,001 822 1,093,723,121 225,488,022,122 823 1,097,716,454 226,585,738,576 824 1,101,719,513 227,687,458,089 825 1,105,732,275 228,793,190,364 826 1,109,754,761 229,902,945,125 827 1,113,786,994 231,016,732,119 828 1,117,828,979 232,134,561,098 829 1,121,880,744 233,256,441,842 830 1,125,942,303 234,382,384,145 831 1,130,013,622 235,512,397,767 832 1,134,094,745 236,646,492,512 833 1,138,185,692 237,784,678,204 834 1,142,286,467 238,926,964,671 835 1,146,397,089 240,073,361,760 836 1,150,517,555 241,223,879,315 837 1,154,647,866 242,378,527,181 838 1,158,788,047 243,537,315,228 839 1,162,938,143 244,700,253,371 840 1,167,098,124 245,867,351,495 841 1,171,268,015 247,038,619,510 842 1,175,447,830 248,214,067,340 843 1,179,637,549 249,393,704,889 844 1,183,837,229 250,577,542,118 845 1,188,046,885 251,765,589,003 846 1,192,266,486 252,957,855,489 847 1,196,496,074 254,154,351,563 848 1,200,735,653 255,355,087,216 849 1,204,985,220 256,560,072,436 850 1,209,244,821 257,769,317,257 851 1,213,514,457 258,982,831,714 852 1,217,794,109 260,200,625,823 853 1,222,083,827 261,422,709,650 854 1,226,383,592 262,649,093,242 855 1,230,693,429 263,879,786,671 856 1,235,013,363 265,114,800,034 857 1,239,343,402 266,354,143,436 858 1,243,683,533 267,597,826,969 859 1,248,033,793 268,845,860,762 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 196.092327 seconds
Crystal
Brute force and slow
Checks if each number upto limit is humble number.
def humble?(i)
return true if (i < 2)
return humble?(i // 2) if (i % 2 == 0)
return humble?(i // 3) if (i % 3 == 0)
return humble?(i // 5) if (i % 5 == 0)
return humble?(i // 7) if (i % 7 == 0)
false
end
count, num = 0, 0_i64
digits = 10 # max digits for humble numbers
limit = 10_i64 ** digits # max numbers to search through
humble = Array.new(digits + 1, 0)
while (num += 1) < limit
if humble?(num)
humble[num.to_s.size] += 1
print num, " " if count < 50
count += 1
end
end
print "\n\nOf the first #{count} humble numbers:\n"
(1..digits).each { |num| printf("%5d have %2d digits\n", humble[num], num) }
- Output:
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
Of the first 7574 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
2381 have 10 digits
Direct Generation: Orders of magnitude faster
Generate humble numbers directly.
require "big"
def humble(digits)
h = [1.to_big_i]
x2, x3, x5, x7 = 2.to_big_i, 3.to_big_i, 5.to_big_i, 7.to_big_i
i, j, k, l = 0, 0, 0, 0
(1..).each do |n|
x = {x2, x3, x5, x7}.min # {} tuple|stack faster [] array|heap
break if x.to_s.size > digits
h << x
x2 = 2 * h[i += 1] if x2 == h[n]
x3 = 3 * h[j += 1] if x3 == h[n]
x5 = 5 * h[k += 1] if x5 == h[n]
x7 = 7 * h[l += 1] if x7 == h[n]
end
h
end
digits = 50 # max digits for humble numbers
h = humble(digits) # humble numbers <= digits size
count = h.size # the total humble numbers count
counts = h.map { |n| n.to_s.size }.tally # hash of digits counts 1..digits
print "First 50 Humble Numbers: \n"; (0...50).each { |i| print "#{h[i]} " }
print "\n\nOf the first #{count} humble numbers:\n"
(1..digits).each { |num| printf("%6d have %2d digits\n", counts[num], num) }
- Output:
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
Of the first 3363713 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
2381 have 10 digits
3113 have 11 digits
3984 have 12 digits
5002 have 13 digits
6187 have 14 digits
7545 have 15 digits
9081 have 16 digits
10815 have 17 digits
12759 have 18 digits
14927 have 19 digits
17323 have 20 digits
19960 have 21 digits
22853 have 22 digits
26015 have 23 digits
29458 have 24 digits
33188 have 25 digits
37222 have 26 digits
41568 have 27 digits
46245 have 28 digits
51254 have 29 digits
56618 have 30 digits
62338 have 31 digits
68437 have 32 digits
74917 have 33 digits
81793 have 34 digits
89083 have 35 digits
96786 have 36 digits
104926 have 37 digits
113511 have 38 digits
122546 have 39 digits
132054 have 40 digits
142038 have 41 digits
152515 have 42 digits
163497 have 43 digits
174986 have 44 digits
187004 have 45 digits
199565 have 46 digits
212675 have 47 digits
226346 have 48 digits
240590 have 49 digits
255415 have 50 digits
D
import std.conv;
import std.stdio;
bool isHumble(int i) {
if (i <= 1) return true;
if (i % 2 == 0) return isHumble(i / 2);
if (i % 3 == 0) return isHumble(i / 3);
if (i % 5 == 0) return isHumble(i / 5);
if (i % 7 == 0) return isHumble(i / 7);
return false;
}
void main() {
auto limit = short.max;
int[int] humble;
auto count = 0;
auto num = 1;
while (count < limit) {
if (isHumble(num)) {
auto str = num.to!string;
auto len = str.length;
humble[len]++;
if (count < 50) write(num, ' ');
count++;
}
num++;
}
writeln('\n');
writeln("Of the first ", count, " humble numbers:");
num = 1;
while (num < humble.length - 1) {
if (num in humble) {
auto c = humble[num];
writefln("%5d have %2d digits", c, num);
num++;
} else {
break;
}
}
}
- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Delphi
See Pascal.
Elm
As discussed in the Elm Hamming numbers contribution and further limited here as to the fastest contribution in the Haskell Direct Generation contribution to this task due to not having an efficient random access read/write data structure such as a linear array, the implementation is limited to using a Minimum Heap binary queue as in the N-smooth Elm contribution, which is reasonably fast using the logarithmic estimates for ordering of the sequence. The following code implements the task (the BigInt package has been installed from "cmditch/elm-bigint" as well as "elm/time"):
module Main exposing (main)
import Browser
import Html exposing (div, pre, text, br)
import Task exposing (Task, succeed, andThen, perform)
import BigInt
import Bitwise exposing (shiftRightBy, and)
import Time exposing (now, posixToMillis)
-- an infinite non-empty non-memoizing Co-Inductive Stream (CIS)...
type CIS a = CIS a (() -> CIS a)
takeCIS2String : Int -> (a -> String) -> CIS a -> String
takeCIS2String n cvtf cis =
let loop i (CIS hd tl) lst =
if i < 1 then List.reverse lst |> String.join ", "
else loop (i - 1) (tl()) (cvtf hd :: lst)
in loop n cis []
-- a Min Heap binary heap Priority Queue...
type PriorityQ comparable v =
Mt
| Br comparable v (PriorityQ comparable v)
(PriorityQ comparable v)
emptyPQ : PriorityQ comparable v
emptyPQ = Mt
peekMinPQ : PriorityQ comparable v -> Maybe (comparable, v)
peekMinPQ pq = case pq of
(Br k v _ _) -> Just (k, v)
Mt -> Nothing
pushPQ : comparable -> v -> PriorityQ comparable v
-> PriorityQ comparable v
pushPQ wk wv pq =
case pq of
Mt -> Br wk wv Mt Mt
(Br vk vv pl pr) ->
if wk <= vk then Br wk wv (pushPQ vk vv pr) pl
else Br vk vv (pushPQ wk wv pr) pl
siftdown : comparable -> v -> PriorityQ comparable v
-> PriorityQ comparable v -> PriorityQ comparable v
siftdown wk wv pql pqr =
case pql of
Mt -> Br wk wv Mt Mt
(Br vkl vvl pll prl) ->
case pqr of
Mt -> if wk <= vkl then Br wk wv pql Mt
else Br vkl vvl (Br wk wv Mt Mt) Mt
(Br vkr vvr plr prr) ->
if wk <= vkl && wk <= vkr then Br wk wv pql pqr
else if vkl <= vkr then Br vkl vvl (siftdown wk wv pll prl) pqr
else Br vkr vvr pql (siftdown wk wv plr prr)
replaceMinPQ : comparable -> v -> PriorityQ comparable v
-> PriorityQ comparable v
replaceMinPQ wk wv pq = case pq of
Mt -> Mt
(Br _ _ pl pr) -> siftdown wk wv pl pr
-- actual humble function implementation...
type alias Mults = { x2 : Int, x3 : Int, x5 : Int, x7 : Int }
type alias LogRep = { lg : Float, mlts : Mults }
oneLogRep : LogRep
oneLogRep = LogRep 0.0 <| Mults 0 0 0 0
lg10 : Float
lg10 = 1.0
lg7 : Float
lg7 = logBase 10 7
lg5 : Float
lg5 = logBase 10.0 5.0
lg3 : Float
lg3 = logBase 10.0 3.0
lg2 : Float
lg2 = lg10 - lg5
multLR2 : LogRep -> LogRep
multLR2 ({ lg, mlts } as lr) =
{ lr | lg = lg + lg2, mlts = { mlts | x2 = mlts.x2 + 1 } }
multLR3 : LogRep -> LogRep
multLR3 ({ lg, mlts } as lr) =
{ lr | lg = lg + lg3, mlts = { mlts | x3 = mlts.x3 + 1 } }
multLR5 : LogRep -> LogRep
multLR5 ({ lg, mlts } as lr) =
{ lr | lg = lg + lg5, mlts = { mlts | x5 = mlts.x5 + 1 } }
multLR7 : LogRep -> LogRep
multLR7 ({ lg, mlts } as lr) =
{ lr | lg = lg + lg7, mlts = { mlts | x7 = mlts.x7 + 1 } }
showLogRep : LogRep -> String
showLogRep lr =
let xpnd x m r =
if x <= 0 then r
else xpnd (shiftRightBy 1 x) (BigInt.mul m m)
(if (and 1 x) /= 0 then BigInt.mul m r else r)
in BigInt.fromInt 1 |> xpnd lr.mlts.x2 (BigInt.fromInt 2)
|> xpnd lr.mlts.x3 (BigInt.fromInt 3) |> xpnd lr.mlts.x5 (BigInt.fromInt 5)
|> xpnd lr.mlts.x7 (BigInt.fromInt 7) |> BigInt.toString
humblesLog : () -> CIS LogRep
humblesLog() =
let prmfs = [multLR7, multLR5, multLR3, multLR2]
fprmf = List.head prmfs |> Maybe.withDefault identity -- never Nothing!
rstps = List.tail prmfs |> Maybe.withDefault [] -- never Nothing!
frstcis =
let nxt lr =
CIS lr <| \ _ -> nxt (fprmf lr)
in nxt (fprmf oneLogRep)
dflt = (0.0, Mults 0 0 0 0)
mkcis lrf cis =
let frst = lrf oneLogRep
scnd = lrf frst
nxt pq (CIS hd tlf as cs) =
let (lgv, v) = peekMinPQ pq |> Maybe.withDefault dflt in
if lgv < hd.lg then let lr = (LogRep lgv v) in CIS lr <| \ _ ->
let { lg, mlts } = lrf lr
in nxt (replaceMinPQ lg mlts pq) cs
else CIS hd <| \ _ ->
let { lg, mlts } = lrf hd
in nxt (pushPQ lg mlts pq) (tlf())
in CIS frst <| \ _ -> nxt (pushPQ scnd.lg scnd.mlts emptyPQ) cis
rest() = List.foldl mkcis frstcis rstps
in CIS oneLogRep <| \ _ -> rest()
-- pretty printing function to add commas every 3 chars from left...
comma3 : String -> String
comma3 s =
let go n lst =
if n < 1 then String.join "," lst else
let nn = max (n - 3) 0
in go nn (String.slice nn n s :: lst)
in go (String.length s) []
humbleDigitCountsTo : Int -> CIS LogRep -> List String
humbleDigitCountsTo n cis =
let go i (CIS hd tlf) cnt cacc lst =
if i >= n then List.reverse lst else
if truncate hd.lg <= i then go i (tlf()) (cnt + 1) cacc lst
else let ni = i + 1
ncacc = cacc + cnt
str =
(String.padLeft 4 ' ' << String.fromInt) ni
++ (String.padLeft 14 ' ' << comma3 << String.fromInt) cnt
++ (String.padLeft 19 ' ' << comma3 << String.fromInt) ncacc
in go ni (tlf()) 1 ncacc (str :: lst) -- always > 1 per dgt
in go 0 cis 0 0 []
-- code to do with testing...
timemillis : () -> Task Never Int -- a side effect function
timemillis() = now |> andThen (\ t -> succeed (posixToMillis t))
test : () -> Cmd Msg
test() =
let numdgt = 100
hdg1 = "The first 50 humble numbers are: "
msg1 = humblesLog() |> takeCIS2String 50 showLogRep
hdg2 = "Count of humble numbers for each digit length 1-"
++ String.fromInt numdgt ++ ":"
msg2 = "Digits Count Accum"
in timemillis()
|> Task.andThen (\ strt ->
let rslt = humblesLog() |> humbleDigitCountsTo numdgt
in timemillis()
|> Task.andThen (\ stop ->
succeed (((hdg1 ++ msg1) :: "" :: hdg2 :: msg2 :: rslt)
++ ["Counting took " ++ String.fromInt (stop - strt)
++ " milliseconds."])))
|> perform Done
-- following code has to do with outputting to a web page using MUV/TEA...
type alias Model = List String
type Msg = Done Model
main : Program () Model Msg
main = Browser.element
{ init = \ _ -> ([], test())
, update = \ (Done mdl) _ -> (mdl, Cmd.none)
, subscriptions = \ _ -> Sub.none
, view = div [] << List.map (\ s ->
if s == "" then br [] []
else pre [] <| List.singleton <| text s)
}
- Output:
The first 50 humble numbers are: 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-100: 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 . . abbreviated with full results checked as per the C++ contribution to this task. . 90 1,463,862 33,914,307 91 1,512,840 35,427,147 92 1,562,897 36,990,044 93 1,614,050 38,604,094 94 1,666,302 40,270,396 95 1,719,669 41,990,065 96 1,774,166 43,764,231 97 1,829,805 45,594,036 98 1,886,590 47,480,626 99 1,944,540 49,425,166 100 2,003,661 51,428,827
This execution time is as run on an Intel i5-6500 at 3.6 GHz boosted for single-threading is perhaps thirty percent slower than if it would be with a loop value holding the minimum "head" value of the queue to avoid the overhead of "peek" operations on the queue when the test doesn't require other access to the queue as used in the Elm contribution to the Hamming numbers task page, but the advantage here is that the code is shorter and easier to read and understand. This code is slower than other languages, even for JavaScript, due to the `O(n log n)` asymptotic execution performance of the persistent priority queue implementation so execution time increases at worse than linear rates with the number of elements processed in the sequence.
F#
The Functions
// Generate humble numbers. Nigel Galloway: June 18th., 2020
let fN g=let mutable n=1UL in (fun()->n<-n*g;n)
let fI (n:uint64) g=let mutable q=n in (fun()->let t=q in q<-n*g();t)
let fG n g=let mutable vn,vg=n(),g() in fun()->match vg<vn with true->let t=vg in vg<-g();t |_->let t=vn in vn<-n();t
let rec fE n=seq{yield n();yield! fE n}
let fL n g=let mutable vn,vg,v=n(),g(),None
fun()->match v with
Some n->v<-None;n
|_->match vg() with
r when r<vn->r
|r->vg<-fG vg (fI vn (g()));vn<-n();v<-Some r;vg()
let humble = seq{yield 1UL;yield! fE(fL (fN 7UL) (fun()->fL (fN 5UL) (fun()->fL (fN 3UL) (fun()->fN 2UL))))}
The Tasks
humble |> Seq.take 50 |> Seq.iter (printf "%d ");printfn ""
- Output:
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
for n in [1..18] do let g=pown 10UL n in printfn "There are %d humble numbers with %d digits" (humble|>Seq.skipWhile(fun n->n<g/10UL)|>Seq.takeWhile(fun n->n<g)|>Seq.length) n
- Output:
There are 9 humble numbers with 1 digits There are 36 humble numbers with 2 digits There are 95 humble numbers with 3 digits There are 197 humble numbers with 4 digits There are 356 humble numbers with 5 digits There are 579 humble numbers with 6 digits There are 882 humble numbers with 7 digits There are 1272 humble numbers with 8 digits There are 1767 humble numbers with 9 digits There are 2381 humble numbers with 10 digits There are 3113 humble numbers with 11 digits There are 3984 humble numbers with 12 digits There are 5002 humble numbers with 13 digits There are 6187 humble numbers with 14 digits There are 7545 humble numbers with 15 digits There are 9081 humble numbers with 16 digits There are 10815 humble numbers with 17 digits There are 12759 humble numbers with 18 digits
Faster by Using Less Seq's and Using Logarithmic Approximations for Ordering
The above code is somewhat of a "toy" implementation in that it is very obscure and difficult to read and understand, even for someone used to F# and functional programming using closures (come now, spell your name with the names of the functions???). The above code is also very slow and therefore limited, even if the minor changes so that it outputs BigInt's is done; this is due to overuse of the very slow Seq's for iteration and the deeply nested closure functions which don't implement memoization other than for the equivalent heads of the "lazy lists" and therefore repeat many operations and thus don't have a linear response with number of elements (which also would not be linear due to the increasing amount of work in doing BigInt computations). The following code doesn't use Seq' but rather a "roll-your-own" Co-Inductive Stream (CIS) and eliminates the need for memoization by keeping track of the required back results in DotNet Queue's; it also uses a logarithmic representation for ordering comparisons to eliminate the BigInt operations:
// a count and logarithmic approximation of the humble value...
type LogRep = struct val lg: uint64; val x2: uint16; val x3: uint16;
val x5: uint16; val x7: uint16
new(lg, x2, x3, x5, x7) =
{lg = lg; x2 = x2; x3 = x3; x5 = x5; x7 = x7 } end
let one: LogRep = LogRep(0UL, 0us, 0us, 0us, 0us)
let logshft = 50
let fac = pown 2.0 logshft
let lg10_10 = 1UL <<< logshft
let lg7_10 = (uint64 << round) <| log 7.0 / log 10.0 * fac
let lg5_10 = (uint64 << round) <| log 5.0 / log 10.0 * fac
let lg3_10 = (uint64 << round) <| log 3.0 / log 10.0 * fac
let lg2_10 = lg10_10 - lg5_10
let inline mul2 (lr: LogRep): LogRep =
LogRep(lr.lg + lg2_10, lr.x2 + 1us, lr.x3, lr.x5, lr.x7)
let inline mul3 (lr: LogRep): LogRep =
LogRep(lr.lg + lg3_10, lr.x2, lr.x3 + 1us, lr.x5, lr.x7)
let inline mul5 (lr: LogRep): LogRep =
LogRep(lr.lg + lg5_10, lr.x2, lr.x3, lr.x5 + 1us, lr.x7)
let inline mul7 (lr: LogRep): LogRep =
LogRep(lr.lg + lg7_10, lr.x2, lr.x3, lr.x5, lr.x7 + 1us)
let lr2BigInt (lr: LogRep) =
let rec xpnd n mlt rslt =
if n <= 0us then rslt
else xpnd (n - 1us) mlt (mlt * rslt)
xpnd lr.x2 2I 1I |> xpnd lr.x3 3I |> xpnd lr.x5 5I |> xpnd lr.x7 7I
type CIS<'a> = CIS of 'a * (Unit -> CIS<'a>) // infinite Co-Inductive Stream...
let cis2Seq cis =
Seq.unfold (fun (CIS(hd, tlf)) -> Some(hd, tlf())) cis
let humblesLog() =
let prmfs = [ mul7; mul5; mul3; mul2 ]
let frstpf = Seq.head prmfs
let rstpfs = Seq.tail prmfs
let frstll =
let rec nxt n = CIS(n, fun () -> nxt (frstpf n))
nxt (frstpf one)
let mkcis cis mf =
let q = Queue<LogRep>(1024)
let fv = mf one
let nv = mf fv
let rec nxt (hdv: LogRep) (CIS(chd: LogRep, ctlf) as cis) =
if hdv.lg < chd.lg then
CIS(hdv, fun () -> q.Enqueue (mf hdv); nxt (q.Dequeue()) cis)
else CIS(chd, fun () -> q.Enqueue (mf chd); nxt hdv (ctlf()))
CIS(fv, fun () -> nxt nv cis)
CIS(one, fun () -> (Seq.fold mkcis frstll rstpfs))
let comma3 v =
let s = string v
let rec loop n lst =
if n < 1 then List.fold (fun s xs ->
s + "," + xs) (List.head lst) <| List.tail lst
else let nn = max (n - 3) 0 in loop nn (s.[nn .. n - 1] :: lst)
loop (String.length s) []
let digitCountTo n ll =
let rec loop i (CIS(hd: LogRep, tlf)) cnt cacc =
if int i <= n then
if hd.lg >>> logshft < i then loop i (tlf()) (cnt + 1) cacc else
let ncacc = cacc + cnt
printfn "%4d%14s%19s" i (comma3 cnt) (comma3 ncacc)
loop (i + 1UL) (tlf()) 1 ncacc
loop 1UL ll 0 0
printfn "The first 50 humble numbers are:"
humblesLog() |> cis2Seq |> Seq.take 50 |> Seq.map lr2BigInt
|> Seq.iter (printf "%A ");printfn ""
printfn ""
let numDigits = 255
printfn "Count of humble numbers for each digit length 1-%d:" numDigits
printfn "Digits Count Accum"
let strt = System.DateTime.Now.Ticks
humblesLog() |> digitCountTo numDigits
let stop = System.DateTime.Now.Ticks
printfn "Counting took %d milliseconds" <| ((stop - strt) / 10000L)
- Output:
The first 50 humble numbers are: 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 85945 milliseconds
This, as run on an Intel i5-6500 at 3.6 GHz when running single-threaded, is about twice as fast as the Haskell version of the same due to the time lost by Haskell in lazy list operations and about twice as slow as the fastest Pascal version due to the Pascal version being completely optimized for the task of counting digits in order, which seems of little point given that ordering before counting digits as in the following code is so much easier and faster.
This takes over three hours to count the digits up to 877.
Even Faster by Using Logarithms but Skipping Ordering Entirely
As per the C++ Direct Generation contribution, there is no need to count the occurrences per digit length in order which saves a lot of code and execution time; as well, there is a slight optimization to do array access via pointer to save about twenty percent of the time used for array bounds checks as implemented in the following code:
open System.Collections.Generic
open Microsoft.FSharp.NativeInterop
// a count and logarithmic approximation of the humble value...
type LogRep = struct val lg: uint64; val x2: uint16; val x3: uint16;
val x5: uint16; val x7: uint16
new(lg, x2, x3, x5, x7) =
{lg = lg; x2 = x2; x3 = x3; x5 = x5; x7 = x7 } end
let one: LogRep = LogRep(0UL, 0us, 0us, 0us, 0us)
let logshft = 50
let fac = pown 2.0 logshft
let lg10_10 = 1UL <<< logshft
let lg7_10 = (uint64 << round) <| log 7.0 / log 10.0 * fac
let lg5_10 = (uint64 << round) <| log 5.0 / log 10.0 * fac
let lg3_10 = (uint64 << round) <| log 3.0 / log 10.0 * fac
let lg2_10 = lg10_10 - lg5_10
let inline mul2 (lr: LogRep): LogRep =
LogRep(lr.lg + lg2_10, lr.x2 + 1us, lr.x3, lr.x5, lr.x7)
let inline mul3 (lr: LogRep): LogRep =
LogRep(lr.lg + lg3_10, lr.x2, lr.x3 + 1us, lr.x5, lr.x7)
let inline mul5 (lr: LogRep): LogRep =
LogRep(lr.lg + lg5_10, lr.x2, lr.x3, lr.x5 + 1us, lr.x7)
let inline mul7 (lr: LogRep): LogRep =
LogRep(lr.lg + lg7_10, lr.x2, lr.x3, lr.x5, lr.x7 + 1us)
let lr2BigInt (lr: LogRep) =
let rec xpnd n mlt rslt =
if n <= 0us then rslt
else xpnd (n - 1us) mlt (mlt * rslt)
xpnd lr.x2 2I 1I |> xpnd lr.x3 3I |> xpnd lr.x5 5I |> xpnd lr.x7 7I
type CIS<'a> = CIS of 'a * (Unit -> CIS<'a>) // infinite Co-Inductive Stream...
let cis2Seq cis =
Seq.unfold (fun (CIS(hd, tlf)) -> Some(hd, tlf())) cis
let humblesLog() =
let prmfs = [ mul7; mul5; mul3; mul2 ]
let frstpf = Seq.head prmfs
let rstpfs = Seq.tail prmfs
let frstll =
let rec nxt n = CIS(n, fun () -> nxt (frstpf n))
nxt (frstpf one)
let mkcis cis mf =
let q = Queue<LogRep>(1024)
let fv = mf one
let nv = mf fv
let rec nxt (hdv: LogRep) (CIS(chd: LogRep, ctlf) as cis) =
if hdv.lg < chd.lg then
CIS(hdv, fun () -> q.Enqueue (mf hdv); nxt (q.Dequeue()) cis)
else CIS(chd, fun () -> q.Enqueue (mf chd); nxt hdv (ctlf()))
CIS(fv, fun () -> nxt nv cis)
CIS(one, fun () -> (Seq.fold mkcis frstll rstpfs))
let comma3 v =
let s = string v
let rec loop n lst =
if n < 1 then List.fold (fun s xs ->
s + "," + xs) (List.head lst) <| List.tail lst
else let nn = max (n - 3) 0 in loop nn (s.[nn .. n - 1] :: lst)
loop (String.length s) []
printfn "The first 50 humble numbers are:"
humblesLog() |> cis2Seq |> Seq.take 50 |> Seq.map lr2BigInt
|> Seq.iter (printf "%A ");printfn ""
printfn ""
let numDigits = 877
printfn "Count of humble numbers for each digit length 1-%d:" numDigits
printfn "Digits Count Accum"
let strt = System.DateTime.Now.Ticks
let bins = Array.zeroCreate numDigits
#nowarn "9" // no warnings for the use of native pointers...
#nowarn "51"
let lmt = uint64 numDigits <<< logshft
let rec loopw w =
if w < lmt then
let rec loopx x =
if x < lmt then
let rec loopy y =
if y < lmt then
let rec loopz z =
if z < lmt then
// let ndx = z >>> logshft |> int
// bins.[ndx] <- bins.[ndx] + 1UL
// use pointers to save array bounds checking...
let ndx = &&bins.[z >>> logshft |> int]
NativePtr.write ndx (NativePtr.read ndx + 1UL)
loopz (z + lg7_10) in loopz y; loopy (y + lg5_10)
loopy x; loopx (x + lg3_10) in loopx w; loopw (w + lg2_10) in loopw 0UL
bins |> Seq.scan (fun (i, _, a) v ->
i + 1, v, a + v) (0, 0UL, 0UL) |> Seq.skip 1
|> Seq.iter (fun (i, c, a) -> printfn "%4d%14s%19s" i (comma3 c) (comma3 a))
let stop = System.DateTime.Now.Ticks
printfn "Counting took %d milliseconds" <| ((stop - strt) / 10000L)
- Output:
The first 50 humble numbers are: 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 316149 milliseconds
This is about twice as slow as the C++ or Haskell versions of the same algorithm due to being run on the DotNet JIT compiled environment.
Factor
USING: accessors assocs combinators deques dlists formatting fry
generalizations io kernel make math math.functions math.order
prettyprint sequences tools.memory.private ;
IN: rosetta-code.humble-numbers
TUPLE: humble-iterator 2s 3s 5s 7s digits
{ #digits initial: 1 } { target initial: 10 } ;
: <humble-iterator> ( -- humble-iterator )
humble-iterator new
1 1dlist >>2s
1 1dlist >>3s
1 1dlist >>5s
1 1dlist >>7s
H{ } clone >>digits ;
: enqueue ( n humble-iterator -- )
{
[ [ 2 * ] [ 2s>> ] ]
[ [ 3 * ] [ 3s>> ] ]
[ [ 5 * ] [ 5s>> ] ]
[ [ 7 * ] [ 7s>> ] ]
} [ bi* push-back ] map-compose 2cleave ;
: count-digits ( humble-iterator n -- )
[ over target>> >=
[ [ 1 + ] change-#digits [ 10 * ] change-target ] when ]
[ drop 1 swap [ #digits>> ] [ digits>> ] bi at+ ] bi ;
: ?pop ( 2s 3s 5s 7s n -- )
'[ dup peek-front _ = [ pop-front* ] [ drop ] if ] 4 napply ;
: next ( humble-iterator -- n )
dup dup { [ 2s>> ] [ 3s>> ] [ 5s>> ] [ 7s>> ] } cleave
4 ndup [ peek-front ] 4 napply min min min
{ [ ?pop ] [ swap enqueue ] [ count-digits ] [ ] } cleave ;
: upto-n-digits ( humble-iterator n -- seq )
1 + swap [ [ 2dup digits>> key? ] [ dup next , ] until ] { }
make [ digits>> delete-at ] dip but-last-slice ;
: .first50 ( seq -- )
"First 50 humble numbers:" print 50 head [ pprint bl ] each
nl ;
: .digit-breakdown ( humble-iterator -- )
"The digit counts of humble numbers:" print digits>> [
commas swap dup 1 = "" "s" ? "%9s have %2d digit%s\n"
printf
] assoc-each ;
: humble-numbers ( -- )
[ <humble-iterator> dup 95 upto-n-digits
[ .first50 nl ] [ drop .digit-breakdown nl ] [
"Total number of humble numbers found: " write length
commas print
] tri ] time ;
MAIN: humble-numbers
- Output:
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
The digit counts of humble numbers:
9 have 1 digit
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1,272 have 8 digits
1,767 have 9 digits
2,381 have 10 digits
3,113 have 11 digits
3,984 have 12 digits
5,002 have 13 digits
6,187 have 14 digits
7,545 have 15 digits
9,081 have 16 digits
10,815 have 17 digits
12,759 have 18 digits
14,927 have 19 digits
17,323 have 20 digits
19,960 have 21 digits
22,853 have 22 digits
26,015 have 23 digits
29,458 have 24 digits
33,188 have 25 digits
37,222 have 26 digits
41,568 have 27 digits
46,245 have 28 digits
51,254 have 29 digits
56,618 have 30 digits
62,338 have 31 digits
68,437 have 32 digits
74,917 have 33 digits
81,793 have 34 digits
89,083 have 35 digits
96,786 have 36 digits
104,926 have 37 digits
113,511 have 38 digits
122,546 have 39 digits
132,054 have 40 digits
142,038 have 41 digits
152,515 have 42 digits
163,497 have 43 digits
174,986 have 44 digits
187,004 have 45 digits
199,565 have 46 digits
212,675 have 47 digits
226,346 have 48 digits
240,590 have 49 digits
255,415 have 50 digits
270,843 have 51 digits
286,880 have 52 digits
303,533 have 53 digits
320,821 have 54 digits
338,750 have 55 digits
357,343 have 56 digits
376,599 have 57 digits
396,533 have 58 digits
417,160 have 59 digits
438,492 have 60 digits
460,533 have 61 digits
483,307 have 62 digits
506,820 have 63 digits
531,076 have 64 digits
556,104 have 65 digits
581,902 have 66 digits
608,483 have 67 digits
635,864 have 68 digits
664,053 have 69 digits
693,065 have 70 digits
722,911 have 71 digits
753,593 have 72 digits
785,141 have 73 digits
817,554 have 74 digits
850,847 have 75 digits
885,037 have 76 digits
920,120 have 77 digits
956,120 have 78 digits
993,058 have 79 digits
1,030,928 have 80 digits
1,069,748 have 81 digits
1,109,528 have 82 digits
1,150,287 have 83 digits
1,192,035 have 84 digits
1,234,774 have 85 digits
1,278,527 have 86 digits
1,323,301 have 87 digits
1,369,106 have 88 digits
1,415,956 have 89 digits
1,463,862 have 90 digits
1,512,840 have 91 digits
1,562,897 have 92 digits
1,614,050 have 93 digits
1,666,302 have 94 digits
1,719,669 have 95 digits
Total number of humble numbers found: 41,990,065
Running time: 335.1803624581294 seconds
FreeBASIC
Function IsHumble(i As Integer) As Boolean
If i <= 1 Then Return True
If i Mod 2 = 0 Then Return IsHumble(i \ 2)
If i Mod 3 = 0 Then Return IsHumble(i \ 3)
If i Mod 5 = 0 Then Return IsHumble(i \ 5)
If i Mod 7 = 0 Then Return IsHumble(i \ 7)
Return False
End Function
Const limiteMax = 7574
Dim As Integer humble(10) 'As New Dictionary(Of Integer, Integer)
Dim As Integer cont = 0, num = 1
Print "Los 50 primeros n£meros de Humble son:";
While cont < limiteMax
If IsHumble(num) Then
Dim As Integer longitud = Len(Str(num))
If longitud > 10 Then
Exit While
End If
If humble(longitud) Then
humble(longitud) += 1
Else
humble(longitud) = 1
End If
If cont < 50 Then
If cont Mod 10 = 0 Then Print
Print Using " ###"; num;
End If
cont += 1
End If
num += 1
Wend
Print !"\n\nDe los primeros"; cont; " n£meros de Humble:"
num = 1
While num < cont
If humble(num) Then
Print Using " #### tienen ## digitos"; humble(num); num
num += 1
Else
Exit While
End If
Wend- Output:
Los 50 primeros números de Humble son:
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
De los primeros 7574 números de Humble:
9 tienen 1 digitos
36 tienen 2 digitos
95 tienen 3 digitos
197 tienen 4 digitos
356 tienen 5 digitos
579 tienen 6 digitos
882 tienen 7 digitos
1272 tienen 8 digitos
1767 tienen 9 digitos
2381 tienen 10 digitos
Go
Not particularly fast and uses a lot of memory but easier to understand than the 'log' based methods for generating 7-smooth numbers.
package main
import (
"fmt"
"math/big"
)
var (
one = new(big.Int).SetUint64(1)
two = new(big.Int).SetUint64(2)
three = new(big.Int).SetUint64(3)
five = new(big.Int).SetUint64(5)
seven = new(big.Int).SetUint64(7)
ten = new(big.Int).SetUint64(10)
)
func min(a, b *big.Int) *big.Int {
if a.Cmp(b) < 0 {
return a
}
return b
}
func humble(n int) []*big.Int {
h := make([]*big.Int, n)
h[0] = new(big.Int).Set(one)
next2, next3 := new(big.Int).Set(two), new(big.Int).Set(three)
next5, next7 := new(big.Int).Set(five), new(big.Int).Set(seven)
var i, j, k, l int
for m := 1; m < len(h); m++ {
h[m] = new(big.Int).Set(min(next2, min(next3, min(next5, next7))))
if h[m].Cmp(next2) == 0 {
i++
next2.Mul(two, h[i])
}
if h[m].Cmp(next3) == 0 {
j++
next3.Mul(three, h[j])
}
if h[m].Cmp(next5) == 0 {
k++
next5.Mul(five, h[k])
}
if h[m].Cmp(next7) == 0 {
l++
next7.Mul(seven, h[l])
}
}
return h
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n)
le := len(s)
for i := le - 3; i >= 1; i -= 3 {
s = s[0:i] + "," + s[i:]
}
return s
}
func main() {
const n = 13 * 1e6 // calculate the first 13 million humble numbers, say
h := humble(n)
fmt.Println("The first 50 humble numbers are:")
fmt.Println(h[0:50])
maxDigits := len(h[len(h)-1].String()) - 1
counts := make([]int, maxDigits+1)
var maxUsed int
digits := 1
pow10 := new(big.Int).Set(ten)
for i := 0; i < len(h); i++ {
for {
if h[i].Cmp(pow10) >= 0 {
pow10.Mul(pow10, ten)
digits++
} else {
break
}
}
if digits > maxDigits {
maxUsed = i
break
}
counts[digits]++
}
fmt.Printf("\nOf the first %s humble numbers:\n", commatize(maxUsed))
for i := 1; i <= maxDigits; i++ {
s := "s"
if i == 1 {
s = ""
}
fmt.Printf("%9s have %2d digit%s\n", commatize(counts[i]), i, s)
}
}
- Output:
The first 50 humble numbers are:
[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]
Of the first 12,591,874 humble numbers:
9 have 1 digit
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1,272 have 8 digits
1,767 have 9 digits
2,381 have 10 digits
3,113 have 11 digits
3,984 have 12 digits
5,002 have 13 digits
6,187 have 14 digits
7,545 have 15 digits
9,081 have 16 digits
10,815 have 17 digits
12,759 have 18 digits
14,927 have 19 digits
17,323 have 20 digits
19,960 have 21 digits
22,853 have 22 digits
26,015 have 23 digits
29,458 have 24 digits
33,188 have 25 digits
37,222 have 26 digits
41,568 have 27 digits
46,245 have 28 digits
51,254 have 29 digits
56,618 have 30 digits
62,338 have 31 digits
68,437 have 32 digits
74,917 have 33 digits
81,793 have 34 digits
89,083 have 35 digits
96,786 have 36 digits
104,926 have 37 digits
113,511 have 38 digits
122,546 have 39 digits
132,054 have 40 digits
142,038 have 41 digits
152,515 have 42 digits
163,497 have 43 digits
174,986 have 44 digits
187,004 have 45 digits
199,565 have 46 digits
212,675 have 47 digits
226,346 have 48 digits
240,590 have 49 digits
255,415 have 50 digits
270,843 have 51 digits
286,880 have 52 digits
303,533 have 53 digits
320,821 have 54 digits
338,750 have 55 digits
357,343 have 56 digits
376,599 have 57 digits
396,533 have 58 digits
417,160 have 59 digits
438,492 have 60 digits
460,533 have 61 digits
483,307 have 62 digits
506,820 have 63 digits
531,076 have 64 digits
556,104 have 65 digits
581,902 have 66 digits
608,483 have 67 digits
635,864 have 68 digits
664,053 have 69 digits
693,065 have 70 digits
Haskell
import Data.Set (deleteFindMin, fromList, union)
import Data.List.Split (chunksOf)
import Data.List (group)
import Data.Bool (bool)
--------------------- HUMBLE NUMBERS ----------------------
humbles :: [Integer]
humbles = go $ fromList [1]
where
go sofar = x : go (union pruned $ fromList ((x *) <$> [2, 3, 5, 7]))
where
(x, pruned) = deleteFindMin sofar
-- humbles = filter (all (< 8) . primeFactors) [1 ..]
-------------------------- TEST ---------------------------
main :: IO ()
main = do
putStrLn "First 50 Humble numbers:"
mapM_ (putStrLn . concat) $
chunksOf 10 $ justifyRight 4 ' ' . show <$> take 50 humbles
putStrLn "\nCount of humble numbers for each digit length 1-25:"
mapM_ print $
take 25 $ ((,) . head <*> length) <$> group (length . show <$> humbles)
------------------------- DISPLAY -------------------------
justifyRight :: Int -> a -> [a] -> [a]
justifyRight n c = (drop . length) <*> (replicate n c ++)
- Output:
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-25: (1,9) (2,36) (3,95) (4,197) (5,356) (6,579) (7,882) (8,1272) (9,1767) (10,2381) (11,3113) (12,3984) (13,5002) (14,6187) (15,7545) (16,9081) (17,10815) (18,12759) (19,14927) (20,17323) (21,19960) (22,22853) (23,26015) (24,29458) (25,33188)
Much Faster Version with Linear Response for Range
The above version is quite "brute force" in that it collects all possible humble numbers in a Set persistent data structure, extracting them from smallest first while adding all of the new humble numbers based on each minimum. The Set takes care of the problem of eliminating repeat values, but at quite a large constant overhead as well as an extra `log n` execution performance cost where `n` is the size of the Set due to all persistent data structures. As well, there is a huge processing time cost due to grouping the values by their length after conversion to a string at the cost of a divide by ten for every digit and multiple levels of list processing.
From this Haskell implementation on the Hamming numbers task page, the method of finding the sequence by non-duplicates as well as using logarithmic approximations for size comparisons costs linear execution time with the length of the sequence, and can be easily adapted for use here just by adding the stage of using the additional prime of seven. In addition, the counting of the digit groups can be greatly simplified and thus faster when accomplished with a single loop based on the logarithmic approximations to reduce the garbage collection overheads of multiple steps of list processing, as per the following code:
{-# OPTIONS_GHC -O2 #-}
import Data.Word (Word16)
import Data.Bits (shiftR, (.&.))
import Data.Function (fix)
import Data.List (group, intercalate)
import Data.Time.Clock.POSIX (getPOSIXTime)
--------------------- HUMBLE NUMBERS ----------------------
data LogRep = LogRep {-# UNPACK #-} !Double
{-# UNPACK #-} !Word16
{-# UNPACK #-} !Word16
{-# UNPACK #-} !Word16
{-# UNPACK #-} !Word16 deriving Show
instance Eq LogRep where
(==) (LogRep la _ _ _ _) (LogRep lb _ _ _ _) = la == lb
instance Ord LogRep where
(<=) (LogRep la _ _ _ _) (LogRep lb _ _ _ _) = la <= lb
logrep2Integer :: LogRep -> Integer
logrep2Integer (LogRep _ w x y z) = xpnd 2 w $ xpnd 3 x $ xpnd 5 y $ xpnd 7 z 1 where
xpnd m v = go v m where
go i mlt acc =
if i <= 0 then acc else
go (i `shiftR` 1) (mlt * mlt) (if i .&. 1 == 0 then acc else acc * mlt)
cOneLR :: LogRep
cOneLR = LogRep 0.0 0 0 0 0
cLgOf2 :: Double
cLgOf2 = logBase 10 2
cLgOf3 :: Double
cLgOf3 = logBase 10 3
cLgOf5 :: Double
cLgOf5 = logBase 10 5
cLgOf7 :: Double
cLgOf7 = logBase 10 7
cLgOf10 :: Double
cLgOf10 = cLgOf2 + cLgOf5
mulLR2 :: LogRep -> LogRep
mulLR2 (LogRep lg w x y z) = LogRep (lg + cLgOf2) (w + 1) x y z
mulLR3 :: LogRep -> LogRep
mulLR3 (LogRep lg w x y z) = LogRep (lg + cLgOf3) w (x + 1) y z
mulLR5 :: LogRep -> LogRep
mulLR5 (LogRep lg w x y z) = LogRep (lg + cLgOf5) w x (y + 1) z
mulLR7 :: LogRep -> LogRep
mulLR7 (LogRep lg w x y z) = LogRep (lg + cLgOf7) w x y (z + 1)
humbleLRs :: () -> [LogRep]
humbleLRs() = cOneLR : foldr u [] [mulLR2, mulLR3, mulLR5, mulLR7] where
u nmf s = fix (merge s . map nmf . (cOneLR:)) where
merge a [] = a
merge [] b = b
merge a@(x:xs) b@(y:ys) | x < y = x : merge xs b
| otherwise = y : merge a ys
-------------------------- TEST ---------------------------
main :: IO ()
main = do
putStrLn "First 50 Humble numbers:"
mapM_ (putStrLn . concat) $
chunksOf 10 $ justifyRight 4 ' ' . show <$> take 50 (map logrep2Integer $ humbleLRs())
strt <- getPOSIXTime
putStrLn "\nCount of humble numbers for each digit length 1-255:"
putStrLn "Digits Count Accum"
mapM_ putStrLn $ take 255 $ groupFormat $ humbleLRs()
stop <- getPOSIXTime
putStrLn $ "Counting took " ++ show (1.0 * (stop - strt)) ++ " seconds."
------------------------- DISPLAY -------------------------
chunksOf :: Int -> [a] -> [[a]]
chunksOf n = go where
go [] = []
go ilst = take n ilst : go (drop n ilst)
justifyRight :: Int -> a -> [a] -> [a]
justifyRight n c = (drop . length) <*> (replicate n c ++)
commaString :: String -> String
commaString =
let grpsOf3 [] = []
grpsOf3 is = let (frst, rest) = splitAt 3 is in frst : grpsOf3 rest
in reverse . intercalate "," . grpsOf3 . reverse
groupFormat :: [LogRep] -> [String]
groupFormat = go (0 :: Int) (0 :: Int) 0 where
go _ _ _ [] = []
go i cnt cacc ((LogRep lg _ _ _ _) : lrtl) =
let nxt = truncate (lg / cLgOf10) :: Int in
if nxt == i then go i (cnt + 1) cacc lrtl else
let ni = i + 1
ncacc = cacc + cnt
str = justifyRight 4 ' ' (show ni) ++
justifyRight 14 ' ' (commaString $ show cnt) ++
justifyRight 19 ' ' (commaString $ show ncacc)
in str : go ni 1 ncacc lrtl
- Output:
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 169.193563058s seconds.
The above abbreviated results match the results from the C++ Direct-Generation contribution as well as the Pascal contribution. The above code, as run on an Intel i5-6500 at 3.6 GHz with single-thread boost, still spends almost two thirds of the execution time doing garbage collection related to lazy list processing, which could almost be eliminated by replacing the inner lists with "deque" first in/first out mutable array-based data structures, which would also reduce memory use by some constant factor, but that hardly seems worth the effort when there are even faster techniques available for this specific set of tasks that almost eliminate the huge memory use and most of the processing time overheads.
Direct Generation of Digit Counts from Logarithms
From the C++ Direct Generation contribution, there is no need to generate the ordered sequence of humble numbers to be able to bin them according to digits length; that technique, which reduces all of the execution time expended in sorting the humble numbers stream in increasing order, works fine in Haskell, as well reducing the huge memory requirements as it only requires space for the binning array. The use of the extended precision logarithmic approximations as in 64-bit rather than the 53-bit precision of 64-bit floats may be faster to compare 64-bit integers rather than `Double`'s. For the trivial exercise of generating the first 50 humble numbers in order, any algorithm including a brute force one may be used, but here a more direct version of the non-duplicates version is used, as per the following code:
{-# OPTIONS_GHC -O2 #-}
{-# LANGUAGE FlexibleContexts #-}
import Data.Word (Word64)
import Data.Bits (shiftR)
import Data.Function (fix)
import Data.Array.Unboxed (UArray, elems)
import Data.Array.Base (MArray(newArray, unsafeRead, unsafeWrite))
import Data.Array.IO (IOUArray)
import Data.List (intercalate)
import Data.Time.Clock.POSIX (getPOSIXTime)
import Data.Array.Unsafe (unsafeFreeze)
cNumDigits :: Int
cNumDigits = 877
cShift :: Int
cShift = 50
cFactor :: Word64
cFactor = 2^cShift
cLogOf10 :: Word64
cLogOf10 = cFactor
cLogOf7 :: Word64
cLogOf7 = round $ (logBase 10 7 :: Double) * fromIntegral cFactor
cLogOf5 :: Word64
cLogOf5 = round $ (logBase 10 5 :: Double) * fromIntegral cFactor
cLogOf3 :: Word64
cLogOf3 = round $ (logBase 10 3 :: Double) * fromIntegral cFactor
cLogOf2 :: Word64
cLogOf2 = cLogOf10 - cLogOf5
cLogLmt :: Word64
cLogLmt = fromIntegral cNumDigits * cLogOf10
humbles :: () -> [Integer]
humbles() = 1 : foldr u [] [2,3,5,7] where
u n s = fix (merge s . map (n*) . (1:)) where
merge a [] = a
merge [] b = b
merge a@(x:xs) b@(y:ys) | x < y = x : merge xs b
| otherwise = y : merge a ys
-------------------------- TEST ---------------------------
main :: IO ()
main = do
putStrLn "First 50 humble numbers:"
mapM_ (putStrLn . concat) $
chunksOf 10 $ justifyRight 4 ' ' . show <$> take 50 (humbles())
putStrLn $ "\nCount of humble numbers for each digit length 1-"
++ show cNumDigits ++ ":"
putStrLn "Digits Count Accum"
strt <- getPOSIXTime
mbins <- newArray (0, cNumDigits - 1) 0 :: IO (IOUArray Int Int)
let loopw w =
if w >= cLogLmt then return () else
let loopx x =
if x >= cLogLmt then loopw (w + cLogOf2) else
let loopy y =
if y >= cLogLmt then loopx (x + cLogOf3) else
let loopz z =
if z >= cLogLmt then loopy (y + cLogOf5) else do
let ndx = fromIntegral (z `shiftR` cShift)
v <- unsafeRead mbins ndx
unsafeWrite mbins ndx (v + 1)
loopz (z + cLogOf7) in loopz y in loopy x in loopx w
loopw 0
stop <- getPOSIXTime
bins <- unsafeFreeze mbins :: IO (UArray Int Int)
mapM_ putStrLn $ format $ elems bins
putStrLn $ "Counting took " ++ show (realToFrac (stop - strt)) ++ " seconds."
------------------------- DISPLAY -------------------------
chunksOf :: Int -> [a] -> [[a]]
chunksOf n = go where
go [] = []
go ilst = take n ilst : go (drop n ilst)
justifyRight :: Int -> a -> [a] -> [a]
justifyRight n c = (drop . length) <*> (replicate n c ++)
commaString :: String -> String
commaString =
let grpsOf3 [] = []
grpsOf3 s = let (frst, rest) = splitAt 3 s in frst : grpsOf3 rest
in reverse . intercalate "," . grpsOf3 . reverse
format :: [Int] -> [String]
format = go (0 :: Int) 0 where
go _ _ [] = []
go i cacc (hd : tl) =
let ni = i + 1
ncacc = cacc + hd
str = justifyRight 4 ' ' (show ni) ++
justifyRight 14 ' ' (commaString $ show hd) ++
justifyRight 19 ' ' (commaString $ show ncacc)
in str : go ni ncacc tl
- Output:
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 177.070331827 seconds.
This result is as run on an Intel i5-6500 at 3.6 GHz for single-threaded boost and the results are the same as for the C++ Direct-Generation contribution as well as the Pascal contribution. The results are perhaps a little faster than the C++ code from which this was translated, likely due to scaling so that the bin index is calculated by a simple bit shift rather than a division.
This code can likely be used to count the number of humble numbers for all digits up to 4000 digits in under a day (untested).
J
Multiply all the humble numbers by all the factors appending the next largest value.
humble=: 4 : 0
NB. x humble y generates x humble numbers based on factors y
result=. , 1
while. x > # result do.
a=. , result */ y
result=. result , <./ a #~ a > {: result
end.
)
p: i.4 2 3 5 7 50 humble p: i. 4 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
humbler tests a number for humbleness by deciding if the prime factors are a subset of the given factors.
humbler=: '' -: (-.~ q:) NB. x humbler y tests whether factors of y are a (improper)subset of x
50 {. (#~ (p:i.4)&humbler&>) >: i. 500
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
Use a class to simulate the python generator. This is a more efficient implementation of the first method.
FACTORS_h_=: p: i. 4
HUMBLE_h_=: 1
next_h_=: 3 : 0
result=. <./ HUMBLE
i=. HUMBLE i. result
HUMBLE=: ~. (((i&{.) , (>:i)&}.) HUMBLE) , result * FACTORS
result
)
reset_h_=: 3 :'0 $ HUMBLE=: 1'
3 :0 [ 50 [ reset_h_'' result=.'' for.i.y do. result=. result,next_h_'' end. ) 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
Tally of humble numbers having up to so many digits. Use up to forty thousand humble numbers.
H=: 3 :0 [ 40000 [ reset_h_''
result=.''
for.i.y do.
result=. result,next_h_''
end.
)
10^.{:H NB. log of tail is number of digits
15.7432
(,. [: +/ H </ 10&^) i.16
0 0
1 9
2 45
3 140
4 337
5 693
6 1272
7 2154
8 3426
9 5193
10 7574
11 10687
12 14671
13 19673
14 25860
15 33405
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class HumbleNumbers {
public static void main(String[] args) {
System.out.println("First 50 humble numbers:");
System.out.println(Arrays.toString(humble(50)));
Map<Integer,Integer> lengthCountMap = new HashMap<>();
BigInteger[] seq = humble(1_000_000);
for ( int i = 0 ; i < seq.length ; i++ ) {
BigInteger humbleNumber = seq[i];
int len = humbleNumber.toString().length();
lengthCountMap.merge(len, 1, (v1, v2) -> v1 + v2);
}
List<Integer> sorted = new ArrayList<>(lengthCountMap.keySet());
Collections.sort(sorted);
System.out.printf("Length Count%n");
for ( Integer len : sorted ) {
System.out.printf(" %2s %5s%n", len, lengthCountMap.get(len));
}
}
private static BigInteger[] humble(int n) {
BigInteger two = BigInteger.valueOf(2);
BigInteger twoTest = two;
BigInteger three = BigInteger.valueOf(3);
BigInteger threeTest = three;
BigInteger five = BigInteger.valueOf(5);
BigInteger fiveTest = five;
BigInteger seven = BigInteger.valueOf(7);
BigInteger sevenTest = seven;
BigInteger[] results = new BigInteger[n];
results[0] = BigInteger.ONE;
int twoIndex = 0, threeIndex = 0, fiveIndex = 0, sevenIndex = 0;
for ( int index = 1 ; index < n ; index++ ) {
results[index] = twoTest.min(threeTest).min(fiveTest).min(sevenTest);
if ( results[index].compareTo(twoTest) == 0 ) {
twoIndex++;
twoTest = two.multiply(results[twoIndex]);
}
if (results[index].compareTo(threeTest) == 0 ) {
threeIndex++;
threeTest = three.multiply(results[threeIndex]);
}
if (results[index].compareTo(fiveTest) == 0 ) {
fiveIndex++;
fiveTest = five.multiply(results[fiveIndex]);
}
if (results[index].compareTo(sevenTest) == 0 ) {
sevenIndex++;
sevenTest = seven.multiply(results[sevenIndex]);
}
}
return results;
}
}
- Output:
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]
Length Count
1 9
2 36
3 95
4 197
5 356
6 579
7 882
8 1272
9 1767
10 2381
11 3113
12 3984
13 5002
14 6187
15 7545
16 9081
17 10815
18 12759
19 14927
20 17323
21 19960
22 22853
23 26015
24 29458
25 33188
26 37222
27 41568
28 46245
29 51254
30 56618
31 62338
32 68437
33 74917
34 81793
35 89083
36 96786
37 63955
JavaScript
(() => {
'use strict';
// ------------------ HUMBLE NUMBERS -------------------
// humbles :: () -> [Int]
function* humbles() {
// A non-finite stream of Humble numbers.
// OEIS A002473
const hs = new Set([1]);
while (true) {
let nxt = Math.min(...hs)
yield nxt;
hs.delete(nxt);
[2, 3, 5, 7].forEach(
x => hs.add(nxt * x)
);
}
};
// ----------------------- TEST ------------------------
// main :: IO ()
const main = () => {
console.log('First 50 humble numbers:\n')
chunksOf(10)(take(50)(humbles())).forEach(
row => console.log(
concat(row.map(compose(justifyRight(4)(' '), str)))
)
);
console.log(
'\nCounts of humble numbers of given digit lengths:'
);
const
counts = map(length)(
group(takeWhileGen(x => 11 > x)(
fmapGen(x => str(x).length)(
humbles()
)
))
);
console.log(
fTable('\n')(str)(str)(
i => counts[i - 1]
)(enumFromTo(1)(10))
);
};
// ----------------- GENERIC FUNCTIONS -----------------
// chunksOf :: Int -> [a] -> [[a]]
const chunksOf = n =>
xs => enumFromThenTo(0)(n)(
xs.length - 1
).reduce(
(a, i) => a.concat([xs.slice(i, (n + i))]),
[]
);
// compose (<<<) :: (b -> c) -> (a -> b) -> a -> c
const compose = (...fs) =>
// A function defined by the right-to-left
// composition of all the functions in fs.
fs.reduce(
(f, g) => x => f(g(x)),
x => x
);
// concat :: [[a]] -> [a]
// concat :: [String] -> String
const concat = xs =>
0 < xs.length ? (() => {
const unit = 'string' !== typeof xs[0] ? (
[]
) : '';
return unit.concat.apply(unit, xs);
})() : [];
// enumFromThenTo :: Int -> Int -> Int -> [Int]
const enumFromThenTo = x1 =>
x2 => y => {
const d = x2 - x1;
return Array.from({
length: Math.floor(y - x2) / d + 2
}, (_, i) => x1 + (d * i));
};
// enumFromTo :: Int -> Int -> [Int]
const enumFromTo = m =>
n => Array.from({
length: 1 + n - m
}, (_, i) => m + i);
// fTable :: String -> (a -> String) -> (b -> String)
// -> (a -> b) -> [a] -> String
const fTable = s => xShow => fxShow => f => xs => {
// Heading -> x display function ->
// fx display function ->
// f -> values -> tabular string
const
ys = xs.map(xShow),
w = Math.max(...ys.map(length));
return s + '\n' + zipWith(
a => b => a.padStart(w, ' ') + ' -> ' + b
)(ys)(
xs.map(x => fxShow(f(x)))
).join('\n');
};
// fmapGen <$> :: (a -> b) -> Gen [a] -> Gen [b]
const fmapGen = f =>
function*(gen) {
let v = take(1)(gen);
while (0 < v.length) {
yield(f(v[0]))
v = take(1)(gen)
}
};
// group :: Eq a => [a] -> [[a]]
const group = xs => {
// A list of lists, each containing only equal elements,
// such that the concatenation of these lists is xs.
const go = xs =>
0 < xs.length ? (() => {
const
h = xs[0],
i = xs.findIndex(x => h !== x);
return i !== -1 ? (
[xs.slice(0, i)].concat(go(xs.slice(i)))
) : [xs];
})() : [];
const v = go(list(xs));
return 'string' === typeof xs ? (
v.map(x => x.join(''))
) : v;
};
// justifyRight :: Int -> Char -> String -> String
const justifyRight = n =>
// The string s, preceded by enough padding (with
// the character c) to reach the string length n.
c => s => n > s.length ? (
s.padStart(n, c)
) : s;
// length :: [a] -> Int
const length = xs =>
// Returns Infinity over objects without finite length.
// This enables zip and zipWith to choose the shorter
// argument when one is non-finite, like cycle, repeat etc
(Array.isArray(xs) || 'string' === typeof xs) ? (
xs.length
) : Infinity;
// list :: StringOrArrayLike b => b -> [a]
const list = xs =>
// xs itself, if it is an Array,
// or an Array derived from xs.
Array.isArray(xs) ? (
xs
) : Array.from(xs);
// map :: (a -> b) -> [a] -> [b]
const map = f =>
// The list obtained by applying f
// to each element of xs.
// (The image of xs under f).
xs => [...xs].map(f);
// str :: a -> String
const str = x =>
x.toString();
// take :: Int -> [a] -> [a]
// take :: Int -> String -> String
const take = n => xs =>
'GeneratorFunction' !== xs.constructor.constructor.name ? (
xs.slice(0, n)
) : [].concat.apply([], Array.from({
length: n
}, () => {
const x = xs.next();
return x.done ? [] : [x.value];
}));
// takeWhileGen :: (a -> Bool) -> Gen [a] -> [a]
const takeWhileGen = p => xs => {
const ys = [];
let
nxt = xs.next(),
v = nxt.value;
while (!nxt.done && p(v)) {
ys.push(v);
nxt = xs.next();
v = nxt.value
}
return ys;
};
// zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
const zipWith = f =>
// Use of `take` and `length` here allows zipping with non-finite lists
// i.e. generators like cycle, repeat, iterate.
xs => ys => {
const n = Math.min(length(xs), length(ys));
return Infinity > n ? (
(([as, bs]) => Array.from({
length: n
}, (_, i) => f(as[i])(
bs[i]
)))([xs, ys].map(
compose(take(n), list)
))
) : zipWithGen(f)(xs)(ys);
};
// MAIN ---
return main();
})();
- Output:
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 Counts of humble numbers of given digit lengths: 1 -> 9 2 -> 36 3 -> 95 4 -> 197 5 -> 356 6 -> 579 7 -> 882 8 -> 1272 9 -> 1767 10 -> 2381
jq
Brute force
First, brute force (because we can) ...
# Input: a positive integer
# Output: true iff the input is humble
def humble:
. as $i
| if ($i < 2) then true
elif ($i % 2 == 0) then ($i / 2 | floor) | humble
elif ($i % 3 == 0) then ($i / 3 | floor) | humble
elif ($i % 5 == 0) then ($i / 5 | floor) | humble
elif ($i % 7 == 0) then ($i / 7 | floor) | humble
else false
end;
def task($digits; $count):
{len: 0, i:0, count: 0}
| until( .len > $digits;
.i += 1
| if (.i|humble)
then .len = (.i | tostring | length)
| if .len <= $digits
then .humble[.len] += 1
| if .count < $count then .count += 1 | .humbles += [.i] else . end
else .
end
else .
end)
| "First \($count):", .humbles,
"Distribution of the number of decimal digits up to \($digits) digits:",
(.humble | range(1;length) as $i | " \($i): \(.[$i])") ;
task(6; 50)
- Output:
First 50: [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] Distribution of the number of decimal digits up to 6 digits: 1: 9 2: 36 3: 95 4: 197 5: 356 6: 579
Fast and economical
(At least up to jq's numerical precision...)
Having already shown one way to display the first few humble numbers, this subsection will focus on the more difficult problem.
# A generator
def humbles($digits):
def update:
.[0] as $h
| ([$h * 2, $h * 3, $h * 5, $h * 7] | map(select(tostring|length <= $digits))) as $next
| if $next == [] then .[1:]
else (.[1:] + $next) | sort
end;
def trim: if length <= 1 then . elif .[0]==.[1] then .[1:]|trim else . end;
{ queue: [1]}
| recurse( select( .queue != [] )
| .h = .queue[0]
| queue |= (update|trim) )
| .h ;
def distribution(stream):
reduce stream as $i ([]; .[$i|tostring|length-1]+=1);
def task($digits):
"Distribution of the number of decimal digits up to \($digits) digits:",
(distribution(humbles($digits)) | range(0;length) as $i | " \($i+1): \(.[$i])") ;
task(16)- Output:
Distribution of the number of decimal digits up to 16 digits: 1: 9 2: 36 3: 95 4: 197 5: 356 6: 579 7: 882 8: 1272 9: 1767 10: 2381 11: 3113 12: 3984 13: 5002 14: 6187 15: 7545 16: 9081
The queue in this case grows to a maximum length of about 18K.
Julia
To spare heap memory, keeps only the last 2 million values found for use in the generation of further values.
function counthumbledigits(maxdigits, returnsequencelength=50)
n, count, adjustindex, maxdiff = BigInt(1), 0, BigInt(0), 0
humble, savesequence = Vector{BigInt}([1]), Vector{BigInt}()
base2, base3, base5, base7 = 1, 1, 1, 1
next2, next3, next5, next7 = BigInt(2), BigInt(3), BigInt(5), BigInt(7)
digitcounts= Dict{Int, Int}(1 => 1)
while n < BigInt(10)^(maxdigits+1)
n = min(next2, next3, next5, next7)
push!(humble, n)
count += 1
if count == returnsequencelength
savesequence = deepcopy(humble[1:returnsequencelength])
elseif count > 2000000
popfirst!(humble)
adjustindex += 1
end
placesbase10 = length(string(n))
if haskey(digitcounts, placesbase10)
digitcounts[placesbase10] += 1
else
digitcounts[placesbase10] = 1
end
maxdiff = max(maxdiff, count - base2, count - base3, count - base5, count - base7)
(next2 <= n) && (next2 = 2 * humble[(base2 += 1) - adjustindex])
(next3 <= n) && (next3 = 3 * humble[(base3 += 1) - adjustindex])
(next5 <= n) && (next5 = 5 * humble[(base5 += 1) - adjustindex])
(next7 <= n) && (next7 = 7 * humble[(base7 += 1) - adjustindex])
end
savesequence, digitcounts, count, maxdiff
end
counthumbledigits(3)
@time first120, digitcounts, count, maxdiff = counthumbledigits(99)
println("\nTotal humble numbers counted: $count")
println("Maximum depth between top of array and a multiplier: $maxdiff\n")
println("The first 50 humble numbers are: $first120\n\nDigit counts of humble numbers:")
for ndigits in sort(collect(keys(digitcounts)))[1:end-1]
println(lpad(digitcounts[ndigits], 10), " have ", lpad(ndigits, 3), " digits.")
end
- Output:
828.693164 seconds (3.61 G allocations: 64.351 GiB, 51.37% gc time)
Total humble numbers counted: 51428827
Maximum depth between top of array and a multiplier: 1697189
The first 50 humble numbers are: BigInt[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]
Digit counts of humble numbers:
9 have 1 digits.
36 have 2 digits.
95 have 3 digits.
197 have 4 digits.
356 have 5 digits.
579 have 6 digits.
882 have 7 digits.
1272 have 8 digits.
1767 have 9 digits.
2381 have 10 digits.
3113 have 11 digits.
3984 have 12 digits.
5002 have 13 digits.
6187 have 14 digits.
7545 have 15 digits.
9081 have 16 digits.
10815 have 17 digits.
12759 have 18 digits.
14927 have 19 digits.
17323 have 20 digits.
19960 have 21 digits.
22853 have 22 digits.
26015 have 23 digits.
29458 have 24 digits.
33188 have 25 digits.
37222 have 26 digits.
41568 have 27 digits.
46245 have 28 digits.
51254 have 29 digits.
56618 have 30 digits.
62338 have 31 digits.
68437 have 32 digits.
74917 have 33 digits.
81793 have 34 digits.
89083 have 35 digits.
96786 have 36 digits.
104926 have 37 digits.
113511 have 38 digits.
122546 have 39 digits.
132054 have 40 digits.
142038 have 41 digits.
152515 have 42 digits.
163497 have 43 digits.
174986 have 44 digits.
187004 have 45 digits.
199565 have 46 digits.
212675 have 47 digits.
226346 have 48 digits.
240590 have 49 digits.
255415 have 50 digits.
270843 have 51 digits.
286880 have 52 digits.
303533 have 53 digits.
320821 have 54 digits.
338750 have 55 digits.
357343 have 56 digits.
376599 have 57 digits.
396533 have 58 digits.
417160 have 59 digits.
438492 have 60 digits.
460533 have 61 digits.
483307 have 62 digits.
506820 have 63 digits.
531076 have 64 digits.
556104 have 65 digits.
581902 have 66 digits.
608483 have 67 digits.
635864 have 68 digits.
664053 have 69 digits.
693065 have 70 digits.
722911 have 71 digits.
753593 have 72 digits.
785141 have 73 digits.
817554 have 74 digits.
850847 have 75 digits.
885037 have 76 digits.
920120 have 77 digits.
956120 have 78 digits.
993058 have 79 digits.
1030928 have 80 digits.
1069748 have 81 digits.
1109528 have 82 digits.
1150287 have 83 digits.
1192035 have 84 digits.
1234774 have 85 digits.
1278527 have 86 digits.
1323301 have 87 digits.
1369106 have 88 digits.
1415956 have 89 digits.
1463862 have 90 digits.
1512840 have 91 digits.
1562897 have 92 digits.
1614050 have 93 digits.
1666302 have 94 digits.
1719669 have 95 digits.
1774166 have 96 digits.
1829805 have 97 digits.
1886590 have 98 digits.
1944540 have 99 digits.
2003661 have 100 digits.
Kotlin
fun isHumble(i: Int): Boolean {
if (i <= 1) return true
if (i % 2 == 0) return isHumble(i / 2)
if (i % 3 == 0) return isHumble(i / 3)
if (i % 5 == 0) return isHumble(i / 5)
if (i % 7 == 0) return isHumble(i / 7)
return false
}
fun main() {
val limit: Int = Short.MAX_VALUE.toInt()
val humble = mutableMapOf<Int, Int>()
var count = 0
var num = 1
while (count < limit) {
if (isHumble(num)) {
val str = num.toString()
val len = str.length
humble.merge(len, 1) { a, b -> a + b }
if (count < 50) print("$num ")
count++
}
num++
}
println("\n")
println("Of the first $count humble numbers:")
num = 1
while (num < humble.size - 1) {
if (humble.containsKey(num)) {
val c = humble[num]
println("%5d have %2d digits".format(c, num))
num++
} else {
break
}
}
}
- Output:
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
Of the first 32767 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Lua
function isHumble(n)
local n2 = math.floor(n)
if n2 <= 1 then
return true
end
if n2 % 2 == 0 then
return isHumble(n2 / 2)
end
if n2 % 3 == 0 then
return isHumble(n2 / 3)
end
if n2 % 5 == 0 then
return isHumble(n2 / 5)
end
if n2 % 7 == 0 then
return isHumble(n2 / 7)
end
return false
end
function main()
local limit = 10000
local humble = {0, 0, 0, 0, 0, 0, 0, 0, 0}
local count = 0
local num = 1
while count < limit do
if isHumble(num) then
local buffer = string.format("%d", num)
local le = string.len(buffer)
if le > 9 then
break
end
humble[le] = humble[le] + 1
if count < 50 then
io.write(num .. " ")
end
count = count + 1
end
num = num + 1
end
print("\n")
print("Of the first " .. count .. " humble numbers:")
for num=1,9 do
print(string.format("%5d have %d digits", humble[num], num))
end
end
main()
- Output:
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
Of the first 5193 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Mathematica/Wolfram Language
Create a simple function which efficiently generates humble numbers up to an inputted max number, then call it twice to generate the output. Finds the number of humble numbers with digits up to 100 in 5 minutes.
HumbleGenerator[max_] :=
Sort[Flatten@ParallelTable[
2^i 3^j 5^k 7^m, {i, 0, Log[2, max]}, {j, 0, Log[3, max/2^i]}, {k,
0, Log[5, max/(2^i 3^j)]}, {m, 0, Log[7, max/(2^i 3^j 5^k)]}]]
{"First 50 Humble Numbers:",
HumbleGenerator[120],
"\nDigits\[Rule]Count",
Rule @@@ Tally[IntegerLength /@ Drop[HumbleGenerator[10^100], -1]] //
Column} // Column
- Output:
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}
Digits->Count
1->9
2->36
3->95
4->197
5->356
6->579
7->882
8->1272
9->1767
10->2381
11->3113
12->3984
13->5002
14->6187
15->7545
16->9081
17->10815
18->12759
19->14927
20->17323
21->19960
22->22853
23->26015
24->29458
25->33188
26->37222
27->41568
28->46245
29->51254
30->56618
31->62338
32->68437
33->74917
34->81793
35->89083
36->96786
37->104926
38->113511
39->122546
40->132054
41->142038
42->152515
43->163497
44->174986
45->187004
46->199565
47->212675
48->226346
49->240590
50->255415
51->270843
52->286880
53->303533
54->320821
55->338750
56->357343
57->376599
58->396533
59->417160
60->438492
61->460533
62->483307
63->506820
64->531076
65->556104
66->581902
67->608483
68->635864
69->664053
70->693065
71->722911
72->753593
73->785141
74->817554
75->850847
76->885037
77->920120
78->956120
79->993058
80->1030928
81->1069748
82->1109528
83->1150287
84->1192035
85->1234774
86->1278527
87->1323301
88->1369106
89->1415956
90->1463862
91->1512840
92->1562897
93->1614050
94->1666302
95->1719669
96->1774166
97->1829805
98->1886590
99->1944540
100->2003661
Nim
A simple algorithm efficient enough to get the number of humble numbers with 18 digits in less than four seconds. To get further, we would have to use big numbers and a more efficient algorithm.
import sets, strformat
proc min[T](s: HashSet[T]): T =
## Return the minimal value in a set.
if s.card == 0:
raise newException(ValueError, "set is empty.")
result = T.high
for n in s:
if n < result: result = n
iterator humbleNumbers(): Positive =
## Yield the successive humble numbers.
var s = [1].toHashSet()
while true:
let m = min(s)
yield m
s.excl(m)
for k in [2, 3, 5, 7]: s.incl(k * m)
proc showHumbleNumbers(maxCount: Positive) =
## Show the first "maxCount" humble numbers.
var currCount = 0
for n in humbleNumbers():
stdout.write n
inc currCount
if currCount == maxCount: break
stdout.write ' '
echo ""
proc showHumbleCount(ndigits: Positive) =
## Show the number of humble numbers with "n <= ndigits" digits.
echo "Digits Count"
echo "------ -----"
var currdigits = 1
var count = 0
var next = 10 # First number with "currdigits + 1" digits.
for n in humbleNumbers():
if n >= next:
# Number of digits has changed.
echo &" {currdigits:2d} {count:5d}"
inc currdigits
if currdigits > ndigits: break
next *= 10
count = 1
else:
inc count
when isMainModule:
echo "First 50 humble numbers:"
showHumbleNumbers(50)
echo ""
echo "Count of humble numbers with n digits:"
showHumbleCount(18)
- Output:
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 with n digits: Digits Count ------ ----- 1 9 2 36 3 95 4 197 5 356 6 579 7 882 8 1272 9 1767 10 2381 11 3113 12 3984 13 5002 14 6187 15 7545 16 9081 17 10815 18 12759
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 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:
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."
- Output:
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.
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:
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."
- Output:
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.
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.
Pascal
modification of hamming-numbers http://rosettacode.org/wiki/Hamming_numbers#a_fast_alternative
Check for the first occurrence of 2^i/5^i -> 10^i
Using float80/extended and float64/double and single/float32 version, to get a possibility to check values.
Up to 877 digits I was able to compare, than float80 run out of memory (13.5 Gb )
float80 and float64 got same values up to 877.float64 is 2,3x faster than float80
float32 get wrong at 37 digits,->37 104925 instead of 104926
runtime: 2 x digits => ~ runtime 2^4
{$IFDEF FPC}
{$MODE DELPHI}
{$OPTIMIZATION ON,ALL}
{$CODEALIGN proc=32,loop=1}
{$ALIGN 16}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils;
const
//PlInit(4); <= maxPrimFakCnt+1
maxPrimFakCnt = 3;//3,7,11 to keep data 8-Byte aligned
minElemCnt = 10;
type
tPrimList = array of NativeUint;
tnumber = extended;
tpNumber= ^tnumber;
tElem = record
n : tnumber;//ln(prime[0]^Pots[0]*...
dummy: array[0..5] of byte;//extend extended to 16 byte
Pots: array[0..maxPrimFakCnt] of word;
end;
tpElem = ^tElem;
tElems = array of tElem;
tElemArr = array [0..0] of tElem;
tpElemArr = ^tElemArr;
tpFaktorRec = ^tFaktorRec;
tFaktorRec = record
frElems : tElems;
frInsElems: tElems;
frAktIdx : NativeUint;
frMaxIdx : NativeUint;
frPotNo : NativeUint;
frActPot : NativeUint;
frNextFr : tpFaktorRec;
frActNumb: tElem;
frLnPrime: tnumber;
end;
tArrFR = array of tFaktorRec;
//LoE == List of Elements
function LoEGetNextElement(pFR :tpFaktorRec):tElem;forward;
var
Pl : tPrimList;
ActIndex : NativeUint;
ArrInsert : tElems;
procedure PlInit(n: integer);
const
cPl : array[0..11] of byte=(2,3,5,7,11,13,17,19,23,29,31,37);
var
i : integer;
Begin
IF n>High(cPl)+1 then
n := High(cPl)
else
IF n < 0 then
n := 1;
IF maxPrimFakCnt+1 < n then
Begin
writeln(' Need to compile with bigger maxPrimFakCnt ');
Halt(0);
end;
setlength(Pl,n);
dec(n);
For i := 0 to n do
Pl[i] := cPl[i];
end;
procedure AusgabeElem(pElem: tElem;ShowPot:Boolean=false);
var
i : integer;
Begin
with pElem do
Begin
IF n < 23 then
write(round(exp(n)),' ')
else
write('ln ',n:13:7);
IF ShowPot then
Begin
For i := 0 to maxPrimFakCnt do
write(' ',PL[i]:2,'^',Pots[i]);
writeln
end;
end;
end;
procedure LoECreate(const Pl: tPrimList;var FA:tArrFR);
var
i : integer;
Begin
setlength(ArrInsert,100);
setlength(FA,Length(PL));
For i := 0 to High(PL) do
with FA[i] do
Begin
//automatic zeroing
IF i < High(PL) then
Begin
setlength(frElems,minElemCnt);
setlength(frInsElems,minElemCnt);
frNextFr := @FA[i+1]
end
else
Begin
setlength(frElems,2);
setlength(frInsElems,0);
frNextFr := NIL;
end;
frPotNo := i;
frLnPrime:= ln(PL[i]);
frMaxIdx := 0;
frAktIdx := 0;
frActPot := 1;
With frElems[0] do
Begin
n := frLnPrime;
Pots[i]:= 1;
end;
frActNumb := frElems[0];
end;
ActIndex := 0;
end;
procedure LoEFree(var FA:tArrFR);
var
i : integer;
Begin
For i := High(FA) downto Low(FA) do
setlength(FA[i].frElems,0);
setLength(FA,0);
end;
function LoEGetActElem(pFr:tpFaktorRec):tElem;inline;
Begin
with pFr^ do
result := frElems[frAktIdx];
end;
function LoEGetActLstNumber(pFr:tpFaktorRec):tpNumber;inline;
Begin
with pFr^ do
result := @frElems[frAktIdx].n;
end;
procedure LoEIncInsArr(var a:tElems);
Begin
setlength(a,Length(a)*8 div 5);
end;
procedure LoEIncreaseElems(pFr:tpFaktorRec;minCnt:NativeUint);
var
newLen: NativeUint;
Begin
with pFR^ do
begin
newLen := Length(frElems);
minCnt := minCnt+frMaxIdx;
repeat
newLen := newLen*8 div 5 +1;
until newLen > minCnt;
setlength(frElems,newLen);
end;
end;
procedure LoEInsertNext(pFr:tpFaktorRec;Limit:tnumber);
var
pNum : tpNumber;
pElems : tpElemArr;
cnt,i,u : NativeInt;
begin
with pFr^ do
Begin
//collect numbers of heigher primes
cnt := 0;
pNum := LoEGetActLstNumber(frNextFr);
while Limit > pNum^ do
Begin
frInsElems[cnt] := LoEGetNextElement(frNextFr);
inc(cnt);
IF cnt > High(frInsElems) then
LoEIncInsArr(frInsElems);
pNum := LoEGetActLstNumber(frNextFr);
end;
if cnt = 0 then
EXIT;
i := frMaxIdx;
u := frMaxIdx+cnt+1;
IF u > High(frElems) then
LoEIncreaseElems(pFr,cnt);
IF frPotNo = 0 then
inc(ActIndex,u);
//Merge
pElems := @frElems[0];
dec(cnt);
dec(u);
frMaxIdx:= u;
repeat
IF pElems^[i].n < frInsElems[cnt].n then
Begin
pElems^[u] := frInsElems[cnt];
dec(cnt);
end
else
Begin
pElems^[u] := pElems^[i];
dec(i);
end;
dec(u);
until (i<0) or (cnt<0);
IF i < 0 then
For u := cnt downto 0 do
pElems^[u] := frInsElems[u];
end;
end;
procedure LoEAppendNext(pFr:tpFaktorRec;Limit:tnumber);
var
pNum : tpNumber;
pElems : tpElemArr;
i : NativeInt;
begin
with pFr^ do
Begin
i := frMaxIdx+1;
pElems := @frElems[0];
pNum := LoEGetActLstNumber(frNextFr);
while Limit > pNum^ do
Begin
IF i > High(frElems) then
Begin
LoEIncreaseElems(pFr,10);
pElems := @frElems[0];
end;
pElems^[i] := LoEGetNextElement(frNextFr);
inc(i);
pNum := LoEGetActLstNumber(frNextFr);
end;
inc(ActIndex,i);
frMaxIdx:= i-1;
end;
end;
procedure LoENextList(pFr:tpFaktorRec);
var
pElems : tpElemArr;
j,PotNum : NativeInt;
lnPrime : tnumber;
begin
with pFR^ do
Begin
//increase all Elements by factor
pElems := @frElems[0];
LnPrime := frLnPrime;
PotNum := frPotNo;
for j := frMaxIdx Downto 0 do
with pElems^[j] do
Begin
n := LnPrime+n;
inc(Pots[PotNum]);
end;
//x^j -> x^(j+1)
j := frActPot+1;
with frActNumb do
begin
n:= j*LnPrime;
Pots[PotNum]:= j;
end;
frActPot := j;
//if something follows
IF frNextFr <> NIL then
LoEInsertNext(pFR,frActNumb.n);
frAktIdx := 0;
end;
end;
function LoEGetNextElementPointer(pFR :tpFaktorRec):tpElem;
Begin
with pFr^ do
Begin
IF frMaxIdx < frAktIdx then
LoENextList(pFr);
result := @frElems[frAktIdx];
inc(frAktIdx);
inc(ActIndex);
end;
end;
function LoEGetNextElement(pFR :tpFaktorRec):tElem;
Begin
with pFr^ do
Begin
result := frElems[frAktIdx];
inc(frAktIdx);
IF frMaxIdx < frAktIdx then
LoENextList(pFr);
inc(ActIndex);
end;
end;
function LoEGetNextNumber(pFR :tpFaktorRec):tNumber;
Begin
with pFr^ do
Begin
result := frElems[frAktIdx].n;
inc(frAktIdx);
IF frMaxIdx < frAktIdx then
LoENextList(pFr);
inc(ActIndex);
end;
end;
procedure LoEGetNumber(pFR :tpFaktorRec;no:NativeUint);
Begin
dec(no);
while ActIndex < no do
LoENextList(pFR);
with pFr^ do
frAktIdx := (no-(ActIndex-frMaxIdx)-1);
end;
procedure first50;
var
FA: tArrFR;
i : integer;
Begin
LoECreate(Pl,FA);
write(1,' ');
For i := 1 to 49 do
AusgabeElem(LoEGetNextElement(@FA[0]));
writeln;
LoEFree(FA);
end;
procedure GetDigitCounts(MaxDigit:Uint32);
var
T1,T0 : TDateTime;
FA: tArrFR;
i,j,LastCnt : NativeUint;
Begin
T0 := now;
inc(MaxDigit);
LoECreate(Pl,FA);
i := 1;
j := 0;
writeln('Digits count total count ');
repeat
LastCnt := j;
repeat
inc(j);
with LoEGetNextElementPointer(@FA[0])^ do
IF (Pots[2]= i) AND (Pots[0]= i) then
break;
until false;
writeln(i:4,j-LastCnt:12,j:15,(now-T0)*86.6e3:10:3,' s');
LastCnt := j;
inc(i);
until i = MaxDigit;
LoEFree(FA);
T1 := now;
writeln('Total number of humble numbers found: ',j);
writeln('Running time: ',(T1-T0)*86.6e6:0:0,' ms');
end;
Begin
//check if PlInit(4); <= maxPrimFakCnt+1
PlInit(4);// 3 -> 2,3,5/ 4 -> 2,3,5,7
first50;
GetDigitCounts(100);
End.
- Output:
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 Digits count total count 1 9 9 0.000 s 2 36 45 0.000 s 3 95 140 0.000 s 4 197 337 0.000 s 5 356 693 0.001 s 6 579 1272 0.001 s 7 882 2154 0.001 s 8 1272 3426 0.001 s 9 1767 5193 0.001 s 10 2381 7574 0.001 s 11 3113 10687 0.001 s 12 3984 14671 0.001 s 13 5002 19673 0.001 s 14 6187 25860 0.002 s 15 7545 33405 0.002 s 16 9081 42486 0.003 s 17 10815 53301 0.003 s 18 12759 66060 0.004 s 19 14927 80987 0.004 s 20 17323 98310 0.005 s 21 19960 118270 0.006 s 22 22853 141123 0.007 s 23 26015 167138 0.008 s 24 29458 196596 0.009 s 25 33188 229784 0.009 s 26 37222 267006 0.010 s 27 41568 308574 0.011 s 28 46245 354819 0.011 s 29 51254 406073 0.012 s 30 56618 462691 0.013 s 31 62338 525029 0.014 s 32 68437 593466 0.015 s 33 74917 668383 0.016 s 34 81793 750176 0.018 s 35 89083 839259 0.019 s 36 96786 936045 0.021 s 37 104926 1040971 0.022 s 38 113511 1154482 0.025 s 39 122546 1277028 0.027 s 40 132054 1409082 0.028 s 41 142038 1551120 0.031 s 42 152515 1703635 0.033 s 43 163497 1867132 0.036 s 44 174986 2042118 0.039 s 45 187004 2229122 0.042 s 46 199565 2428687 0.045 s 47 212675 2641362 0.050 s 48 226346 2867708 0.053 s 49 240590 3108298 0.056 s 50 255415 3363713 0.061 s 51 270843 3634556 0.065 s 52 286880 3921436 0.070 s 53 303533 4224969 0.076 s 54 320821 4545790 0.081 s 55 338750 4884540 0.087 s 56 357343 5241883 0.094 s 57 376599 5618482 0.100 s 58 396533 6015015 0.106 s 59 417160 6432175 0.113 s 60 438492 6870667 0.122 s 61 460533 7331200 0.130 s 62 483307 7814507 0.137 s 63 506820 8321327 0.148 s 64 531076 8852403 0.156 s 65 556104 9408507 0.165 s 66 581902 9990409 0.177 s 67 608483 10598892 0.186 s 68 635864 11234756 0.197 s 69 664053 11898809 0.210 s 70 693065 12591874 0.222 s 71 722911 13314785 0.235 s 72 753593 14068378 0.250 s 73 785141 14853519 0.262 s 74 817554 15671073 0.275 s 75 850847 16521920 0.292 s 76 885037 17406957 0.305 s 77 920120 18327077 0.320 s 78 956120 19283197 0.338 s 79 993058 20276255 0.354 s 80 1030928 21307183 0.370 s 81 1069748 22376931 0.391 s 82 1109528 23486459 0.409 s 83 1150287 24636746 0.428 s 84 1192035 25828781 0.451 s 85 1234774 27063555 0.470 s 86 1278527 28342082 0.490 s 87 1323301 29665383 0.516 s 88 1369106 31034489 0.538 s 89 1415956 32450445 0.559 s 90 1463862 33914307 0.582 s 91 1512840 35427147 0.612 s 92 1562897 36990044 0.636 s 93 1614050 38604094 0.662 s 94 1666302 40270396 0.694 s 95 1719669 41990065 0.721 s 96 1774166 43764231 0.748 s 97 1829805 45594036 0.785 s 98 1886590 47480626 0.815 s 99 1944540 49425166 0.845 s 100 2003661 51428827 0.884 s 101 2063972 53492799 0.916 s 102 2125486 55618285 0.948 s 103 2188204 57806489 0.991 s 104 2252146 60058635 1.026 s 105 2317319 62375954 1.062 s 106 2383733 64759687 1.110 s 107 2451413 67211100 1.147 s 108 2520360 69731460 1.186 s 109 2590584 72322044 1.237 s 110 2662102 74984146 1.278 s .. shortened 120 3450981 105831533 1.795 s 130 4382079 145340170 2.450 s 140 5467187 194997130 3.298 s 150 6718091 256407329 4.318 s 170 9764417 421496245 7.055 s 180 11583420 528974086 8.822 s 190 13615367 655803459 10.972 s 200 15872045 804178604 13.424 s 300 53391941 4039954757 66.882 s 400 126350163 12719142480 207.622 s 500 246533493 30980806733 503.269 s 600 425728730 64142692268 1039.928 s 700 675722681 118701223590 1920.325 s 800 1008302151 202331504969 3265.469 s .. 876 1323548095 290737888948 4690.230 s 877 1328082553 292065971501 4709.931 s
Perl
use strict;
use warnings;
use List::Util 'min';
#use bigint # works, but slow
use Math::GMPz; # this module gives roughly 16x speed-up
sub humble_gen {
my @s = ([1], [1], [1], [1]);
my @m = (2, 3, 5, 7);
@m = map { Math::GMPz->new($_) } @m; # comment out to NOT use Math::GMPz
return sub {
my $n = min $s[0][0], $s[1][0], $s[2][0], $s[3][0];
for (0..3) {
shift @{$s[$_]} if $s[$_][0] == $n;
push @{$s[$_]}, $n * $m[$_]
}
return $n
}
}
my $h = humble_gen;
my $i = 0;
my $upto = 50;
my $list;
++$i, $list .= $h->(). " " until $i == $upto;
print "$list\n";
$h = humble_gen; # from the top...
my $count = 0;
my $digits = 1;
while ($digits <= $upto) {
++$count and next if $digits == length $h->();
printf "Digits: %2d - Count: %s\n", $digits++, $count;
$count = 1;
}
- Output:
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 Digits: 1 - Count: 9 Digits: 2 - Count: 36 Digits: 3 - Count: 95 Digits: 4 - Count: 197 Digits: 5 - Count: 356 Digits: 6 - Count: 579 Digits: 7 - Count: 882 Digits: 8 - Count: 1272 Digits: 9 - Count: 1767 Digits: 10 - Count: 2381 Digits: 11 - Count: 3113 Digits: 12 - Count: 3984 Digits: 13 - Count: 5002 Digits: 14 - Count: 6187 Digits: 15 - Count: 7545 Digits: 16 - Count: 9081 Digits: 17 - Count: 10815 Digits: 18 - Count: 12759 Digits: 19 - Count: 14927 Digits: 20 - Count: 17323 Digits: 21 - Count: 19960 Digits: 22 - Count: 22853 Digits: 23 - Count: 26015 Digits: 24 - Count: 29458 Digits: 25 - Count: 33188 Digits: 26 - Count: 37222 Digits: 27 - Count: 41568 Digits: 28 - Count: 46245 Digits: 29 - Count: 51254 Digits: 30 - Count: 56618 Digits: 31 - Count: 62338 Digits: 32 - Count: 68437 Digits: 33 - Count: 74917 Digits: 34 - Count: 81793 Digits: 35 - Count: 89083 Digits: 36 - Count: 96786 Digits: 37 - Count: 104926 Digits: 38 - Count: 113511 Digits: 39 - Count: 122546 Digits: 40 - Count: 132054 Digits: 41 - Count: 142038 Digits: 42 - Count: 152515 Digits: 43 - Count: 163497 Digits: 44 - Count: 174986 Digits: 45 - Count: 187004 Digits: 46 - Count: 199565 Digits: 47 - Count: 212675 Digits: 48 - Count: 226346 Digits: 49 - Count: 240590 Digits: 50 - Count: 255415
Phix
I felt pretty good about the performance of this, until I ran the Go version - humbled indeed!
It will go all the way to 100 digits if you give it time (18 mins, on 64bit - 32bit runs out of memory after printing the 99th line)
I also tried a log version (similar to Hamming_numbers) but inaccuracies with floor(h[n][LOG]) crept in quite early, at just 10 digits.
-- demo/rosetta/humble.exw with javascript_semantics include mpfr.e procedure humble(integer n, bool countdigits=false) -- if countdigits is false: show first n humble numbers, -- if countdigits is true: count them up to n digits. sequence humble = {mpz_init(1)}, nexts = {2,3,5,7}, indices = repeat(1,4) for i=1 to 4 do nexts[i] = mpz_init(nexts[i]) end for integer digits = 1, count = 1, dead = 1, tc = 0 atom t0 = time() mpz p10 = mpz_init(10) while ((not countdigits) and length(humble)<n) or (countdigits and digits<=n) do mpz x = mpz_init_set(mpz_min(nexts)) humble = append(humble,x) if countdigits then if mpz_cmp(x,p10)>=0 then mpz_mul_si(p10,p10,10) integer d = min(indices) for k=dead to d-1 do humble[k] = mpz_free(humble[k]) end for dead = d string s = iff(digits=1?"":"s"), e = elapsed(time()-t0) tc += count -- e &= sprintf(", %,d dead",{dead-1}) e &= sprintf(", total:%,d",{tc}) printf(1,"%,12d humble numbers have %d digit%s (%s)\n",{count,digits,s,e}) digits += 1 count = 1 else count += 1 end if end if for j=1 to 4 do if mpz_cmp(nexts[j],x)<=0 then indices[j] += 1 mpz_mul_si(nexts[j],humble[indices[j]],get_prime(j)) end if end for end while if not countdigits then for i=1 to length(humble) do humble[i] = shorten(mpz_get_str(humble[i]),ml:=10) end for printf(1,"First %d humble numbers: %s\n\n",{n,join(humble," ")}) end if end procedure humble(50) humble(42,true)
- Output:
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
9 humble numbers have 1 digit (0s, total:9)
36 humble numbers have 2 digits (0s, total:45)
95 humble numbers have 3 digits (0s, total:140)
197 humble numbers have 4 digits (0s, total:337)
356 humble numbers have 5 digits (0s, total:693)
579 humble numbers have 6 digits (0.0s, total:1,272)
882 humble numbers have 7 digits (0.0s, total:2,154)
1,272 humble numbers have 8 digits (0.0s, total:3,426)
1,767 humble numbers have 9 digits (0.1s, total:5,193)
2,381 humble numbers have 10 digits (0.1s, total:7,574)
3,113 humble numbers have 11 digits (0.2s, total:10,687)
3,984 humble numbers have 12 digits (0.2s, total:14,671)
5,002 humble numbers have 13 digits (0.3s, total:19,673)
6,187 humble numbers have 14 digits (0.4s, total:25,860)
7,545 humble numbers have 15 digits (0.5s, total:33,405)
9,081 humble numbers have 16 digits (0.6s, total:42,486)
10,815 humble numbers have 17 digits (0.8s, total:53,301)
12,759 humble numbers have 18 digits (0.9s, total:66,060)
14,927 humble numbers have 19 digits (1.2s, total:80,987)
17,323 humble numbers have 20 digits (1.4s, total:98,310)
19,960 humble numbers have 21 digits (1.7s, total:118,270)
22,853 humble numbers have 22 digits (2.0s, total:141,123)
26,015 humble numbers have 23 digits (2.3s, total:167,138)
29,458 humble numbers have 24 digits (2.7s, total:196,596)
33,188 humble numbers have 25 digits (3.2s, total:229,784)
37,222 humble numbers have 26 digits (3.7s, total:267,006)
41,568 humble numbers have 27 digits (4.3s, total:308,574)
46,245 humble numbers have 28 digits (5.0s, total:354,819)
51,254 humble numbers have 29 digits (5.7s, total:406,073)
56,618 humble numbers have 30 digits (6.5s, total:462,691)
62,338 humble numbers have 31 digits (7.3s, total:525,029)
68,437 humble numbers have 32 digits (8.3s, total:593,466)
74,917 humble numbers have 33 digits (9.4s, total:668,383)
81,793 humble numbers have 34 digits (10.5s, total:750,176)
89,083 humble numbers have 35 digits (11.7s, total:839,259)
96,786 humble numbers have 36 digits (13.1s, total:936,045)
104,926 humble numbers have 37 digits (14.6s, total:1,040,971)
113,511 humble numbers have 38 digits (16.2s, total:1,154,482)
122,546 humble numbers have 39 digits (17.9s, total:1,277,028)
132,054 humble numbers have 40 digits (19.7s, total:1,409,082)
142,038 humble numbers have 41 digits (21.7s, total:1,551,120)
152,515 humble numbers have 42 digits (23.9s, total:1,703,635)
PL/I
PL/M
This can be compiled with the original 8080 PL/M compiler and run under CP/M or an emulator or clone.Only handles Humble numbers with up to 4 digits as 8080 PL/M only has unsigned 8 and 16 bit integers.
100H: /* FIND SOME HUMBLE NUMBERS - NUMBERS WITH NO PRIME FACTORS ABOVE 7 */
BDOS: PROCEDURE( FN, ARG ); /* CP/M BDOS SYSTEM CALL */
DECLARE FN BYTE, ARG ADDRESS;
GOTO 5;
END BDOS;
PRINT$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PRINT$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PRINT$NL: PROCEDURE; CALL PRINT$STRING( .( 0DH, 0AH, '$' ) ); END;
PRINT$NUMBER: PROCEDURE( N );
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR( 6 ) BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PRINT$STRING( .N$STR( W ) );
END PRINT$NUMBER;
MIN: PROCEDURE( A, B ) ADDRESS;
DECLARE ( A, B ) ADDRESS;
IF A < B THEN RETURN( A ); ELSE RETURN( B );
END MIN;
/* DISPLAY A STATISTIC ABOUT HUMBLE NUMBERS */
PRINT$H$STAT: PROCEDURE( S, D );
DECLARE ( S, D ) ADDRESS;
CALL PRINT$STRING( .'THERE ARE $' );
IF S < 10 THEN CALL PRINT$CHAR( ' ' );
IF S < 100 THEN CALL PRINT$CHAR( ' ' );
CALL PRINT$NUMBER( S );
CALL PRINT$STRING( .' HUMBLE NUMBERS WITH $' );
CALL PRINT$NUMBER( D );
CALL PRINT$STRING( .' DIGIT$' );
IF D > 1 THEN CALL PRINT$CHAR( 'S' );
CALL PRINT$NL;
END PRINT$H$STAT;
/* FIND AND PRINT HUMBLE NUMBERS */
DECLARE MAX$HUMBLE LITERALLY '400';
DECLARE H( MAX$HUMBLE ) ADDRESS;
DECLARE ( P2, P3, P5, P7, M
, LAST2, LAST3, LAST5, LAST7
, H1, H2, H3, H4, H5, H6, HPOS
) ADDRESS;
/* 1 IS THE FIRST HUMBLE NUMBER */
H( 0 ) = 1;
H1 = 0; H2 = 0; H3 = 0; H4 = 0; H5 = 0; H6 = 0;
LAST2 = 0; LAST3 = 0; LAST5 = 0; LAST7 = 0;
P2 = 2; P3 = 3; P5 = 5; P7 = 7;
DO HPOS = 1 TO MAX$HUMBLE - 1;
/* THE NEXT HUMBLE NUMBER IS THE LOWEST OF THE NEXT MULTIPLES OF */
/* 2, 3, 5, 7 */
M = MIN( MIN( MIN( P2, P3 ), P5 ), P7 );
H( HPOS ) = M;
IF M = P2 THEN P2 = 2 * H( LAST2 := LAST2 + 1 );
IF M = P3 THEN P3 = 3 * H( LAST3 := LAST3 + 1 );
IF M = P5 THEN P5 = 5 * H( LAST5 := LAST5 + 1 );
IF M = P7 THEN P7 = 7 * H( LAST7 := LAST7 + 1 );
END;
DO HPOS = 0 TO 49;
CALL PRINT$CHAR( ' ' );
CALL PRINT$NUMBER( H( HPOS ) );
END;
CALL PRINT$NL;
DO HPOS = 0 TO MAX$HUMBLE - 1;
M = H( HPOS );
IF M < 10 THEN H1 = H1 + 1;
ELSE IF M < 100 THEN H2 = H2 + 1;
ELSE IF M < 1000 THEN H3 = H3 + 1;
ELSE IF M < 10000 THEN H4 = H4 + 1;
END;
CALL PRINT$H$STAT( H1, 1 );
CALL PRINT$H$STAT( H2, 2 );
CALL PRINT$H$STAT( H3, 3 );
CALL PRINT$H$STAT( H4, 4 );
EOF- Output:
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 THERE ARE 9 HUMBLE NUMBERS WITH 1 DIGIT THERE ARE 36 HUMBLE NUMBERS WITH 2 DIGITS THERE ARE 95 HUMBLE NUMBERS WITH 3 DIGITS THERE ARE 197 HUMBLE NUMBERS WITH 4 DIGITS
See also #Polyglot:PL/I and PL/M
Polyglot:PL/I and PL/M
... under CP/M (or an emulator)Should work with many PL/I implementations.
The PL/I include file "pg.inc" can be found on the Polyglot:PL/I and PL/M page.
Note the use of text in column 81 onwards to hide the PL/I specifics from the PL/M compiler.
Based on the PL/M version, note PL/I does not have the "walrus operator" (:=) which allows assignments to be nested in expressions, so it can't be used in the non-PL/M specific parts of this.
/* FIND SOME HUMBLE NUMBERS - NUMBERS WITH NO PRIME FACTORS ABOVE 7 */
humble_100H: procedure options (main);
/* PROGRAM-SPECIFIC %REPLACE STATEMENTS MUST APPEAR BEFORE THE %INCLUDE AS */
/* E.G. THE CP/M PL/I COMPILER DOESN'T LIKE THEM TO FOLLOW PROCEDURES */
/* PL/I */
%replace dclhumble by 400;
/* PL/M */ /*
DECLARE DCLHUMBLE LITERALLY '401';
/* */
/* PL/I DEFINITIONS */
%include 'pg.inc';
/* PL/M DEFINITIONS: CP/M BDOS SYSTEM CALL AND CONSOLE I/O ROUTINES, ETC. */ /*
DECLARE BINARY LITERALLY 'ADDRESS', CHARACTER LITERALLY 'BYTE';
DECLARE FIXED LITERALLY ' ', BIT LITERALLY 'BYTE';
DECLARE STATIC LITERALLY ' ', RETURNS LITERALLY ' ';
DECLARE FALSE LITERALLY '0', TRUE LITERALLY '1';
DECLARE HBOUND LITERALLY 'LAST', SADDR LITERALLY '.';
BDOSF: PROCEDURE( FN, ARG )BYTE;
DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PRCHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PRSTRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PRNL: PROCEDURE; CALL PRCHAR( 0DH ); CALL PRCHAR( 0AH ); END;
PRNUMBER: PROCEDURE( N );
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR( 6 ) BYTE, W BYTE;
N$STR( W := LAST( N$STR ) ) = '$';
N$STR( W := W - 1 ) = '0' + ( ( V := N ) MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL BDOS( 9, .N$STR( W ) );
END PRNUMBER;
MODF: PROCEDURE( A, B )ADDRESS;
DECLARE ( A, B )ADDRESS;
RETURN( A MOD B );
END MODF;
MIN: PROCEDURE( A, B ) ADDRESS;
DECLARE ( A, B ) ADDRESS;
IF A < B THEN RETURN( A ); ELSE RETURN( B );
END MIN;
/* END LANGUAGE DEFINITIONS */
/* TASK */
/* DISPLAY A STATISTIC ABOUT HUMBLE NUMBERS */
PRHUMBLESTAT: PROCEDURE( S, D );
DECLARE ( S, D ) FIXED BINARY;
CALL PRSTRING( SADDR( 'THERE ARE $' ) );
IF S < 10 THEN CALL PRCHAR( ' ' );
IF S < 100 THEN CALL PRCHAR( ' ' );
CALL PRNUMBER( S );
CALL PRSTRING( SADDR( ' HUMBLE NUMBERS WITH $' ) );
CALL PRNUMBER( D );
CALL PRSTRING( SADDR( ' DIGIT$' ) );
IF D > 1 THEN CALL PRCHAR( 'S' );
CALL PRNL;
END PRHUMBLESTAT;
/* FIND AND PRINT HUMBLE NUMBERS */
DECLARE H( DCLHUMBLE ) FIXED BINARY;
DECLARE ( P2, P3, P5, P7, M
, LAST2, LAST3, LAST5, LAST7
, H1, H2, H3, H4, HPOS
, MAXHUMBLE
) FIXED BINARY;
MAXHUMBLE = HBOUND( H ,1
);
/* 1 IS THE FIRST HUMBLE NUMBER */
H( 1 ) = 1;
LAST2 = 1; LAST3 = 1; LAST5 = 1; LAST7 = 1;
P2 = 2; P3 = 3; P5 = 5; P7 = 7;
DO HPOS = 2 TO MAXHUMBLE;
/* THE NEXT HUMBLE NUMBER IS THE LOWEST OF THE NEXT MULTIPLES OF */
/* 2, 3, 5, 7 */
M = MIN( MIN( MIN( P2, P3 ), P5 ), P7 );
H( HPOS ) = M;
IF M = P2 THEN DO; LAST2 = LAST2 + 1; P2 = 2 * H( LAST2 ); END;
IF M = P3 THEN DO; LAST3 = LAST3 + 1; P3 = 3 * H( LAST3 ); END;
IF M = P5 THEN DO; LAST5 = LAST5 + 1; P5 = 5 * H( LAST5 ); END;
IF M = P7 THEN DO; LAST7 = LAST7 + 1; P7 = 7 * H( LAST7 ); END;
END;
/* SHOW THE FIRST 50 HUMBLE NUMBERS */
DO HPOS = 1 TO 50;
CALL PRCHAR( ' ' );
M = H( HPOS );
IF M < 10 THEN CALL PRCHAR( ' ' );
IF M < 100 THEN CALL PRCHAR( ' ' );
CALL PRNUMBER( H( HPOS ) );
IF MODF( HPOS, 16 ) = 0 THEN CALL PRNL;
END;
CALL PRNL;
/* SHOW THE NUMBER OF HUMBLE NUMBERS UP TO VARIOUS POWERS OF 10 */
H1 = 0; H2 = 0; H3 = 0; H4 = 0;
DO HPOS = 1 TO MAXHUMBLE;
M = H( HPOS );
IF M < 10 THEN H1 = H1 + 1;
ELSE IF M < 100 THEN H2 = H2 + 1;
ELSE IF M < 1000 THEN H3 = H3 + 1;
ELSE IF M < 10000 THEN H4 = H4 + 1;
END;
CALL PRHUMBLESTAT( H1, 1 );
CALL PRHUMBLESTAT( H2, 2 );
CALL PRHUMBLESTAT( H3, 3 );
CALL PRHUMBLESTAT( H4, 4 );
EOF: end humble_100H;- Output:
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 THERE ARE 9 HUMBLE NUMBERS WITH 1 DIGIT THERE ARE 36 HUMBLE NUMBERS WITH 2 DIGITS THERE ARE 95 HUMBLE NUMBERS WITH 3 DIGITS THERE ARE 197 HUMBLE NUMBERS WITH 4 DIGITS
Python
'''Humble numbers'''
from itertools import groupby, islice
from functools import reduce
# humbles :: () -> [Int]
def humbles():
'''A non-finite stream of Humble numbers.
OEIS A002473
'''
hs = set([1])
while True:
nxt = min(hs)
yield nxt
hs.remove(nxt)
hs.update(nxt * x for x in [2, 3, 5, 7])
# TEST ----------------------------------------------------
# main :: IO ()
def main():
'''First 50, and counts with N digits'''
print('First 50 Humble numbers:\n')
for row in chunksOf(10)(
take(50)(humbles())
):
print(' '.join(map(
lambda x: str(x).rjust(3),
row
)))
print('\nCounts of Humble numbers with n digits:\n')
for tpl in take(10)(
(k, len(list(g))) for k, g in
groupby(len(str(x)) for x in humbles())
):
print(tpl)
# GENERIC -------------------------------------------------
# chunksOf :: Int -> [a] -> [[a]]
def chunksOf(n):
'''A series of lists of length n, subdividing the
contents of xs. Where the length of xs is not evenly
divible, the final list will be shorter than n.
'''
return lambda xs: reduce(
lambda a, i: a + [xs[i:n + i]],
range(0, len(xs), n), []
) if 0 < n else []
# take :: Int -> [a] -> [a]
# take :: Int -> String -> String
def take(n):
'''The prefix of xs of length n,
or xs itself if n > length xs.
'''
return lambda xs: (
list(islice(xs, n))
)
# MAIN ---
if __name__ == '__main__':
main()
- Output:
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 Counts of Humble numbers with n digits: (1, 9) (2, 36) (3, 95) (4, 197) (5, 356) (6, 579) (7, 882) (8, 1272) (9, 1767) (10, 2381)
Quackery
Uses smoothwith from N-smooth numbers#Quackery, and searchwith from Binary search#Quackery.
' [ 2 3 5 7 ] smoothwith
[ -1 peek [ 10 12 ** ] constant = ]
-1 split drop
dup 50 split drop echo cr cr
12 times
[ dup 0 over size
rot 10 i^ 1+ ** searchwith <
drop split swap size echo
sp i^ 1+ echo
say "-digit humble numbers" cr ]
drop- Output:
[ 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 ] 9 1-digit humble numbers 36 2-digit humble numbers 95 3-digit humble numbers 197 4-digit humble numbers 356 5-digit humble numbers 579 6-digit humble numbers 882 7-digit humble numbers 1272 8-digit humble numbers 1767 9-digit humble numbers 2381 10-digit humble numbers 3113 11-digit humble numbers 3984 12-digit humble numbers
Racket
#lang racket
(define (gen-humble-numbers N (kons #f) (k0 (void)))
(define rv (make-vector N 1))
(define (loop n 2-idx 3-idx 5-idx 7-idx next-2 next-3 next-5 next-7 k)
(if (= n N)
rv
(let ((mn (min next-2 next-3 next-5 next-7)))
(vector-set! rv n mn)
(define (add-1-if-min n x) (if (= mn n) (add1 x) x))
(define (*vr.i-if-min n m i) (if (= mn n) (* m (vector-ref rv i)) n))
(let* ((2-idx (add-1-if-min next-2 2-idx))
(next-2 (*vr.i-if-min next-2 2 2-idx))
(3-idx (add-1-if-min next-3 3-idx))
(next-3 (*vr.i-if-min next-3 3 3-idx))
(5-idx (add-1-if-min next-5 5-idx))
(next-5 (*vr.i-if-min next-5 5 5-idx))
(7-idx (add-1-if-min next-7 7-idx))
(next-7 (*vr.i-if-min next-7 7 7-idx))
(k (and kons (kons mn k))))
(loop (add1 n) 2-idx 3-idx 5-idx 7-idx next-2 next-3 next-5 next-7 k)))))
(loop 1 0 0 0 0 2 3 5 7 (and kons (kons 1 k0))))
(define ((digit-tracker breaker) h last-ten.count)
(let ((last-ten (car last-ten.count)))
(if (< h last-ten)
(cons last-ten (add1 (cdr last-ten.count)))
(begin
(printf "~a humble numbers with ~a digits~%" (cdr last-ten.count) (order-of-magnitude last-ten))
(cons (breaker (* 10 last-ten)) 1)))))
(define (Humble-numbers)
(displayln (gen-humble-numbers 50))
(time
(let/ec break
(void (gen-humble-numbers
100000000
(digit-tracker (λ (o) (if (> (order-of-magnitude o) 100) (break) o)))
'(10 . 0))))))
(module+ main
(Humble-numbers))
- Output:
output has been elided manually, to avoid repetition with the numbers you've already seen elsewhere:
#(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) 9 humble numbers with 1 digits 36 humble numbers with 2 digits 95 humble numbers with 3 digits 197 humble numbers with 4 digits 356 humble numbers with 5 digits 579 humble numbers with 6 digits 882 humble numbers with 7 digits 1272 humble numbers with 8 digits 1767 humble numbers with 9 digits 2381 humble numbers with 10 digits ... 17323 humble numbers with 20 digits ... 56618 humble numbers with 30 digits ... 132054 humble numbers with 40 digits ... 255415 humble numbers with 50 digits ... 438492 humble numbers with 60 digits ... 693065 humble numbers with 70 digits ... 1030928 humble numbers with 80 digits ... 1463862 humble numbers with 90 digits ... 2003661 humble numbers with 100 digits cpu time: 234970 real time: 235489 gc time: 189187
Raku
(formerly Perl 6)
sub smooth-numbers (*@list) {
cache my \Smooth := gather {
my %i = (flat @list) Z=> (Smooth.iterator for ^@list);
my %n = (flat @list) Z=> 1 xx *;
loop {
take my $n := %n{*}.min;
for @list -> \k {
%n{k} = %i{k}.pull-one * k if %n{k} == $n;
}
}
}
}
my $humble := smooth-numbers(2,3,5,7);
put $humble[^50];
say '';
my $upto = 50;
my $digits = 1;
my $count;
$humble.map: -> \h {
++$count and next if h.chars == $digits;
printf "Digits: %2d - Count: %s\n", $digits++, $count;
$count = 1;
last if $digits > $upto;
}
- Output:
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 Digits: 1 - Count: 9 Digits: 2 - Count: 36 Digits: 3 - Count: 95 Digits: 4 - Count: 197 Digits: 5 - Count: 356 Digits: 6 - Count: 579 Digits: 7 - Count: 882 Digits: 8 - Count: 1272 Digits: 9 - Count: 1767 Digits: 10 - Count: 2381 Digits: 11 - Count: 3113 Digits: 12 - Count: 3984 Digits: 13 - Count: 5002 Digits: 14 - Count: 6187 Digits: 15 - Count: 7545 Digits: 16 - Count: 9081 Digits: 17 - Count: 10815 Digits: 18 - Count: 12759 Digits: 19 - Count: 14927 Digits: 20 - Count: 17323 Digits: 21 - Count: 19960 Digits: 22 - Count: 22853 Digits: 23 - Count: 26015 Digits: 24 - Count: 29458 Digits: 25 - Count: 33188 Digits: 26 - Count: 37222 Digits: 27 - Count: 41568 Digits: 28 - Count: 46245 Digits: 29 - Count: 51254 Digits: 30 - Count: 56618 Digits: 31 - Count: 62338 Digits: 32 - Count: 68437 Digits: 33 - Count: 74917 Digits: 34 - Count: 81793 Digits: 35 - Count: 89083 Digits: 36 - Count: 96786 Digits: 37 - Count: 104926 Digits: 38 - Count: 113511 Digits: 39 - Count: 122546 Digits: 40 - Count: 132054 Digits: 41 - Count: 142038 Digits: 42 - Count: 152515 Digits: 43 - Count: 163497 Digits: 44 - Count: 174986 Digits: 45 - Count: 187004 Digits: 46 - Count: 199565 Digits: 47 - Count: 212675 Digits: 48 - Count: 226346 Digits: 49 - Count: 240590 Digits: 50 - Count: 255415
REXX
/*REXX program computes and displays humble numbers, also will display counts of sizes.*/
parse arg n m . /*obtain optional arguments from the CL*/
if n=='' | n=="," then n= 50 /*Not specified? Then use the default.*/
if m=='' | m=="," then m= 60 /* " " " " " " */
numeric digits 1 + max(20, m) /*be able to handle some big numbers. */
$.= 0 /*a count array for X digit humble #s*/
call humble n; list= /*call HUMBLE sub; initialize the list.*/
do j=1 for n; list= list @.j /*append a humble number to the list.*/
end /*j*/
if list\='' then do; say "A list of the first " n ' humble numbers are:'
say strip(list) /*elide the leading blank in the list. */
end
say
call humble -m /*invoke subroutine for counting nums. */
if $.1==0 then exit /*if no counts, then we're all finished*/
total= 0 /*initialize count of humble numbers. */
$.1= $.1 + 1 /*adjust count for absent 1st humble #.*/
say ' The digit counts of humble numbers:'
say ' ═════════════════════════════════════════'
do c=1 while $.c>0; s= left('s', length($.c)>1) /*count needs pluralization?*/
say right( commas($.c), 30) ' have ' right(c, 2) " digit"s
total= total + $.c /* ◄─────────────────────────────────┐ */
end /*k*/ /*bump humble number count (so far)──┘ */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: procedure; arg _; do i=length(_)-3 to 1 by -3; _=insert(',', _, i); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
humble: procedure expose @. $.; parse arg x; if x==0 then return
y= abs(x); a= y; noCount= x>0; if x<0 then y= 999999999
#2= 1; #3= 1; #5= 1; #7= 1 /*define the initial humble constants. */
$.= 0; @.= 0; @.1= 1 /*initialize counts and humble numbers.*/
do h=2 for y-1
@.h= min(2*@.#2,3*@.#3,5*@.#5,7*@.#7) /*pick the minimum of 4 humble numbers.*/
m= @.h /*M: " " " " " " */
if 2*@.#2 == m then #2 = #2 + 1 /*Is number already defined? Use next #*/
if 3*@.#3 == m then #3 = #3 + 1 /* " " " " " " "*/
if 5*@.#5 == m then #5 = #5 + 1 /* " " " " " " "*/
if 7*@.#7 == m then #7 = #7 + 1 /* " " " " " " "*/
if noCount then iterate /*Not counting digits? Then iterate. */
L= length(m); if L>a then leave /*Are we done with counting? Then quit*/
$.L= $.L + 1 /*bump the digit count for this number.*/
end /*h*/ /*the humble numbers are in the @ array*/
return /* " count results " " " $ " */
- output when using the default inputs:
(Shown at 7/8 size.)
A list of the first 50 humble numbers are:
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
The digit counts of humble numbers:
═════════════════════════════════════════
9 have 1 digit
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1,272 have 8 digits
1,767 have 9 digits
2,381 have 10 digits
3,113 have 11 digits
3,984 have 12 digits
5,002 have 13 digits
6,187 have 14 digits
7,545 have 15 digits
9,081 have 16 digits
10,815 have 17 digits
12,759 have 18 digits
14,927 have 19 digits
17,323 have 20 digits
19,960 have 21 digits
22,853 have 22 digits
26,015 have 23 digits
29,458 have 24 digits
33,188 have 25 digits
37,222 have 26 digits
41,568 have 27 digits
46,245 have 28 digits
51,254 have 29 digits
56,618 have 30 digits
62,338 have 31 digits
68,437 have 32 digits
74,917 have 33 digits
81,793 have 34 digits
89,083 have 35 digits
96,786 have 36 digits
104,926 have 37 digits
113,511 have 38 digits
122,546 have 39 digits
132,054 have 40 digits
142,038 have 41 digits
152,515 have 42 digits
163,497 have 43 digits
174,986 have 44 digits
187,004 have 45 digits
199,565 have 46 digits
212,675 have 47 digits
226,346 have 48 digits
240,590 have 49 digits
255,415 have 50 digits
270,843 have 51 digits
286,880 have 52 digits
303,533 have 53 digits
320,821 have 54 digits
338,750 have 55 digits
357,343 have 56 digits
376,599 have 57 digits
396,533 have 58 digits
417,160 have 59 digits
438,492 have 60 digits
total number of humble numbers found: 6,870,667
Ring
load "stdlib.ring"
limit = 10
numList = []
for n2 = 0 to limit
for n3 = 0 to limit
for n5 = 0 to limit
for n7 = 0 to limit
num = pow(2,n2) * pow(3,n3) * pow(5,n5) * pow(7,n7)
add(numList,num)
next
next
next
next
numList = sort(numList)
see "The first 50 Humble numbers: " + nl
for n = 1 to 50
see "" + numList[n] + " "
nextOutput:
The 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
RPL
Brute force is summoned to generate the list of the first humble numbers, but direct humble numbers generation and immediate digit counting allow to obtain the distribution of humble numbers below a user-defined value in a few seconds with an emulator, requiring a few stack levels and a vector sized at the maximum number of digits only.
On a 4-bit calculator with 32 Mb RAM, it takes 25 seconds to get the first 50 numbers and 9 minutes to count all the humble numbers under 1,000,000
≪ CASE
DUP 1 == THEN END
DUP 2 MOD NOT THEN 2 / H? END
DUP 3 MOD NOT THEN 3 / H? END
DUP 5 MOD NOT THEN 5 / H? END
DUP 7 MOD NOT THEN 7 / H? END
NOT END
≫ 'H?' STO
≪ { } 0 → j
≪ WHILE DUP2 SIZE > REPEAT
'j' INCR
IF H? THEN j + END
END SWAP DROP
≫ ≫ 'HLIST' STO
≪ 1 4 FOR j
DUP LN { 2 3 5 7 } j GET LN / CEIL SWAP NEXT @ max exponent for n factor is ceil(log(max value)/log(n))
DUP LOG CEIL { } + 0 CON @ create vector of counters
→ m2 m3 m5 m7 max cnt
≪ 1
0 m2 FOR j
j 2 1 IFTE * DUP
0 m3 FOR k
k 3 1 IFTE * DUP
0 m5 FOR m
m 5 1 IFTE * DUP
0 m7 FOR n
n 7 1 IFTE *
IF DUP max ≤ THEN
cnt OVER XPON 1 + DUP2 GET 1 + PUT 'cnt' STO
ELSE m7 'n' STO END @ only way to break a FOR..NEXT loop in RPL
NEXT DROP
NEXT DROP
NEXT DROP
NEXT DROP
cnt
≫ ≫ 'HCNT' STO
50 HLIST 999999999999 HCNT
- Output:
2: { 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 }
1: [ 9 36 95 197 356 579 882 1272 1767 2381 3113 3984 ]
Ruby
Brute force and slow
Checks if each number upto limit is humble number.
def humble?(i)
while i % 2 == 0; i /= 2 end
while i % 3 == 0; i /= 3 end
while i % 5 == 0; i /= 5 end
while i % 7 == 0; i /= 7 end
i == 1
end
count, num = 0, 0
digits = 10 # max digits for humble numbers
limit = 10 ** digits # max numbers to search through
humble = Array.new(digits + 1, 0)
while (num += 1) < limit
if humble?(num)
humble[num.to_s.size] += 1
print num, " " if count < 50
count += 1
end
end
print "\n\nOf the first #{count} humble numbers:\n"
(1..digits).each { |num| printf("%5d have %2d digits\n", humble[num], num) }
- Output:
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
Of the first 7574 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
2381 have 10 digits
Direct Generation: Orders of magnitude faster
Generate humble numbers directly.
def humble(digits)
h = [1]
x2, x3, x5, x7 = 2, 3, 5, 7
i, j, k, l = 0, 0, 0, 0
n = 0
while n += 1 # ruby => 2.6: (1..).each do |n|
x = [x2, x3, x5, x7].min
break if x.to_s.size > digits
h[n] = x
x2 = 2 * h[i += 1] if x2 == h[n]
x3 = 3 * h[j += 1] if x3 == h[n]
x5 = 5 * h[k += 1] if x5 == h[n]
x7 = 7 * h[l += 1] if x7 == h[n]
end
h
end
digits = 50 # max digits for humble numbers
h = humble(digits) # humble numbers <= digits size
count = h.size # the total humble numbers count
#counts = h.map { |n| n.to_s.size }.tally # hash of digits counts 1..digits: Ruby => 2.7
counts = h.map { |n| n.to_s.size }.group_by(&:itself).transform_values(&:size) # Ruby => 2.4
print "First 50 Humble Numbers: \n"; (0...50).each { |i| print "#{h[i]} " }
print "\n\nOf the first #{count} humble numbers:\n"
(1..digits).each { |num| printf("%6d have %2d digits\n", counts[num], num) }
- Output:
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
Of the first 3363713 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
2381 have 10 digits
3113 have 11 digits
3984 have 12 digits
5002 have 13 digits
6187 have 14 digits
7545 have 15 digits
9081 have 16 digits
10815 have 17 digits
12759 have 18 digits
14927 have 19 digits
17323 have 20 digits
19960 have 21 digits
22853 have 22 digits
26015 have 23 digits
29458 have 24 digits
33188 have 25 digits
37222 have 26 digits
41568 have 27 digits
46245 have 28 digits
51254 have 29 digits
56618 have 30 digits
62338 have 31 digits
68437 have 32 digits
74917 have 33 digits
81793 have 34 digits
89083 have 35 digits
96786 have 36 digits
104926 have 37 digits
113511 have 38 digits
122546 have 39 digits
132054 have 40 digits
142038 have 41 digits
152515 have 42 digits
163497 have 43 digits
174986 have 44 digits
187004 have 45 digits
199565 have 46 digits
212675 have 47 digits
226346 have 48 digits
240590 have 49 digits
255415 have 50 digits
Sidef
func smooth_generator(primes) {
var s = primes.len.of { [1] }
{
var n = s.map { .first }.min
{ |i|
s[i].shift if (s[i][0] == n)
s[i] << (n * primes[i])
} * primes.len
n
}
}
with (smooth_generator([2,3,5,7])) {|g|
say 50.of { g.run }.join(' ')
}
say "\nThe digit counts of humble numbers"
say '═'*35
with (smooth_generator([2,3,5,7])) {|g|
for (var(d=1,c=0); d <= 20; ++c) {
var n = g.run
n.len > d || next
say "#{'%10s'%c.commify} have #{'%2d'%d} digit#{[:s,''][d==1]}"
(c, d) = (0, n.len)
}
}
- Output:
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
The digit counts of humble numbers
═══════════════════════════════════
9 have 1 digit
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1,272 have 8 digits
1,767 have 9 digits
2,381 have 10 digits
3,113 have 11 digits
3,984 have 12 digits
5,002 have 13 digits
6,187 have 14 digits
7,545 have 15 digits
9,081 have 16 digits
10,815 have 17 digits
12,759 have 18 digits
14,927 have 19 digits
17,323 have 20 digits
Tcl
proc humble? x {
foreach f {2 3 5 7} {
while {$x % $f == 0} {set x [expr {$x / $f}]}
}
return [expr {$x == 1}]
}
set t1 {}
for {set i 1} {[llength $t1] < 50} {incr i} {
if [humble? $i] {lappend t1 $i}
}
puts $t1
Task 1:
- Output:
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
Task 2, took a long while due to brute force:
proc task2 {nmax} {
puts "Distribution of digit length for the first $nmax humble numbers"
set nHumble 0
for {set i 1} {$nHumble < $nmax} {incr i} {
if {[humble? $i]} {
incr nHumble
incr N([string length $i])
}
}
parray N
}
task2 4096
- Output:
~ $ time ./humble.tcl Distribution of digit length for the first 4096 humble numbers N(1) = 9 N(2) = 36 N(3) = 95 N(4) = 197 N(5) = 356 N(6) = 579 N(7) = 882 N(8) = 1272 N(9) = 670 real 38m11.922s user 0m0.000s sys 0m0.093s
Visual Basic .NET
Module Module1
Function IsHumble(i As Long) As Boolean
If i <= 1 Then
Return True
End If
If i Mod 2 = 0 Then
Return IsHumble(i \ 2)
End If
If i Mod 3 = 0 Then
Return IsHumble(i \ 3)
End If
If i Mod 5 = 0 Then
Return IsHumble(i \ 5)
End If
If i Mod 7 = 0 Then
Return IsHumble(i \ 7)
End If
Return False
End Function
Sub Main()
Dim LIMIT = Short.MaxValue
Dim humble As New Dictionary(Of Integer, Integer)
Dim count = 0L
Dim num = 1L
While count < LIMIT
If (IsHumble(num)) Then
Dim str = num.ToString
Dim len = str.Length
If len > 10 Then
Exit While
End If
If humble.ContainsKey(len) Then
humble(len) += 1
Else
humble(len) = 1
End If
If count < 50 Then
Console.Write("{0} ", num)
End If
count += 1
End If
num += 1
End While
Console.WriteLine(vbNewLine)
Console.WriteLine("Of the first {0} humble numbers:", count)
num = 1
While num < humble.Count
If humble.ContainsKey(num) Then
Dim c = humble(num)
Console.WriteLine("{0,5} have {1,2} digits", c, num)
num += 1
Else
Exit While
End If
End While
End Sub
End Module
- Output:
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
Of the first 7574 humble numbers:
9 have 1 digits
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1272 have 8 digits
1767 have 9 digits
Wren
Wren doesn't have arbitrary precision arithmetic and 'safe' integer operations are limited to a maximum absolute value of 2^53-1 (a 16 digit number). So there is no point in trying to generate humble numbers beyond that.
import "/fmt" for Fmt
import "/math" for Int, Nums
import "/sort" for Find
var humble = Fn.new { |n|
var h = List.filled(n, 0)
h[0] = 1
var next2 = 2
var next3 = 3
var next5 = 5
var next7 = 7
var i = 0
var j = 0
var k = 0
var l = 0
for (m in 1...n) {
h[m] = Nums.min([next2, next3, next5, next7])
if (h[m] == next2) {
i = i + 1
next2 = 2 * h[i]
}
if (h[m] == next3) {
j = j + 1
next3 = 3 * h[j]
}
if (h[m] == next5) {
k = k + 1
next5 = 5 * h[k]
}
if (h[m] == next7) {
l = l + 1
next7 = 7 * h[l]
}
}
return h
}
var n = 43000 // say
var h = humble.call(n)
System.print("The first 50 humble numbers are:")
System.print(h[0..49])
var f = Find.all(h, Int.maxSafe) // binary search
var maxUsed = f[0] ? f[2].min + 1 : f[2].min
var maxDigits = 16 // Int.maxSafe (2^53 -1) has 16 digits
var counts = List.filled(maxDigits + 1, 0)
var digits = 1
var pow10 = 10
for (i in 0...maxUsed) {
while (true) {
if (h[i] >= pow10) {
pow10 = pow10 * 10
digits = digits + 1
} else break
}
counts[digits] = counts[digits] + 1
}
System.print("\nOf the first %(Fmt.dc(0, maxUsed)) humble numbers:")
for (i in 1..maxDigits) {
var s = (i != 1) ? "s" : ""
System.print("%(Fmt.dc(9, counts[i])) have %(Fmt.d(2, i)) digit%(s)")
}- Output:
The first 50 humble numbers are:
[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]
Of the first 42,037 humble numbers:
9 have 1 digit
36 have 2 digits
95 have 3 digits
197 have 4 digits
356 have 5 digits
579 have 6 digits
882 have 7 digits
1,272 have 8 digits
1,767 have 9 digits
2,381 have 10 digits
3,113 have 11 digits
3,984 have 12 digits
5,002 have 13 digits
6,187 have 14 digits
7,545 have 15 digits
8,632 have 16 digits
XPL0
func Humble(N); \Return 'true' if N is a humble number
int N;
[if N = 1 then return true;
if rem(N/2) = 0 then return Humble(N/2);
if rem(N/3) = 0 then return Humble(N/3);
if rem(N/5) = 0 then return Humble(N/5);
if rem(N/7) = 0 then return Humble(N/7);
return false;
];
int N, C, D, P;
[N:= 1; C:= 0;
loop [if Humble(N) then
[C:= C+1;
IntOut(0, N); ChOut(0, ^ );
if C >= 50 then quit;
];
N:= N+1;
];
CrLf(0);
D:= 1; P:= 10; N:= 1; C:= 0;
loop [if Humble(N) then
[if N >= P then
[IntOut(0, D);
Text(0, ": ");
IntOut(0, C);
CrLf(0);
C:= 0;
D:= D+1;
if D > 9 then quit;
P:= P*10;
];
C:= C+1;
];
N:= N+1;
];
]- Output:
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 1: 9 2: 36 3: 95 4: 197 5: 356 6: 579 7: 882 8: 1272 9: 1767
zkl
GNU Multiple Precision Arithmetic Libraryvar [const] BI=Import("zklBigNum"); // libGMP
var one = BI(1), two = BI(2), three = BI(3),
five = BI(5), seven = BI(7);
fcn humble(n){ // --> List of BigInt Humble numbers
h:=List.createLong(n); h.append(one);
next2,next3 := two.copy(), three.copy();
next5,next7 := five.copy(), seven.copy();
reg i=0,j=0,k=0,l=0;
do(n-1){
h.append( hm:=BI(next2.min(next3.min(next5.min(next7)))) );
if(hm==next2) next2.set(two) .mul(h[i+=1]);
if(hm==next3) next3.set(three).mul(h[j+=1]);
if(hm==next5) next5.set(five) .mul(h[k+=1]);
if(hm==next7) next7.set(seven).mul(h[l+=1]);
}
h
}fcn __main__{
const N = 5 * 1e6; // calculate the first 1 million humble numbers, say
h:=humble(N);
println("The first 50 humble numbers are:\n ",h[0,50].concat(" "));
counts:=Dictionary(); // tally the number of digits in each number
h.apply2('wrap(n){ counts.incV(n.numDigits) });
println("\nOf the first %,d humble numbers:".fmt(h.len()));
println("Digits Count");
foreach n in (counts.keys.apply("toInt").sort()){
println("%2d %,9d".fmt(n,counts[n], n));
}
}- Output:
The first 50 humble numbers are: 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 Of the first 5,000,000 humble numbers: Digits Count 1 9 2 36 3 95 4 197 5 356 6 579 7 882 8 1,272 9 1,767 10 2,381 11 3,113 12 3,984 13 5,002 14 6,187 15 7,545 16 9,081 17 10,815 18 12,759 19 14,927 20 17,323 21 19,960 22 22,853 23 26,015 24 29,458 25 33,188 26 37,222 27 41,568 28 46,245 29 51,254 30 56,618 31 62,338 32 68,437 33 74,917 34 81,793 35 89,083 36 96,786 37 104,926 38 113,511 39 122,546 40 132,054 41 142,038 42 152,515 43 163,497 44 174,986 45 187,004 46 199,565 47 212,675 48 226,346 49 240,590 50 255,415 51 270,843 52 286,880 53 303,533 54 320,821 55 338,750 56 115,460
Categories:
- Pages with syntax highlighting errors
- Programming Tasks
- Prime Numbers
- 11l
- Action!
- Ada
- ALGOL 68
- ALGOL W
- AutoHotkey
- AWK
- C
- C sharp
- System.Numerics
- C++
- Crystal
- D
- Delphi
- Elm
- F Sharp
- Factor
- FreeBASIC
- Go
- Haskell
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- Lua
- Mathematica
- Wolfram Language
- Nim
- Pascal
- Perl
- Phix
- Phix/mpfr
- PL/I
- PL/M
- Polyglot:PL/I and PL/M
- Python
- Quackery
- Racket
- Raku
- REXX
- Ring
- Ring examples needing attention
- Examples needing attention
- RPL
- Ruby
- Sidef
- Tcl
- Visual Basic .NET
- Wren
- Wren-fmt
- Wren-math
- Wren-sort
- XPL0
- Zkl
- GMP