Integer Division Q15 Result

[Sebastian]

Hi,

I am working on the software for an eCompass and I am using an appnote from Freescale Semiconductor which includes c# source code. (http://cache.freescale.com/files/sensors/doc/app_note/AN4248.pdf) I am working with Microchip xc8 which is c so I had to change some things and it does not work. I am not even sure if the original code is correct because I tried to follow it with a calculator and it did not work either. The part I am working on at the moment is some integer division thing. The whole software uses 16bit integers only which is nice because the µC I have to use is a pretty ancient (PIC16F690 8bit) so I can't just use floating point.
The problem is I am not exactly a math genius and I am also not that good with c and i never heard about or even used c# before.
That's what is in the appnote:

Code: [Select]

const UInt16 MINDELTADIV = 1; /* final step size for iDivide */
/* function to calculate ir = iy / ix with iy <= ix, and ix, iy both > 0 */
static Int16 iDivide(Int16 iy, Int16 ix)
{
Int16 itmp; /* scratch */
Int16 ir; /* result = iy / ix range 0., 1. returned in range 0 to 32767 */
Int16 idelta; /* delta on candidate result dividing each stage by factor of 2 */
/* set result r to zero and binary search step to 16384 = 0.5 */
ir = 0;
idelta = 16384; /* set as 2^14 = 0.5 */
/* to reduce quantization effects, boost ix and iy to the maximum signed 16 bit value */
while ((ix < 16384) && (iy < 16384))
{
ix = (Int16)(ix + ix);
iy = (Int16)(iy + iy);
}
/* loop over binary sub-division algorithm solving for ir*ix = iy */
do
{
/* generate new candidate solution for ir and test if we are too high or too low */
itmp = (Int16)(ir + idelta); /* itmp=ir+delta, the candidate solution */
itmp = (Int16)((itmp * ix) >> 15);
if (itmp <= iy) ir += idelta;
idelta = (Int16)(idelta >> 1); /* divide by 2 using right shift one bit */
} while (idelta >= MINDELTADIV); /* last loop is performed for idelta=MINDELTADIV */
return (ir);
}
That's my code:

Code: [Select]

int16_t Divide(int16_t y, int16_t x) {

    int16_t Tmp, Result, Delta;

    Result = 0;
    Delta = 16384;

    while ((x < 16384) && (y < 16384)) {
        x = (x + x);
        y = (y + y);
    }

    do {
        Tmp = (Result + Delta);
        Tmp = ((Tmp * x) >> 15);
        if (Tmp <= y) {
            Result += Delta;
        }
        Delta = (Delta >> 1);
    } while (Delta >= 1);


    return Result;
}
Does anybody see any mistakes or is there an easier way of doing that?

[dannyf]

I can confirm for you that the Freescale code you posted is wrong - not sure if it is wrong originally or you transcribed it wrong. But as posted, it doesn't work. It is also very slow: I timed it to be 10k - 20k cycles.

A few options:

1) Find a Q15 library: lots of them over the web so it shouldn't be difficult.
2) Find a AD9850 library: it should have a routine converting the desired frequency to tuning word. With minimum modification, it works for you here - I timed it to be 500 - 600 cycles.
3) Port the ARM DSP library here: or just the Q15 portion of it. Daunting but doable.

[dannyf]

Code: [Select]

itmp = (Int16)((itmp * ix) >> 15);
That's why the Freescale code doesn't work on XC8. Once fixed, it works just fine.

[krivx]

Code: [Select]

itmp = (Int16)((itmp * ix) >> 15);
That's why the Freescale code doesn't work on XC8. Once fixed, it works just fine.

Yep, you need to be familiar with the integer promotion rules of your compiler.

[Sebastian]

Do you know what I need to change in order to make it work?

[dannyf]

Trust me, you will learn a ton more if you figure out the problem. Enough hints have been provided here.

[westfw]

