Humble numbers: Difference between revisions
→{{header|F_Sharp|F#}}: add faster implementations at the end...
GordonBGood (talk | contribs) (→{{header|F_Sharp|F#}}: added Elm contribution above...) |
GordonBGood (talk | contribs) (→{{header|F_Sharp|F#}}: add faster implementations at the end...) |
||
Line 2,561:
There are 12759 humble numbers with 18 digits
</pre>
===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:
<syntaxhighlight lang="fsharp">// 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)</syntaxhighlight>
{{out}}
<pre>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</pre>
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:
<syntaxhighlight lang="fsharp">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)</syntaxhighlight>
{{out}}
<pre>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</pre>
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.
=={{header|Factor}}==
| |||