Knuth's power tree: Difference between revisions
m Describe data structure |
J: bugfix - this kpt is much slower than before but it was already orders of magnitude slower than the builtin |
||
Line 86: | Line 86: | ||
for_n. 1+i.y do. |
for_n. 1+i.y do. |
||
findpath=: [: |. {&P^:a: |
findpath=: [: |. {&P^:a: |
||
lev=.>./L-._ |
|||
for_n.(/: findpath&>)I.L=lev do. |
|||
⚫ | |||
for_a. findpath n do. |
|||
j=. n+a |
|||
l=. 1+n{L |
|||
if. j>y do. break. end. |
|||
if. l>:j{ L do. continue. end. |
|||
L=: l j} L |
|||
⚫ | |||
end. |
|||
end. |
end. |
||
end. |
end. |
||
Line 111: | Line 114: | ||
<lang J> kpt 191 NB. generate sufficiently large tree |
<lang J> kpt 191 NB. generate sufficiently large tree |
||
0 1 1 2 2 3 3 5 4 6 5 |
0 1 1 2 2 3 3 5 4 6 5 10 6 10 7 10 8 16 9 14 10 14 11 13 12 15 13 18 14 28 15 28 16 17 17 21 18 36 19 26 20 40 21 40 22 30 23 42 24 48 25 48 26 52 27 44 28 38 29 31 30 56 31 42 32 64 33 66 34 46 35 57 36 37 37 50 38 76 39 76 40 41 41 43 42 80 43 84 44 47 45 70 46 62 47 57 48 49 49 51 50 100 51 100 52 70 53 104 54 104 55 108 56 112 57 112 58 61 59 112 60 120 61 120 62 75 63 126 64 65 65 129 66 67 67 90 68 136 69 138 70 140 71 140 72 144 73 144 74 132 75 138 76 144 77 79 78 152 79 152 80 160 81 160 82 85 83 162 84 168 85 114 86 168 87 105 88 118 89 176 90 176 91 122 92 184 93 176 94 126 95 190 |
||
findpath 0 |
findpath 0 |
||
0 |
0 |
||
Line 157: | Line 161: | ||
1024 |
1024 |
||
findpath 11 |
findpath 11 |
||
1 2 3 |
1 2 3 5 10 11 |
||
2 usepath 11 |
2 usepath 11 |
||
2048 |
2048 |
||
Line 173: | Line 177: | ||
16384 |
16384 |
||
findpath 15 |
findpath 15 |
||
1 2 3 |
1 2 3 5 10 15 |
||
2 usepath 15 |
2 usepath 15 |
||
32768 |
32768 |
||
Line 185: | Line 189: | ||
131072 |
131072 |
||
findpath 191 |
findpath 191 |
||
1 2 |
1 2 3 5 7 14 19 38 57 95 190 191 |
||
3x usepath 191 |
3x usepath 191 |
||
13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 |
13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 |
||
findpath 81 |
|||
1 2 3 5 10 20 40 41 81 |
|||
(x:1.1) usepath 81 |
(x:1.1) usepath 81 |
||
2253240236044012487937308538033349567966729852481170503814810577345406584190098644811r1000000000000000000000000000000000000000000000000000000000000000000000000000000000 |
2253240236044012487937308538033349567966729852481170503814810577345406584190098644811r1000000000000000000000000000000000000000000000000000000000000000000000000000000000 |
Revision as of 07:33, 3 June 2015
(Used for computing xn efficiently using Knuth's power tree.)
The requirements of this draft task are to compute and show the list of Knuth's power tree integers necessary for the computation of Xn for any real X and any non-negative integer n.
Then, using those integers, calculate and show the exact (not approximate) value of (at least) the integer powers below:
- 2n where n ranges from 0 ──► 17 (inclusive)
- 3191
- 1.181
- 2n where n ranges from 0 ──► 17 (inclusive)
A zero power is often handled separately as a special case.
Optionally, support negative integers (for the power).
An example of a small power tree for some low integers:
1 \ 2 ___________________________________________/ \ / \ 3 4 / \____________________________________ \ / \ \ 5 6 8 / \____________ / \ \ / \ / \ \ 7 10 9 12 16 / //\\ │ │ /\ / _____// \\________ │ │ / \ 14 / / \ \ │ │ / \ /│ \ 11 13 15 20 18 24 17 32 / │ \ │ /\ /\ │ /\ │ /\ │ / │ \ │ / \ / \ │ / \ │ / \ │ 19 21 28 22 23 26 25 30 40 27 36 48 33 34 64 │ /\ /│\ │ │ /\ │ /\ /│\ │ /\ /│\ │ │ /\ │ / \ / │ \ │ │ / \ │ / \ / │ \ │ / \ / │ \ │ │ / \ 38 35 42 29 31 56 44 46 39 52 50 45 60 41 43 80 54 37 72 49 51 96 66 68 65 128
Where, for the power 43, following the tree "downwards" from 1:
- (for 2) compute square of X, store X2
- (for 3) compute X * X2, store X3
- (for 5) compute X3 * X2, store X5
- (for 10) compute square of X5, store X10
- (for 20) compute square of X10, store X20
- (for 40) compute square of X20, store X40
- (for 43) compute X40 * X3 (result).
Note that for every even integer (in the power tree), one just squares the previous value.
For an odd integer, multiple the previous value with an appropriate odd power of X (which was previously calculated).
For the last multiplication in the above example, it would be (43-40), or 3.
According to Dr. Knuth (see below), computer tests have shown that this power tree gives optimum results for all of the n listed above in the graph.
For n ≤ 100,000, the power tree method:
- bests the factor method 88,803 times,
- ties 11,191 times,
- loses 6 times.
See: Donald E. Knuth's book: The Art of Computer Programming, Vol. 2, Second Edition, Seminumerical Algorithms, section 4.6.3: Evaluation of Powers.
See: link codegolf.stackexchange.com/questions/3177/knuths-power-tree It shows a Haskel, Python, and a Ruby computer program example (but they are mostly code golf).
See: link comeoncodeon.wordpress.com/tag/knuth/ (See the section on Knuth's Power Tree.) It shows a C++ computer program example.
See: link to Rosetta Code addition-chain exponentiation.
J
This task is a bit verbose...
We can represent the tree as a list of indices. Each entry in the list gives the value of n for the index a+n. (We can find a using subtraction.)
<lang J>kpt=:3 :0
L=: %P=: (1+y){.0 1 for_n. 1+i.y do. findpath=: [: |. {&P^:a: lev=.>./L-._ for_n.(/: findpath&>)I.L=lev do. for_a. findpath n do. j=. n+a l=. 1+n{L if. j>y do. break. end. if. l>:j{ L do. continue. end. L=: l j} L P=: n j} P end. end. end. P
)
usepath=:4 :0
path=. findpath y exp=. 1,({:path)#x for_ex.(,.~2 -~/\"1])2 ,\path do. 'ea eb ec'=. ex exp=.((ea{exp)*eb{exp) ec} exp end. {:exp
)</lang>
Task examples:
<lang J> kpt 191 NB. generate sufficiently large tree 0 1 1 2 2 3 3 5 4 6 5 10 6 10 7 10 8 16 9 14 10 14 11 13 12 15 13 18 14 28 15 28 16 17 17 21 18 36 19 26 20 40 21 40 22 30 23 42 24 48 25 48 26 52 27 44 28 38 29 31 30 56 31 42 32 64 33 66 34 46 35 57 36 37 37 50 38 76 39 76 40 41 41 43 42 80 43 84 44 47 45 70 46 62 47 57 48 49 49 51 50 100 51 100 52 70 53 104 54 104 55 108 56 112 57 112 58 61 59 112 60 120 61 120 62 75 63 126 64 65 65 129 66 67 67 90 68 136 69 138 70 140 71 140 72 144 73 144 74 132 75 138 76 144 77 79 78 152 79 152 80 160 81 160 82 85 83 162 84 168 85 114 86 168 87 105 88 118 89 176 90 176 91 122 92 184 93 176 94 126 95 190
findpath 0
0
2 usepath 0
1
findpath 1
1
2 usepath 1
2
findpath 2
1 2
2 usepath 2
4
findpath 3
1 2 3
2 usepath 3
8
findpath 4
1 2 4
2 usepath 4
16
findpath 5
1 2 3 5
2 usepath 5
32
findpath 6
1 2 3 6
2 usepath 6
64
findpath 7
1 2 3 5 7
2 usepath 7
128
findpath 8
1 2 4 8
2 usepath 8
256
findpath 9
1 2 3 6 9
2 usepath 9
512
findpath 10
1 2 3 5 10
2 usepath 10
1024
findpath 11
1 2 3 5 10 11
2 usepath 11
2048
findpath 12
1 2 3 6 12
2 usepath 12
4096
findpath 13
1 2 3 5 10 13
2 usepath 13
8192
findpath 14
1 2 3 5 7 14
2 usepath 14
16384
findpath 15
1 2 3 5 10 15
2 usepath 15
32768
findpath 16
1 2 4 8 16
2 usepath 16
65536
findpath 17
1 2 4 8 16 17
2 usepath 17
131072
findpath 191
1 2 3 5 7 14 19 38 57 95 190 191
3x usepath 191
13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
findpath 81
1 2 3 5 10 20 40 41 81
(x:1.1) usepath 81
2253240236044012487937308538033349567966729852481170503814810577345406584190098644811r1000000000000000000000000000000000000000000000000000000000000000000000000000000000 </lang>
Note that an 'r' in a number indicates a rational number with the numerator to the left of the r and the denominator to the right of the r. We could instead use decimal notation by indicating how many characters of result we want to see, as well as how many characters to the right of the decimal point we want to see.
Thus, for example:
<lang J> 90j83 ": (x:1.1) usepath 81
2253.24023604401248793730853803334956796672985248117050381481057734540658419009864481100</lang>
Python
<lang python>from __future__ import print_function
- remember the tree generation state and expand on demand
def path(n, p = {1:0}, lvl=1): if not n: return [] while n not in p: q = [] for x,y in ((x, x+y) for x in lvl[0] for y in path(x) if not x+y in p): p[y] = x q.append(y) lvl[0] = q
return path(p[n]) + [n]
def tree_pow(x, n):
r, p = {0:1, 1:x}, 0 for i in path(n): r[i] = r[i-p] * r[p] p = i return r[n]
def show_pow(x, n):
fmt = "%d: %s\n" + ["%g^%d = %f", "%d^%d = %d"][x==int(x)] + "\n" print(fmt % (n, repr(path(n)), x, n, tree_pow(x, n)))
for x in range(18): show_pow(2, x) show_pow(3, 191) show_pow(1.1, 81)</lang>
- Output:
0: [] 2^0 = 1 1: [1] 2^1 = 2 2: [1, 2] 2^2 = 4 <... snipped ...> 17: [1, 2, 4, 8, 16, 17] 2^17 = 131072 191: [1, 2, 3, 5, 7, 14, 19, 38, 57, 95, 190, 191] 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 81: [1, 2, 3, 5, 10, 20, 40, 41, 81] 1.1^81 = 2253.240236
REXX
This REXX version supports results up to 1,000 decimal digits (which can be expanded with the numeric digits nnn REXX statement. Also, negative powers are supported. <lang rexx>/*REXX program produces & shows a power tree for P, and calculates & shows X^P*/ numeric digits 1000 /*be able to handle some large numbers.*/ parse arg nums /*get set stuff: sets of three numbers.*/ if nums= then nums='2 -4 17 3 191 191 1.1 81' /*Not specified? Use defaults*/
/*X lowPow highPow ··· repeat*/ do until nums= parse var nums x pL pH nums; x=x/1 /*get X, lowP, highP; and normalize X. */ if pH= then pH=pL /*No highPower? Then assume lowPower. */
do e=pL to pH; p=abs(e)/1 /*use a range of powers; use │E│ */ $=powerTree(p); w=length(pH) /*construct the power tree, (pow list).*/ /* [↑] W≡length for show*/ do i=1 for words($); @.i=word($,i); end /*build a fast pow array.*/
if p==0 then do; z=1; call show; iterate; end /*handle case of 0 power.*/ !.=.; z=x; !.1=z; prv=z /*define the first power of X. */
do k=2 to words($); n=@.k /*obtain the power (number) to be used.*/ prev=k-1; diff=n-@.prev /*these are used for the odd powers. */ if n//2==0 then z=prv**2 /*Even power? Then square the number.*/ else z=z*!.diff /* Odd " " mult. by pow diff.*/ !.n=z /*remember for other multiplications. */ prv=z /*remember for squaring the numbers. */ end /*k*/ call show /*display the expression and its value.*/ end /*e*/ end /*until nums ···*/
exit /*stick a fork in it, we're all done. */ /*────────────────────────────────POWERTREE subroutine────────────────────────*/ powerTree: arg y 1 oy; $= /*Z is the result; $ is the power tree.*/ if y=0 | y=1 then return y /*handle special cases for zero & unity*/
- .=0; @.=0; #.0=1 /*define default & initial array values*/
/* [↓] add blank "flag" thingy──►list.*/ do while \(y//2); $=$ ' ' /*reduce "front" even power #s to odd #*/ if y\==oy then $=y $ /*(only) ignore the first power number*/ y=y%2 /*integer divide the power (it's even).*/ end /*while ···*/
if $\== then $=y $ /*re─introduce the last power number. */ $=$ oy /*insert last power number 1st in list.*/ if y>1 then do while @.y==0; n=#.0; m=0
do while n\==0; q=0; s=n do while s\==0; _=n+s if @._==0 then do; if q==0 then m_=_; #._=q; @._=n; q=_ end s=@.s end /*while s¬==0*/ if q\==0 then do; #.m=q; m=m_; end n=#.n end /*while n¬==0*/ #.m=0 end /*while @.y==0*/
z=@.y
do while z\==0; $=z $; z=@.z; end /*build power list.*/
return space($) /*────────────────────────────────SHOW subroutine─────────────────────────────*/ show: if e<0 then z=format(1/z,,40)/1; _=right(e,w) /*use reciprocal ? */ say left('power tree for ' _ " is: " $,60) '═══' x"^"_ ' is: ' z; return</lang> output when using the default inputs:
power tree for -4 is: 1 2 4 ═══ 2^-4 is: 0.0625 power tree for -3 is: 1 2 3 ═══ 2^-3 is: 0.125 power tree for -2 is: 1 2 ═══ 2^-2 is: 0.25 power tree for -1 is: 1 ═══ 2^-1 is: 0.5 power tree for 0 is: 0 ═══ 2^ 0 is: 1 power tree for 1 is: 1 ═══ 2^ 1 is: 2 power tree for 2 is: 1 2 ═══ 2^ 2 is: 4 power tree for 3 is: 1 2 3 ═══ 2^ 3 is: 8 power tree for 4 is: 1 2 4 ═══ 2^ 4 is: 16 power tree for 5 is: 1 2 3 5 ═══ 2^ 5 is: 32 power tree for 6 is: 1 2 3 6 ═══ 2^ 6 is: 64 power tree for 7 is: 1 2 3 5 7 ═══ 2^ 7 is: 128 power tree for 8 is: 1 2 4 8 ═══ 2^ 8 is: 256 power tree for 9 is: 1 2 3 6 9 ═══ 2^ 9 is: 512 power tree for 10 is: 1 2 3 5 10 ═══ 2^10 is: 1024 power tree for 11 is: 1 2 3 5 10 11 ═══ 2^11 is: 2048 power tree for 12 is: 1 2 3 6 12 ═══ 2^12 is: 4096 power tree for 13 is: 1 2 3 5 10 13 ═══ 2^13 is: 8192 power tree for 14 is: 1 2 3 5 7 14 ═══ 2^14 is: 16384 power tree for 15 is: 1 2 3 5 10 15 ═══ 2^15 is: 32768 power tree for 16 is: 1 2 4 8 16 ═══ 2^16 is: 65536 power tree for 17 is: 1 2 4 8 16 17 ═══ 2^17 is: 131072 power tree for 191 is: 1 2 3 5 7 14 19 38 57 95 190 191 ═══ 3^191 is: 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 power tree for 81 is: 1 2 3 5 10 20 40 41 81 ═══ 1.1^81 is: 2253.240236044012487937308538033349567966729852481170503814810577345406584190098644811