How to add function argument for nthroot() ?

Hi

I want to calculate nthroot() for 'embedded.fi', following is the code.

I don't understand , how to fix rem(n(:),2) ? / reminder after division

There is claimed fix https://in.mathworks.com/matlabcentral/answers/100386-why-do-i-receive-an-error-when-i-use-the-nthroot-function-within-matlab#answer_109734 from MathWorks Staff

x = fi(1,0,32,31);
n = fi(1,0,32,31);
if isreal(x) || isreal(n)
    nroot = nthroot(x, n);
end

Check for incorrect argument data type or missing argument in call to function 'rem'.

Error in nthroot (line 27)
if any((x(:) < 0) & (n(:) ~= fix(n(:)) | rem(n(:),2)==0))

rem(1,1) is 0 

Answer 1

[Walter Roberson]

That fix you link to has no relationship to your current difficulty.

https://www.mathworks.com/help/matlab/ref/nthroot.html#buchym2-X

Data Types: single | double

Notice that fi is not listed. The nthroot() function is not defined for fixed point objects.

x = fi(1,0,32,31);
n = fi(1,0,32,31);
if isreal(x) && isreal(n)
   nroot = sign(x) .* abs(x).^(1./n);
else
   nroot = x .^ (1./n);
end
nroot

nroot = 
     1

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 40
        FractionLength: 31

[Life is Wonderful]

@Walter Roberson, Thanks, Wondering what makes wordlength and fractional length goes for a toss.

I think , it is incorrect WordLength: 40 FractionLength: 31

please see this below numeric calculation. Can please help me to understand what kind of cast is needed to make your proposal work i.e. without fi object ?

Thank you!!

nthroot(1024,2)

ans = 32

x = fi(1,0,32,31);
n = fi(1,0,32,31);
if isreal(x) && isreal(n)
    nroot = sign(x) .* abs(x).^(1./n);
else
    nroot = x .^ (1./n);
end
nroot

nroot = 
     1

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 40
        FractionLength: 31

x = fi(1024,0,32,21);
n = fi(2,0,32,2);

srt_val = fi(sqrt(x),0,32,26)

srt_val = 
    32

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Unsigned
            WordLength: 32
        FractionLength: 26

nrt_val = fi(nthroot(1024,2),0,32,26)

nrt_val = 
    32

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Unsigned
            WordLength: 32
        FractionLength: 26

x = 1024,

x = 1024

n = 2;

if isreal(x) && isreal(n)
    nroot = sign(x) .* abs(x).^(1./n);
else
    nroot = x .^ (1./n);
end
nroot

nroot = 32

n = fi(2,0,32,31);
x = fi(1024,0,32,22);
if isreal(x) && isreal(n)
    nroot = sign(x) .* abs(x).^(1./n);
else
    nroot = x .^ (1./n);
end

Error using fixed.internal.power.validateK (line 13)
Exponent input to 'power' must be a real scalar and the value must be a non-negative integer.

Error in fixed.internal.power.powImpl>validateInputs (line 39)
fixed.internal.power.validateK(k);

Error in fixed.internal.power.powImpl (line 9)
validateInputs(a, k);

Error in .^ (line 32)
y = fixed.internal.power.powImpl(a, k);

nroot

[Walter Roberson]

I see what you mean.

Is the root to be taken fixed (like always sqrt()), or is it dynamic? If dynamic, is it restricted to one out of a small list?

[Walter Roberson]

If the root to be taken is not fixed, and is not one out of a small list, then:

If you know your maximum range of x values, and your maximum root, you can use a taylor expansion on x^(1/n) to get any desired (fixed) degree of accuracy, in the form of a polynomial. For example,

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10)

T = 

Tx = collect(simplify(T), x)

Tx = 

misfit = subs(Tx, n, 10) - x^(1/sym(10))

misfit = 

fplot(misfit, [0 1024])

vpa(subs(misfit,x,1))

ans = 

0.3975946820153549595981366392819

You can see that degree 10 expansion with n = 10 is not all that accurate for x = 1 (and is worse at 0) so you might need to use a higher order polynomial in the taylor expansion; also a higher order would probably be needed if your maximum n is above 10, or if your range of values is more than 1000.

