Knuth's power tree
You are encouraged to solve this task according to the task description, using any language you may know.
(Knuth's power tree is used for computing xn efficiently.)
- Task
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 values 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 integer powers.
- Example
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, multiply 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.
- References
-
- Donald E. Knuth's book: The Art of Computer Programming, Vol. 2, Second Edition, Seminumerical Algorithms, section 4.6.3: Evaluation of Powers.
- link codegolf.stackexchange.com/questions/3177/knuths-power-tree It shows a Haskell, Python, and a Ruby computer program example (but they are mostly code golf).
- link comeoncodeon.wordpress.com/tag/knuth/ (See the section on Knuth's Power Tree.) It shows a C++ computer program example.
- link to Rosetta Code addition-chain exponentiation.
V p = [1 = 0]
V lvl = [[1]]
F path(n)
I !n
R [Int]()
L n !C :p
[Int] q
L(x) :lvl[0]
L(y) path(x)
I !(x + y C :p)
:p[x + y] = x
q.append(x + y)
:lvl[0] = q
R path(:p[n]) [+] [n]
F tree_pow(x, n)
T Ty = T(x)
V (r, p) = ([0 = Ty(1), 1 = x], 0)
L(i) path(n)
r[i] = r[i - p] * r[p]
p = i
R r[n]
F show_pow_i(x, n)
print("#.: #.\n#.^#. = #.\n".format(n, path(n), x, n, tree_pow(BigInt(x), n)))
F show_pow_f(x, n)
print("#.: #.\n#.^#. = #.6\n".format(n, path(n), x, n, tree_pow(x, n)))
L(x) 0.<18
show_pow_i(2, x)
show_pow_i(3, 191)
show_pow_f(1.1, 81)- Output:
0: [] 2^0 = 1 1: [1] 2^1 = 2 2: [1, 2] 2^2 = 4 3: [1, 2, 3] 2^3 = 8 4: [1, 2, 4] 2^4 = 16 5: [1, 2, 3, 5] 2^5 = 32 6: [1, 2, 3, 6] 2^6 = 64 7: [1, 2, 3, 5, 7] 2^7 = 128 8: [1, 2, 4, 8] 2^8 = 256 9: [1, 2, 3, 6, 9] 2^9 = 512 10: [1, 2, 3, 5, 10] 2^10 = 1024 11: [1, 2, 3, 5, 10, 11] 2^11 = 2048 12: [1, 2, 3, 6, 12] 2^12 = 4096 13: [1, 2, 3, 5, 10, 13] 2^13 = 8192 14: [1, 2, 3, 5, 7, 14] 2^14 = 16384 15: [1, 2, 3, 5, 10, 15] 2^15 = 32768 16: [1, 2, 4, 8, 16] 2^16 = 65536 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
In this program the power tree is in a generic package which is instantiated for each base value.
with Ada.Containers.Ordered_Maps;
with Ada.Numerics.Big_Numbers.Big_Integers;
with Ada.Text_IO; use Ada.Text_IO;
procedure Power_Tree is
Debug_On : constant Boolean := False;
generic
type Result_Type is private;
Base : Result_Type;
Identity : Result_Type;
with function "*" (Left, Right : Result_Type) return Result_Type is <>;
package Knuth_Power_Tree is
subtype Exponent is Natural;
function Power (Exp : Exponent) return Result_Type;
end Knuth_Power_Tree;
package body Knuth_Power_Tree is
package Power_Trees is
new Ada.Containers.Ordered_Maps (Key_Type => Exponent,
Element_Type => Result_Type);
Tree : Power_Trees.Map;
procedure Debug (Item : String) is
begin
if Debug_On then
Put_Line (Standard_Error, Item);
end if;
end Debug;
function Power (Exp : Exponent) return Result_Type is
Pow : Result_Type;
begin
if Tree.Contains (Exp) then
return Tree.Element (Exp);
else
Debug ("lookup failed of " & Exp'Image);
end if;
if Exp mod 2 = 0 then
Debug ("lookup half " & Exponent'(Exp / 2)'Image);
Pow := Power (Exp / 2);
Pow := Pow * Pow;
else
Debug ("lookup one less " & Exponent'(Exp - 1)'Image);
Pow := Power (Exp - 1);
Pow := Result_Type (Base) * Pow;
end if;
Debug ("insert " & Exp'Image);
Tree.Insert (Key => Exp, New_Item => Pow);
return Pow;
end Power;
begin
Tree.Insert (Key => 0, New_Item => Identity);
end Knuth_Power_Tree;
procedure Part_1
is
package Power_2 is new Knuth_Power_Tree (Result_Type => Long_Integer,
Base => 2,
Identity => 1);
R : Long_Integer;
begin
Put_Line ("=== Part 1 ===");
for N in 0 .. 25 loop
R := Power_2.Power (N);
Put ("2 **"); Put (N'Image);
Put (" ="); Put (R'Image);
New_Line;
end loop;
end Part_1;
procedure Part_2
is
use Ada.Numerics.Big_Numbers.Big_Integers;
package Power_3 is new Knuth_Power_Tree (Result_Type => Big_Integer,
Base => 3,
Identity => 1);
R : Big_Integer;
begin
Put_Line ("=== Part 2 ===");
for E in 190 .. 192 loop
R := Power_3.Power (E);
Put ("3 **" & E'Image & " ="); Put (R'Image); New_Line;
end loop;
end Part_2;
procedure Part_3
is
subtype Real is Long_Long_Float;
package Real_IO is new Ada.Text_IO.Float_IO (Real);
package Power_1_1 is new Knuth_Power_Tree (Result_Type => Real,
Base => 1.1,
Identity => 1.0);
R : Real;
begin
Put_Line ("=== Part 3 ===");
for E in 81 .. 84 loop
R := Power_1_1.Power (E);
Put ("1.1 **" & E'Image & " = ");
Real_IO.Put (R, Exp => 0, Aft => 6);
New_Line;
end loop;
end Part_3;
begin
Part_1;
Part_2;
Part_3;
end Power_Tree;
- Output:
=== Part 1 === 2 ** 0 = 1 2 ** 1 = 2 2 ** 2 = 4 2 ** 3 = 8 2 ** 4 = 16 2 ** 5 = 32 2 ** 6 = 64 2 ** 7 = 128 2 ** 8 = 256 2 ** 9 = 512 2 ** 10 = 1024 2 ** 11 = 2048 2 ** 12 = 4096 2 ** 13 = 8192 2 ** 14 = 16384 2 ** 15 = 32768 2 ** 16 = 65536 2 ** 17 = 131072 2 ** 18 = 262144 2 ** 19 = 524288 2 ** 20 = 1048576 2 ** 21 = 2097152 2 ** 22 = 4194304 2 ** 23 = 8388608 2 ** 24 = 16777216 2 ** 25 = 33554432 === Part 2 === 3 ** 190 = 4498196224760364601242719132174628305800834098010033971355568455673974002968757862019419449 3 ** 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 3 ** 192 = 40483766022843281411184472189571654752207506882090305742200116101065766026718820758174775041 === Part 3 === 1.1 ** 81 = 2253.240236 1.1 ** 82 = 2478.564260 1.1 ** 83 = 2726.420686 1.1 ** 84 = 2999.062754
Power tree
We build the tree using tree.lib, adding leaves until the target n is found.
(lib 'tree)
;; displays a chain hit
(define (power-hit target chain)
(vector-push chain target)
(printf "L(%d) = %d - chain:%a "
target (1- (vector-length chain)) chain)
(vector-pop chain))
;; build the power-tree : add 1 level of leaf nodes
;; display all chains which lead to target
(define (add-level node chain target nums (new))
(vector-push chain (node-datum node))
(cond
[(node-leaf? node)
;; add leaves by summing this node to all nodes in chain
;; do not add leaf if number already known
(for [(prev chain)]
(set! new (+ prev (node-datum node)))
(when (= new target) (power-hit target chain ))
#:continue (vector-search* new nums)
(node-add-leaf node new)
(vector-insert* nums new)
)]
[else ;; not leaf node -> recurse
(for [(son (node-sons node))]
(add-level son chain target nums )) ])
(vector-pop chain))
;; add levels in tree until target found
;; return (number of nodes . upper-bound for L(target))
(define (power-tree target)
(define nums (make-vector 1 1)) ;; known nums = 1
(define T (make-tree 1)) ;; root node has value 1
(printf "Looking for %d in %a." target T)
(while #t
#:break (vector-search* target nums) => (tree-count T)
(add-level T init-chain: (make-vector 0) target nums)
))
- Output:
(for ((n (in-range 2 18))) (power-tree n)) L(2) = 1 - chain:#( 1 2) L(3) = 2 - chain:#( 1 2 3) [ ... ] (power-tree 17) Looking for 17 in (🌴 1). L(17) = 5 - chain:#( 1 2 4 8 16 17) (power-tree 81) Looking for 81 in (🌴 1). L(81) = 8 - chain:#( 1 2 3 5 10 20 40 41 81) L(81) = 8 - chain:#( 1 2 3 5 10 20 40 80 81) L(81) = 8 - chain:#( 1 2 3 6 9 18 27 54 81) L(81) = 8 - chain:#( 1 2 3 6 9 18 36 72 81) L(81) = 8 - chain:#( 1 2 4 8 16 32 64 65 81) (power-tree 191) Looking for 191 in (🌴 1). L(191) = 11 - chain:#( 1 2 3 5 7 14 19 38 57 95 190 191) L(191) = 11 - chain:#( 1 2 3 5 7 14 21 42 47 94 188 191) L(191) = 11 - chain:#( 1 2 3 5 7 14 21 42 63 126 189 191) L(191) = 11 - chain:#( 1 2 3 5 7 14 28 31 59 118 177 191) L(191) = 11 - chain:#( 1 2 3 5 7 14 28 31 62 93 186 191) L(191) = 11 - chain:#( 1 2 3 5 10 11 22 44 88 176 181 191) (power-tree 12509) ;; not optimal Looking for 12509 in (🌴 1). L(12509) = 18 - chain:#( 1 2 3 5 10 13 26 39 78 156 312 624 1248 2496 2509 3757 6253 12506 12509) L(12509) = 18 - chain:#( 1 2 3 5 10 15 25 50 75 125 250 500 1000 2000 2003 4003 8006 12009 12509) (power-tree 222222) Looking for 222222 in (🌴 1). L(222222) = 22 - chain:#( 1 2 3 5 7 14 21 35 70 105 210 420 840 1680 1687 3367 6734 13468 26936 53872 57239 111111 222222)
Exponentiation
;; j such as chain[i] = chain[i-1] + chain[j]
(define (adder chain i)
(for ((j i)) #:break (= [chain i] (+ [chain(1- i)] [chain j])) => j ))
(define (power-exp x chain)
(define lg (vector-length chain))
(define pow (make-vector lg x))
(for ((i (in-range 1 lg)))
(vector-set! pow i ( * [pow [1- i]] [pow (adder chain i)])))
[pow (1- lg)])
- Output:
(power-exp 2 #( 1 2 4 8 16 17) )
→ 131072
(power-exp 1.1 #( 1 2 3 5 10 20 40 41 81) )
→ 2253.2402360440283
(lib 'bigint)
bigint.lib v1.4 ® EchoLisp
Lib: bigint.lib loaded.
(power-exp 3 #( 1 2 3 5 7 14 19 38 57 95 190 191) )
→ 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
Integer Exponentiation
// Integer exponentiation using Knuth power tree. Nigel Galloway: October 29th., 2020
let kT α=let n=Array.zeroCreate<int*int list>((pown 2 (α+1))+1)
let fN g=let rec fN p=[yield g+p; if p>0 then let g,_=n.[p] in yield! fN g] in (g+g)::(fN (fst n.[g]))|>List.rev
let fG g=[for α,β in g do for g in β do let p,_=n.[g] in n.[g]<-(p,fN g|>List.filter(fun β->if box n.[β]=null then n.[β]<-(g,[]); true else false)); yield n.[g]]
let rec kT n g=match g with 0->() |_->let n=fG n in kT n (g-1)
let fE X g=let α=let rec fE g=[|yield g; if g>1 then yield! fE (fst n.[g])|] in fE g
let β=Array.zeroCreate<bigint>α.Length
let l=β.Length-1
β.[l]<-bigint (X+0)
for e in l-1.. -1..0 do β.[e]<-match α.[e]%2 with 0->β.[e+1]*β.[e+1] |_->let l=α.[e+1] in β.[e+1]*β.[α|>Array.findIndex(fun n->l+n=α.[e])]
β.[0]
n.[1]<-(0,[2]); n.[2]<-(1,[]); kT [n.[1]] α; (fun n g->if g=0 then 1I else fE n g)
let xp=kT 11
[0..17]|>List.iter(fun n->printfn "2**%d=%A\n" n (xp 2 n))
printfn "3**191=%A" (xp 3 191)
{out}
2**0=1 2**1=2 2**2=4 2**3=8 2**4=16 2**5=32 2**6=64 2**7=128 2**8=256 2**9=512 2**10=1024 2**11=2048 2**12=4096 2**13=8192 2**14=16384 2**15=32768 2**16=65536 2**17=131072 3**191=13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
Float Exponentiation
// Float exponentiation using Knuth power tree. Nigel Galloway: October 29th., 2020
let kTf α=let n=Array.zeroCreate<int*int list>((pown 2 (α+1))+1)
let fN g=let rec fN p=[yield g+p; if p>0 then let g,_=n.[p] in yield! fN g] in (g+g)::(fN (fst n.[g]))|>List.rev
let fG g=[for α,β in g do for g in β do let p,_=n.[g] in n.[g]<-(p,fN g|>List.filter(fun β->if box n.[β]=null then n.[β]<-(g,[]); true else false)); yield n.[g]]
let rec kT n g=match g with 0->() |_->let n=fG n in kT n (g-1)
let fE X g=let α=let rec fE g=[|yield g; if g>1 then yield! fE (fst n.[g])|] in fE g
let β=Array.zeroCreate<float>α.Length
let l=β.Length-1
β.[l]<-X
for e in l-1.. -1..0 do β.[e]<-match α.[e]%2 with 0->β.[e+1]*β.[e+1] |_->let l=α.[e+1] in β.[e+1]*β.[α|>Array.findIndex(fun n->l+n=α.[e])]
β.[0]
n.[1]<-(0,[2]); n.[2]<-(1,[]); kT [n.[1]] α; (fun n g->if g=0 then 1.0 else fE n g)
let xpf=kTf 11
printfn "1.1**81=%f" (xpf 1.1 81)
- Output:
1.1**81=2253.240236
module hash_map_module
use iso_fortran_env
use qdmodule
implicit none
private
public :: HashEntry, HashMap
type HashEntry
integer(int64) :: key, value
type(HashEntry), pointer :: sgte => null()
end type HashEntry
! Container for pointer to HashEntry
type HashEntryContainer
type(HashEntry), pointer :: ptr => null()
end type HashEntryContainer
type HashMap
type(HashEntryContainer), allocatable :: buckets(:)
contains
procedure :: put => hash_map_put
procedure :: get => hash_map_get
procedure :: contains => hash_map_contains
end type HashMap
contains
subroutine hash_map_put(this, key, value)
class(HashMap), intent(inout) :: this
integer(int64), intent(in) :: key, value
integer(int64) :: bucket
type(HashEntry), pointer :: entry
bucket = iand(key, int(255, int64))
entry => this % buckets(bucket) % ptr
! Check for existing key
do while (associated(entry))
if (entry % key == key) then
entry % value = value
return
end if
entry => entry % sgte
end do
! Insert new entry at head
allocate (entry)
entry % key = key
entry % value = value
entry % sgte => this % buckets(bucket) % ptr ! Correct pointer assignment
this % buckets(bucket) % ptr => entry ! Update head pointer
end subroutine hash_map_put
function hash_map_get(this, key, defaultValue) result(res)
class(HashMap), intent(in) :: this
integer(int64), intent(in) :: key
integer(int64), optional, intent(in) :: defaultValue
integer(int64) :: res, bucket
type(HashEntry), pointer :: entry
res = 0
if (present(defaultValue)) res = defaultValue
bucket = iand(key, int(255, int64))
entry => this % buckets(bucket) % ptr
do while (associated(entry))
if (entry % key == key) then
res = entry % value
return
end if
entry => entry % sgte
end do
end function hash_map_get
function hash_map_contains(this, key) result(res)
class(HashMap), intent(in) :: this
integer(int64), intent(in) :: key
logical :: res
integer(int64) :: bucket
type(HashEntry), pointer :: entry
bucket = iand(key, int(255, int64))
entry => this % buckets(bucket) % ptr
res = .false.
do while (associated(entry))
if (entry % key == key) then
res = .true.
return
end if
entry => entry % sgte
end do
end function hash_map_contains
end module hash_map_module
module int_array_module
use iso_fortran_env
implicit none
private
public :: IntArray
type IntArray
integer(int64), allocatable :: dato(:)
integer(int64) :: size = 0
contains
procedure :: add => int_array_add
procedure :: get => int_array_get
procedure :: getSize => int_array_size
procedure :: clear => int_array_clear
end type
contains
subroutine int_array_add(this, value)
use iso_fortran_env
class(IntArray), intent(inout) :: this
integer(int64), intent(in) :: value
integer(int64), allocatable :: temp(:)
integer(int64) :: current_capacity
if (.not. allocated(this % dato)) then
allocate (this % dato(30))
current_capacity = 30_8
else
current_capacity = size(this % dato)
end if
if (this % size >= current_capacity) then
call move_alloc(this % dato, temp)
allocate (this % dato(2 * current_capacity))
this % dato(1:this % size) = temp(1:this % size)
end if
this % size = this % size + 1
this % dato(this % size) = value
end subroutine int_array_add
function int_array_get(this, index) result(res)
class(IntArray), intent(in) :: this
integer(int64), intent(in) :: index
integer(int64) :: res
res = 0
if (index >= 1 .and. index <= this % size) res = this % dato(index)
end function
function int_array_size(this) result(res)
use iso_fortran_env
class(IntArray), intent(in) :: this
integer(int64) :: res
res = this % size
end function
subroutine int_array_clear(this)
class(IntArray), intent(inout) :: this
this % size = 0
end subroutine
end module
module power_module
use iso_fortran_env
use hash_map_module
use int_array_module
implicit none
private
public :: showPow, initialize_power ! Add initialization routine
type(HashMap) :: p
type(IntArray) :: levels(1)
contains
subroutine initialize_power
integer :: i
if (allocated(p % buckets)) deallocate (p % buckets)
allocate (p % buckets(0:255))
do i = 0, 255
nullify (p % buckets(i) % ptr) ! Explicitly nullify all bucket pointers
end do
call levels(1) % clear()
call p % put(1_8, 0_8)
call levels(1) % add(1_8)
end subroutine
recursive function path(n) result(res)
use iso_fortran_env
use hash_map_module
use int_array_module
implicit none
integer(int64), intent(in) :: n
type(IntArray) :: res
integer(int64) :: i, j, x, y, sum, curr
type(IntArray) :: q, tempPath
integer(int64) :: max_iter, iter
! Safety limit to prevent infinite loops
max_iter = 100000000
iter = 0
if (n == 0) return
do while (.not. p % contains(n))
q = IntArray()
do i = 1, levels(1) % getSize()
x = levels(1) % get(i)
associate (xPath => path(x))
do j = 1, xPath % getSize()
y = xPath % get(j)
sum = x + y
if (p % contains(sum)) cycle
call p % put(sum, x)
call q % add(sum)
end do
end associate
end do
call levels(1) % clear()
do i = 1, q % getSize()
call levels(1) % add(q % get(i))
end do
if (q % getSize() == 0) exit
iter = iter + 1
if (iter > max_iter) then
print *, "Path search exceeded maximum iterations for n=", n
exit
end if
end do
curr = n
tempPath = IntArray()
do while (curr /= 0)
call tempPath % add(curr)
curr = p % get(curr, 0_8)
end do
do i = tempPath % getSize(), 1, -1
call res % add(tempPath % get(i))
end do
return
end function path
function treePowBig(x, n) result(res)
use iso_fortran_env
integer(int64), intent(in) :: x, n
integer(kind=16) :: res
integer(kind=16), allocatable :: r(:)
integer(int64) :: i, curr, p
type(IntArray) :: pathToN
! Allocate array to hold intermediate results
allocate (r(0:n))
r(0) = 1_16
if (n >= 1) r(1) = x
p = 0
pathToN = path(n)
do i = 1, pathToN % getSize()
curr = pathToN % get(i)
r(curr) = r(curr - p) * r(p)
p = curr
end do
res = r(n)
deallocate (r)
end function treePowBig
function treePowDecimal(x, n) result(res)
use qdmodule
type(qd_real), intent(in) :: x
integer(int64), intent(in) :: n
type(qd_real) :: res
type(qd_real), allocatable :: r(:)
integer(int64) :: i, curr, p
type(IntArray) :: pathToN
allocate (r(0:n))
r(0) = "1.0" ! r(0) = 1.0 in quad-double
if (n >= 1) r(1) = x
p = 0
pathToN = path(n)
do i = 1, pathToN % getSize()
curr = pathToN % get(i)
r(curr) = r(curr - p) * r(p)
p = curr
end do
res = r(n)
deallocate (r)
end function treePowDecimal
subroutine showPow(x, n, xx)
use iso_fortran_env
use qdmodule
use int_array_module
implicit none
character(len=*) :: xx
type(qd_real), intent(in) :: x
integer(int64), intent(in) :: n
type(IntArray) :: pathToN
character(len=1024) :: pathStr
integer(int64) :: i
type(qd_real) :: result
character(len=1024) :: str
character(len=5) :: str2
double precision :: try
integer :: unit
!
pathToN = path(n)
pathStr = ""
do i = 1, pathToN % getSize()
if (i > 1) pathStr = trim(pathStr)//", "
write (pathStr(len_trim(pathStr) + 1:), '(I0)') pathToN % get(i)
end do
print '(I0,": [",A,"]" )', n, trim(pathStr)
result = treePowDecimal(x, n)
! call qdwrite(6, result) ! Writes result to stdout (unit 6)
! Create a temporary internal file (string)
open (newunit=unit, status='scratch', action='readwrite', form='formatted')
! Write qd_real to the internal file
call qdwrite(unit, result)
! Read back into the string
rewind (unit)
read (unit, '(A)') str
close (unit)
write (str2, '(i0)') n
write (*, '(a)', advance="no") 'Result: '//xx//'^'//trim(adjustl(str2))//' = '
read (str, *) try
if ((mod(try, 1.0) == 0.0) .and. (try < real(huge(1_16), kind=real64))) then
write (str, '(I0)') int(try, kind=16)
write (*, '(A)') trim(adjustl(str))
else
if (try < 1.0d10) then
write (str, '(f20.13)') try
write (*, '(a)') trim(adjustl(str))
else
write (*, '(a)') trim(adjustl(str))
end if
end if
print *
end subroutine showPow
end module
program main
use iso_fortran_env
use qdmodule
use power_module
implicit none
integer(int64) :: i
type(qd_real) :: x
call initialize_power()
x = "2.0"
do i = 0, 17
call showPow(x, i, "2.0")
end do
x = "1.1"
call showPow(x, 81_8, "1.1")
call initialize_power()
x = "3.0"
call showPow(x, 191_8, "3.0")
stop 'Normal Finish'
end program
- Output:
0: [] Result: 2.0^0 = 1 1: [1] Result: 2.0^1 = 2 2: [1,2] Result: 2.0^2 = 4 3: [1,2,3] Result: 2.0^3 = 8 4: [1,2,4] Result: 2.0^4 = 16 5: [1,2,3,5] Result: 2.0^5 = 32 6: [1,2,3,6] Result: 2.0^6 = 64 7: [1,2,3,5,7] Result: 2.0^7 = 128 8: [1,2,4,8] Result: 2.0^8 = 256 9: [1,2,3,6,9] Result: 2.0^9 = 512 10: [1,2,3,5,10] Result: 2.0^10 = 1024 11: [1,2,3,5,10,11] Result: 2.0^11 = 2048 12: [1,2,3,6,12] Result: 2.0^12 = 4096 13: [1,2,3,5,10,13] Result: 2.0^13 = 8192 14: [1,2,3,5,7,14] Result: 2.0^14 = 16384 15: [1,2,3,5,10,15] Result: 2.0^15 = 32768 16: [1,2,4,8,16] Result: 2.0^16 = 65536 17: [1,2,4,8,16,17] Result: 2.0^17 = 131072 81: [1,2,3,5,10,20,40,41,81] Result: 1.1^81 = 2253.2402360440124 191: [1,2,3,5,7,14,19,38,57,95,190,191] Result: 3.0^191 = 1.34945886742810938037281573965238849174025022940301019140667054E+91
#include once "big_int\big_integer.bi"
Const NULL As Any Ptr = 0
Type HashEntry
key As Integer
value As Integer
sgte As HashEntry Ptr
End Type
Type HashMap
buckets(255) As HashEntry Ptr
Declare Sub put(key As Integer, value As Integer)
Declare Function get(key As Integer, defaultValue As Integer = 0) As Integer
Declare Function contains(key As Integer) As Boolean
End Type
Sub HashMap.put(key As Integer, value As Integer)
Dim As Integer bucket = key And 255
Dim As HashEntry Ptr entry = buckets(bucket)
' Check if key already exists
While entry <> NULL
If entry->key = key Then
entry->value = value
Exit Sub
End If
entry = entry->sgte
Wend
' Create new entry
entry = New HashEntry
entry->key = key
entry->value = value
entry->sgte = buckets(bucket)
buckets(bucket) = entry
End Sub
Function HashMap.get(key As Integer, defaultValue As Integer = 0) As Integer
Dim As Integer bucket = key And 255
Dim As HashEntry Ptr entry = buckets(bucket)
While entry <> NULL
If entry->key = key Then Return entry->value
entry = entry->sgte
Wend
Return defaultValue
End Function
Function HashMap.contains(key As Integer) As Boolean
Dim As Integer bucket = key And 255
Dim As HashEntry Ptr entry = buckets(bucket)
While entry <> NULL
If entry->key = key Then Return True
entry = entry->sgte
Wend
Return False
End Function
Type IntArray
dato As Integer Ptr
size As Integer
capacity As Integer
Declare Sub add(value As Integer)
Declare Function get(index As Integer) As Integer
Declare Function getSize() As Integer
Declare Sub clear()
End Type
Sub IntArray.add(value As Integer)
If size >= capacity Then
capacity = Iif(capacity = 0, 8, capacity * 2)
Dim As Integer Ptr newdato = New Integer[capacity]
If dato <> NULL Then
For i As Integer = 0 To size - 1
newdato[i] = dato[i]
Next
Delete[] dato
End If
dato = newdato
End If
dato[size] = value
size += 1
End Sub
Function IntArray.get(index As Integer) As Integer
If index >= 0 And index < size Then Return dato[index]
Return 0
End Function
Function IntArray.getSize() As Integer
Return size
End Function
Sub IntArray.clear()
size = 0
End Sub
' Global variables
Dim Shared As HashMap p
Dim Shared As IntArray levels(0) ' Array of levels
' Function to compute the path to n
Function path(n As Integer) As IntArray
Dim As Integer i, j, x, y, sum
Dim As IntArray result
' Base case
If n = 0 Then Return result
' Check if we already have a path to n
While Not p.contains(n)
Dim q As IntArray
' Process current level
For i = 0 To levels(0).getSize() - 1
x = levels(0).get(i)
' Get path to x
Dim As IntArray xPath = path(x)
' Try to extend the path
For j = 0 To xPath.getSize() - 1
y = xPath.get(j)
sum = x + y
' Check if we already have a path to sum
If p.contains(sum) Then Exit For
' Record the path
p.put(sum, x)
q.add(sum)
Next j
Next i
' Update level
levels(0).clear()
For i = 0 To q.getSize() - 1
levels(0).add(q.get(i))
Next i
' If we can't make progress, break
If q.getSize() = 0 Then Exit While
Wend
' Reconstruct the path
Dim As Integer curr = n
Dim As IntArray tempPath
While curr <> 0
tempPath.add(curr)
curr = p.get(curr)
Wend
' Reverse the path
For i = tempPath.getSize() - 1 To 0 Step -1
result.add(tempPath.get(i))
Next i
Return result
End Function
' Function to compute x^n using the tree method with BigInteger for integers
Function treePowBig(x As Integer, n As Integer) As BigInt
Dim As BigInt r(0 To n)
' Initialize r[0] = 1 and r[1] = x
r(0) = 1
r(1) = x
Dim As Integer i, curr, p = 0
Dim As IntArray pathToN = path(n)
For i = 0 To pathToN.getSize() - 1
curr = pathToN.get(i)
r(curr) = r(curr - p) * r(p)
p = curr
Next i
Return r(n)
End Function
' Function to compute x^n using the tree method with high precision for decimals
Function treePowDecimal(x As Double, n As Integer) As Double
Dim As Double r(0 To n)
' Initialize r[0] = 1 and r[1] = x
r(0) = 1
r(1) = x
Dim As Integer i, curr, p = 0
Dim As IntArray pathToN = path(n)
For i = 0 To pathToN.getSize() - 1
curr = pathToN.get(i)
r(curr) = r(curr - p) * r(p)
p = curr
Next i
Return r(n)
End Function
' Function to display the power calculation
Sub showPow(x As Double, n As Integer)
Dim As IntArray pathToN = path(n)
Dim As String pathStr = ""
For i As Integer = 0 To pathToN.getSize() - 1
If i > 0 Then pathStr &= ", "
pathStr &= Str(pathToN.get(i))
Next i
Print n & ": [" & pathStr & "]"
If Int(x) = x And n > 20 Then
' Use BigInt for large integer powers
Dim As BigInt result = treePowBig(Int(x), n)
Print x & "^" & n & " = " & result
Elseif Int(x) <> x Then
' Use high precision for decimal numbers
Dim As Double result = treePowDecimal(x, n)
Print Using "##.#^## = ####.######"; x; n; result
Else
' Use standard calculation for small integer powers
Dim As Double result = treePowDecimal(x, n)
Print Using "##^## = &"; Int(x); n; result
End If
Print
End Sub
' Initialize
p.put(1, 0)
levels(0).add(1)
' Main program
For i As Integer = 0 To 17
showPow(2, i)
Next i
showPow(1.1, 81)
showPow(3, 191)
Sleep
- Output:
0: [] 2^ 0 = 1 1: [1] 2^ 1 = 2 2: [1, 2] 2^ 2 = 4 3: [1, 2, 3] 2^ 3 = 8 4: [1, 2, 4] 2^ 4 = 16 5: [1, 2, 4, 5] 2^ 5 = 32 6: [1, 2, 4, 6] 2^ 6 = 64 7: [1, 2, 4, 6, 7] 2^ 7 = 128 8: [1, 2, 4, 8] 2^ 8 = 256 9: [1, 2, 4, 8, 9] 2^ 9 = 512 10: [1, 2, 4, 8, 10] 2^10 = 1024 11: [1, 2, 4, 8, 10, 11] 2^11 = 2048 12: [1, 2, 4, 8, 12] 2^12 = 4096 13: [1, 2, 4, 8, 12, 13] 2^13 = 8192 14: [1, 2, 4, 8, 12, 14] 2^14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2^15 = 32768 16: [1, 2, 4, 8, 16] 2^16 = 65536 17: [1, 2, 4, 8, 16, 17] 2^17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.1^81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3^191 = +13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
package main
import (
"fmt"
"math/big"
)
var (
p = map[int]int{1: 0}
lvl = [][]int{[]int{1}}
)
func path(n int) []int {
if n == 0 {
return []int{}
}
for {
if _, ok := p[n]; ok {
break
}
var q []int
for _, x := range lvl[0] {
for _, y := range path(x) {
z := x + y
if _, ok := p[z]; ok {
break
}
p[z] = x
q = append(q, z)
}
}
lvl[0] = q
}
r := path(p[n])
r = append(r, n)
return r
}
func treePow(x float64, n int) *big.Float {
r := map[int]*big.Float{0: big.NewFloat(1), 1: big.NewFloat(x)}
p := 0
for _, i := range path(n) {
temp := new(big.Float).SetPrec(320)
temp.Mul(r[i-p], r[p])
r[i] = temp
p = i
}
return r[n]
}
func showPow(x float64, n int, isIntegral bool) {
fmt.Printf("%d: %v\n", n, path(n))
f := "%f"
if isIntegral {
f = "%.0f"
}
fmt.Printf(f, x)
fmt.Printf(" ^ %d = ", n)
fmt.Printf(f+"\n\n", treePow(x, n))
}
func main() {
for n := 0; n <= 17; n++ {
showPow(2, n, true)
}
showPow(1.1, 81, false)
showPow(3, 191, true)
}
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1 2] 2 ^ 2 = 4 3: [1 2 3] 2 ^ 3 = 8 4: [1 2 4] 2 ^ 4 = 16 5: [1 2 4 5] 2 ^ 5 = 32 6: [1 2 4 6] 2 ^ 6 = 64 7: [1 2 4 6 7] 2 ^ 7 = 128 8: [1 2 4 8] 2 ^ 8 = 256 9: [1 2 4 8 9] 2 ^ 9 = 512 10: [1 2 4 8 10] 2 ^ 10 = 1024 11: [1 2 4 8 10 11] 2 ^ 11 = 2048 12: [1 2 4 8 12] 2 ^ 12 = 4096 13: [1 2 4 8 12 13] 2 ^ 13 = 8192 14: [1 2 4 8 12 14] 2 ^ 14 = 16384 15: [1 2 4 8 12 14 15] 2 ^ 15 = 32768 16: [1 2 4 8 16] 2 ^ 16 = 65536 17: [1 2 4 8 16 17] 2 ^ 17 = 131072 81: [1 2 4 8 16 32 64 80 81] 1.100000 ^ 81 = 2253.240236 191: [1 2 4 8 16 32 64 128 160 176 184 188 190 191] 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
class PowerTree {
private static Map<Integer, Integer> p = new HashMap<>()
private static List<List<Integer>> lvl = new ArrayList<>()
static {
p[1] = 0
List<Integer> temp = new ArrayList<Integer>()
temp.add 1
lvl.add temp
}
private static List<Integer> path(int n) {
if (n == 0) return new ArrayList<Integer>()
while (!p.containsKey(n)) {
List<Integer> q = new ArrayList<>()
for (Integer x in lvl.get(0)) {
for (Integer y in path(x)) {
if (p.containsKey(x + y)) break
p[x + y] = x
q.add x + y
}
}
lvl[0].clear()
lvl[0].addAll q
}
List<Integer> temp = path p[n]
temp.add n
temp
}
private static BigDecimal treePow(double x, int n) {
Map<Integer, BigDecimal> r = new HashMap<>()
r[0] = BigDecimal.ONE
r[1] = BigDecimal.valueOf(x)
int p = 0
for (Integer i in path(n)) {
r[i] = r[i - p] * r[p]
p = i
}
r[n]
}
private static void showPos(double x, int n, boolean isIntegral) {
printf("%d: %s\n", n, path(n))
String f = isIntegral ? "%.0f" : "%f"
printf(f, x)
printf(" ^ %d = ", n)
printf(f, treePow(x, n))
println()
println()
}
static void main(String[] args) {
for (int n = 0; n <= 17; ++n) {
showPos 2.0, n, true
}
showPos 1.1, 81, false
showPos 3.0, 191, true
}
}
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1, 2] 2 ^ 2 = 4 3: [1, 2, 3] 2 ^ 3 = 8 4: [1, 2, 4] 2 ^ 4 = 16 5: [1, 2, 4, 5] 2 ^ 5 = 32 6: [1, 2, 4, 6] 2 ^ 6 = 64 7: [1, 2, 4, 6, 7] 2 ^ 7 = 128 8: [1, 2, 4, 8] 2 ^ 8 = 256 9: [1, 2, 4, 8, 9] 2 ^ 9 = 512 10: [1, 2, 4, 8, 10] 2 ^ 10 = 1024 11: [1, 2, 4, 8, 10, 11] 2 ^ 11 = 2048 12: [1, 2, 4, 8, 12] 2 ^ 12 = 4096 13: [1, 2, 4, 8, 12, 13] 2 ^ 13 = 8192 14: [1, 2, 4, 8, 12, 14] 2 ^ 14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2 ^ 15 = 32768 16: [1, 2, 4, 8, 16] 2 ^ 16 = 65536 17: [1, 2, 4, 8, 16, 17] 2 ^ 17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.100000 ^ 81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
{-# LANGUAGE ScopedTypeVariables #-}
module Rosetta.PowerTree
( Natural
, powerTree
, power
) where
import Data.Foldable (toList)
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import Data.Maybe (fromMaybe)
import Data.List (foldl')
import Data.Sequence (Seq (..), (|>))
import qualified Data.Sequence as Seq
import Numeric.Natural (Natural)
type M = Map Natural S
type S = Seq Natural
levels :: [M]
levels = let s = Seq.singleton 1 in fst <$> iterate step (Map.singleton 1 s, s)
step :: (M, S) -> (M, S)
step (m, xs) = foldl' f (m, Empty) xs
where
f :: (M, S) -> Natural -> (M, S)
f (m', ys) n = foldl' g (m', ys) ns
where
ns :: S
ns = m' Map.! n
g :: (M, S) -> Natural -> (M, S)
g (m'', zs) k =
let l = n + k
in case Map.lookup l m'' of
Nothing -> (Map.insert l (ns |> l) m'', zs |> l)
Just _ -> (m'', zs)
powerTree :: Natural -> [Natural]
powerTree n
| n <= 0 = []
| otherwise = go levels
where
go :: [M] -> [Natural]
go [] = error "impossible branch"
go (m : ms) = fromMaybe (go ms) $ toList <$> Map.lookup n m
power :: forall a. Num a => a -> Natural -> a
power _ 0 = 1
power a n = go a 1 (Map.singleton 1 a) $ tail $ powerTree n
where
go :: a -> Natural -> Map Natural a -> [Natural] -> a
go b _ _ [] = b
go b k m (l : ls) =
let b' = b * m Map.! (l - k)
m' = Map.insert l b' m
in go b' l m' ls
- Output:
(The CReal type from package numbers is used to get the exact result for the last example.)
powerTree 0 = [], power 2 0 = 1 powerTree 1 = [1], power 2 1 = 2 powerTree 2 = [1,2], power 2 2 = 4 powerTree 3 = [1,2,3], power 2 3 = 8 powerTree 4 = [1,2,4], power 2 4 = 16 powerTree 5 = [1,2,3,5], power 2 5 = 32 powerTree 6 = [1,2,3,6], power 2 6 = 64 powerTree 7 = [1,2,3,5,7], power 2 7 = 128 powerTree 8 = [1,2,4,8], power 2 8 = 256 powerTree 9 = [1,2,3,6,9], power 2 9 = 512 powerTree 10 = [1,2,3,5,10], power 2 10 = 1024 powerTree 11 = [1,2,3,5,10,11], power 2 11 = 2048 powerTree 12 = [1,2,3,6,12], power 2 12 = 4096 powerTree 13 = [1,2,3,5,10,13], power 2 13 = 8192 powerTree 14 = [1,2,3,5,7,14], power 2 14 = 16384 powerTree 15 = [1,2,3,5,10,15], power 2 15 = 32768 powerTree 16 = [1,2,4,8,16], power 2 16 = 65536 powerTree 17 = [1,2,4,8,16,17], power 2 17 = 131072 powerTree 191 = [1,2,3,5,7,14,19,38,57,95,190,191], power 3 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 powerTree 81 = [1,2,3,5,10,20,40,41,81], power 1.1 81 = 2253.240236044012487937308538033349567966729852481170503814810577345406584190098644811
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.)
knuth_power_tree=:3 :0
L=: P=: %(1+y){._ 1
findpath=: ]
while. _ e.P do.
for_n.(/: findpath&>)I.L=>./L-._ 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.
findpath=: [: |. {&P^:a:
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
)
Task examples:
knuth_power_tree 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
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:
90j83 ": (x:1.1) usepath 81
2253.24023604401248793730853803334956796672985248117050381481057734540658419009864481100
import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class PowerTree {
private static Map<Integer, Integer> p = new HashMap<>();
private static List<List<Integer>> lvl = new ArrayList<>();
static {
p.put(1, 0);
ArrayList<Integer> temp = new ArrayList<>();
temp.add(1);
lvl.add(temp);
}
private static List<Integer> path(int n) {
if (n == 0) return new ArrayList<>();
while (!p.containsKey(n)) {
List<Integer> q = new ArrayList<>();
for (Integer x : lvl.get(0)) {
for (Integer y : path(x)) {
if (p.containsKey(x + y)) break;
p.put(x + y, x);
q.add(x + y);
}
}
lvl.get(0).clear();
lvl.get(0).addAll(q);
}
List<Integer> temp = path(p.get(n));
temp.add(n);
return temp;
}
private static BigDecimal treePow(double x, int n) {
Map<Integer, BigDecimal> r = new HashMap<>();
r.put(0, BigDecimal.ONE);
r.put(1, BigDecimal.valueOf(x));
int p = 0;
for (Integer i : path(n)) {
r.put(i, r.get(i - p).multiply(r.get(p)));
p = i;
}
return r.get(n);
}
private static void showPow(double x, int n, boolean isIntegral) {
System.out.printf("%d: %s\n", n, path(n));
String f = isIntegral ? "%.0f" : "%f";
System.out.printf(f, x);
System.out.printf(" ^ %d = ", n);
System.out.printf(f, treePow(x, n));
System.out.println("\n");
}
public static void main(String[] args) {
for (int n = 0; n <= 17; ++n) {
showPow(2.0, n, true);
}
showPow(1.1, 81, false);
showPow(3.0, 191, true);
}
}
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1, 2] 2 ^ 2 = 4 3: [1, 2, 3] 2 ^ 3 = 8 4: [1, 2, 4] 2 ^ 4 = 16 5: [1, 2, 4, 5] 2 ^ 5 = 32 6: [1, 2, 4, 6] 2 ^ 6 = 64 7: [1, 2, 4, 6, 7] 2 ^ 7 = 128 8: [1, 2, 4, 8] 2 ^ 8 = 256 9: [1, 2, 4, 8, 9] 2 ^ 9 = 512 10: [1, 2, 4, 8, 10] 2 ^ 10 = 1024 11: [1, 2, 4, 8, 10, 11] 2 ^ 11 = 2048 12: [1, 2, 4, 8, 12] 2 ^ 12 = 4096 13: [1, 2, 4, 8, 12, 13] 2 ^ 13 = 8192 14: [1, 2, 4, 8, 12, 14] 2 ^ 14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2 ^ 15 = 32768 16: [1, 2, 4, 8, 16] 2 ^ 16 = 65536 17: [1, 2, 4, 8, 16, 17] 2 ^ 17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.100000 ^ 81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
class PowerTree {
static p = new Map();
static lvl = [];
static {
PowerTree.p.set(1, 0);
let temp = [1];
PowerTree.lvl.push(temp);
}
static path(n) {
if (n === 0) return [];
while (!PowerTree.p.has(n)) {
let q = [];
for (let x of PowerTree.lvl[0]) {
for (let y of PowerTree.path(x)) {
if (PowerTree.p.has(x + y)) break;
PowerTree.p.set(x + y, x);
q.push(x + y);
}
}
PowerTree.lvl[0] = [...q]; // Important: create a new array! Otherwise, the original array is mutated.
}
let temp = PowerTree.path(PowerTree.p.get(n));
temp.push(n);
return temp;
}
static treePow(x, n) {
const r = new Map();
r.set(0, 1n); // Changed to bigint for accurate integer results
r.set(1, BigInt(Math.round(x))); // Convert to BigInt
if (Number.isInteger(x) && Number.isInteger(n) && n>30){
r.set(0, 1n);
r.set(1, BigInt(x));
} else {
r.set(0, 1);
r.set(1, x);
}
let p = 0;
const path_result = PowerTree.path(n);
for (let i of path_result) {
let a = r.get(i - p);
let b = r.get(p);
if (Number.isInteger(x) && Number.isInteger(n) && n>30){
r.set(i, a * b);
} else {
if (typeof a === 'bigint' && typeof b === 'bigint'){
r.set(i, Number(a) * Number(b));
} else if (typeof a === 'bigint'){
r.set(i, Number(a) * b);
} else if (typeof b === 'bigint'){
r.set(i, a * Number(b));
} else {
r.set(i, a * b);
}
}
p = i;
}
if (Number.isInteger(x) && Number.isInteger(n) && n>30){
return Number(r.get(n));
}
return r.get(n);
}
static showPow(x, n, isIntegral) {
console.log(`${n}: ${PowerTree.path(n)}`);
let f = isIntegral ? "%.0f" : "%f";
//console.log(f, x); // Javascript doesn't have printf-style formatting built in
let formatted_x;
if (isIntegral) {
formatted_x = Math.round(x);
} else {
formatted_x = x.toFixed(6); // Or whatever precision you need.
}
console.log(`${formatted_x} ^ ${n} = ${PowerTree.treePow(x, n)}`);
console.log("\n");
}
static main() {
for (let n = 0; n <= 17; ++n) {
PowerTree.showPow(2.0, n, true);
}
PowerTree.showPow(1.1, 81, false);
PowerTree.showPow(3.0, 191, true);
}
}
PowerTree.main();
- Output:
0: 2 ^ 0 = 1 1: 1 2 ^ 1 = 2 2: 1,2 2 ^ 2 = 4 3: 1,2,3 2 ^ 3 = 8 4: 1,2,4 2 ^ 4 = 16 5: 1,2,4,5 2 ^ 5 = 32 6: 1,2,4,6 2 ^ 6 = 64 7: 1,2,4,6,7 2 ^ 7 = 128 8: 1,2,4,8 2 ^ 8 = 256 9: 1,2,4,8,9 2 ^ 9 = 512 10: 1,2,4,8,10 2 ^ 10 = 1024 11: 1,2,4,8,10,11 2 ^ 11 = 2048 12: 1,2,4,8,12 2 ^ 12 = 4096 13: 1,2,4,8,12,13 2 ^ 13 = 8192 14: 1,2,4,8,12,14 2 ^ 14 = 16384 15: 1,2,4,8,12,14,15 2 ^ 15 = 32768 16: 1,2,4,8,16 2 ^ 16 = 65536 17: 1,2,4,8,16,17 2 ^ 17 = 131072 81: 1,2,4,8,16,32,64,80,81 1.100000 ^ 81 = 2253.240236044025 191: 1,2,4,8,16,32,64,128,160,176,184,188,190,191 3 ^ 191 = 1.3494588674281094e+91
Adapted from #Wren
Works with jq (*)
Works with gojq, the Go implementation of jq
(*) Use gojq for infinite-precision integer arithmetic.
# Input: {p, lvl, path}
def kpath($n):
if $n == 0 then .path=[]
else until( .p[$n|tostring];
.q = []
| reduce .lvl[0][] as $x (.;
kpath($x)
| label $out
| foreach (.path[], null) as $y (.;
if $y == null then .
else (($x + $y)|tostring) as $xy
| if .p[$xy] then .
else .p[$xy] = $x
| .q += [$x + $y]
end
end;
select($y == null) ) )
| .lvl[0] = .q )
| kpath(.p[$n|tostring])
| .path += [$n]
end ;
# Input: as for kpath
def treePow($x; $n):
reduce kpath($n).path[] as $i (
{r: { "0": 1, "1": $x }, pp: 0 };
.r[$i|tostring] = .r[($i - .pp)|tostring] * .r[.pp|tostring]
| .pp = $i )
| .r[$n|tostring] ;
def showPow($x; $n):
{ p: {"1": 0},
lvl: [[1]],
path: []}
| "\($n): \(kpath($n).path)",
"\($x) ^ \($n) = \(treePow($x; $n))";
(range(0;18) as $n | showPow(2; $n)),
showPow(1.1; 81),
showPow(3; 191)- Output:
Using gojq, the output is the same as for Wren.
Module:
module KnuthPowerTree
const p = Dict(1 => 0)
const lvl = [[1]]
function path(n)
global p, lvl
iszero(n) && return Int[]
while n ∉ keys(p)
q = Int[]
for x in lvl[1], y in path(x)
if (x + y) ∉ keys(p)
p[x + y] = x
push!(q, x + y)
end
end
lvl[1] = q
end
return push!(path(p[n]), n)
end
function pow(x::Number, n::Integer)
r = Dict{typeof(n), typeof(x)}(0 => 1, 1 => x)
p = 0
for i in path(n)
r[i] = r[i - p] * r[p]
p = i
end
return r[n]
end
end # module KnuthPowerTree
Main:
using .KnuthPowerTree: path, pow
for n in 0:17
println("2 ^ $n:\n - path: ", join(path(n), ", "), "\n - result: ", pow(2, n))
end
for (x, n) in ((big(3), 191), (1.1, 81))
println("$x ^ $n:\n - path: ", join(path(n), ", "), "\n - result: ", pow(x, n))
end
- Output:
2 ^ 0: - path: - result: 1 2 ^ 1: - path: 1 - result: 2 2 ^ 2: - path: 1, 2 - result: 4 2 ^ 3: - path: 1, 2, 3 - result: 8 2 ^ 4: - path: 1, 2, 4 - result: 16 2 ^ 5: - path: 1, 2, 3, 5 - result: 32 2 ^ 6: - path: 1, 2, 3, 6 - result: 64 2 ^ 7: - path: 1, 2, 3, 5, 7 - result: 128 2 ^ 8: - path: 1, 2, 4, 8 - result: 256 2 ^ 9: - path: 1, 2, 3, 6, 9 - result: 512 2 ^ 10: - path: 1, 2, 3, 5, 10 - result: 1024 2 ^ 11: - path: 1, 2, 3, 5, 10, 11 - result: 2048 2 ^ 12: - path: 1, 2, 3, 6, 12 - result: 4096 2 ^ 13: - path: 1, 2, 3, 5, 10, 13 - result: 8192 2 ^ 14: - path: 1, 2, 3, 5, 7, 14 - result: 16384 2 ^ 15: - path: 1, 2, 3, 5, 10, 15 - result: 32768 2 ^ 16: - path: 1, 2, 4, 8, 16 - result: 65536 2 ^ 17: - path: 1, 2, 4, 8, 16, 17 - result: 131072 3 ^ 191: - path: 1, 2, 3, 5, 7, 14, 19, 38, 57, 95, 190, 191 - result: 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 1.1 ^ 81: - path: 1, 2, 3, 5, 10, 20, 40, 41, 81 - result: 2253.2402360440283
// version 1.1.3
import java.math.BigDecimal
var p = mutableMapOf(1 to 0)
var lvl = mutableListOf(listOf(1))
fun path(n: Int): List<Int> {
if (n == 0) return emptyList<Int>()
while (n !in p) {
val q = mutableListOf<Int>()
for (x in lvl[0]) {
for (y in path(x)) {
if ((x + y) in p) break
p[x + y] = x
q.add(x + y)
}
}
lvl[0] = q
}
return path(p[n]!!) + n
}
fun treePow(x: Double, n: Int): BigDecimal {
val r = mutableMapOf(0 to BigDecimal.ONE, 1 to BigDecimal(x.toString()))
var p = 0
for (i in path(n)) {
r[i] = r[i - p]!! * r[p]!!
p = i
}
return r[n]!!
}
fun showPow(x: Double, n: Int, isIntegral: Boolean = true) {
println("$n: ${path(n)}")
val f = if (isIntegral) "%.0f" else "%f"
println("${f.format(x)} ^ $n = ${f.format(treePow(x, n))}\n")
}
fun main(args: Array<String>) {
for (n in 0..17) showPow(2.0, n)
showPow(1.1, 81, false)
showPow(3.0, 191)
}
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1, 2] 2 ^ 2 = 4 3: [1, 2, 3] 2 ^ 3 = 8 4: [1, 2, 4] 2 ^ 4 = 16 5: [1, 2, 4, 5] 2 ^ 5 = 32 6: [1, 2, 4, 6] 2 ^ 6 = 64 7: [1, 2, 4, 6, 7] 2 ^ 7 = 128 8: [1, 2, 4, 8] 2 ^ 8 = 256 9: [1, 2, 4, 8, 9] 2 ^ 9 = 512 10: [1, 2, 4, 8, 10] 2 ^ 10 = 1024 11: [1, 2, 4, 8, 10, 11] 2 ^ 11 = 2048 12: [1, 2, 4, 8, 12] 2 ^ 12 = 4096 13: [1, 2, 4, 8, 12, 13] 2 ^ 13 = 8192 14: [1, 2, 4, 8, 12, 14] 2 ^ 14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2 ^ 15 = 32768 16: [1, 2, 4, 8, 16] 2 ^ 16 = 65536 17: [1, 2, 4, 8, 16, 17] 2 ^ 17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.100000 ^ 81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
ClearAll[NextStep, TreePow]
NextStep[pows_List] := Module[{maxlen, sel, new, vals, knows},
maxlen = Max[Length /@ pows[[All, "Path"]]];
sel = Select[pows, Length[#["Path"]] == maxlen &];
knows = pows[[All, "P"]];
new = {};
Do[
vals = s["P"] + s["Path"];
vals = DeleteCases[vals, Alternatives @@ Join[s["Path"], knows]];
new =
Join[
new, <|"Path" -> Append[s["Path"], #], "P" -> #|> & /@ vals];
,
{s, sel}
];
new //= DeleteDuplicatesBy[#["P"] &];
SortBy[Join[pows, new], #["P"] &]
]
TreePow[path_List, base_] := Module[{db, tups},
db = <|1 -> base|>;
Do[
tups = Tuples[Keys[db], 2];
tups = Select[tups, #[[2]] >= #[[1]] &];
tups = Select[tups, Total[#] == next &];
If[Length[tups] < 1, Abort[]];
tups //= First;
AssociateTo[db, Total[tups] -> (Times @@ (db /@ tups))]
,
{next, Rest[path]}
];
db[Last[path]]
]
pows = {<|"Path" -> {1}, "P" -> 1|>};
steps = Nest[NextStep, pows, 7];
LayeredGraphPlot[DirectedEdge @@@ steps[[2 ;;, "Path", -2 ;;]], VertexLabels -> Automatic]
pows = {<|"Path" -> {1}, "P" -> 1|>};
steps = Nest[NextStep, pows, 5];
assoc = Association[#["P"] -> #["Path"] & /@ steps];
Dataset[assoc]
TreePow[assoc[#], 2] & /@ Range[1, 17]
pows = {<|"Path" -> {1}, "P" -> 1|>};
steps = NestWhile[NextStep, pows, Not[MemberQ[#[[All, "P"]], 191]] &];
SelectFirst[steps, #["P"] == 191 &]["Path"];
TreePow[%, 3]
pows = {<|"Path" -> {1}, "P" -> 1|>};
steps = NestWhile[NextStep, pows, Not[MemberQ[#[[All, "P"]], 81]] &];
SelectFirst[steps, #["P"] == 81 &]["Path"];
TreePow[%, 1.1]
- Output:
[Graphics object showing the tree]
1 {1}
2 {1,2}
3 {1,2,3}
4 {1,2,4}
5 {1,2,3,5}
6 {1,2,3,6}
7 {1,2,3,5,7}
8 {1,2,4,8}
9 {1,2,3,6,9}
10 {1,2,3,5,10}
11 {1,2,3,6,9,11}
12 {1,2,3,6,12}
13 {1,2,3,5,10,13}
14 {1,2,3,5,7,14}
15 {1,2,3,6,9,15}
16 {1,2,4,8,16}
17 {1,2,4,8,16,17}
18 {1,2,3,6,9,18}
20 {1,2,3,5,10,20}
24 {1,2,3,6,12,24}
... ...
{2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072}
13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
2253.24
This is a translation of the Kotlin/Python program. The algorithm is the same but there are some changes: replacement of global variables by a tree object, use of longer names, replacement of the "lvl" list of lists by a simple list, use of generics to handle the case of big integers...
import tables
import bignum
type
Id = Natural # Node identifier (= exponent value).
Tree = object
parents: Table[Id, Id] # Mapping node id -> parent node id.
lastLevel: seq[Id] # List of node ids in current last level.
func initTree(): Tree =
## Return an initialized tree.
const Root = Id(1)
Tree(parents: {Root: Id(0)}.toTable, lastLevel: @[Root])
func path(tree: var Tree; id: Id): seq[Id] =
## Return the path to node with given id.
if id == 0: return
while id notin tree.parents:
# Node "id" not yet present in the tree: build a new level.
var newLevel: seq[Id]
for x in tree.lastLevel:
for y in tree.path(x):
let newId = x + y
if newId in tree.parents: break # Node already created.
# Create a new node.
tree.parents[newId] = x
newLevel.add newId
tree.lastLevel = move(newLevel)
# Node path is the concatenation of parent node path and node id.
result = tree.path(tree.parents[id]) & id
func treePow[T: SomeNumber | Int](tree: var Tree; x: T; n: Natural): T =
## Compute x^n using the power tree.
let one = when T is Int: newInt(1) else: T(1)
var results = {0: one, 1: x}.toTable # Intermediate and last results.
var k = 0
for i in tree.path(n):
results[i] = results[i - k] * results[k]
k = i
return results[n]
proc showPow[T: SomeNumber | Int](tree: var Tree; x: T; n: Natural) =
echo n, " → ", ($tree.path(n))[1..^1]
let result = tree.treePow(x, n)
echo x, "^", n, " = ", result
when isMainModule:
var tree = initTree()
for n in 0..17: tree.showPow(2, n)
echo ""
tree.showPow(1.1, 81)
echo ""
tree.showPow(newInt(3), 191)
- Output:
0 → [] 2^0 = 1 1 → [1] 2^1 = 2 2 → [1, 2] 2^2 = 4 3 → [1, 2, 3] 2^3 = 8 4 → [1, 2, 4] 2^4 = 16 5 → [1, 2, 4, 5] 2^5 = 32 6 → [1, 2, 4, 6] 2^6 = 64 7 → [1, 2, 4, 6, 7] 2^7 = 128 8 → [1, 2, 4, 8] 2^8 = 256 9 → [1, 2, 4, 8, 9] 2^9 = 512 10 → [1, 2, 4, 8, 10] 2^10 = 1024 11 → [1, 2, 4, 8, 10, 11] 2^11 = 2048 12 → [1, 2, 4, 8, 12] 2^12 = 4096 13 → [1, 2, 4, 8, 12, 13] 2^13 = 8192 14 → [1, 2, 4, 8, 12, 14] 2^14 = 16384 15 → [1, 2, 4, 8, 12, 14, 15] 2^15 = 32768 16 → [1, 2, 4, 8, 16] 2^16 = 65536 17 → [1, 2, 4, 8, 16, 17] 2^17 = 131072 81 → [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.1^81 = 2253.240236044025 191 → [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
uses NumLibABC;
var
p := Dict((1, 0));
lvl := Lst(Lst(word(1)));
function tofraction(s: string): fraction;
begin
var period := s.IndexOf('.');
if period = -1 then result := Frc(s.ToBigInteger)
else result := Frc(s.Remove(period, 1).ToBigInteger, Power(10bi, s.Length - period - 1));
end;
function path(n: word): List<word>;
begin
result := new List<word>;
if n = 0 then exit;
while not p.ContainsKey(n) do
begin
var q := new List<word>;
foreach var x in lvl[0] do
foreach var y in path(x) do
begin
if (x + y) in p then break;
p[x + y] := x;
q.add(x + y);
end;
lvl[0] := q;
end;
result := path(p[n]) + Lst(n);
end;
function treePow(x: real; n: word): Fraction;
begin
var frac := tofraction(FloatToStr(x));
var r := Dict((0, Frc(1)), (1, frac));
var p := 0;
foreach var i in path(n) do
begin
r[i] := r[i - p] * r[p];
p := i;
end;
result := r[n];
end;
procedure showPow(x: real; n: word);
begin
writeln(n, ': ', path(n));
writeln(x, '^', n, ' = ', treePow(x, n), #10);
end;
begin
for var n := 0 to 17 do showPow(2, n);
showPow(1.1, 81);
showPow(3, 191);
end.
- Output:
0: [] 2^0 = 1 1: [1] 2^1 = 2 2: [1,2] 2^2 = 4 3: [1,2,3] 2^3 = 8 4: [1,2,4] 2^4 = 16 5: [1,2,4,5] 2^5 = 32 6: [1,2,4,6] 2^6 = 64 7: [1,2,4,6,7] 2^7 = 128 8: [1,2,4,8] 2^8 = 256 9: [1,2,4,8,9] 2^9 = 512 10: [1,2,4,8,10] 2^10 = 1024 11: [1,2,4,8,10,11] 2^11 = 2048 12: [1,2,4,8,12] 2^12 = 4096 13: [1,2,4,8,12,13] 2^13 = 8192 14: [1,2,4,8,12,14] 2^14 = 16384 15: [1,2,4,8,12,14,15] 2^15 = 32768 16: [1,2,4,8,16] 2^16 = 65536 17: [1,2,4,8,16,17] 2^17 = 131072 81: [1,2,4,8,16,32,64,80,81] 1.1^81 = 2253240236044012487937308538033349567966729852481170503814810577345406584190098644811/1000000000000000000000000000000000000000000000000000000000000000000000000000000000 191: [1,2,4,8,16,32,64,128,160,176,184,188,190,191] 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
my @lvl = [1];
my %p = (1 => 0);
sub path {
my ($n) = @_;
return () if ($n == 0);
until (exists $p{$n}) {
my @q;
foreach my $x (@{$lvl[0]}) {
foreach my $y (path($x)) {
my $z = $x + $y;
last if exists($p{$z});
$p{$z} = $x;
push @q, $z;
}
}
$lvl[0] = \@q;
}
(path($p{$n}), $n);
}
sub tree_pow {
my ($x, $n) = @_;
my %r = (0 => 1, 1 => $x);
my $p = 0;
foreach my $i (path($n)) {
$r{$i} = $r{$i - $p} * $r{$p};
$p = $i;
}
$r{$n};
}
sub show_pow {
my ($x, $n) = @_;
my $fmt = "%d: %s\n" . ("%g^%s = %f", "%s^%s = %s")[$x == int($x)] . "\n";
printf($fmt, $n, "(" . join(" ", path($n)) . ")", $x, $n, tree_pow($x, $n));
}
show_pow(2, $_) for 0 .. 17;
show_pow(1.1, 81);
{
use bigint (try => 'GMP');
show_pow(3, 191);
}
- Output:
0: () 2^0 = 1 1: (1) 2^1 = 2 2: (1 2) 2^2 = 4 3: (1 2 3) 2^3 = 8 4: (1 2 4) 2^4 = 16 5: (1 2 4 5) 2^5 = 32 6: (1 2 4 6) 2^6 = 64 7: (1 2 4 6 7) 2^7 = 128 8: (1 2 4 8) 2^8 = 256 9: (1 2 4 8 9) 2^9 = 512 10: (1 2 4 8 10) 2^10 = 1024 11: (1 2 4 8 10 11) 2^11 = 2048 12: (1 2 4 8 12) 2^12 = 4096 13: (1 2 4 8 12 13) 2^13 = 8192 14: (1 2 4 8 12 14) 2^14 = 16384 15: (1 2 4 8 12 14 15) 2^15 = 32768 16: (1 2 4 8 16) 2^16 = 65536 17: (1 2 4 8 16 17) 2^17 = 131072 81: (1 2 4 8 16 32 64 80 81) 1.1^81 = 2253.240236 191: (1 2 4 8 16 32 64 128 160 176 184 188 190 191) 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
with javascript_semantics constant p = new_dict({{1,0}}) sequence lvl = {1} function path(integer n) if n=0 then return {} end if while getd_index(n,p)=NULL do sequence q = {} for i=1 to length(lvl) do integer x = lvl[i] sequence px = path(x) for j=1 to length(px) do integer y = x+px[j] if getd_index(y,p)!=NULL then exit end if setd(y,x,p) q &= y end for end for lvl = q end while return path(getd(n,p))&n end function include mpfr.e mpfr_set_default_precision(500) function treepow(object x, integer n, sequence pn = {}) -- x can be atom or string (but not mpfr) -- (asides: sequence r uses out-by-1 indexing, ie r[1] is for 0. -- sequence c is used to double-check we are not trying -- to use something which has not yet been calculated.) if pn={} then pn=path(n) end if sequence r = {mpfr_init(1),mpfr_init(x)}, c = {1,1}&repeat(0,max(0,n-1)) for i=1 to max(0,n-1) do r &= mpfr_init() end for integer p = 0 for i=1 to length(pn) do integer pi = pn[i] if c[pi-p+1]=0 then ?9/0 end if if c[p+1]=0 then ?9/0 end if mpfr_mul(r[pi+1],r[pi-p+1],r[p+1]) c[pi+1] = 1 p = pi end for -- string res = shorten(trim_tail(mpfr_sprintf("%.83Rf",r[n+1]),".0")) string res = shorten(trim_tail(mpfr_get_fixed(r[n+1],83),".0")) r = mpfr_free(r) return res end function procedure showpow(object x, integer n) sequence pn = path(n) string xs = iff(string(x)?x:sprintf("%3g",x)) printf(1,"%48v : %3s ^ %d = %s\n", {pn,xs,n,treepow(x,n,pn)}) end procedure for n=0 to 17 do showpow(2,n) end for showpow("1.1",81) showpow(3,191)
- Output:
{} : 2 ^ 0 = 1
{1} : 2 ^ 1 = 2
{1,2} : 2 ^ 2 = 4
{1,2,3} : 2 ^ 3 = 8
{1,2,4} : 2 ^ 4 = 16
{1,2,4,5} : 2 ^ 5 = 32
{1,2,4,6} : 2 ^ 6 = 64
{1,2,4,6,7} : 2 ^ 7 = 128
{1,2,4,8} : 2 ^ 8 = 256
{1,2,4,8,9} : 2 ^ 9 = 512
{1,2,4,8,10} : 2 ^ 10 = 1024
{1,2,4,8,10,11} : 2 ^ 11 = 2048
{1,2,4,8,12} : 2 ^ 12 = 4096
{1,2,4,8,12,13} : 2 ^ 13 = 8192
{1,2,4,8,12,14} : 2 ^ 14 = 16384
{1,2,4,8,12,14,15} : 2 ^ 15 = 32768
{1,2,4,8,16} : 2 ^ 16 = 65536
{1,2,4,8,16,17} : 2 ^ 17 = 131072
{1,2,4,8,16,32,64,80,81} : 1.1 ^ 81 = 2253.240236044012487937308538033349567966729852481170503814810577345406584190098644811
{1,2,4,8,16,32,64,128,160,176,184,188,190,191} : 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
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)
- 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
Using the mpmath library
""" https://rosettacode.org/wiki/Knuth%27s_power_tree """
from typing import List, Dict
from mpmath import mp
# Set decimal precision for mpmath calculations (like Go example 320 bits)
mp.dps = 140
# Global variables
p: Dict[int, int] = {1: 0}
lvl: List[List[int]] = [[1]]
def path(n: int) -> List[int]:
""" return path to the node with id ''n'' """
if n == 0:
return []
while n not in p:
q: List[int] = []
for x in lvl[0]:
for y in path(x):
z = x + y
if z in p:
break
p[z] = x
q.append(z)
lvl[0] = q
r = path(p[n]) + [n]
return r
def tree_pow(x: float, n: int) -> mp.mpf:
""" use power tree to compute x ** n """
r: Dict[int, mp.mpf] = {0: mp.mpf('1'), 1: mp.mpf(str(x))}
prev = 0
for i in path(n):
r[i] = mp.mpf(r[i - prev] * r[prev])
prev = i
return r[n]
def show_pow(x: float, n: int) -> None:
""" print the power calculation """
print(f"{n}: {path(n)}")
y = tree_pow(x, n)
ans = mp.nstr(y, n=100, strip_zeros=True)[
:-2] if mp.isint(y) else mp.nstr(y, n=10)
print(f"{mp.nstr(x)} ^ {n} = {ans}\n")
if __name__ == "__main__":
for expo in range(18):
show_pow(2, expo)
show_pow(1.1, 81)
show_pow(3.0, 191)
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1, 2] 2 ^ 2 = 4 3: [1, 2, 3] 2 ^ 3 = 8 4: [1, 2, 4] 2 ^ 4 = 16 5: [1, 2, 4, 5] 2 ^ 5 = 32 6: [1, 2, 4, 6] 2 ^ 6 = 64 7: [1, 2, 4, 6, 7] 2 ^ 7 = 128 8: [1, 2, 4, 8] 2 ^ 8 = 256 9: [1, 2, 4, 8, 9] 2 ^ 9 = 512 10: [1, 2, 4, 8, 10] 2 ^ 10 = 1024 11: [1, 2, 4, 8, 10, 11] 2 ^ 11 = 2048 12: [1, 2, 4, 8, 12] 2 ^ 12 = 4096 13: [1, 2, 4, 8, 12, 13] 2 ^ 13 = 8192 14: [1, 2, 4, 8, 12, 14] 2 ^ 14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2 ^ 15 = 32768 16: [1, 2, 4, 8, 16] 2 ^ 16 = 65536 17: [1, 2, 4, 8, 16, 17] 2 ^ 17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.1 ^ 81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3.0 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
#lang racket
(define pow-path-cache (make-hash '((0 . (0)) (1 . (0 1)))))
(define pow-path-level '(1))
(define (pow-path-extend!)
(define next-level
(for*/fold ([next-level '()])
([x (in-list pow-path-level)]
[y (in-list (pow-path x))]
[s (in-value (+ x y))]
#:when (not (hash-has-key? pow-path-cache s)))
(hash-set! pow-path-cache s (append (hash-ref pow-path-cache x) (list s)))
(cons s next-level)))
(set! pow-path-level (reverse next-level)))
(define (pow-path n)
(let loop ()
(unless (hash-has-key? pow-path-cache n)
(pow-path-extend!)
(loop)))
(hash-ref pow-path-cache n))
(define (pow-tree x n)
(define pows (make-hash `((0 . 1) (1 . ,x))))
(for/fold ([prev 0])
([i (in-list (pow-path n))])
(hash-set! pows i (* (hash-ref pows (- i prev)) (hash-ref pows prev)))
i)
(hash-ref pows n))
(define (show-pow x n)
(printf "~a: ~a\n" n (cdr (pow-path n)))
(printf "~a^~a = ~a\n" x n (pow-tree x n)))
(for ([x (in-range 18)])
(show-pow 2 x))
(show-pow 3 191)
(show-pow 1.1 81)
- Output:
0: () 2^0 = 1 1: (1) 2^1 = 2 2: (1 2) 2^2 = 4 3: (1 2 3) 2^3 = 8 4: (1 2 4) 2^4 = 16 5: (1 2 3 5) 2^5 = 32 6: (1 2 3 6) 2^6 = 64 7: (1 2 3 5 7) 2^7 = 128 8: (1 2 4 8) 2^8 = 256 9: (1 2 3 6 9) 2^9 = 512 10: (1 2 3 5 10) 2^10 = 1024 11: (1 2 3 5 10 11) 2^11 = 2048 12: (1 2 3 6 12) 2^12 = 4096 13: (1 2 3 5 10 13) 2^13 = 8192 14: (1 2 3 5 7 14) 2^14 = 16384 15: (1 2 3 5 10 15) 2^15 = 32768 16: (1 2 4 8 16) 2^16 = 65536 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.2402360440283
(formerly Perl 6)
Paths are random. It is possible replace .pick(*) with .reverse if you want paths as in Perl, or remove it for Python paths.
use v6;
sub power-path ($n ) {
state @unused_nodes = (2,);
state @power-tree = (False,0,1);
until @power-tree[$n].defined {
my $node = @unused_nodes.shift;
for $node X+ power-path($node).pick(*) {
next if @power-tree[$_].defined;
@unused_nodes.push($_);
@power-tree[$_]= $node;
}
}
( $n, { @power-tree[$_] } ...^ 0 ).reverse;
}
multi power ( $, 0 ) { 1 };
multi power ( $n, $exponent ) {
state %p;
my %r = %p{$n} // ( 0 => 1, 1 => $n ) ;
for power-path( $exponent ).rotor( 2 => -1 ) -> ( $p, $c ) {
%r{ $c } = %r{ $p } * %r{ $c - $p }
}
%p{$n} := %r ;
%r{ $exponent }
}
say 'Power paths: ', pairs map *.&power-path, ^18;
say '2 ** key = value: ', pairs map { 2.&power($_) }, ^18;
say 'Path for 191: ', power-path 191;
say '3 ** 191 = ', power 3, 191;
say 'Path for 81: ', power-path 81;
say '1.1 ** 81 = ', power 1.1, 81;
- Output:
Power paths: (0 => () 1 => (1) 2 => (1 2) 3 => (1 2 3) 4 => (1 2 4) 5 => (1 2 3 5) 6 => (1 2 3 6) 7 => (1 2 3 6 7) 8 => (1 2 4 8) 9 => (1 2 3 6 9) 10 => (1 2 3 5 10) 11 => (1 2 3 6 9 11) 12 => (1 2 3 6 12) 13 => (1 2 3 6 12 13) 14 => (1 2 3 6 12 14) 15 => (1 2 3 6 9 15) 16 => (1 2 4 8 16) 17 => (1 2 4 8 16 17)) 2 ** key = value: (0 => 1 1 => 2 2 => 4 3 => 8 4 => 16 5 => 32 6 => 64 7 => 128 8 => 256 9 => 512 10 => 1024 11 => 2048 12 => 4096 13 => 8192 14 => 16384 15 => 32768 16 => 65536 17 => 131072) Path for 191: (1 2 3 6 9 18 27 54 108 162 189 191) 3 ** 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347 Path for 81: (1 2 3 6 9 18 27 54 81) 1.1 ** 81 = 2253.24023604401
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.
/*REXX program produces & displays a power tree for P, and calculates & displays X^P.*/
numeric digits 1000 /*be able to handle some large numbers.*/
parse arg XP /*get sets: X, low power, high power.*/
if XP='' then XP='2 -4 17 3 191 191 1.1 81' /*Not specified? Then use the default.*/
/*────── X LP HP X LP HP X LP ◄── X, low power, high power ··· repeat*/
do until XP=''
parse var XP x pL pH XP; 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 an aligned display*/
do i=1 for words($); @.i= word($, i) /*build a fast Knuth's power tree array*/
end /*i*/
if p==0 then do; z= 1; call show; iterate; end /*handle case of zero power.*/
!.= .; z= x; !.1= z; prv= z /*define/construct 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 XP ···*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
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($) /*del extra blanks*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
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
- 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
var lvl = [[1]]
var p = Hash(1 => 0)
func path(n) is cached {
n || return []
while (n !~ p) {
var q = []
for x in lvl[0] {
for y in path(x) {
break if (x+y ~~ p)
y = x+y
p{y} = x
q << y
}
}
lvl[0] = q
}
path(p{n}) + [n]
}
func tree_pow(x, n) {
var r = Hash(0 => 1, 1 => x)
var p = 0
for i in path(n) {
r{i} = (r{i-p} * r{p})
p = i
}
r{n}
}
func show_pow(x, n) {
var fmt = ("%d: %s\n" + ["%g^%s = %f", "%s^%s = %s"][x.is_int] + "\n")
print(fmt % (n, path(n), x, n, tree_pow(x, n)))
}
for x in ^18 { show_pow(2, x) }
show_pow(1.1, 81)
show_pow(3, 191)
- Output:
0: [] 2^0 = 1 1: [1] 2^1 = 2 2: [1, 2] 2^2 = 4 3: [1, 2, 3] 2^3 = 8 4: [1, 2, 4] 2^4 = 16 5: [1, 2, 4, 5] 2^5 = 32 6: [1, 2, 4, 6] 2^6 = 64 7: [1, 2, 4, 6, 7] 2^7 = 128 8: [1, 2, 4, 8] 2^8 = 256 9: [1, 2, 4, 8, 9] 2^9 = 512 10: [1, 2, 4, 8, 10] 2^10 = 1024 11: [1, 2, 4, 8, 10, 11] 2^11 = 2048 12: [1, 2, 4, 8, 12] 2^12 = 4096 13: [1, 2, 4, 8, 12, 13] 2^13 = 8192 14: [1, 2, 4, 8, 12, 14] 2^14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2^15 = 32768 16: [1, 2, 4, 8, 16] 2^16 = 65536 17: [1, 2, 4, 8, 16, 17] 2^17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.1^81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
import "./big" for BigRat
import "./fmt" for Fmt
var p = { 1: 0 }
var lvl = [[1]]
var path // recursive
path = Fn.new { |n|
if (n == 0) return []
while (!p.containsKey(n)) {
var q = []
for (x in lvl[0]) {
System.write("") // guard against VM recursion bug
for (y in path.call(x)) {
if (p.containsKey(x + y)) break
p[x + y] = x
q.add(x + y)
}
}
lvl[0] = q
}
System.write("") // guard against VM recursion bug
var l = path.call(p[n])
l.add(n)
return l
}
var treePow = Fn.new { |x, n|
var r = { 0: BigRat.one, 1: BigRat.fromDecimal(x) }
var p = 0
for (i in path.call(n)) {
r[i] = r[i-p] * r[p]
p = i
}
return r[n]
}
var showPow = Fn.new { |x, n|
System.print("%(n): %(path.call(n))")
Fmt.print("$s ^ $d = $s\n", x, n, treePow.call(x, n).toDecimal(6))
}
for (n in 0..17) showPow.call(2, n)
showPow.call(1.1, 81)
showPow.call(3, 191)
- Output:
0: [] 2 ^ 0 = 1 1: [1] 2 ^ 1 = 2 2: [1, 2] 2 ^ 2 = 4 3: [1, 2, 3] 2 ^ 3 = 8 4: [1, 2, 4] 2 ^ 4 = 16 5: [1, 2, 4, 5] 2 ^ 5 = 32 6: [1, 2, 4, 6] 2 ^ 6 = 64 7: [1, 2, 4, 6, 7] 2 ^ 7 = 128 8: [1, 2, 4, 8] 2 ^ 8 = 256 9: [1, 2, 4, 8, 9] 2 ^ 9 = 512 10: [1, 2, 4, 8, 10] 2 ^ 10 = 1024 11: [1, 2, 4, 8, 10, 11] 2 ^ 11 = 2048 12: [1, 2, 4, 8, 12] 2 ^ 12 = 4096 13: [1, 2, 4, 8, 12, 13] 2 ^ 13 = 8192 14: [1, 2, 4, 8, 12, 14] 2 ^ 14 = 16384 15: [1, 2, 4, 8, 12, 14, 15] 2 ^ 15 = 32768 16: [1, 2, 4, 8, 16] 2 ^ 16 = 65536 17: [1, 2, 4, 8, 16, 17] 2 ^ 17 = 131072 81: [1, 2, 4, 8, 16, 32, 64, 80, 81] 1.1 ^ 81 = 2253.240236 191: [1, 2, 4, 8, 16, 32, 64, 128, 160, 176, 184, 188, 190, 191] 3 ^ 191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347
# remember the tree generation state and expand on demand
fcn path(n,p=Dictionary(1,0),lvl=List(List(1))){
if(n==0) return(T);
while(not p.holds(n)){
q:=List();
foreach x,y in (lvl[0],path(x,p,lvl)){
if(p.holds(x+y)) break; // not this y
y=x+y; p[y]=x;
q.append(y);
}
lvl[0]=q
}
path(p[n],p,lvl) + n
}
fcn tree_pow(x,n,path){
r,p:=Dictionary(0,1, 1,x), 0;
foreach i in (path){ r[i]=r[i-p]*r[p]; p=i; }
r[n]
}
fcn show_pow(x,n){
fmt:="%d: %s\n" + T("%g^%d = %f", "%d^%d = %d")[x==Int(x)] + "\n";
println(fmt.fmt(n,p:=path(n),x,n,tree_pow(x,n,p)))
}foreach x in (18){ show_pow(2,x) }
show_pow(1.1,81);
var [const] BN=Import("zklBigNum"); // GNU GMP big ints
show_pow(BN(3),191);- Output:
0: L() 2^0 = 1 1: L(1) 2^1 = 2 2: L(1,2) 2^2 = 4 3: L(1,2,3) 2^3 = 8 4: L(1,2,4) 2^4 = 16 5: L(1,2,4,5) 2^5 = 32 6: L(1,2,4,6) 2^6 = 64 7: L(1,2,4,6,7) 2^7 = 128 8: L(1,2,4,8) 2^8 = 256 9: L(1,2,4,8,9) 2^9 = 512 10: L(1,2,4,8,10) 2^10 = 1024 11: L(1,2,4,8,10,11) 2^11 = 2048 12: L(1,2,4,8,12) 2^12 = 4096 13: L(1,2,4,8,12,13) 2^13 = 8192 14: L(1,2,4,8,12,14) 2^14 = 16384 15: L(1,2,4,8,12,14,15) 2^15 = 32768 16: L(1,2,4,8,16) 2^16 = 65536 17: L(1,2,4,8,16,17) 2^17 = 131072 81: L(1,2,4,8,16,32,64,80,81) 1.1^81 = 2253.240236 191: L(1,2,4,8,16,32,64,128,160,176,184,188,190,191) 3^191 = 13494588674281093803728157396523884917402502294030101914066705367021922008906273586058258347