Algorithm for Fractional power calculation

Hi

1. fractional base is no problem with current implementation.

2. To support fractional exponents, get the n-th root for any given number b. How to implement algorithm to get a numerical approximation ?

3. Current approach is inefficient, because it loops e times

b = [-32:32]; % example input values
e = [-3:3];   % example input values but doesn't support fraction's
power_function(b,e)
p = 1;
if e < 0
    e = abs(e);
    multiplier = 1/b;
else
    multiplier = b;
end

for k = 1:e
    p(:) = p * multiplier; %  n-th root for any given number 
end

Answer 1

[Alan Stevens]

How about using the Newton-Raphson algorithm. Here's the basic idea:

% x^n = b
% Let f(x) = x^n - b
% dfdx(x) = n*x^(n-1)
%
% Use Newton-Raphson iteration
% x1 = initial guess
% err = 1;
% while err > tol
%       x = x1 - f(x1)/dfdx(x1);
%       err = abs(x - x1);
%       x1 = x
% end

[Life is Wonderful]

Can you please elaborate a bit more. Input the values with base, exponent and return.

I see Newton's method as nthroot()

x = -5

n = -5

I see Newton's method as below

y = (sign(x) + (x==0)) .* (abs(x).^(1./n));
m = x ~= 0 & (abs(x) < (1/eps(class(y)))) & isfinite(n);
if isscalar(n)
    y(m) = y(m) - (y(m).^n - x(m)) ./ (n .* y(m).^(n-1));
else
    y(m) = y(m) - (y(m).^n(m) - x(m)) ./ (n(m) .* y(m).^(n(m)-1));
end

Thank you!!

[Alan Stevens]

Here's a simple example

% x = b^(1/n)   The n'th root of b
% This can be written as x^n = b  or x^n - b = 0
%
% Let f(x) = x^n - b  and dfdx(x) = n*x^(n-1)
%
% Then Newton-Raphson algorithm is
% x(n) = x(n-1) - f(x(n-1))/dfdx(x(n-1))

% Simple example:
b = 5;
n = 3;
% i.e. we want the cube root of 5

f = @(x) x^n - b;
dfdx = @(x) n*x^(n-1);
tol = 10^-4;
err = 1;

x = 1;  % initial guess
while err > tol
    xn = x - f(x)/dfdx(x);
    err = abs(xn - x);
    x = xn;
end
disp(x)

    1.7100

% Check
disp([x b^(1/n)])

    1.7100    1.7100

Your original values of b include negative values, for which the roots will, in general, involve complex roots. You will need to modify the routine to account for the real and imaginary parts for these cases.

[Life is Wonderful]

It seems have issues with embedded.fi.

% x = b^(1/n) The n'th root of b

% This can be written as x^n = b or x^n - b = 0

%

% Let f(x) = x^n - b and dfdx(x) = n*x^(n-1)

%

% Then Newton-Raphson algorithm is

% x(n) = x(n-1) - f(x(n-1))/dfdx(x(n-1))

% Simple example:

b = fi([5],0,32, 28);
n = fi([3],0,32, 29);
% i.e. we want the cube root of 5

f = @(x) x^n - b;
dfdx = @(x) n*x^(n-1);
tol = 10^-4;
err = 1;

x = 1;  % initial guess
while err > tol
    xn = x - f(x)/dfdx(x);
    err = abs(xn - x);
    x = xn;
end

Error using ^ (line 31)
The first input argument must be a fi object.

Error in solution (line 5)
f = @(x) x^n - b;

disp(x)

% Check
disp([x b^(1/n)])

[Life is Wonderful]

Your original values of b include negative values, for which the roots will, in general, involve complex roots. You will need to modify the routine to account for the real and imaginary parts for these cases.

Yes, my initial implementation works fine for negative base and exponent, At the moment, I am looking for the implementation which supports 32 bit size DSP. which means p = power(base ,exp)

[Walter Roberson]

For fixed point work with non-integer exponents, you will find that it is typically easiest to take the fixed-point log, multiply by the exponent, and then take the fixed-point exp() .

There are particular integer exponents and particular roots where those particular powers can be expressed in better ways, but CORDIC approaches are the way to go in most other cases.

[Life is Wonderful]

You mean,

function = power(base, exponent)
p = exp(exponent * log(base))
end

@Walter Roberson, Thanks, Do you have an exp() code for fixed point implementation ?

Pls, I don't need lookup table based exp() and cordicexp().

Thank you!!

[Walter Roberson]

https://ieeexplore.ieee.org/document/8738918

[Life is Wonderful]

I don't need cordic , paper talks about generilzed cordic algorithm. Any other suggestion for exponential function would be help.

Thank you!!

[Walter Roberson]

You asked for Fixed Point implementations, but you say that you do not need CORDIC, and you also say that you do not need lookup tables.

When you say that you do not need lookup tables or CORDIC, do you mean that you only want a loose approximation of the value? Do you mean that you have a lot of memory restrictions so you cannot afford to store the range lookup tables? Do you mean that you are dealing with patent / trade secret issues and cannot use CORDIC techniques for legal reasons?

It would help if I were to understand why you do not want to use well-established CORDIC methods ???

https://stackoverflow.com/questions/4657468/fast-fixed-point-pow-log-exp-and-sqrt

[Life is Wonderful]

[WR] :Do you mean that you have a lot of memory restrictions so you cannot afford to store the range lookup tables?

Yes,

[WR] do you mean that you only want a loose approximation of the value?

Ofcourse not, but expectation would be not less than -90dBc ( THDN) ( which most of the standard accepts)

pow(sqrt(x), 2*frac(y)

is a good point . I'II try out .

Thanks a lot for helping me

[Walter Roberson]

https://cradpdf.drdc-rddc.gc.ca/PDFS/unc82/p531366.pdf talks about a reduced-complexity lookup table method.

[Life is Wonderful]

Hi Walter,

Please see below code which is computing exp() and taylor series. At the 1 iteration final exp() value is calculated/converge. Please find the attachment

I could see with one multiplication and one addition, I could get the exp() THDN = -108dBc with x = +23 & -145 with x = -23 dBc which is very much efficient. The same works for fixed point as well . One need to be careful with division operator ( this takes lots time ).

I have added Excel spread sheet ( to replicate my below implementation ) for offline calculation. It has upto 10 fraction place of accuracy ( which pretty good for any standard as acceptance criteria). Please find the attachment.

Next step - With below code, do you have a proposal where number iteration is reduced [ I think one can optimize 90 iteration ] and If one do a loop unroll [smaller than 5 iteration ] , it can be done one max 5 steps. Any proposal/ suggestion ?

Thanks !!

n = 1e2;
x = 23.0; % for  taylor series [ -22:22] for 32 bit DSP
inpval = 23.0;  % for built-in exp()
% create a function and pass the args
expval = 1.0; % initialize as double
fprintf('%20s|%30s|%26s|%25s|%25s|\n--------------------+------------------------------+--------------------------+-------------------------+-------------------------+\n',"indx","expval","exp(x)",'diff','Nbits');
for i = n-1:-1:1
    expval = 1 + x * expval / i;
    fprintf('%20d|%30.10f|%26.15f|%25.10f|%25.10f|\n',i,expval,exp(inpval),abs(exp(inpval) - expval), log2(exp(inpval)));
end

[Walter Roberson]

Sorry, I do not know.