Evaluate

Powers from square roots

Created 11 April 2018 by Ronnie Mainieri.
Last updated 15 April 2018

This notebook has a method for computing
a
c
numerically using the four operations and the square root by expanding
c
into a binary float. There are two algorithms: a simpler one that handles only the case of
0≤c<1
and another that handles
c>0
(by re-using functions names, so use the Clear).

The algorithms use the binary expansion of the exponent
c
to apply in sequence one of two functions to an initial value of 1. The integer part is known as the exponentiation by squaring algorithm. It uses two functions
h
1
(x)=ax
2
and
h
0
(x)=x
2
that are applied in sequence to 1. For example, to compute
3
6
, one sets
a:=3
, expands 6 into binary notation, 110, and then computes
h
0
◦h
1
◦h
1
(1)
to get 729. Notice that the functions are applied based on their index starting from the left of the binary number and proceeding to right.

For the fractional part of the exponent the algorithm is similar. One uses the functions
g
0
(x)=
x
and
g
1
(x)=
ax
. To compute
a
c
, one expands
c
into a binary float and then applies
g
0
and
g
1
in sequence to 1. One starts with the function corresponding to the rightmost index and proceeds to the left. For example, to compute
5
23/64
, one expands 23/64 into its binary float
.010111
and computes
g
0
◦g
1
◦g
0
◦g
1
◦g
1
◦g
1
(1)
to obtain 1.7817.

For exponents smaller than 1

Given a positive number
x
smaller than one, computes its floating expansion up to a certain number of binary digits.

In[]:=
binaryDigits[x_,opts___Rule]:=Module[​
​{a,shift,len},​
​​
​len="Digits"/.{opts}/."Digits"20;​
​--len;​
​​
​{a,shift}=RealDigits[x,2,len];​
​ConstantArray[0,-shift]~Join~a​
​]/;0≤x<1

The function accepts an option specifying how many binary digits to return

In[]:=
binaryDigits[1/3,"Digits"10]
Out[]=
{0,1,0,1,0,1,0,1,0,1}

The iter function implements the repeated square root algorithm to compute a power of a. It returns the values of applying the square root
h
0
or multiply and square root
h
1
operations right to left on the input list.

In[]:=
iter[a_][list_]:=Module[{x=1},N@FoldList[Sqrt@If[#21,a#1,#1]&,1,Reverse@list]]

To compute the cube root of 8, call

In[]:=
iter[8]@binaryDigits[1/3]
Out[]=
{1.,2.82843,1.68179,3.66802,1.91521,3.91429,1.97846,3.9784,1.99459,3.99459,1.99865,3.99865,1.99966,
3.99966,1.99992,3.99992,1.99998,3.99998,1.99999,3.99999,2.}

And we can check that it handles the edge case

In[]:=
iter[8]@{0,0,0,0,0,0}
Out[]=
{1.,1.,1.,1.,1.,1.,1.}

For all positive exponents

Re-write of the binaryDigits and iter functions for positive powers

In[]:=
Clear[binaryDigits,iter];

The function now returns two lists: one for the integer part and another for the fractional part. The option “Digits” is an integer controlling the total number of digits returned among both lists.

In[]:=
binaryDigits[x_,opts___Rule]:=Module[​
​{a,shift,len},​
​​
​len="Digits"/.{opts}/."Digits"20;​
​​
​{a,shift}=RealDigits[x,2,len];​
​If[shift≤0,​
​{{0},ConstantArray[0,-shift]~Join~a}​
​(*else*),​
​{a〚1;;shift〛,a〚(shift+1);;〛}​
​]​
​]/;x>0
In[]:=
binaryDigits[10.76514]
Out[]=
{{1,0,1,0},{1,1,0,0,0,0,1,1,1,1,1,0,0,0,0,0}}

The iter function now accepts a double list and applies the square root algorithm to the fractional part and a powers of 2 algorithm to the integer part.

In[]:=
iter[a_][{i_,f_}]:=Module​
​{x=1},​
​​
​FoldList[If[#21,a#1
2
,#1
2
]&,1,i]*FoldListIf[#21,a#1,#1]
1/2
&,1,Reverse@f〚-1〛​
​//N​
​
In[]:=
binaryDigits[]
Out[]=
{{1,0},{1,0,1,1,0,1,1,1,1,1,1,0,0,0,0,1,0,1}}
In[]:=
8^//N
Out[]=
285.005