The main benefit of this approach is that it fixes the amount of work to be done, which could be important for VHDL purposes.

You could probably get a more accurate result faster using some kind of loop.

[Life is Wonderful]

Sorry, I didn't understand fixed vs dynamic!!

I consider root as square ,cube and it goes on from Max to Min number . Return root needs fixed word length otherwise we have truncated value

WL = 32;
fprintf('%10s|%30s|\n----------+------------------------------+\n','power', 'root');

     power|                          root|
----------+------------------------------+

for i = 1:WL % Num bits in Wordlength
        fprintf('%10d|%30.10f|\n',i, nthroot(2^WL,i))
end

         1|         4294967296.0000000000|
         2|              65536.0000000000|
         3|               1625.4986772154|
         4|                256.0000000000|
         5|                 84.4485062895|
         6|                 40.3174735966|
         7|                 23.7759086262|
         8|                 16.0000000000|
         9|                 11.7578759382|
        10|                  9.1895868400|
        11|                  7.5114472853|
        12|                  6.3496042079|
        13|                  5.5080378048|
        14|                  4.8760546168|
        15|                  4.3872999188|
        16|                  4.0000000000|
        17|                  3.6867585638|
        18|                  3.4289759314|
        19|                  3.2136449063|
        20|                  3.0314331330|
        21|                  2.8754933949|
        22|                  2.7407019694|
        23|                  2.6231573457|
        24|                  2.5198420998|
        25|                  2.4283897688|
        26|                  2.3469209200|
        27|                  2.2739258102|
        28|                  2.2081790273|
        29|                  2.1486764788|
        30|                  2.0945882456|
        31|                  2.0452228712|
        32|                  2.0000000000|

Another question is , for negative base value

example nthroot(-1 * (2^32),i) - what is the formula

[Walter Roberson]

I think at this point, your best approach would be to use a cordic algorithm to take logarithm, divide by n, and use a cordic exp().

Remember that nthroot() of a negative number is only defined for odd integer n, for which it is -1*nthroot(-x,n).

nthroot is not just another way to write .^ and is instead specific to finding a real root... which is not possible for complex x or for negative x with an even root.

[Life is Wonderful]

Thank @Walter Roberson. At the moment, I tested small dataset and it looks promising . Thanks a lot for your efforts.

[Life is Wonderful]

@Walter Roberson, could you please share example snippet , Error message

Function input at args{1} does not have a valid type.

Caused by:

Class sym is not supported by coder.Type as it is a handle class.

Use help codegen for more information on using this command.

Error using codegen

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10)

how can I use matlab coder [command line ] to generate taylor series here ?

coder.extrinsic('sym');
coder.extrinsic('taylor');
x = sym('x');
n = sym('n');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10)

T = 

codegen  taylor -args {x(1), n(1)} ...
    -o taylor_fixedpoint_mex -config:mex

Unable to run the 'fevalJSON' function because it calls the 'codegen' function, which is not supported for this product offering.

[Walter Roberson]

MATLAB Coder can never be used with Symbolic Toolbox.

What you would do is

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);

Tfun = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor.m', 'optimize', false);

Tfun_opt = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

You can use code generation with nthrootTaylor.m or nthrootTaylor_opt.m

You do not need both of the .m files. I show generating both of them because it is important that you test the optimized version for accuracy compared to the non-optimized version. There have been serious bugs in the optimization process in some of the recent releases.

Reminder: the order 10 approximation is not a good one as you get further away from the expansion point (here 512). A cordic log() and exp() would probably be significantly more accurate.

[Life is Wonderful]

Thanks Walter,

code generation with nthrootTaylor.m or nthrootTaylor_opt.m

Unfortunately , I didn't get the steps - how to use codegen for

Tfun = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor.m', 'optimize', false);

&

Tfun_opt = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

One more thing, can you please insight on cordic log() and exp().

[Walter Roberson]

codegen nthrootTaylor.m -args {x(1), n(1)} ...
    -o nthrootTaylor_fixedpoint_mex -config:mex
codegen nthrootTaylor_opt.m -args {x(1), n(1)} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex

Remember, in practice you will only need one of these two. The reason I show both of them is that in recent versions 'optimize', true for matlabFunction() has produced incorrect code, so you need to test the two versions against each other before you can trust the version with 'optimize', true set on. The version with 'optimize', true should be notably more efficient if it works properly .

