Carmichael lambda function: Difference between revisions
→{{header|jq}}: software engineering |
|||
(12 intermediate revisions by 7 users not shown) | |||
Line 1: | Line 1: | ||
{{draft task}} |
{{draft task|Prime Numbers}} |
||
{{task|Prime Numbers}} |
|||
;Background |
;Background |
||
Line 14: | Line 13: | ||
;Task |
;Task |
||
* Write a function to obtain the Carmichael |
* Write a function to obtain the Carmichael lambda of a positive integer. If the function is supplied by a core library of the language, this also may be used. |
||
* Write a function to count the number of iterations '''k''' of the '''k'''-iterated lambda function needed to get to a value of 1. Show the value of {{math | ''λ''(''n'')}} and the number of '''k'''-iterations for integers from 1 to 25. |
* Write a function to count the number of iterations '''k''' of the '''k'''-iterated lambda function needed to get to a value of 1. Show the value of {{math | ''λ''(''n'')}} and the number of '''k'''-iterations for integers from 1 to 25. |
||
Line 27: | Line 26: | ||
:* [[wp:Carmichael_function|Wikipedia entry for Carmichael function]] |
:* [[wp:Carmichael_function|Wikipedia entry for Carmichael function]] |
||
:* [[oeis:A181776|OEIS lambda( |
:* [[oeis:A181776|OEIS lambda(lambda(n)) function sequence]] |
||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
'''Works with jq, the C implementation of jq''' (*) |
|||
'''Works with gojq, the Go implementation of jq''' (**) |
|||
'''Works with jaq, the Rust implementation of jq''' (***) |
|||
A point of interest in the following implementation is the |
|||
way in which the filter `iteratedToOne/1` prevents contamination of the input object |
|||
by encapsulating it using the pattern: `def F: {"global": .} |.... | .global;` |
|||
(*) The program seemingly made no progress after i == 12 in the second |
|||
table and was terminated. |
|||
(**) gojq's space requirements for this task beyond i == 10 in the |
|||
second table are prohibitive. |
|||
(***) The program shown below is compatible with jaq except that jaq does not |
|||
allow arrays to be grown implicitly by assignment; a simple workaround |
|||
would be to change .cache to be a JSON object. However the space |
|||
requirements effectively prevent computation of the rows in the second |
|||
table beyond the row with i == 11. |
|||
<syntaxhighlight lang="jq"> |
|||
### Generic functions |
|||
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .; |
|||
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in); |
|||
# Emit an array of the prime factors of 'n' in order using a wheel with basis [2, 3, 5] |
|||
# e.g. 44 | primeFactors => [2,2,11] |
|||
# From https://rosettacode.org/wiki/Product_of_min_and_max_prime_factors#jq |
|||
def primeFactors: |
|||
def out($i): until (.n % $i != 0; .factors += [$i] | .n = ((.n/$i)|floor) ); |
|||
if . < 2 then [] |
|||
else [4, 2, 4, 2, 4, 6, 2, 6] as $inc |
|||
| { n: ., |
|||
factors: [] } |
|||
| out(2) |
|||
| out(3) |
|||
| out(5) |
|||
| .k = 7 |
|||
| .i = 0 |
|||
| until(.k * .k > .n; |
|||
if .n % .k == 0 |
|||
then .factors += [.k] |
|||
| .n = ((.n/.k)|floor) |
|||
else .k += $inc[.i] |
|||
| .i = ((.i + 1) % 8) |
|||
end) |
|||
| if .n > 1 then .factors += [ .n ] else . end |
|||
| .factors |
|||
end; |
|||
# Define the helper function to take advantage of jq's tail-recursion optimization |
|||
def lcm($m; $n): |
|||
def _lcm: |
|||
# state is [m, n, i] |
|||
if (.[2] % .[1]) == 0 then .[2] else (.[0:2] + [.[2] + $m]) | _lcm end; |
|||
[$m, $n, $m] | _lcm; |
|||
def lcm(stream): |
|||
reduce stream as $s (1; lcm(.; $s)); |
|||
### Carmichael Lambda Function |
|||
def primePows: |
|||
. as $n |
|||
| { factPows: [] } |
|||
| ($n | primeFactors) as $pf |
|||
| .currFact = $pf[0] |
|||
| .count = 1 |
|||
| reduce $pf[1:][] as $fact (.; |
|||
if $fact != .currFact |
|||
then .factPows += [[.currFact, .count]] |
|||
| .currFact = $fact |
|||
| .count = 1 |
|||
else .count += 1 |
|||
end) |
|||
| .factPows + [[.currFact, .count]] ; |
|||
def phi($p; $r): |
|||
($p | power($r-1)) * ($p - 1); |
|||
# Initial cache values |
|||
def cache: [0, 1, 1, null, 2]; # [_, 1, phi(2; 1), _, phi(2; 2)] |
|||
# input: {cache} |
|||
# output: {cache, result} |
|||
def CarmichaelHelper($p; $r): |
|||
($p|power($r)) as $n |
|||
| .cache[$n] as $cached |
|||
| if $cached then .result=$cached |
|||
else |
|||
if $p > 2 |
|||
then .cache[$n] = phi($p; $r) |
|||
else .cache[$n] = phi($p; $r - 1) |
|||
end |
|||
| .result = .cache[$n] |
|||
end; |
|||
# input: {cache} |
|||
# output: {cache, result, a} |
|||
def CarmichaelLambda($n): |
|||
.cache |= (. // cache) |
|||
| .result = null |
|||
| if $n < 1 then "CarmichaelLambda: argument must be a positive integer (vs\($n))" | error |
|||
else |
|||
if .cache[$n] then .result = .cache[$n] |
|||
else ($n|primePows) as $pps |
|||
| if ($pps|length) == 1 |
|||
then $pps[0][0] as $p |
|||
| $pps[0][1] as $r |
|||
| if $p > 2 then .cache[$n] = phi($p; $r) |
|||
else .cache[$n] = phi($p; $r - 1) |
|||
end |
|||
| .result = .cache[$n] |
|||
else .a = [] |
|||
| reduce $pps[] as $pp (.; |
|||
CarmichaelHelper($pp[0]; $pp[1]) |
|||
| .a += [ .result ] ) |
|||
| .cache[$n] = lcm(.a[]) |
|||
| .result = .cache[$n] |
|||
end |
|||
end |
|||
end ; |
|||
# input: {cache} |
|||
# output: {cache, result} |
|||
def iteratedToOne($j): |
|||
{global: ., |
|||
k: 0, |
|||
$j} |
|||
| until( .j <= 1; |
|||
.j as $j |
|||
| .global |= CarmichaelLambda($j) |
|||
| .j = .global.result |
|||
| .k += 1 ) |
|||
| .global.result = .k |
|||
| .global; |
|||
def task1($m): |
|||
" n λ k", |
|||
"----------", |
|||
(foreach range(1; 1+$m) as $n ({}; |
|||
CarmichaelLambda($n) |
|||
| .result as $lambda |
|||
| iteratedToOne($n) |
|||
| .result as $k |
|||
| .emit = "\($n|lpad(2)) \($lambda|lpad(2)) \($k|lpad(2))" ) |
|||
| .emit ) ; |
|||
def task2($maxI; $maxN): |
|||
"\nIterations to 1 i lambda(i)", |
|||
"=====================================", |
|||
" 0 1 1", |
|||
( . // {} |
|||
| .cache |= (. // cache) |
|||
| .found = [true, (range(0; $maxN)|false)] |
|||
| .i=1 |
|||
| .n=null |
|||
| while (.i <= $maxI and .n != $maxN; |
|||
.emit = null |
|||
| iteratedToOne(.i) |
|||
| .n = .result |
|||
| if (.found[.n]|not) |
|||
then .found[.n] = true |
|||
| CarmichaelLambda(.i) |
|||
| .result as $lambda |
|||
| .emit = "\(.n|lpad(4)) \(.i|lpad(18)) \($lambda|lpad(12))" |
|||
end |
|||
| .i += 1 ) |
|||
| select(.emit).emit ); |
|||
task1(25), |
|||
"", |
|||
task2(5*1e6; 15) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
n λ k |
|||
---------- |
|||
1 1 0 |
|||
2 1 1 |
|||
3 2 2 |
|||
4 2 2 |
|||
5 4 3 |
|||
6 2 2 |
|||
7 6 3 |
|||
8 2 2 |
|||
9 6 3 |
|||
10 4 3 |
|||
11 10 4 |
|||
12 2 2 |
|||
13 12 3 |
|||
14 6 3 |
|||
15 4 3 |
|||
16 4 3 |
|||
17 16 4 |
|||
18 6 3 |
|||
19 18 4 |
|||
20 4 3 |
|||
21 6 3 |
|||
22 10 4 |
|||
23 22 5 |
|||
24 2 2 |
|||
25 20 4 |
|||
Iterations to 1 i lambda(i) |
|||
===================================== |
|||
0 1 1 |
|||
1 2 1 |
|||
2 3 2 |
|||
3 5 4 |
|||
4 11 10 |
|||
5 23 22 |
|||
6 47 46 |
|||
7 283 282 |
|||
8 719 718 |
|||
9 1439 1438 |
|||
10 2879 2878 |
|||
11 34549 34548 |
|||
12 138197 138196 |
|||
[...terminated...] |
|||
</pre> |
|||
=={{header|Julia}}== |
|||
<syntaxhighlight lang="julia">using DelimitedFiles |
|||
using Primes |
|||
const maxcached = typemax(Int32) |
|||
const carmicache = zeros(Int32, maxcached) # NB: 8 gigabytes memory used for this size cache |
|||
const rlock = ReentrantLock() |
|||
""" Carmichael reduced totient function lambda(n) """ |
|||
function lambda(n::Integer) |
|||
@assert n > 0 |
|||
n <= maxcached && carmicache[n] > 0 && return Int(carmicache[n]) |
|||
lam = 1 |
|||
for (p, e) in factor(n) |
|||
if p == 2 && e > 2 |
|||
lam = lcm(lam, 2^(e - 2)) |
|||
else |
|||
lam = lcm(lam, (p - 1) * p^(e - 1)) |
|||
end |
|||
end |
|||
if n <= maxcached |
|||
lock(rlock) |
|||
carmicache[n] = lam |
|||
unlock(rlock) |
|||
end |
|||
return lam |
|||
end |
|||
""" |
|||
Return k for the k-fold iterated lambda function, where k is count of iterations to 1 |
|||
""" |
|||
function iterated_lambdas_to_one(i) |
|||
k = 0 |
|||
while i > 1 |
|||
i = lambda(i) |
|||
k += 1 |
|||
end |
|||
return k |
|||
end |
|||
""" |
|||
Calculate values for the task. Due to a runtime of several days for a sequence length |
|||
of 25, the function switches to multithreading for values over about 2 billion. Note |
|||
that (empirically) adjacent pairs of values within the sequence above 100000 are less |
|||
than a factor of 6 apart, so each next large value is sought within that range. |
|||
""" |
|||
function print_iterated_lambdas(upto = 25) |
|||
println("Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25:") |
|||
firsts = zeros(Int, upto) |
|||
if isfile("lambdas.txt") |
|||
precalc = readdlm("lambdas.txt", Int64) |
|||
if length(precalc) == upto |
|||
firsts .= precalc # restore work done if resterted (since it takes days) |
|||
end |
|||
end |
|||
for i = 1:25 |
|||
lam = lambda(i) |
|||
k = iterated_lambdas_to_one(i) |
|||
print("($i, $lam, $k)", i % 5 == 0 ? "\n" : " ") |
|||
end |
|||
target = something(findlast(!iszero, firsts), 1) |
|||
if firsts[target] < maxcached ÷ 7 # get first ones with single thread |
|||
for i = 2:maxcached |
|||
n = iterated_lambdas_to_one(i) |
|||
if n <= upto && (firsts[n] == 0 || firsts[n] > i) |
|||
firsts[n] = i |
|||
end |
|||
end |
|||
end |
|||
target = something(findlast(!iszero, firsts), 1) + 1 |
|||
while target <= upto # multithreaded for larger targets |
|||
startval = firsts[target - 1] + 1 |
|||
for j in startval:startval:startval*6 |
|||
j += iseven(j) # make odd |
|||
@Threads.threads for i in j:2:j+startval-1 # should not be even if i > 2 |
|||
n = iterated_lambdas_to_one(i) |
|||
if n == target |
|||
if firsts[n] == 0 || firsts[n] > i |
|||
try |
|||
lock(rlock) |
|||
firsts[n] = i |
|||
writedlm("lambdas.txt", firsts) # save work done for restarts |
|||
unlock(rlock) |
|||
catch y |
|||
unlock(rlock) |
|||
@warn y |
|||
rethrow() |
|||
end |
|||
end |
|||
break |
|||
end |
|||
firsts[target] > 0 && firsts[target] < i && break |
|||
end |
|||
end |
|||
target += 1 |
|||
end |
|||
println("\nIterations to 1 n lambda(n)\n=====================================") |
|||
println(" 0 1 1") |
|||
for n = 1:upto |
|||
println(lpad(n, 4) * lpad(firsts[n], 17) * lpad(lambda(firsts[n]), 14)) |
|||
end |
|||
end |
|||
print_iterated_lambdas() |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25: |
|||
(1, 1, 0) (2, 1, 1) (3, 2, 2) (4, 2, 2) (5, 4, 3) |
|||
(6, 2, 2) (7, 6, 3) (8, 2, 2) (9, 6, 3) (10, 4, 3) |
|||
(11, 10, 4) (12, 2, 2) (13, 12, 3) (14, 6, 3) (15, 4, 3) |
|||
(16, 4, 3) (17, 16, 4) (18, 6, 3) (19, 18, 4) (20, 4, 3) |
|||
(21, 6, 3) (22, 10, 4) (23, 22, 5) (24, 2, 2) (25, 20, 4) |
|||
Iterations to 1 n lambda(n) |
|||
===================================== |
|||
0 1 1 |
|||
1 2 1 |
|||
2 3 2 |
|||
3 5 4 |
|||
4 11 10 |
|||
5 23 22 |
|||
6 47 46 |
|||
7 283 282 |
|||
8 719 718 |
|||
9 1439 1438 |
|||
10 2879 2878 |
|||
11 34549 34548 |
|||
12 138197 138196 |
|||
13 531441 354294 |
|||
14 1594323 1062882 |
|||
15 4782969 3188646 |
|||
16 14348907 9565938 |
|||
17 43046721 28697814 |
|||
18 86093443 86093442 |
|||
19 258280327 258280326 |
|||
20 688747547 688747546 |
|||
21 3486784401 2324522934 |
|||
22 10460353203 6973568802 |
|||
23 31381059609 20920706406 |
|||
24 94143178827 62762119218 |
|||
25 282429536481 188286357654 |
|||
</pre> |
|||
=={{header|Mathematica}}|{{header|Wolfram Language}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="Mathematica"> |
|||
IteratedToOne[i_] := Module[{k = 0, iter = i}, |
|||
While[iter > 1, iter = CarmichaelLambda[iter]; |
|||
k++;]; |
|||
k] |
|||
Print["Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25:"]; |
|||
Do[lam = CarmichaelLambda[i]; |
|||
k = IteratedToOne[i]; |
|||
If[Mod[i, 5] == 0, Print[{i, lam, k}], Print[{i, lam, k}, " "]], {i, 1, 25}] |
|||
upTo = 20; |
|||
maxToTest = 100000000000; |
|||
firsts = Table[0, {upTo + 1}]; |
|||
firsts[[1]] = 1; |
|||
Print["\nIterations to 1 n lambda(n)\n=================================="]; |
|||
Print[StringForm["`` `` ``", 0, 1, 1]]; |
|||
Do[n = IteratedToOne[i]; |
|||
If[0 < n <= upTo && firsts[[n + 1]] == 0, firsts[[n + 1]] = i; |
|||
j = If[n < 1, 0, CarmichaelLambda[i]]; |
|||
Print[StringForm["`` `` ``", n, i, j]]; |
|||
If[n >= upTo, Break[]];], {i, 2, maxToTest}] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25: |
|||
{1, 1, 0} |
|||
{2, 1, 1} |
|||
{3, 2, 2} |
|||
{4, 2, 2} |
|||
{5, 4, 3} |
|||
{6, 2, 2} |
|||
{7, 6, 3} |
|||
{8, 2, 2} |
|||
{9, 6, 3} |
|||
{10, 4, 3} |
|||
{11, 10, 4} |
|||
{12, 2, 2} |
|||
{13, 12, 3} |
|||
{14, 6, 3} |
|||
{15, 4, 3} |
|||
{16, 4, 3} |
|||
{17, 16, 4} |
|||
{18, 6, 3} |
|||
{19, 18, 4} |
|||
{20, 4, 3} |
|||
{21, 6, 3} |
|||
{22, 10, 4} |
|||
{23, 22, 5} |
|||
{24, 2, 2} |
|||
{25, 20, 4} |
|||
Iterations to 1 n lambda(n) |
|||
================================== |
|||
0 1 1 |
|||
1 2 1 |
|||
2 3 2 |
|||
3 5 4 |
|||
4 11 10 |
|||
5 23 22 |
|||
6 47 46 |
|||
7 283 282 |
|||
8 719 718 |
|||
9 1439 1438 |
|||
10 2879 2878 |
|||
11 34549 34548 |
|||
12 138197 138196 |
|||
13 531441 354294 |
|||
14 1594323 1062882 |
|||
15 4782969 3188646 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Line 101: | Line 544: | ||
<pre> |
<pre> |
||
16 14,348,907 9,565,938 |
16 14,348,907 9,565,938 |
||
17 43,046,721 28,697,814 |
|||
tbc... |
|||
</pre> |
</pre> |
||
I estimate that took ~47mins, at which point I killed it.<br> |
|||
The Python code took 13mins 58 secs to get to 15, on the same box. |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
Line 173: | Line 618: | ||
19 258280327 258280326 |
19 258280327 258280326 |
||
20 688747547 688747546 |
20 688747547 688747546 |
||
</pre> |
|||
=={{header|Raku}}== |
|||
{{trans|Wren}} |
|||
<syntaxhighlight lang="raku" line># 20240228 Raku programming solution |
|||
use Prime::Factor; |
|||
sub phi(Int $p, Int $r) { return $p**($r - 1) * ($p - 1) } |
|||
sub CarmichaelLambda(Int $n) { |
|||
state %cache = 1 => 1, 2 => 1, 4 => 2; |
|||
sub CarmichaelHelper(Int $p, Int $r) { |
|||
my Int $n = $p ** $r; |
|||
return %cache{$n} if %cache{$n}:exists; |
|||
return %cache{$n} = $p > 2 ?? phi($p, $r) !! phi($p, $r - 1) |
|||
} |
|||
if $n < 1 { die "'n' must be a positive integer." } |
|||
return %cache{$n} if %cache{$n}:exists; |
|||
if ( my %pps = prime-factors($n).Bag ).elems == 1 { |
|||
my ($p, $r) = %pps.kv>>.Int; |
|||
return %cache{$n} = $p > 2 ?? phi($p, $r) !! phi($p, $r - 1) |
|||
} |
|||
return [lcm] %pps.kv.map: -> $k, $v { CarmichaelHelper($k.Int, $v) } |
|||
} |
|||
sub iteratedToOne($i is copy) { |
|||
my $k = 0; |
|||
while $i > 1 { $i = CarmichaelLambda($i) andthen $k++ } |
|||
return $k |
|||
} |
|||
say " n λ k"; |
|||
say "----------"; |
|||
for 1..25 -> $n { |
|||
printf "%2d %2d %2d\n", $n, CarmichaelLambda($n), iteratedToOne($n) |
|||
} |
|||
say "\nIterations to 1 i lambda(i)"; |
|||
say "====================================="; |
|||
say " 0 1 1"; |
|||
my ($maxI, $maxN) = 5e6, 10; # for N=15, takes around 47 minutes with an i5-10500T |
|||
my @found = True, |( False xx $maxN ); |
|||
for 1 .. $maxI -> $i { |
|||
unless @found[ my $n = iteratedToOne($i) ] { |
|||
printf "%4d %18d %12d\n", $n, $i, CarmichaelLambda($i); |
|||
@found[$n] = True andthen ( last if $n == $maxN ) |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} Same as Wren example .. and it is probably 30X times slower 😅 |
|||
=={{header|Sidef}}== |
|||
Built-in as '''Number#carmichael_lambda''': |
|||
<syntaxhighlight lang="ruby">func iteratedToOne(n) { |
|||
var k = 0 |
|||
while (n > 1) { |
|||
n.carmichael_lambda! |
|||
++k |
|||
} |
|||
return k |
|||
} |
|||
say " n λ k" |
|||
say "----------" |
|||
for n in (1..25) { |
|||
printf("%2d %2d %2d\n", n, n.carmichael_lambda, iteratedToOne(n)) |
|||
} |
|||
say "\nIterations to 1 i lambda(i)" |
|||
say "=====================================" |
|||
var table = [] |
|||
1..Inf -> each {|k| |
|||
var n = iteratedToOne(k) |
|||
if (!table[n]) { |
|||
table[n] = true |
|||
printf("%4s %18s %12s\n", n, k, k.carmichael_lambda) |
|||
break if (n == 15) |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
n λ k |
|||
---------- |
|||
1 1 0 |
|||
2 1 1 |
|||
3 2 2 |
|||
4 2 2 |
|||
5 4 3 |
|||
6 2 2 |
|||
7 6 3 |
|||
8 2 2 |
|||
9 6 3 |
|||
10 4 3 |
|||
11 10 4 |
|||
12 2 2 |
|||
13 12 3 |
|||
14 6 3 |
|||
15 4 3 |
|||
16 4 3 |
|||
17 16 4 |
|||
18 6 3 |
|||
19 18 4 |
|||
20 4 3 |
|||
21 6 3 |
|||
22 10 4 |
|||
23 22 5 |
|||
24 2 2 |
|||
25 20 4 |
|||
Iterations to 1 i lambda(i) |
|||
===================================== |
|||
0 1 1 |
|||
1 2 1 |
|||
2 3 2 |
|||
3 5 4 |
|||
4 11 10 |
|||
5 23 22 |
|||
6 47 46 |
|||
7 283 282 |
|||
8 719 718 |
|||
9 1439 1438 |
|||
10 2879 2878 |
|||
11 34549 34548 |
|||
12 138197 138196 |
|||
13 531441 354294 |
|||
14 1594323 1062882 |
|||
15 4782969 3188646 |
|||
</pre> |
</pre> |
||
Revision as of 19:07, 11 July 2024
- Background
The Carmichael function, or Carmichael lambda function, is a function in number theory. The Carmichael lambda λ(n) of a positive integer n is the smallest positive integer m such that
holds for every integer coprime to n.
The Carmichael lambda function can be iterated, that is, called repeatedly on its result. If this iteration is
performed k times, this is considered the k-iterated Carmichael lambda function. Thus, λ(λ(n))
would be the 2-iterated lambda function. With repeated, sufficiently large but finite k-iterations, the iterated Carmichael function will eventually compute as 1 for all positive integers n.
- Task
- Write a function to obtain the Carmichael lambda of a positive integer. If the function is supplied by a core library of the language, this also may be used.
- Write a function to count the number of iterations k of the k-iterated lambda function needed to get to a value of 1. Show the value of λ(n) and the number of k-iterations for integers from 1 to 25.
- Find the lowest integer for which the value of iterated k is i, for i from 1 to 15.
- Stretch task (an extension of the third task above)
- Find, additionally, for i from 16 to 25, the lowest integer n for which the number of iterations k of the Carmichael lambda function from n to get to 1 is i.
- See also
jq
Adapted from Wren
Works with jq, the C implementation of jq (*)
Works with gojq, the Go implementation of jq (**)
Works with jaq, the Rust implementation of jq (***)
A point of interest in the following implementation is the way in which the filter `iteratedToOne/1` prevents contamination of the input object by encapsulating it using the pattern: `def F: {"global": .} |.... | .global;`
(*) The program seemingly made no progress after i == 12 in the second table and was terminated.
(**) gojq's space requirements for this task beyond i == 10 in the second table are prohibitive.
(***) The program shown below is compatible with jaq except that jaq does not allow arrays to be grown implicitly by assignment; a simple workaround would be to change .cache to be a JSON object. However the space requirements effectively prevent computation of the rows in the second table beyond the row with i == 11.
### Generic functions
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
# Emit an array of the prime factors of 'n' in order using a wheel with basis [2, 3, 5]
# e.g. 44 | primeFactors => [2,2,11]
# From https://rosettacode.org/wiki/Product_of_min_and_max_prime_factors#jq
def primeFactors:
def out($i): until (.n % $i != 0; .factors += [$i] | .n = ((.n/$i)|floor) );
if . < 2 then []
else [4, 2, 4, 2, 4, 6, 2, 6] as $inc
| { n: .,
factors: [] }
| out(2)
| out(3)
| out(5)
| .k = 7
| .i = 0
| until(.k * .k > .n;
if .n % .k == 0
then .factors += [.k]
| .n = ((.n/.k)|floor)
else .k += $inc[.i]
| .i = ((.i + 1) % 8)
end)
| if .n > 1 then .factors += [ .n ] else . end
| .factors
end;
# Define the helper function to take advantage of jq's tail-recursion optimization
def lcm($m; $n):
def _lcm:
# state is [m, n, i]
if (.[2] % .[1]) == 0 then .[2] else (.[0:2] + [.[2] + $m]) | _lcm end;
[$m, $n, $m] | _lcm;
def lcm(stream):
reduce stream as $s (1; lcm(.; $s));
### Carmichael Lambda Function
def primePows:
. as $n
| { factPows: [] }
| ($n | primeFactors) as $pf
| .currFact = $pf[0]
| .count = 1
| reduce $pf[1:][] as $fact (.;
if $fact != .currFact
then .factPows += [[.currFact, .count]]
| .currFact = $fact
| .count = 1
else .count += 1
end)
| .factPows + [[.currFact, .count]] ;
def phi($p; $r):
($p | power($r-1)) * ($p - 1);
# Initial cache values
def cache: [0, 1, 1, null, 2]; # [_, 1, phi(2; 1), _, phi(2; 2)]
# input: {cache}
# output: {cache, result}
def CarmichaelHelper($p; $r):
($p|power($r)) as $n
| .cache[$n] as $cached
| if $cached then .result=$cached
else
if $p > 2
then .cache[$n] = phi($p; $r)
else .cache[$n] = phi($p; $r - 1)
end
| .result = .cache[$n]
end;
# input: {cache}
# output: {cache, result, a}
def CarmichaelLambda($n):
.cache |= (. // cache)
| .result = null
| if $n < 1 then "CarmichaelLambda: argument must be a positive integer (vs\($n))" | error
else
if .cache[$n] then .result = .cache[$n]
else ($n|primePows) as $pps
| if ($pps|length) == 1
then $pps[0][0] as $p
| $pps[0][1] as $r
| if $p > 2 then .cache[$n] = phi($p; $r)
else .cache[$n] = phi($p; $r - 1)
end
| .result = .cache[$n]
else .a = []
| reduce $pps[] as $pp (.;
CarmichaelHelper($pp[0]; $pp[1])
| .a += [ .result ] )
| .cache[$n] = lcm(.a[])
| .result = .cache[$n]
end
end
end ;
# input: {cache}
# output: {cache, result}
def iteratedToOne($j):
{global: .,
k: 0,
$j}
| until( .j <= 1;
.j as $j
| .global |= CarmichaelLambda($j)
| .j = .global.result
| .k += 1 )
| .global.result = .k
| .global;
def task1($m):
" n λ k",
"----------",
(foreach range(1; 1+$m) as $n ({};
CarmichaelLambda($n)
| .result as $lambda
| iteratedToOne($n)
| .result as $k
| .emit = "\($n|lpad(2)) \($lambda|lpad(2)) \($k|lpad(2))" )
| .emit ) ;
def task2($maxI; $maxN):
"\nIterations to 1 i lambda(i)",
"=====================================",
" 0 1 1",
( . // {}
| .cache |= (. // cache)
| .found = [true, (range(0; $maxN)|false)]
| .i=1
| .n=null
| while (.i <= $maxI and .n != $maxN;
.emit = null
| iteratedToOne(.i)
| .n = .result
| if (.found[.n]|not)
then .found[.n] = true
| CarmichaelLambda(.i)
| .result as $lambda
| .emit = "\(.n|lpad(4)) \(.i|lpad(18)) \($lambda|lpad(12))"
end
| .i += 1 )
| select(.emit).emit );
task1(25),
"",
task2(5*1e6; 15)
- Output:
n λ k ---------- 1 1 0 2 1 1 3 2 2 4 2 2 5 4 3 6 2 2 7 6 3 8 2 2 9 6 3 10 4 3 11 10 4 12 2 2 13 12 3 14 6 3 15 4 3 16 4 3 17 16 4 18 6 3 19 18 4 20 4 3 21 6 3 22 10 4 23 22 5 24 2 2 25 20 4 Iterations to 1 i lambda(i) ===================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1439 1438 10 2879 2878 11 34549 34548 12 138197 138196 [...terminated...]
Julia
using DelimitedFiles
using Primes
const maxcached = typemax(Int32)
const carmicache = zeros(Int32, maxcached) # NB: 8 gigabytes memory used for this size cache
const rlock = ReentrantLock()
""" Carmichael reduced totient function lambda(n) """
function lambda(n::Integer)
@assert n > 0
n <= maxcached && carmicache[n] > 0 && return Int(carmicache[n])
lam = 1
for (p, e) in factor(n)
if p == 2 && e > 2
lam = lcm(lam, 2^(e - 2))
else
lam = lcm(lam, (p - 1) * p^(e - 1))
end
end
if n <= maxcached
lock(rlock)
carmicache[n] = lam
unlock(rlock)
end
return lam
end
"""
Return k for the k-fold iterated lambda function, where k is count of iterations to 1
"""
function iterated_lambdas_to_one(i)
k = 0
while i > 1
i = lambda(i)
k += 1
end
return k
end
"""
Calculate values for the task. Due to a runtime of several days for a sequence length
of 25, the function switches to multithreading for values over about 2 billion. Note
that (empirically) adjacent pairs of values within the sequence above 100000 are less
than a factor of 6 apart, so each next large value is sought within that range.
"""
function print_iterated_lambdas(upto = 25)
println("Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25:")
firsts = zeros(Int, upto)
if isfile("lambdas.txt")
precalc = readdlm("lambdas.txt", Int64)
if length(precalc) == upto
firsts .= precalc # restore work done if resterted (since it takes days)
end
end
for i = 1:25
lam = lambda(i)
k = iterated_lambdas_to_one(i)
print("($i, $lam, $k)", i % 5 == 0 ? "\n" : " ")
end
target = something(findlast(!iszero, firsts), 1)
if firsts[target] < maxcached ÷ 7 # get first ones with single thread
for i = 2:maxcached
n = iterated_lambdas_to_one(i)
if n <= upto && (firsts[n] == 0 || firsts[n] > i)
firsts[n] = i
end
end
end
target = something(findlast(!iszero, firsts), 1) + 1
while target <= upto # multithreaded for larger targets
startval = firsts[target - 1] + 1
for j in startval:startval:startval*6
j += iseven(j) # make odd
@Threads.threads for i in j:2:j+startval-1 # should not be even if i > 2
n = iterated_lambdas_to_one(i)
if n == target
if firsts[n] == 0 || firsts[n] > i
try
lock(rlock)
firsts[n] = i
writedlm("lambdas.txt", firsts) # save work done for restarts
unlock(rlock)
catch y
unlock(rlock)
@warn y
rethrow()
end
end
break
end
firsts[target] > 0 && firsts[target] < i && break
end
end
target += 1
end
println("\nIterations to 1 n lambda(n)\n=====================================")
println(" 0 1 1")
for n = 1:upto
println(lpad(n, 4) * lpad(firsts[n], 17) * lpad(lambda(firsts[n]), 14))
end
end
print_iterated_lambdas()
- Output:
Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25: (1, 1, 0) (2, 1, 1) (3, 2, 2) (4, 2, 2) (5, 4, 3) (6, 2, 2) (7, 6, 3) (8, 2, 2) (9, 6, 3) (10, 4, 3) (11, 10, 4) (12, 2, 2) (13, 12, 3) (14, 6, 3) (15, 4, 3) (16, 4, 3) (17, 16, 4) (18, 6, 3) (19, 18, 4) (20, 4, 3) (21, 6, 3) (22, 10, 4) (23, 22, 5) (24, 2, 2) (25, 20, 4) Iterations to 1 n lambda(n) ===================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1439 1438 10 2879 2878 11 34549 34548 12 138197 138196 13 531441 354294 14 1594323 1062882 15 4782969 3188646 16 14348907 9565938 17 43046721 28697814 18 86093443 86093442 19 258280327 258280326 20 688747547 688747546 21 3486784401 2324522934 22 10460353203 6973568802 23 31381059609 20920706406 24 94143178827 62762119218 25 282429536481 188286357654
Mathematica |Wolfram Language
IteratedToOne[i_] := Module[{k = 0, iter = i},
While[iter > 1, iter = CarmichaelLambda[iter];
k++;];
k]
Print["Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25:"];
Do[lam = CarmichaelLambda[i];
k = IteratedToOne[i];
If[Mod[i, 5] == 0, Print[{i, lam, k}], Print[{i, lam, k}, " "]], {i, 1, 25}]
upTo = 20;
maxToTest = 100000000000;
firsts = Table[0, {upTo + 1}];
firsts[[1]] = 1;
Print["\nIterations to 1 n lambda(n)\n=================================="];
Print[StringForm["`` `` ``", 0, 1, 1]];
Do[n = IteratedToOne[i];
If[0 < n <= upTo && firsts[[n + 1]] == 0, firsts[[n + 1]] = i;
j = If[n < 1, 0, CarmichaelLambda[i]];
Print[StringForm["`` `` ``", n, i, j]];
If[n >= upTo, Break[]];], {i, 2, maxToTest}]
- Output:
Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25: {1, 1, 0} {2, 1, 1} {3, 2, 2} {4, 2, 2} {5, 4, 3} {6, 2, 2} {7, 6, 3} {8, 2, 2} {9, 6, 3} {10, 4, 3} {11, 10, 4} {12, 2, 2} {13, 12, 3} {14, 6, 3} {15, 4, 3} {16, 4, 3} {17, 16, 4} {18, 6, 3} {19, 18, 4} {20, 4, 3} {21, 6, 3} {22, 10, 4} {23, 22, 5} {24, 2, 2} {25, 20, 4} Iterations to 1 n lambda(n) ================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1439 1438 10 2879 2878 11 34549 34548 12 138197 138196 13 531441 354294 14 1594323 1062882 15 4782969 3188646
Phix
with javascript_semantics
function lambda(sequence pk)
integer {p,k} = pk, res = phi(power(p,k))
return iff(p=2 and k>=3 ? res/2 : res)
end function
function reduced_totient(integer n)
sequence p = prime_powers(n)
for i,pk in p do
p[i] = lambda(pk)
end for
return lcm(p)
-- alt: neater but ~2* slower
-- return lcm(apply(prime_powers(n),lambda))
end function
function k_iter(integer i)
return iff(i>1?1+k_iter(reduced_totient(i)):0)
end function
function ident(integer i) return i end function
sequence n = tagset(25)
for i,d in {"n","eulers_totient","reduced_totient","k"} do
sequence r = apply(n,{ident,phi,reduced_totient,k_iter}[i])
printf(1,"%17s: %s\n",{d,join(r," ",fmt:="%2d")})
end for
atom t0 = time()
printf(1,"\n k i lambda(i)\n")
printf(1, "==========================\n")
integer i = 1
for k=0 to iff(platform()=JS?12:15) do -- (keep JS < 3s)
while k_iter(i)!=k do i += 1 end while
printf(1,"%2d %,9d %,9d\n",{k,i,reduced_totient(i)})
end for
?elapsed(time()-t0)
- Output:
n: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 eulers_totient: 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 reduced_totient: 1 1 2 2 4 2 6 2 6 4 10 2 12 6 4 4 16 6 18 4 6 10 22 2 20 k: 0 1 2 2 3 2 3 2 3 3 4 2 3 3 3 3 4 3 4 3 3 4 5 2 4 k i lambda(i) ========================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1,439 1,438 10 2,879 2,878 11 34,549 34,548 12 138,197 138,196 13 531,441 354,294 14 1,594,323 1,062,882 15 4,782,969 3,188,646 "2 minutes and 27s"
stretch
16 14,348,907 9,565,938 17 43,046,721 28,697,814
I estimate that took ~47mins, at which point I killed it.
The Python code took 13mins 58 secs to get to 15, on the same box.
Python
Python has the Carmichael function in the SymPy library, where it is called reduced_totient.
from sympy import reduced_totient
def iterated_to_one(i):
""" return k for the k-fold iterated lambda function where k is the first time iteration reaches 1 """
k = 0
while i > 1:
i = reduced_totient(i)
k += 1
return k
if __name__ == "__main__":
print("Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25:")
for i in range(1, 26):
lam = reduced_totient(i)
k = iterated_to_one(i)
print(f'({i}, {lam}, {k})', end = '\n' if i % 5 == 0 else ' ')
UP_TO = 20
MAX_TO_TEST = 100_000_000_000
FIRSTS = [0] * (UP_TO + 1)
FIRSTS[0] = 1
print('\nIterations to 1 n lambda(n)\n==================================')
print(' 0 1 1')
for i in range(2, MAX_TO_TEST):
n = iterated_to_one(i)
if 0 < n <= UP_TO and FIRSTS[n] == 0:
FIRSTS[n] = i
j = 0 if n < 1 else int(reduced_totient(i))
print(f"{n:>4}{i:>17}{j:>11}")
if n >= UP_TO:
break
- Output:
Listing of (n, lambda(n), k for iteration to 1) for integers from 1 to 25: (1, 1, 0) (2, 1, 1) (3, 2, 2) (4, 2, 2) (5, 4, 3) (6, 2, 2) (7, 6, 3) (8, 2, 2) (9, 6, 3) (10, 4, 3) (11, 10, 4) (12, 2, 2) (13, 12, 3) (14, 6, 3) (15, 4, 3) (16, 4, 3) (17, 16, 4) (18, 6, 3) (19, 18, 4) (20, 4, 3) (21, 6, 3) (22, 10, 4) (23, 22, 5) (24, 2, 2) (25, 20, 4) Iterations to 1 n lambda(n) ================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1439 1438 10 2879 2878 11 34549 34548 12 138197 138196 13 531441 354294 14 1594323 1062882 15 4782969 3188646 16 14348907 9565938 17 43046721 28697814 18 86093443 86093442 19 258280327 258280326 20 688747547 688747546
Raku
# 20240228 Raku programming solution
use Prime::Factor;
sub phi(Int $p, Int $r) { return $p**($r - 1) * ($p - 1) }
sub CarmichaelLambda(Int $n) {
state %cache = 1 => 1, 2 => 1, 4 => 2;
sub CarmichaelHelper(Int $p, Int $r) {
my Int $n = $p ** $r;
return %cache{$n} if %cache{$n}:exists;
return %cache{$n} = $p > 2 ?? phi($p, $r) !! phi($p, $r - 1)
}
if $n < 1 { die "'n' must be a positive integer." }
return %cache{$n} if %cache{$n}:exists;
if ( my %pps = prime-factors($n).Bag ).elems == 1 {
my ($p, $r) = %pps.kv>>.Int;
return %cache{$n} = $p > 2 ?? phi($p, $r) !! phi($p, $r - 1)
}
return [lcm] %pps.kv.map: -> $k, $v { CarmichaelHelper($k.Int, $v) }
}
sub iteratedToOne($i is copy) {
my $k = 0;
while $i > 1 { $i = CarmichaelLambda($i) andthen $k++ }
return $k
}
say " n λ k";
say "----------";
for 1..25 -> $n {
printf "%2d %2d %2d\n", $n, CarmichaelLambda($n), iteratedToOne($n)
}
say "\nIterations to 1 i lambda(i)";
say "=====================================";
say " 0 1 1";
my ($maxI, $maxN) = 5e6, 10; # for N=15, takes around 47 minutes with an i5-10500T
my @found = True, |( False xx $maxN );
for 1 .. $maxI -> $i {
unless @found[ my $n = iteratedToOne($i) ] {
printf "%4d %18d %12d\n", $n, $i, CarmichaelLambda($i);
@found[$n] = True andthen ( last if $n == $maxN )
}
}
- Output:
Same as Wren example .. and it is probably 30X times slower 😅
Sidef
Built-in as Number#carmichael_lambda:
func iteratedToOne(n) {
var k = 0
while (n > 1) {
n.carmichael_lambda!
++k
}
return k
}
say " n λ k"
say "----------"
for n in (1..25) {
printf("%2d %2d %2d\n", n, n.carmichael_lambda, iteratedToOne(n))
}
say "\nIterations to 1 i lambda(i)"
say "====================================="
var table = []
1..Inf -> each {|k|
var n = iteratedToOne(k)
if (!table[n]) {
table[n] = true
printf("%4s %18s %12s\n", n, k, k.carmichael_lambda)
break if (n == 15)
}
}
- Output:
n λ k ---------- 1 1 0 2 1 1 3 2 2 4 2 2 5 4 3 6 2 2 7 6 3 8 2 2 9 6 3 10 4 3 11 10 4 12 2 2 13 12 3 14 6 3 15 4 3 16 4 3 17 16 4 18 6 3 19 18 4 20 4 3 21 6 3 22 10 4 23 22 5 24 2 2 25 20 4 Iterations to 1 i lambda(i) ===================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1439 1438 10 2879 2878 11 34549 34548 12 138197 138196 13 531441 354294 14 1594323 1062882 15 4782969 3188646
Wren
Takes about 88 seconds on my machine (Core i7) to get up to n = 15 for the third part of the task, so haven't attempted the stretch goal.
import "./math" for Int
import "./fmt" for Fmt
var primePows = Fn.new { |n|
var factPows = []
var pf = Int.primeFactors(n)
var currFact = pf[0]
var count = 1
for (fact in pf.skip(1)) {
if (fact != currFact) {
factPows.add([currFact, count])
currFact = fact
count = 1
} else {
count = count + 1
}
}
factPows.add([currFact, count])
return factPows
}
var phi = Fn.new { |p, r| p.pow(r-1) * (p - 1) }
var cache = { 1: 1, 2: phi.call(2, 1), 4: phi.call(2, 2) }
var CarmichaelHelper = Fn.new { |p, r|
var n = p.pow(r)
if (cache.containsKey(n)) return cache[n]
if (p > 2) return cache[n] = phi.call(p, r)
return cache[n] = phi.call(p, r - 1)
}
var CarmichaelLambda = Fn.new { |n|
if (n < 1) Fiber.abort("'n' must be a positive integer.")
if (cache.containsKey(n)) return cache[n]
var pps = primePows.call(n)
if (pps.count == 1) {
var p = pps[0][0]
var r = pps[0][1]
if (p > 2) return cache[n] = phi.call(p, r)
return cache[n] = phi.call(p, r - 1)
}
var a = []
for (pp in pps) a.add(CarmichaelHelper.call(pp[0], pp[1]))
return cache[n] = Int.lcm(a)
}
var iteratedToOne = Fn.new { |i|
var k = 0
while (i > 1) {
i = CarmichaelLambda.call(i)
k = k + 1
}
return k
}
System.print(" n λ k")
System.print("----------")
for (n in 1..25) {
var lambda = CarmichaelLambda.call(n)
var k = iteratedToOne.call(n)
Fmt.print("$2d $2d $2d", n, lambda, k)
}
System.print("\nIterations to 1 i lambda(i)")
System.print("=====================================")
System.print(" 0 1 1")
var maxI = 5 * 1e6
var maxN = 15
var found = List.filled(maxN + 1, false)
found[0] = true
var i = 1
while (i <= maxI) {
var n = iteratedToOne.call(i)
if (!found[n]) {
found[n] = true
var lambda = CarmichaelLambda.call(i)
Fmt.print("$4d $,18d $,12d", n, i, lambda)
if (n == maxN) break
}
i = i + 1
}
- Output:
n λ k ---------- 1 1 0 2 1 1 3 2 2 4 2 2 5 4 3 6 2 2 7 6 3 8 2 2 9 6 3 10 4 3 11 10 4 12 2 2 13 12 3 14 6 3 15 4 3 16 4 3 17 16 4 18 6 3 19 18 4 20 4 3 21 6 3 22 10 4 23 22 5 24 2 2 25 20 4 Iterations to 1 i lambda(i) ===================================== 0 1 1 1 2 1 2 3 2 3 5 4 4 11 10 5 23 22 6 47 46 7 283 282 8 719 718 9 1,439 1,438 10 2,879 2,878 11 34,549 34,548 12 138,197 138,196 13 531,441 354,294 14 1,594,323 1,062,882 15 4,782,969 3,188,646