If your compiler has 32bit math that you trust (and/or are already using), try

Code: [Select]

uint16_t idiv(uint16_t iy, uint16_t ix)
{
    uint32_t t;

    t = iy*32768;
    t += ix/2;
    return (uint16_t)(t/ix);
}

[Sebastian]

If your compiler has 32bit math that you trust (and/or are already using), try

Code: [Select]

uint16_t idiv(uint16_t iy, uint16_t ix)
{
    uint32_t t;

    t = iy*32768;
    t += ix/2;
    return (uint16_t)(t/ix);
}

That works awesome, thank you!

[dannyf]

That works awesome,

It also negates the purpose of a fixed point math format, like Q15.

[Sebastian]

I just want to get that thing running. If I had the choice I would have used a better µC which could do it in hardware, is cheaper and would have made the hardware smaller and cheaper because I wouldn't have to use a 5V to 3V3 level shifter which would make other stuff easier as well. But my Teacher insisted in using that one and now I have to program stuff I never learned. I will try to get the rest of the software working on my own but I will probably have to ask you again.

[westfw]

It also negates the purpose of a fixed point math format, like Q15.

Huh?  It still produces a Q15 fixed point number (which the trig functions from the app note go on to use, I believe.)  Compared to the code from the app note, it adds exactly one 32bit division, while removing multiple (~13?) 32bit multiplications.