See https://www.mathworks.com/help/fixedpoint/ref/cordiccexp.html for cordic exp()

The page at https://www.quinapalus.com/efunc.html shows the general algorithm that could be used for log(), and see also https://pdfs.semanticscholar.org/fde2/099428b378dbffa2cb6445a070aafbd6c915.pdf

Ah, there is a more direct route:

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

[Life is Wonderful]

Thanks a lot for your efforts !! But I am unable to make progress here ....

Can you please help in correcting the mistake

x = fi(linspace(0,31.99999,10),0,32,27);    % x can be negative as well, but start with positive numbers
n = fi(linspace(0,3,10),0,32,30);           % n can be negative as well, but start with positive numbers
[T] = taylor_fun(x,n);
Tfun_opt = nthrootTaylor_opt(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

function [T] = taylor_fun(x,n)
syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);

end

function [Tfun_opt] = nthrootTaylor_opt(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true)

Invalid expression. Check for missing multiplication operator, missing or unbalanced delimiters, or other syntax error. To construct matrices, use brackets instead of parentheses.

end


codegen nthrootTaylor_opt.m -args {x(1), n(1)} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex

[Walter Roberson]

What? No!

Preparation, not to be deployed

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);
Tfun = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor.m', 'optimize', false);
Tfun_opt = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

x_example = fi(0,0,32,27);    % x can be negative as well, but start with positive numbers
n_example = fi(0,0,32,30);           % n can be negative as well, but start with positive numbers

codegen nthrootTaylor.m -args {x_example, n_example} ...
    -o nthrootTaylor_fixedpoint_mex -config:mex
codegen nthrootTaylor_opt.m -args {x_example, n_example} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex

Once code is generated then for deployment

x = fi(linspace(0,31.99999,10),0,32,27);    % x can be negative as well, but start with positive numbers
n = fi(linspace(0,3,10),0,32,30);           % n can be negative as well, but start with positive numbers

num_x = length(x);
num_n = length(n);

results = fi(zeros(num_x, num_n),0,32,30);
results_opt = fi(zeros(num_x, num_n),0,32,30);

for nidx = 1 : num_n
    for xidx = 1 : num_x
        results(xidx,nidx) = nthrootTaylor(x(xidx), n(nidx) );
        results_opt(xidx,nidx) = nthrootTaylor_opt(x(xidx), n(nidx) );
    end
end

now compare results to results_opt to verify whether the optimized version produces results that are acceptable quality compared to the non-optimized version.

At the moment, I do not know whether the code that would be generated would have looped over vector n and vector x -- and I didn't want to have to look up how to signal to codegen to expect arrays.

At the moment i do not know what fi would be used for the results, so I had to fill in something.

[Life is Wonderful]

Still no success in codegen !! Also I have think fi object is not working

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);
Tfun = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor.m', 'optimize', false);

Error using sym/matlabFunction>checkVarsSubset (line 253)
Free variable 'n' must be included in 'Vars' value.

Error in sym/matlabFunction>checkVars (line 240)
checkVarsSubset(vexpanded,funvars);

Error in sym/matlabFunction (line 158)
vars = checkVars(funvars,opts);

