Addition-chain exponentiation: Difference between revisions
m
Python example
m (typo) |
m (Python example) |
||
Line 1,751:
--1.00003 ^ 27182 (20) = 1.9999999999727968238298855
--1.00003 ^ 27182 (19) = 1.9999999999727968238298527</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">''' Rosetta Code task Addition-chain_exponentiation '''
from math import sqrt
from numpy import array
from mpmath import mpf
class AdditionChains:
''' two methods of calculating addition chains '''
def __init__(self):
''' memoization for knuth_path '''
self.chains, self.idx, self.pos = [[1]], 0, 0
self.pat, self.lvl = {1: 0}, [[1]]
def add_chain(self):
''' method 1: add chains depth then breadth first until done '''
newchain = self.chains[self.idx].copy()
newchain.append(self.chains[self.idx][-1] +
self.chains[self.idx][self.pos])
self.chains.append(newchain)
if self.pos == len(self.chains[self.idx])-1:
self.idx += 1
self.pos = 0
else:
self.pos += 1
return newchain
def find_chain(self, nexp):
''' method 1 interface: search for chain ending with n, adding more as needed '''
assert nexp > 0
if nexp == 1:
return [1]
chn = next((a for a in self.chains if a[-1] == nexp), None)
if chn is None:
while True:
chn = self.add_chain()
if chn[-1] == nexp:
break
return chn
def knuth_path(self, ngoal):
''' method 2: knuth method, uses memoization to search for a shorter chain '''
if ngoal < 1:
return []
while not ngoal in self.pat:
new_lvl = []
for i in self.lvl[0]:
for j in self.knuth_path(i):
if not i + j in self.pat:
self.pat[i + j] = i
new_lvl.append(i + j)
self.lvl[0] = new_lvl
returnpath = self.knuth_path(self.pat[ngoal])
returnpath.append(ngoal)
return returnpath
def cpow(xbase, chain):
''' raise xbase by an addition exponentiation chain for what becomes x**chain[-1] '''
pows, products = 0, {0: 1, 1: xbase}
for i in chain:
products[i] = products[pows] * products[i - pows]
pows = i
return products[chain[-1]]
if __name__ == '__main__':
# test both addition chain methods
acs = AdditionChains()
print('First one hundred addition chain lengths:')
for k in range(1, 101):
print(f'{len(acs.find_chain(k))-1:3}', end='\n'if k % 10 == 0 else '')
print('\nKnuth chains for addition chains of 31415 and 27182:')
chns = {m: acs.knuth_path(m) for m in [31415, 27182]}
for (num, cha) in chns.items():
print(f'Exponent: {num:10}\n Addition Chain: {cha[:-1]}')
print('\n1.00002206445416^31415 =', cpow(1.00002206445416, chns[31415]))
print('1.00002550055251^27182 =', cpow(1.00002550055251, chns[27182]))
print('1.000025 + 0.000058i)^27182 =',
cpow(complex(1.000025, 0.000058), chns[27182]))
print('1.000022 + 0.000050i)^31415 =',
cpow(complex(1.000022, 0.000050), chns[31415]))
sq05 = mpf(sqrt(0.5))
mat = array([[sq05, 0, sq05, 0, 0, 0], [0, sq05, 0, sq05, 0, 0], [0, sq05, 0, -sq05, 0, 0],
[-sq05, 0, sq05, 0, 0, 0], [0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 1, 0]])
print('matrix A ^ 27182 =')
print(cpow(mat, chns[27182]))
print('matrix A ^ 31415 =')
print(cpow(mat, chns[31415]))
print('(matrix A ** 27182) ** 31415 =')
print(cpow(cpow(mat, chns[27182]), chns[31415]))
</syntaxhighlight>{{out}}
<pre>
First one hundred addition chain lengths:
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: 31415
Addition Chain: [1, 2, 3, 5, 7, 14, 28, 56, 61, 122, 244, 488, 976, 1952, 3904, 7808, 15616, 15677, 31293]
Exponent: 27182
Addition Chain: [1, 2, 3, 5, 7, 14, 21, 35, 70, 140, 143, 283, 566, 849, 1698, 3396, 6792, 6799, 13591]
1.00002206445416^31415 = 1.9999999998924638
1.00002550055251^27182 = 1.9999999999792688
1.000025 + 0.000058i)^27182 = (-0.01128636963542673+1.9730308496660347j)
1.000022 + 0.000050i)^31415 = (0.00016144681325535107+1.9960329014194498j)
matrix A ^ 27182 =
[[mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0 0]
[0 mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0]
[0 mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0]
[mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0 0]
[0 0 0 0 0 1]
[0 0 0 0 1 0]]
matrix A ^ 31415 =
[[mpf('3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0 0]
[0 mpf('3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0]
[0 mpf('3.7268602513250562e-4729') 0 mpf('-3.7268602513250562e-4729') 0
0]
[mpf('-3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0
0]
[0 0 0 0 0 1]
[0 0 0 0 1 0]]
(matrix A ** 27182) ** 31415 =
[[mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148')
0 0 0]
[0 mpf('1.77158544475749e-128528148') 0
mpf('1.77158544475749e-128528148') 0 0]
[0 mpf('1.77158544475749e-128528148') 0
mpf('1.77158544475749e-128528148') 0 0]
[mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148')
0 0 0]
[0 0 0 0 0 1]
[0 0 0 0 1 0]]
</pre>
=={{header|Racket}}==
| |||