Powers from square roots
Powers from square roots
Created 11 April 2018 by Ronnie Mainieri.
Last updated 15 April 2018
Last updated 15 April 2018
This notebook has a method for computing numerically using the four operations and the square root by expanding into a binary float. There are two algorithms: a simpler one that handles only the case of and another that handles (by re-using functions names, so use the Clear).
a
c
c
0≤c<1
c>0
The algorithms use the binary expansion of the exponent 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 and that are applied in sequence to 1. For example, to compute , one sets , expands 6 into binary notation, 110, and then computes 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.
c
h(x)=ax
1
2
h(x)=x
0
2
3
6
a:=3
h◦h◦h(1)
0
1
1
For the fractional part of the exponent the algorithm is similar. One uses the functions and . To compute , one expands into a binary float and then applies and in sequence to 1. One starts with the function corresponding to the rightmost index and proceeds to the left. For example, to compute , one expands 23/64 into its binary float and computes to obtain 1.7817.
g(x)=
0
x
g(x)=
1
ax
a
c
c
g
0
g
1
5
23/64
.010111
g◦g◦g◦g◦g◦g(1)
0
1
0
1
1
1
For exponents smaller than 1
For exponents smaller than 1
Given a positive number smaller than one, computes its floating expansion up to a certain number of binary digits.
x
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 or multiply and square root operations right to left on the input list.
h
0
h
1
In[]:=
iter[a_][list_]:=Module[{x=1},N@FoldList[Sqrt@If[#21,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
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[#21,a#1,#1]&,1,i]*FoldListIf[#21,a#1,#1]&,1,Reverse@f〚-1〛
2
2
1/2
//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