Addition-chain exponentiation: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) m (Automated syntax highlighting fixup (second round - minor fixes)) |
m (julia example) |
||
Line 1,051: | Line 1,051: | ||
[0.0,0.0,0.0,0.0,1.0,0.0], |
[0.0,0.0,0.0,0.0,1.0,0.0], |
||
[0.0,0.0,0.0,0.0,0.0,1.0]]</pre> |
[0.0,0.0,0.0,0.0,0.0,1.0]]</pre> |
||
=={{header|Julia}}== |
|||
Uses an iterator to generate A003313 (the first 100 values). Uses the Knuth path method for the larger integers. This gives a chain length of 18 for both 27182 and 31415. |
|||
<syntaxhighlight lang="julia">import Base.iterate |
|||
using LinearAlgebra |
|||
mutable struct AdditionChains{T} |
|||
chains::Vector{Vector{T}} |
|||
work_chain::Int |
|||
work_element::Int |
|||
AdditionChains{T}() where T = new{T}([[one(T)]], 1, 1) |
|||
end |
|||
function Base.iterate(acs::AdditionChains, state = 1) |
|||
i, j = acs.work_chain, acs.work_element |
|||
newchain = [acs.chains[i]; acs.chains[i][end] + acs.chains[i][j]] |
|||
push!(acs.chains, newchain) |
|||
if j == length(acs.chains[i]) |
|||
acs.work_chain += 1 |
|||
acs.work_element = 1 |
|||
else |
|||
acs.work_element += 1 |
|||
end |
|||
return newchain, state + 1 |
|||
end |
|||
function findchain!(acs::AdditionChains, n) |
|||
n == 0 && return zero(eltype(first(acs.chains))) |
|||
n == 1 && return [one(eltype(first(acs.chains)))] |
|||
idx = findfirst(a -> a[end] == n, acs.chains) |
|||
if idx == nothing |
|||
for (i, chain) in enumerate(acs) |
|||
chain[end] == n && return chain |
|||
end |
|||
end |
|||
return acs.chains[idx] |
|||
end |
|||
""" memoization for knuth_path """ |
|||
const knuth_path_p, knuth_path_lvl = Dict(1 => 0), [[1]] |
|||
""" knuth path method for addition chains """ |
|||
function knuth_path(n)::Vector{Int} |
|||
iszero(n) && return Int[] |
|||
while !haskey(knuth_path_p, n) |
|||
q = Int[] |
|||
for x in first(knuth_path_lvl), y in knuth_path(x) |
|||
if !haskey(knuth_path_p, x + y) |
|||
knuth_path_p[x + y] = x |
|||
push!(q, x + y) |
|||
end |
|||
end |
|||
knuth_path_lvl[begin] = q |
|||
end |
|||
return push!(knuth_path(knuth_path_p[n]), n) |
|||
end |
|||
function pow(x, chain) |
|||
p, products = 0, Dict{Int, typeof(x)}(0 => one(x), 1 => x) |
|||
for i in chain |
|||
products[i] = products[p] * products[i - p] |
|||
p = i |
|||
end |
|||
return products[chain[end]] |
|||
end |
|||
function test_addition_chains() |
|||
additionchain = AdditionChains{Int}() |
|||
println("First one hundred addition chain lengths:") |
|||
for i in 1:100 |
|||
print(rpad(length(findchain!(additionchain, i)) -1, 3), i % 10 == 0 ? "\n" : "") |
|||
end |
|||
println("\nKnuth chains for addition chains of 31415 and 27182:") |
|||
expchains = Dict(i => knuth_path(i) for i in [31415, 27182]) |
|||
for (n, chn) in expchains |
|||
println("Exponent: ", rpad(n, 10), "\n Addition Chain: $(chn[begin:end-1]))") |
|||
end |
|||
println("\n1.00002206445416^31415 = ", pow(1.00002206445416, expchains[31415])) |
|||
println("1.00002550055251^27182 = ", pow(1.00002550055251, expchains[27182])) |
|||
println("1.00002550055251^(27182 * 31415) = ", pow(BigFloat(pow(1.00002550055251, expchains[27182])), expchains[31415])) |
|||
println("1.000025 + 0.000058i)^27182 = ", pow(Complex(1.000025, 0.000058), expchains[27182])) |
|||
println("1.000022 + 0.000050i)^31415 = ", pow(Complex(1.000022, 0.000050), expchains[31415])) |
|||
end |
|||
test_addition_chains() |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
0 1 2 2 3 3 4 3 4 4 |
|||
5 4 5 5 5 4 5 5 6 5 |
|||
6 6 6 5 6 6 6 6 7 6 |
|||
7 5 6 6 7 6 7 7 7 6 |
|||
7 7 7 7 7 7 8 6 7 7 |
|||
7 7 8 7 8 7 8 8 8 7 |
|||
8 8 8 6 7 7 8 7 8 8 |
|||
9 7 8 8 8 8 8 8 9 7 |
|||
8 8 8 8 8 8 9 8 9 8 |
|||
9 8 9 9 9 7 8 8 8 8 |
|||
Knuth chains for addition chains of 31415 and 27182: |
|||
Exponent: 27182 |
|||
Addition Chain: [1, 2, 3, 5, 7, 14, 21, 35, 70, 140, 143, 283, 566, 849, 1698, 3396, 6792, 6799, 13591]) |
|||
Exponent: 31415 |
|||
Addition Chain: [1, 2, 3, 5, 7, 14, 28, 56, 61, 122, 244, 488, 976, 1952, 3904, 7808, 15616, 15677, 31293]) |
|||
1.00002206445416^31415 = 1.9999999998924638 |
|||
1.00002550055251^27182 = 1.9999999999792688 |
|||
1.00002550055251^(27182 * 31415) = 7.199687435551025768365237017391630520475805934721292725697031724530209692195819e+9456 |
|||
1.000025 + 0.000058i)^27182 = -0.01128636963542673 + 1.9730308496660347im |
|||
1.000022 + 0.000050i)^31415 = 0.00016144681325535107 + 1.9960329014194498im |
|||
</pre> |
|||
=={{header|Nim}}== |
=={{header|Nim}}== |