Hello!
There is a recip16 (void recip16 (DATA *x, DATA *r, DATA *rexp, ushort nx)) function.
I'd like to implement it in C.
I've taken DSP_recip16 function (from C64x DSPLib) algorithm as a basement.
My function:
void DSP_recip16_c(DATA *x, DATA *rfrac, DATA *rexp, ushort nx)
{
LDATA a, b;
DATA neg = 0,
normal,
i,
j;
for(i = nx; i > 0; i--)
{
a = *(x++);
if(a < 0) /* take absolute value */
{
a = -a;
neg = 1;
}
normal = _norm(a); /* normalize */
a = a << 16 + (1<<normal);
*(rexp++) = 2 << normal;
b = 0x80000000;
for(j = 15; j > 0; j--) {
b = subc(b, a);
}
b = b & 0x7FFF;
if(neg) b = -b; /* if negative, negate */
*(rfrac++) = b; /* store fraction */
}
}
I've implemented subc function in this way:
static inline LDATA subc(LDATA src1, LDATA src2) {
if ((src1 - src2) >= 0)
return ((src1 - src2) << 1) + 1;
else
return src1 << 1;
}
The exponent value is right, but the fraction is not.
Where is the problem?
Thank you!