The reliance of the original code on (hidden) 32bit multiplication is annoying to the point where I would call it broken.  If you're going to implement a bit-wise Q15 division algorithm for "improved efficiency", you shouldn't be relying on 32bit multiplication for it to work.  And especially not on 32bit multiplications resulting from implicitly-promoted shorter datatypes.
(And "example portable code" written in C#?  WTF?)

[Sebastian]

I am now trying to get this HundredAtan function to work. It should give a value between 0° to 90° *100, so 0 to 9000. There is another function called HundredAtan2 that inverts and adds 18000 or something like that to put in in the right quadrant. But I understand that. The problem is the HundredAtan function and I think it is the same problem as before with the shifts or something like that. I tried to make my own function using the integrated atan2 function of the xc8 compiler like westfw did it with the divide function but it did not work.

That is the original c# code:

Code: [Select]

/* fifth order of polynomial approximation giving 0.05 deg max error */
const Int16 K1 = 5701;
const Int16 K2 = -1645;
const Int16 K3 = 446;
/* calculates 100*atan(iy/ix) range 0 to 9000 for all ix, iy positive in range 0 to 32767 */
static Int16 iHundredAtanDeg(Int16 iy, Int16 ix)
{
Int32 iAngle; /* angle in degrees times 100 */
Int16 iRatio; /* ratio of iy / ix or vice versa */
Int32 iTmp; /* temporary variable */
/* check for pathological cases */
if ((ix == 0) && (iy == 0)) return (0);
if ((ix == 0) && (iy != 0)) return (9000);
/* check for non-pathological cases */
if (iy <= ix)
iRatio = iDivide(iy, ix); /* return a fraction in range 0. to 32767 = 0. to 1. */
else
iRatio = iDivide(ix, iy); /* return a fraction in range 0. to 32767 = 0. to 1. */
/* first, third and fifth order polynomial approximation */
iAngle = (Int32) K1 * (Int32) iRatio;
iTmp = ((Int32) iRatio >> 5) * ((Int32) iRatio >> 5) * ((Int32) iRatio >> 5);
iAngle += (iTmp >> 15) * (Int32) K2;
iTmp = (iTmp >> 20) * ((Int32) iRatio >> 5) * ((Int32) iRatio >> 5)
iAngle += (iTmp >> 15) * (Int32) K3;
iAngle = iAngle >> 15;
/* check if above 45 degrees */
if (iy > ix) iAngle = (Int16)(9000 - iAngle);
/* for tidiness, limit result to range 0 to 9000 equals 0.0 to 90.0 degrees */
if (iAngle < 0) iAngle = 0;
if (iAngle > 9000) iAngle = 9000;
return ((Int16) iAngle);
}
That is my code so far (btw , what does (Int16) or in xc8 (int16_t) do?):

Code: [Select]

int16_t HundredAtan(int16_t y, int16_t x) {

    const int16_t K1 = 5701;
    const int16_t K2 = -1645;
    const int16_t K3 = 446;
    int32_t Angle;
    int16_t Ratio;
    int32_t Tmp;

    if ((x == 0) && (y == 0)) {
        return 0;
    }
    if ((x == 0) && (y != 0)) {
        return 9000;
    }

    if (y <= x) {
        Ratio = Divide(y, x);
    } else {
        Ratio = Divide(x, y);
    }

    Angle = (int16_t) K1 * (int16_t) Ratio;
    Tmp = ((int16_t) Ratio >> 5) * ((int16_t) Ratio >> 5) * ((int16_t) Ratio >> 5);
    Angle += ((Ratio >> 15) * (int16_t) K2);
    Tmp = (Tmp >> 20) * ((int16_t) Ratio >> 5) * ((int16_t) Ratio >> 5);
    Angle += (Tmp >> 15) * (int16_t) K3;
    Angle = Angle >> 15;

    if (y > x) {
        Angle = (9000 - Angle);
    }

    if (Angle < 0) {
        Angle = 0;
    }
    if (Angle > 9000) {
        Angle = 9000;
    }

    return Angle;
}
I would be really thankful if you could help me. In general I would be interested in how to solve it on my own but unfortunately I don't have the time to learn it properly, but if someone could describe the problem and how to solve it in a way so I can understand it without in-depth knowledge of that kind of math and c, that would be awesome and I wouldn't have to let you do my work all the time .

[westfw]

what does (Int16) or in xc8 (int16_t) do?

You really do need to understand that, at least somewhat, to port this code, since it's using it all over the place.  An expression like:

Code: [Select]

(int32_t) a is called a "cast"; it means that regardless of the type of variable "a", it should be treated as an 32bit integer (int32_t) "at this point in the code."  Sometimes this involves format conversion - casting a float to int32_t actually converts to an integer, and casting a int16 to a int32 does 16bits of sign extension.  (other times it's more of a "re-interpret" the bits you have; for example casting an int to a pointer.)

In this code, the inputs and outputs are all 16bit integers (and sometimes Q15 scaled integers), so when the code does something like:

Code: [Select]

iAngle = (Int32) K1 * (Int32) iRatio; it means "even though K1 and iRation are 16bit numbers, I'm doing some calculations that are going to need 32bits to get enough range or precision, so you need to convert everything to 32bits first, and use 32bit math."  That means that your attempted conversion:

Code: [Select]

iAngle = (Int32) K1 * (Int32) iRatio;
//   to
Angle = (int16_t) K1 * (int16_t) Ratio;is particularly wrong.  The original says to do calculations in 32bits, and you're saying to only use 16bits.  It looks like the original code can be mostly converted simply by changing "Int16" to "int16_t" and "Int32" to "int32_t"

The original idivide code had an additional subtlety (bug, IMO) because it left out one of these casts.  In

Code: [Select]

itmp = (Int16)((itmp * ix) >> 15);the calculation of itmp*ix needs to be done with 32bits, even though the inputs are 16 and the output of the whole expression is also 16bits eventually.  Here, there are some default "promotion" rules that come into play.  C (and I guess C# as well), says that intermediate results in an expression are of type "int."  On an x86 or ARM (where C#runs) that's 32bits (and the code works), and on an AVR it's 16bits and the code fails.  It could be fixed

Code: [Select]

itmp = (Int16)(((Int32)itmp * ix) >> 15);

[Sebastian]

Thank you very much for that explanation. The funny thing is I didn't even changed the Int32s to Int16s on purpose, but I just copy and pasted them wrong because I tried different things and somehow overlooked it. I changed them now and it works. But the good thing is at least I kinda know what I am doing now. There is another function for sine and cosine and I think I will get them working on my own and again, thanks for the help.

[Sebastian]

Damn it, stuck again.
I am pretty sure the problem is similar to the one with the division.
There is a do while loop which does some kind of division. i tried to fix it but it does not work.

Original code:

Code: [Select]

const UInt16 MINDELTATRIG = 1; /* final step size for iTrig */
/* function to calculate ir = ix / sqrt(ix*ix+iy*iy) using binary division */
static Int16 iTrig(Int16 ix, Int16 iy)
{
UInt32 itmp; /* scratch */
UInt32 ixsq; /* ix * ix */
Int16 isignx; /* storage for sign of x. algorithm assumes x >= 0 then corrects later */
UInt32 ihypsq; /* (ix * ix) + (iy * iy) */
Int16 ir; /* result = ix / sqrt(ix*ix+iy*iy) range -1, 1 returned as signed Int16 */
Int16 idelta; /* delta on candidate result dividing each stage by factor of 2 */
/* stack variables */
/* ix, iy: signed 16 bit integers representing sensor reading in range -32768 to 32767 */
/* function returns signed Int16 as signed fraction (ie +32767=0.99997, -32768=-1.0000) */
/* algorithm solves for ir*ir*(ix*ix+iy*iy)=ix*ix */
/* correct for pathological case: ix==iy==0 */
if ((ix == 0) && (iy == 0)) ix = iy = 1;
/* check for -32768 which is not handled correctly */
if (ix == -32768) ix = -32767;
if (iy == -32768) iy = -32767;
/* store the sign for later use. algorithm assumes x is positive for convenience */
isignx = 1;
if (ix < 0)
{
ix = (Int16)-ix;
isignx = -1;
}
/* for convenience in the boosting set iy to be positive as well as ix */
iy = (Int16)Math.Abs(iy);
/* to reduce quantization effects, boost ix and iy but keep below maximum signed 16 bit */
while ((ix < 16384) && (iy < 16384))
{
ix = (Int16)(ix + ix);
iy = (Int16)(iy + iy);
}
/* calculate ix*ix and the hypotenuse squared */
ixsq = (UInt32)(ix * ix); /* ixsq=ix*ix: 0 to 32767^2 = 1073676289 */
ihypsq = (UInt32)(ixsq + iy * iy); /* ihypsq=(ix*ix+iy*iy) 0 to 2*32767*32767=2147352578 */
/* set result r to zero and binary search step to 16384 = 0.5 */
ir = 0;
idelta = 16384; /* set as 2^14 = 0.5 */
/* loop over binary sub-division algorithm */
do
{
/* generate new candidate solution for ir and test if we are too high or too low */
/* itmp=(ir+delta)^2, range 0 to 32767*32767 = 2^30 = 1073676289 */
itmp = (UInt32)((ir + idelta) * (ir + idelta));
/* itmp=(ir+delta)^2*(ix*ix+iy*iy), range 0 to 2^31 = 2147221516 */
itmp = (itmp >> 15) * (ihypsq >> 15);
if (itmp <= ixsq) ir += idelta;
idelta = (Int16)(idelta >> 1); /* divide by 2 using right shift one bit */
} while (idelta >= MINDELTATRIG); /* last loop is performed for idelta=MINDELTATRIG */
/* correct the sign before returning */
return (Int16)(ir * isignx);
}
My code:

Code: [Select]

int16_t Trig(int16_t x, int16_t y) {
    uint32_t Tmp;
    uint32_t xsq;
    int16_t signx;
    uint32_t hypsq;
    int16_t Result;
    int16_t Delta;

    if ((x == 0) && (y == 0)) {
        x = 1;
        y = 1;
    }
    if (x == -32768) {
        x = -32767;
    }
    if (y == -32768) {
        y = -32767;
    }

    signx = 1;
    if (x < 0) {
        x = (int16_t) - x;
        signx = -1;
    }

    y = (int16_t) myAbs(y);

    while ((x < 16384) && (y < 16384)) {
        x = (int16_t) (x + x);
        y = (int16_t) (y + y);
    }

    xsq = (uint32_t) (x * x);
    hypsq = (uint32_t) (xsq + y * y);

    Result = 0;
    Delta = 16384;

    do {
        Tmp = (uint32_t) ((Result + Delta) * (Result + Delta));
        Tmp = (uint32_t) (Tmp >> 15) * (hypsq >> 15);
        if (Tmp <= xsq) {
            Result += Delta;
        }
        Delta = (int16_t) (Delta >> 1);
    } while (Delta >= 1);

    return (int16_t) (Result * signx);
}
I hope I don't bother you to much with my stupid questions.

[dannyf]

Some of the calculations overflow for int16_t types.

[Sebastian]

I thought so. I tried to find them and change them to uint32_t but it did not work. I tried different combinations but without fully understanding it it's hopeless.

[dannyf]

I am sure you will figure it out: just look at what you did in the prior example.

It takes about 27K cycles per iTrig() invocation.

[Sebastian]

I tried to change the calculations that I think do overflow to uint32_t but it did not work.

Code: [Select]

int16_t Trig(int16_t x, int16_t y) {
    uint32_t Tmp;
    uint32_t xsq;
    int16_t signx;
    uint32_t hypsq;
    int16_t Result;
    int16_t Delta;

    if ((x == 0) && (y == 0)) {
        x = 1;
        y = 1;
    }
    if (x == -32768) {
        x = -32767;
    }
    if (y == -32768) {
        y = -32767;
    }

    signx = 1;
    if (x < 0) {
        x = (int16_t) - x;
        signx = -1;
    }

    y = (int16_t) myAbs(y);

    while ((x < 16384) && (y < 16384)) {
        x = (int16_t) (x + x);
        y = (int16_t) (y + y);
    }



    xsq = (uint32_t) (x * x);
    hypsq = (uint32_t) (xsq + (uint32_t) y * y);

    Result = 0;
    Delta = 16384;

    do {
        Tmp = (uint32_t) ((Result + Delta) * (Result + Delta));
        Tmp = (uint32_t) ( Tmp >> 15) * (hypsq >> 15);
        if (Tmp <= xsq) {
            Result += Delta;
        }
        Delta = (int16_t) (Delta >> 1);
    } while (Delta >= 1);

    return (int16_t) (Result * signx);
}

Then I tried to get rid of that complicated stuff and replace it with something I understand better, but it did not work either and it also takes up a lot of ROM which is bad because I don't have a lot.

Code: [Select]

int16_t Trig(int16_t x, int16_t y) {
    uint32_t Tmp;
    uint32_t xsq;
    int16_t signx;
    uint32_t hypsq;
    int16_t Result;
    int16_t Delta;

    if ((x == 0) && (y == 0)) {
        x = 1;
        y = 1;
    }
    if (x == -32768) {
        x = -32767;
    }
    if (y == -32768) {
        y = -32767;
    }

    signx = 1;
    if (x < 0) {
        x = (int16_t) - x;
        signx = -1;
    }

    y = (int16_t) myAbs(y);

    while ((x < 16384) && (y < 16384)) {
        x = (int16_t) (x + x);
        y = (int16_t) (y + y);
    }

    xsq = (uint32_t) (x * x);
    hypsq = (uint32_t) (xsq + (uint32_t) y * y);

    Result = (double) ((x / (sqrt(hypsq))) * 32768);

    return (int16_t) (Result * signx);
}
I have no idea what else I could try.

[Sebastian]

Could anyone please help me with that, I need that thing working tomorrow and I really don't know what to do anymore.