Tfun_opt = matlabFunction(T, 'vars', x, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

x_example = fi(0,0,32,27);    % x can be negative as well, but start with positive numbers
n_example = fi(0,0,32,30);           % n can be negative as well, but start with positive numbers

codegen nthrootTaylor.m -args {x_example, n_example} ...
    -o nthrootTaylor_fixedpoint_mex -config:mex
codegen nthrootTaylor_opt.m -args {x_example, n_example} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex

x = fi(linspace(0,31.99999,10),0,32,27);    % x can be negative as well, but start with positive numbers
n = fi(linspace(0,3,10),0,32,30);           % n can be negative as well, but start with positive numbers
num_x = length(x);
num_n = length(n);
results = fi(zeros(num_x, num_n),0,32,30);
results_opt = fi(zeros(num_x, num_n),0,32,30);
for nidx = 1 : num_n
    for xidx = 1 : num_x
        results(xidx,nidx) = nthrootTaylor(x(xidx), n(nidx) );
        results_opt(xidx,nidx) = nthrootTaylor_opt(x(xidx), n(nidx) );
    end
end

I tried with and without fi object verification . I think there is an issue

Output argument "Tfun" (and possibly others) not assigned a value in the execution with "pow_taylor_v1>nthrootTaylor" function.

Error in pow_taylor_v1 (line 14)

results(xidx,nidx) = nthrootTaylor(x(xidx), n(nidx) );

clearvars;clc;
x = fi(linspace(0,31.99999,10),0,32,27);    % x can be negative as well, but start with positive numbers
n = fi(linspace(0,3,10),0,32,30);           % n can be negative as well, but start with positive numbers
% x = linspace(0,31.99999,10)
% n = linspace(0,3,10)
num_x = length(x);
num_n = length(n);

results = fi(zeros(num_x, num_n),0,32,30);
results_opt = fi(zeros(num_x, num_n),0,32,30);

for nidx = 1 : num_n
    for xidx = 1 : num_x
        results(xidx,nidx) = nthrootTaylor(x(xidx), n(nidx) );
%         results(xidx,nidx) = nthroot(x(xidx), n(nidx) )
%         results_opt(xidx,nidx) = nthrootTaylor_opt(x(xidx), n(nidx) );
    end
end

function [Tfun] = nthrootTaylor(x,n)
% syms x real
% syms n integer
% assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);
end

[Walter Roberson]

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 10);
Tfun = matlabFunction(T, 'vars', {x, n}, 'File', 'nthrootTaylor.m', 'optimize', false);
Tfun_opt = matlabFunction(T, 'vars', {x, n}, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

x_example = fi(0,0,32,27);    % x can be negative as well, but start with positive numbers
n_example = fi(0,0,32,30);           % n can be negative as well, but start with positive numbers

codegen nthrootTaylor.m -args {x_example, n_example} ...
    -o nthrootTaylor_fixedpoint_mex -config:mex
codegen nthrootTaylor_opt.m -args {x_example, n_example} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex

[Life is Wonderful]

still there is an issue!! codegen is unsuccessful

Please ignore the editor error since it is unsupported

syms x real
syms n integer
assumeAlso(n, 'positive');
T = taylor(x^(1/n), x, 'ExpansionPoint', 512, 'Order', 12);
Tfun = matlabFunction(T, 'vars', {x, n}, 'File', 'nthrootTaylor.m', 'optimize', false);
Tfun_opt = matlabFunction(T, 'vars', {x, n}, 'File', 'nthrootTaylor_opt.m', 'optimize', true);

x_example = fi(0,0,32,27);    % x can be negative as well, but start with positive numbers
n_example = fi(0,0,32,30);           % n can be negative as well, but start with positive numbers

codegen nthrootTaylor.m -args {x_example, n_example} ...
    -o nthrootTaylor_fixedpoint_mex -config:mex
??? The computed word length of the result is 142 bits. This exceeds the maximum supported wordlength of 128 bits.

Unable to run the 'fevalJSON' function because it calls the 'codegen' function, which is not supported for this product offering.

codegen nthrootTaylor_opt.m -args {x_example, n_example} ...
    -o nthrootTaylor_opt_fixedpoint_mex -config:mex
??? Function 'exp' is not defined for values of class 'embedded.fi'.

[Walter Roberson]

This does not astonish me. When you have a series of expressions that are not being stored into a variable with defined fixed point characteristics, then the Fixed Point Toolbox dynamically expands the word length in order to try to preserve precision.

I suspect it might be possible to attach behaviour to the fi() so that the output results of expressions does not automatically expand for precision, but that is not something I have researched.

One approach would be to break up the expression into smaller segments that you assign into variables that you have set the precision for.

For example if you were to take children(T) then you would get a cell array of terms that are to be added together. You could then vertcat() the cell expansion to turn it into a vector, and you could compute the vector in the generated code, and fi() the results and sum() them.

You would still have some difficulty -- the individual term for the derivative to the 9th power would want to blow up the precision.

Probably it would be more rigourous to calculate the terms in a loop, applying fi() at each step.

But... really I think this approach is a dead end. Look at the precision loss near 0 or 1 !!!