Iterated digits squaring: Difference between revisions
(J, with a perhaps ironic note on timing) |
|||
| Line 254: | Line 254: | ||
{{out}} |
{{out}} |
||
<pre>85744333</pre> |
<pre>85744333</pre> |
||
=={{header|J}}== |
|||
Here's an expression to turn a number into digits: |
|||
<lang J>digits=: 10&#.inv</lang> |
|||
And here's an expression to square them and find their sum: |
|||
<lang J>sumdigsq=: +/"1@:*:@digits</lang> |
|||
But note that while the task description claims "you always end with either 1 or 89", that claim is somewhat arbitrary. |
|||
<lang J> sumdigsq^:(i.16) 15 |
|||
15 26 40 16 37 58 89 145 42 20 4 16 37 58 89 145</lang> |
|||
You could just as easily claim that you always end with either 1 or 4. So here's a routine which repeats the sum-square process until the sequence converges, or until it reaches the value 4: |
|||
<lang J>itdigsq4=:4 = sumdigsq^:(0=e.&4)^:_"0</lang> |
|||
But we do not actually need to iterate. The largest value after the first iteration would be: |
|||
<lang J> sumdigsq 99999999 |
|||
648</lang> |
|||
So we could write a routine which works for the intended range, and stops after the first iteration: |
|||
<lang J>digsq1e8=:(I.itdigsq1 i.649) e.~ sumdigsq<lang> |
|||
And this is sufficient to find our result. We don't want to compute the entire batch of values in one pass, however, so let's break this up in to 100 batches of one million each: |
|||
<lang J> +/+/@:-.@digsq1e8"1(1+i.100 1e6) |
|||
85744333</lang> |
|||
Of course, there are faster ways of obtaining that result. The fastest is probably this: |
|||
<lang J> 85744333 |
|||
85744333</lang> |
|||
This might be thought of as representing the behavior of a highly optimized compiled program. We could abstract this further by using the previous expression at compile time, so we would not have to hard code it. |
|||
=={{header|Java}}== |
=={{header|Java}}== |
||
Revision as of 23:18, 15 September 2014
You are encouraged to solve this task according to the task description, using any language you may know.
If you add the square of the digits of a Natural number (an integer bigger than zero), you always end with either 1 or 89:
15 -> 26 -> 40 -> 16 -> 37 -> 58 -> 89 7 -> 49 -> 97 -> 130 -> 10 -> 1
An example in Python:
<lang python>>>> step = lambda x: sum(int(d) ** 2 for d in str(x)) >>> iterate = lambda x: x if x in [1, 89] else iterate(step(x)) >>> [iterate(x) for x in xrange(1, 20)] [1, 89, 89, 89, 89, 89, 1, 89, 89, 1, 89, 89, 1, 89, 89, 89, 89, 89, 1]</lang>
- Task
- Count how many number chains for integers 1 <= n < 100_000_000 end with a value 89.
Or, for much less credit - (showing that your algorithm and/or language is slow):
- Count how many number chains for integers 1 <= n < 1_000_000 end with a value 89.
This problem derives from the Project Euler problem 92.
For a quick algorithm for this task see the talk page
- Cf
D
A simple memoizing partially-imperative brute-force solution: <lang d>import std.stdio, std.algorithm, std.range, std.functional;
uint step(uint x) pure nothrow @safe @nogc {
uint total = 0;
while (x) {
total += (x % 10) ^^ 2;
x /= 10;
}
return total;
}
uint iterate(in uint x) nothrow @safe {
return (x == 89 || x == 1) ? x : x.step.memoize!iterate;
}
void main() {
iota(1, 100_000_000).filter!(x => x.iterate == 89).count.writeln;
}</lang>
- Output:
85744333
The run-time is about 10 seconds compiled with ldc2.
A fast imperative brute-force solution: <lang d>void main() nothrow @nogc {
import core.stdc.stdio: printf;
enum uint magic = 89; enum uint limit = 100_000_000; uint[(9 ^^ 2) * 8 + 1] lookup = void;
uint[10] squares;
foreach (immutable i, ref x; squares)
x = i ^^ 2;
foreach (immutable uint i; 1 .. lookup.length) {
uint x = i;
while (x != magic && x != 1) {
uint total = 0;
while (x) {
total += squares[(x % 10)];
x /= 10;
}
x = total;
}
lookup[i] = x == magic; }
uint magicCount = 0;
foreach (immutable uint i; 1 .. limit) {
uint x = i;
uint total = 0;
while (x) {
total += squares[(x % 10)];
x /= 10;
}
magicCount += lookup[total]; }
printf("%u\n", magicCount);
}</lang> The output is the same. The run-time is less than 3 seconds compiled with ldc2.
A more efficient solution: <lang d>import core.stdc.stdio, std.algorithm, std.range;
enum factorial = (in uint n) pure nothrow @safe @nogc
=> reduce!q{a * b}(1u, iota(1u, n + 1));
uint iLog10(in uint x) pure nothrow @safe @nogc in {
assert(x > 0);
} body {
return (x >= 1_000_000_000) ? 9 :
(x >= 100_000_000) ? 8 :
(x >= 10_000_000) ? 7 :
(x >= 1_000_000) ? 6 :
(x >= 100_000) ? 5 :
(x >= 10_000) ? 4 :
(x >= 1_000) ? 3 :
(x >= 100) ? 2 :
(x >= 10) ? 1 : 0;
}
uint nextStep(uint x) pure nothrow @safe @nogc {
typeof(return) result = 0;
while (x > 0) {
result += (x % 10) ^^ 2;
x /= 10;
}
return result;
}
uint check(in uint[] number) pure nothrow @safe @nogc {
uint candidate = reduce!((tot, n) => tot * 10 + n)(0, number);
while (candidate != 89 && candidate != 1)
candidate = candidate.nextStep;
if (candidate == 89) {
uint[10] digitsCount;
foreach (immutable d; number)
digitsCount[d]++;
return reduce!((r, c) => r / c.factorial)
(number.length.factorial, digitsCount);
}
return 0;
}
void main() nothrow @nogc {
enum uint limit = 100_000_000; immutable uint cacheSize = limit.iLog10;
uint[cacheSize] number; uint result = 0; uint i = cacheSize - 1;
while (true) {
if (i == 0 && number[i] == 9)
break;
if (i == cacheSize - 1 && number[i] < 9) {
number[i]++;
result += number.check;
} else if (number[i] == 9) {
i--;
} else {
number[i]++;
number[i + 1 .. $] = number[i];
i = cacheSize - 1;
result += number.check;
}
}
printf("%u\n", result);
}</lang> The output is the same. The run-time is about 0.04 seconds or less. This third version was ported to D and improved from: mathblog.dk/project-euler-92-square-digits-number-chain/
A purely functional version, from the Haskell code. It includes two functions currently missing in Phobos used in the Haskell code.
<lang d>import std.stdio, std.typecons, std.traits, std.typetuple, std.range, std.algorithm;
auto divMod(T)(T x, T y) pure nothrow @safe @nogc {
return tuple(x / y, x % y);
}
auto expand(alias F, B)(B x) pure nothrow @safe @nogc if (isCallable!F &&
is(ParameterTypeTuple!F == TypeTuple!B)
&& __traits(isSame, TemplateOf!(ReturnType!F), Nullable)
&& isTuple!(TemplateArgsOf!(ReturnType!F)[0])
&& is(TemplateArgsOf!(TemplateArgsOf!(ReturnType!F)[0])[1] == B)) {
alias NAB = ReturnType!F; alias AB = TemplateArgsOf!NAB[0]; alias A = AB.Types[0];
struct Expand {
bool first;
NAB last;
@property bool empty() pure nothrow @safe @nogc {
if (first) {
first = false;
popFront;
}
return last.isNull;
}
@property A front() pure nothrow @safe @nogc {
if (first) {
first = false;
popFront;
}
return last.get[0];
}
void popFront() pure nothrow @safe @nogc { last = F(last.get[1]); }
}
return Expand(true, NAB(AB(A.init, x)));
}
//------------------------------------------------
uint step(uint x) pure nothrow @safe @nogc {
Nullable!(Tuple!(uint, uint)) f(uint n) pure nothrow @safe @nogc {
return (n == 0) ? typeof(return)() : typeof(return)(divMod(n, 10u).reverse);
}
return expand!f(x).map!(x => x ^^ 2).sum;
}
uint iter(uint x) pure nothrow @nogc {
return x.recurrence!((a, n) => step(a[n - 1])).filter!(y => y.among!(1, 89)).front;
}
void main() {
iota(1u, 100_000u).filter!(n => n.iter == 89).count.writeln;
}</lang> With a small back-porting (to run it with the Phobos of LDC2 2.065) it runs in about 15.5 seconds.
Haskell
Basic solution that contains just a little more than the essence of this computation. This runs in less than eight minutes: <lang haskell>import Data.List (unfoldr) import Data.Tuple (swap)
step :: Int -> Int step = sum . map (^ 2) . unfoldr f where
f 0 = Nothing f n = Just . swap $ n `divMod` 10
iter :: Int -> Int iter = head . filter (`elem` [1, 89]) . iterate step
main = do
print $ length $ filter ((== 89) . iter) [1 .. 99999999]</lang>
- Output:
85744333
J
Here's an expression to turn a number into digits:
<lang J>digits=: 10&#.inv</lang>
And here's an expression to square them and find their sum: <lang J>sumdigsq=: +/"1@:*:@digits</lang>
But note that while the task description claims "you always end with either 1 or 89", that claim is somewhat arbitrary.
<lang J> sumdigsq^:(i.16) 15 15 26 40 16 37 58 89 145 42 20 4 16 37 58 89 145</lang>
You could just as easily claim that you always end with either 1 or 4. So here's a routine which repeats the sum-square process until the sequence converges, or until it reaches the value 4:
<lang J>itdigsq4=:4 = sumdigsq^:(0=e.&4)^:_"0</lang>
But we do not actually need to iterate. The largest value after the first iteration would be:
<lang J> sumdigsq 99999999 648</lang>
So we could write a routine which works for the intended range, and stops after the first iteration: <lang J>digsq1e8=:(I.itdigsq1 i.649) e.~ sumdigsq<lang>
And this is sufficient to find our result. We don't want to compute the entire batch of values in one pass, however, so let's break this up in to 100 batches of one million each:
<lang J> +/+/@:-.@digsq1e8"1(1+i.100 1e6) 85744333</lang>
Of course, there are faster ways of obtaining that result. The fastest is probably this: <lang J> 85744333 85744333</lang>
This might be thought of as representing the behavior of a highly optimized compiled program. We could abstract this further by using the previous expression at compile time, so we would not have to hard code it.
Java
<lang java>import java.util.stream.IntStream;
public class IteratedDigitsSquaring {
public static void main(String[] args) {
long r = IntStream.range(1, 100_000_000)
.parallel()
.filter(n -> calc(n) == 89)
.count();
System.out.println(r);
}
private static int calc(int n) {
while (n != 89 && n != 1) {
int total = 0;
while (n > 0) {
total += Math.pow(n % 10, 2);
n /= 10;
}
n = total;
}
return n;
}
}</lang>
85744333
Perl
<lang perl>use warnings; use strict;
my @sq = map { $_ ** 2 } 0 .. 9; my %cache; my $cnt = 0;
sub Euler92 {
my $n = 0 + join( , sort split( , shift ) );
$cache{$n} //= ($n == 1 || $n == 89) ? $n :
Euler92( sum( @sq[ split , $n ] ) )
}
sub sum {
my $sum; $sum += shift while @_; $sum;
}
for (1 .. 100_000_000) {
++$cnt if Euler92( $_ ) == 89;
}
print $cnt;</lang>
85744333
Perl 6
This fairly abstract version does caching and filtering to reduce the number of values it needs to check and moves calculations out of the hot loop, but is still interminably slow... even for just up to 1,000,000.
<lang perl6>constant @sq = ^10 X** 2; my $cnt = 0; my %cache;
sub Euler92($n) {
%cache{$n} //=
$n == any(1,89) ?? $n !!
Euler92( [+] @sq[$n.comb] )
}
for 1 .. 1_000_000 -> $n {
my $i = +$n.comb.sort.join; ++$cnt if Euler92($i) == 89;
}
say $cnt;</lang>
- Output:
856929
All is not lost, however. Through the use of gradual typing, Perl 6 scales down as well as up, so this jit-friendly version is performant enough to brute force the larger calculation: <lang perl6>my @cache; @cache[1] = 1; @cache[89] = 89;
sub Euler92(int $n) {
$n < 649 # 99,999,999 sums to 648, so no point remembering more
?? (@cache.at_pos($n) //= ids($n))
!! ids($n)
}
sub ids(int $num --> int) {
my int $n = $num;
my int $ten = 10;
my int $sum = 0;
my int $t;
my int $c;
repeat until $n == 89 or $n == 1 {
$sum = 0;
repeat {
$t = $n div $ten;
$c = $n - $t * $ten;
$sum = $sum + $c * $c;
} while $n = $t;
$n = @cache.at_pos($sum) // $sum;
}
$n;
}
my int $cnt = 0; for 1 .. 100_000_000 -> int $n {
$cnt = $cnt + 1 if Euler92($n) == 89;
} say $cnt;</lang>
- Output:
85744333
This runs in under ten minutes. We can reduce this to 4 minutes by writing in the NQP (Not Quite Perl6) subset of the language: <lang perl6>my $cache := nqp::list_i(); nqp::bindpos_i($cache, 650, 0); nqp::bindpos_i($cache, 1, 1); nqp::bindpos_i($cache, 89, 89);
sub Euler92(int $n) {
$n < 650
?? nqp::bindpos_i($cache,$n,ids($n)) !! ids($n) }
sub ids(int $num --> int) {
my int $n = $num;
my int $ten = 10;
my int $sum = 0;
my int $t;
my int $c;
repeat until $n == 89 or $n == 1 {
$sum = 0; repeat { $t = nqp::div_i($n, $ten); $c = $n - $t * $ten; $sum = $sum + $c * $c; } while $n = $t; $n = nqp::atpos_i($cache,$sum) || $sum;
} $n;
}
my int $cnt = 0; for 1 .. 100_000_000 -> int $n {
$cnt = $cnt + 1 if Euler92($n) == 89;
} say $cnt;</lang>
Python
Combinatorics
Translation of D
<lang python>from math import ceil, log10, factorial
def next_step(x):
result = 0
while x > 0:
result += (x % 10) ** 2
x /= 10
return result
def check(number):
candidate = 0
for n in number:
candidate = candidate * 10 + n
while candidate != 89 and candidate != 1:
candidate = next_step(candidate)
if candidate == 89:
digits_count = [0] * 10
for d in number:
digits_count[d] += 1
result = factorial(len(number))
for c in digits_count:
result /= factorial(c)
return result
return 0
def main():
limit = 100000000 cache_size = int(ceil(log10(limit))) assert 10 ** cache_size == limit
number = [0] * cache_size result = 0 i = cache_size - 1
while True:
if i == 0 and number[i] == 9:
break
if i == cache_size - 1 and number[i] < 9:
number[i] += 1
result += check(number)
elif number[i] == 9:
i -= 1
else:
number[i] += 1
for j in xrange(i + 1, cache_size):
number[j] = number[i]
i = cache_size - 1
result += check(number)
print result
main()</lang>
- Output:
85744333
The run-time is less than half a second.
Translation of Ruby
<lang ruby> from itertools import combinations_with_replacement from array import array from time import clock D = 8 F = [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000] def b(n):
yield 1
for g in range(1,n+1):
gn = g
res = 0
while gn > 0:
gn,rem = divmod(gn,10)
res += rem**2
if res==89:
yield 0
else:
yield res
N = array('I',b(81*D)) for n in range(2,len(N)):
q = N[n]
while q>1:
q = N[q]
N[n] = q
es = clock() z = 0 for n in combinations_with_replacement(range(10),D):
t = 0
for g in n:
t += g*g
if N[t] == 0:
continue
t = [0,0,0,0,0,0,0,0,0,0]
for g in n:
t[g] += 1
t1 = F[D]
for g in t:
t1 /= F[g]
z += t1
ee = clock() - es print "\nD==" + str(D) + "\n " + str(z) + " numbers produce 1 and " + str(10**D-z) + " numbers produce 89" print "Time ~= " + str(ee) + " secs" </lang>
- Output:
D==8 14255667 numbers produce 1 and 85744333 numbers produce 89 Time ~= 0.14 secs D==11 15091199357 numbers produce 1 and 84908800643 numbers produce 89 Time ~= 1.12 secs D==14 13770853279685 numbers produce 1 and 86229146720315 numbers produce 89 Time ~= 7.46 secs D==17 12024696404768025 numbers produce 1 and 87975303595231975 numbers produce 89 Time ~= 34.16 secs
Python: Simple caching
<lang python>>>> from functools import lru_cache >>> @lru_cache(maxsize=1024) def ids(n): if n in {1, 89}: return n else: return ids(sum(int(d) ** 2 for d in str(n)))
>>> ids(15)
89
>>> [ids(x) for x in range(1, 21)]
[1, 89, 89, 89, 89, 89, 1, 89, 89, 1, 89, 89, 1, 89, 89, 89, 89, 89, 1, 89]
>>> sum(ids(x) == 89 for x in range(1, 100000000))
85744333
>>> </lang>
This took a much longer time, in the order of hours.
Python: Enhanced caching
Notes that the order of digits in a number does not affect the final result so caches the digits of the number in sorted order for more hits. <lang python>>>> from functools import lru_cache >>> @lru_cache(maxsize=1024) def _ids(nt): if nt in {('1',), ('8', '9')}: return nt else: return _ids(tuple(sorted(str(sum(int(d) ** 2 for d in nt)))))
>>> def ids(n):
return int(.join(_ids(tuple(sorted(str(n))))))
>>> ids(1), ids(15) (1, 89) >>> [ids(x) for x in range(1, 21)] [1, 89, 89, 89, 89, 89, 1, 89, 89, 1, 89, 89, 1, 89, 89, 89, 89, 89, 1, 89] >>> sum(ids(x) == 89 for x in range(1, 100000000)) 85744333 >>> _ids.cache_info() CacheInfo(hits=99991418, misses=5867462, maxsize=1024, currsize=1024) >>> </lang>
This took tens of minutes to run.
Racket
This contains two versions (in one go). The naive version which can (and should, probably) be used for investigating a single number. The second version can count the IDSes leading to 89 for powers of 10.
<lang racket>#lang racket
- Tim-brown 2014-09-11
- The basic definition.
- It is possible to memoise this or use fixnum (native) arithmetic, but frankly iterating over a
- hundred million, billion, trillion numbers will be slow. No matter how you do it.
(define (digit^2-sum n)
(let loop ((n n) (s 0)) (if (= 0 n) s (let-values ([(q r) (quotient/remainder n 10)]) (loop q (+ s (sqr r)))))))
(define (iterated-digit^2-sum n)
(match (digit^2-sum n) [0 0] [1 1] [89 89] [(app iterated-digit^2-sum rv) rv]))
- Note that
- ids(345) = ids(354) = ids(435) = ids(453) = ids(534) = ids(543) = 50 --> 89
- One calculation does for 6 candidates.
- The plan
- - get all the ordered combinations of digits including 0's which can be used both as digits and
- "padding" digits in the most significant digits. (n.b. all-zeros is not in the range to be
- tested and should be dropped)
- - find the digit sets that have an IDS of 89
- - find out how many combinations there are of these digits
- output
- a list of n-digits long lists containing all of the digit combinations.
- a smart bunny would figure out the sums of the digits as they're generated but I'll plod
- along step-by-step. a truly smart bunny would also count the combinations. that said, I
- don't think I do much unnecessary computation here.
(define (all-digit-lists n-digits)
(define (inner remain acc least-digit)
(cond
[(zero? remain) (list (list))]
[(= least-digit 10) null]
[else
(for*/list
((ld+ (in-range least-digit 10))
(rgt (in-list (inner (sub1 remain) empty ld+))))
(append acc (cons ld+ rgt)))]))
(inner n-digits '() 0))
- We calculate IDS differently since we're presented with a list of digits rather than a number
(define (digit-list-IDS c)
(define (digit-combo-IDS c) (apply + (map sqr c))) (iterated-digit^2-sum (digit-combo-IDS c)))
- ! (factiorial) -- everyone's favourite combinatorial function! (that's just an exclamation mark)
- there's one in (require math/number-theory) for any heavy lifting, but we're not or I could import
- it from math/number-theory -- but this is about all I need. A lookup table is going to be faster
- than a more general function.
(define (! n)
(case n [(0 1) 1] [(2) 2] [(3) 6] [(4) 24] [(5) 120] [(6) 720] [(7) 5040] [(8) 40320] [(9) 362880] [else (* n (! (sub1 n)))] ; I expect this clause'll never be called ))
- We need to count the permutations -- digits are in order so we can use the tail (cdr) function for
- determining my various k's. See
- https://en.wikipedia.org/wiki/Combination
(define (count-digit-list-permutations c #:length (l (length c)) #:length! (l! (! l)))
(let loop ((c c) (i 0) (prev -1 #;"never a digit") (p l!))
(match c
[(list) (/ p (! i))]
[(cons (== prev) d) (loop d (+ i 1) prev p)]
[(cons a d) (loop d 1 a (/ p (! i)))])))
- Wrap it all up in a neat function
(define (count-89s-in-100... n)
(define n-digits (order-of-magnitude n)) (define combos (drop (all-digit-lists n-digits) 1)) ; don't want first one which is "all-zeros" (for/sum ((c (in-list combos)) #:when (= 89 (digit-list-IDS c))) (count-digit-list-permutations c #:length n-digits)))
(displayln "Testing permutations:") (time (printf "1000000:\t~a~%" (count-89s-in-100... 1000000))) (time (printf "100000000:\t~a~%" (count-89s-in-100... 100000000))) (time (printf "1000000000:\t~a~%" (count-89s-in-100... 1000000000))) (time (printf "1000000000000:\t~a~%" (count-89s-in-100... 1000000000000))) (newline)
- Do these last, since the 10^8 takes longer than my ADHD can cope with
(displayln "Testing one number at a time (somewhat slower):") (time (printf "1000000:\t~a~%" (for/sum ((n (in-range 1 1000000))
#:when (= 89 (iterated-digit^2-sum n))) 1)))
(time (printf "100000000:\t~a~%" (for/sum ((n (in-range 1 100000000))
#:when (= 89 (iterated-digit^2-sum n))) 1)))
{module+ test
(require tests/eli-tester) [test (iterated-digit^2-sum 15) => 89 (iterated-digit^2-sum 7) => 1 (digit-combo-perms '()) => 1 (digit-combo-perms '(1 2 3)) => 6 (digit-combo-perms '(1 1 3)) => 3 (for/sum ((n (in-range 1 1000000)) #:when (= 89 (iterated-digit^2-sum n))) 1) => 856929 (all-digit-lists 1) => '((0) (1) (2) (3) (4) (5) (6) (7) (8) (9)) (length (all-digit-lists 2)) => 55 (length (all-digit-lists 3)) => 220 (count-89s-in-100... 1000000) => 856929] }</lang>
- Output:
Testing permutations: 1000000: 856929 cpu time: 8 real time: 8 gc time: 0 100000000: 85744333 cpu time: 44 real time: 43 gc time: 0 1000000000: 854325192 cpu time: 112 real time: 110 gc time: 20 1000000000000: 850878696414 cpu time: 1108 real time: 1110 gc time: 472 Testing one number at a time (somewhat slower): 1000000: 856929 cpu time: 1168 real time: 1171 gc time: 0 100000000: 85744333 cpu time: 130720 real time: 130951 gc time: 0
Ok, so maybe 131 seconds is not so flattering -- but I have not memoised or anything fancy like that, because even doing that isn't going to come anywhere near competing with 44ms.
REXX
version 1
This REXX version is essentially a brute force approach and is slow. <lang rexx>/*REXX program to perform iterated digits squaring ('til sum=1 | sum=89)*/ parse arg n . /*get optional N from the C.L. */ if n== then n=1000000 /*Was N given? No, use default.*/ !.=0; do m=1 to 9; !.m=m**2; end /*build a short-cut for squares. */
- .=0 /*count of 1 & 89 results so far.*/
do j=1 for n /* [↓] process each num in range*/
x=j /*use X for a proxy for the J var*/
do until s==1 | s==89 /*add the squared digits of X.*/
s=0 /*set the sum to zero initially. */
do k=1 for length(x) /*process each of the digits in X*/
_=substr(x,k,1) /*pick off a particular X digit*/
s=s+!._ /*do a fast squaring of it & sum.*/
end /*k*/ /* [↑] S≡is sum of squared digs.*/
x=s /*subsitute the sum for "new" X. */
end /*until*/ /* [↑] keep looping 'til S=1|89.*/
#.s=#.s+1 /*bump the counter for 1's | 89's*/
end /*j*/
do i=1 by 89-1 for 2; c=' 'right('"'i'"',4)' chains '
say 'count of' c 'for all natural numbers up to ' n " is " #.i
end /*i*/
/*stick a fork in it, we're done.*/</lang>
output when using the default input:
count of "1" chains for all natural numbers up to 1000000 is 143071 count of "89" chains for all natural numbers up to 1000000 is 856929
version 2
This Rexx version is optimized (by processing the numbers in four-digit chunks), but it's still on the slow side. <lang rexx>/*REXX program to perform iterated digits squaring ('til sum=1 | sum=89)*/ parse arg n . /*get optional N from the C.L. */ if n== then n=100000000 /*Was N given? No, use default.*/ !.=0; do m=1 to 9; !.m=m**2; end /*build a short-cut for squares. */ $.=.; $.0=0; $.00=0; $.000=0; $.0000=0; @.=. /*short-cuts for some sums*/ @.=. /*placeholder for computed sums. */
- .=0 /*count of 1 & 89 results so far.*/
do j=1 for n /* [↓] process each num in range*/
s=sumds(j) /*get the sum of squared digits. */
#.s=#.s+1 /*bump the counter for 1's | 89's*/
end /*j*/
do i=1 by 89-1 for 2; c=' 'right('"'i'"',4)' chains '
say 'count of' c 'for all natural numbers up to ' n " is " #.i
end /*i*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────SUMDS subroutine────────────────────*/ sumds: parse arg z; p=0
do m=1 by 4 to length(z)
p=p + summer(substr(z,m,4))
end /*m*/
if $.p\==. then return $.p /*if computed before, use the val*/ y=p
do until s==1 | s==89 /*add the squared digits of P.*/
s=0 /*set the sum to zero initially. */
do k=1 for length(y) /*process each of the digits in X*/
_=substr(y,k,1) /*pick off a particular X digit*/
s=s+!._ /*do a fast squaring of it & sum.*/
end /*k*/ /* [↑] S≡is sum of squared digs.*/
y=s /*subsitute the sum for "new" X. */
end /*until*/ /* [↑] keep looping 'til S=1|89.*/
$.p=s return s /*──────────────────────────────────SUMMER subroutine───────────────────*/ summer: parse arg y . 1 oy .; if @.y\==. then return @.y /*use old val*/ a=0
do k=1 for length(y) /*process each of the digits in X*/
_=substr(y,k,1) /*pick off a particular X digit*/
a=a+!._ /*do a fast squaring of it & sum.*/
end /*k*/ /* [↑] S≡is sum of squared digs.*/
@.oy=a return a</lang> output when using the default input:
count of "1" chains for all natural numbers up to 1000000 is 14255667 count of "89" chains for all natural numbers up to 1000000 is 85744333
Ruby
<lang ruby>
- Count how many number chains for Natural Numbers <= 100,000,000 end with a value 89.
- Nigel_Galloway
- August 26th., 2014.
require 'benchmark' D = 8 #Calculate from 1 to 10**D (8 for task) F = Array.new(D+1){|n| n==0?1:(1..n).inject(:*)} #Some small factorials g = -> n, gn=[n,0], res=0 { while gn[0]>0
gn = gn[0].divmod(10)
res += gn[1]**2
end
return res==89?0:res
- An array: N[n]==0 means that n translates to 89 and 1 means that n translates to 1 }
N = (G=Array.new(D*81+1){|n| n==0? 1:(i=g.call(n))==89 ? 0:i}).collect{|n| while n>1 do n = G[n] end; n }
z = 0 #Running count of numbers translating to 1 t = Benchmark.measure do (0..9).to_a.repeated_combination(D).each{|n| #Iterate over unique digit combinations
next if N[n.inject(0){|g,n| g+n*n}] == 0 #Count only ones
nn = [0,0,0,0,0,0,0,0,0,0] #Determine how many numbers this digit combination corresponds to
n.each{|n| nn[n] += 1}
z += nn.inject(F[D]){|gn,n| gn/F[n]} #Add to the count of numbers terminating in 1
} end puts "\n\n#{z} numbers produce 1 and #{10**D-z} numbers produce 89"
puts "\n\nTiming\n#{t}" </lang>
- Output:
100 Million (D==8)
14255667 numbers produce 1 and 85744333 numbers produce 89 Timing 0.090000 0.000000 0.090000 ( 0.098350)
100 Billion (D==11)
15091199357 numbers produce 1 and 84908800643 numbers produce 89 Timing 0.720000 0.000000 0.720000 ( 0.715757)
100 Trillion (D==14)
13770853279685 numbers produce 1 and 86229146720315 numbers produce 89 Timing 4.920000 0.000000 4.920000 ( 4.916912)
Anyone know what's after Trillion? Well here's 100 of whatever (D==17)
12024696404768025 numbers produce 1 and 87975303595231975 numbers produce 89 Timing 23.520000 0.010000 23.530000 ( 23.559256)
Tcl
All three versions below produce identical output (85744333), but the third is fastest and the first is the slowest, both by substantial margins.
Very Naïve Version
<lang tcl>proc ids n {
while {$n != 1 && $n != 89} {
set n [tcl::mathop::+ {*}[lmap x [split $n ""] {expr {$x**2}}]]
} return $n
} for {set i 1} {$i <= 100000000} {incr i} {
incr count [expr {[ids $i] == 89}]
} puts $count</lang>
Intelligent Version
Conversion back and forth between numbers and strings is slow. Using math operations directly is much faster (around 4 times in informal testing). <lang tcl>proc ids n {
while {$n != 1 && $n != 89} {
for {set m 0} {$n} {set n [expr {$n / 10}]} { incr m [expr {($n%10)**2}] } set n $m
} return $n
} for {set i 1} {$i <= 100000000} {incr i} {
incr count [expr {[ids $i] == 89}]
} puts $count</lang>
Substantially More Intelligent Version
Using the observation that the maximum value after 1 step is obtained for 999999999, which is . Thus, running one step of the reduction and then using a lookup table (which we can construct quickly at the start of the run, and which has excellent performance) is much faster overall, approximately 3–4 times than the second version above (and over 12 times faster than the first version).
- Donald, you have 1 too many 9's the value after step 1 is 81*8 = 648. Not that that is the problem here, you can not afford to go around this loop 100 million times. Notice that IDS[21] == IDS[12], IDS[123] == IDS[132] == IDS[213} ... etc, etc. The Ruby version takes about a tenth of a second.--Nigel Galloway (talk) 12:47, 31 August 2014 (UTC)
<lang tcl># Basic implementation proc ids n {
while {$n != 1 && $n != 89} {
for {set m 0} {$n} {set n [expr {$n / 10}]} { incr m [expr {($n%10)**2}] } set n $m
} return $n
}
- Build the optimised version
set body {
# Microoptimisation to avoid an unnecessary alloc in the loop
for {set m 0} {$n} {set n [expr {"$n[unset n]" / 10}]} {
incr m [expr {($n%10)**2}]
}
} set map 0 for {set i 1} {$i <= 729} {incr i} {
lappend map [ids $i]
} proc ids2 n [append body "return \[lindex [list $map] \$m\]"]
- Put this in a lambda context for a little extra speed.
apply {{} {
set count 0
for {set i 1} {$i <= 100000000} {incr i} {
incr count [expr {[ids2 $i] == 89}]
} puts $count
}}</lang>
zkl
Using brute force is a never ending process so need to be clever, which takes under a second.
<lang zkl>fcn check(number){ // a list of digits: 13 is L(0,0,0,0,0,0,1,3)
candidate:=number.reduce(fcn(sum,n){ sum*10 + n },0); // digits to int
while(candidate != 89 and candidate != 1) // repeatedly sum squares of digits
{ candidate = candidate.split().reduce(fcn(sum,c){ sum + c*c },0); }
if(candidate == 89){ // count permutations of these digits, they all sum to 89
digitsCount:=List(0,0,0,0,0,0,0,0,0,0);
foreach d in (number){ digitsCount[d] += 1; }
return(digitsCount.reduce(fcn(r,c){ r/factorial(c) },cacheBang)); // cacheBang==number.len()!
}
0 // this number doesn't sum to 89 (ie sums to 1)
} fcn factorial(n) { (1).reduce(n,fcn(N,n){ N*n },1) }
limit:=0d100_000_000; cacheSize:=limit.toFloat().log10().ceil().toInt(); number:=(0).pump(cacheSize,List().write,0); // list of zeros result:=0; i:=cacheSize - 1; var cacheBang=factorial(cacheSize); //== number.len()!
while(True){ // create numbers s.t. no set of digits is repeated
if(i == 0 and number[i] == 9) break;
if(i == cacheSize - 1 and number[i] < 9){ number[i] += 1; result += check(number); }
else if(number[i] == 9) i -= 1;
else{
number[i] += 1;
foreach j in ([i + 1 .. cacheSize - 1]){ number[j] = number[i]; }
i = cacheSize - 1;
result += check(number);
}
} println(result);</lang>
- Output:
85744333