32 bit div

[GoogleCodeExporter]

Hi,


32 bit lib seems very expensive in terms of number of cycles and accuracy. 

Is there any one working on it ?

Original issue reported on code.google.com by birdla...@gmail.com on 3 Mar 2011 at 3:23

[GoogleCodeExporter]

There are only 2 people on this project at the minute since it's very new. I've 
implemented all the 32-bit stuff so far, but I only have so much time to work 
on it amongst full-time work many other projects and commitments. If you want 
better/faster algorithms and are interested in implementing them then I'd be 
happy to give you committer rights.

As the division goes, division will always be a slow operation, there aren't 
that many optimizations I can think to make to the algorithm I used (not 
without assembler anyway) but that doesn't mean that there isn't a 
fundamentally better algorithm out there that could be used.

The inaccuracy must be due to a small error somewhere in the code, the fact 
that it's such a small error indicates that it's likely to be a rounding error 
or an error in some extreme case where the result overflows.

Original comment by Flatmush@googlemail.com on 3 Mar 2011 at 3:35

• Added labels: Performance

[GoogleCodeExporter]

I'd be interested to hear what you're using the library for birdland (if you 
can tell me), I might be able to help you achieve what you want better that way 
with the limited time I have to work on this project.

Original comment by Flatmush@googlemail.com on 3 Mar 2011 at 3:45

[GoogleCodeExporter]

hi 
I need these fixed point functions in the baseband processing of some wireless 
communication protocol such as WLAN channel estimation. Timing rquirements are 
very strict. Please run the following code and see how many cycles division is 
taking. Please which processor you are runing on. Please tell me also what 
value of d you get ? 


#include <stdio.h>
#include <stdlib.h>


typedef __int32_t fix16_t;

fix16_t fix16_div(fix16_t inArg0, fix16_t inArg1);
fix16_t fix16_mul(fix16_t inArg0, fix16_t inArg1);
volatile int d;

int main(){

    int x = 5063265; //77.2593;
    int y = -6338360;//-96.7157;

    int a;
    d = fix16_div(x, y);
    printf("%d \n",d);
    //printf("%d \n",a);

}

fix16_t fix16_div(fix16_t inArg0, fix16_t inArg1) {
int neg = ((inArg0 < 0) != (inArg1 < 0));
        inArg0 = (inArg0 < 0 ? -inArg0 : inArg0);
        inArg1 = (inArg1 < 0 ? -inArg1 : inArg1);

        while(((inArg0 | inArg1) & 1) == 0) {
                inArg0 >>= 1;
                inArg1 >>= 1;
        }

        __uint32_t r_hi = (inArg0 / inArg1);

        __uint32_t n_lo = (inArg0 % inArg1);
        __uint32_t n_hi = (n_lo >> 16);
        n_lo <<= 16;

        __uint32_t i, arg;
        for(i = 1, arg = inArg1; ((n_lo | arg) & 1) == 0; i <<= 1) {
                n_lo = ((n_lo >> 1) | (n_hi << 31));
                n_hi =  (n_hi >> 1);
                arg >>= 1;
        }

        __uint32_t res = 0;
        if(n_hi) {
                __uint32_t arg_lo, arg_hi;
                for(arg_lo = inArg1; (arg_lo >> 31) == 0; arg_lo <<= 1, i <<= 1);
                for(arg_hi = (arg_lo >> 31), arg_lo <<= 1, i <<= 1; arg_hi < n_hi; arg_hi = (arg_hi << 1) | (arg_lo >> 31), arg_lo <<= 1, i <<= 1);

                do {
                        arg_lo = (arg_lo >> 1) | (arg_hi << 31);
                        arg_hi = (arg_hi >> 1);
                        i >>= 1;
                        if(arg_hi < n_hi) {
                                n_hi -= arg_hi;
                                if(arg_lo > n_lo)
                                        n_hi--;
                                n_lo -= arg_lo;
                                res += i;
                        } else if((arg_hi == n_hi) && (arg_lo <= n_lo)) {
                                n_hi -= arg_hi;
                                n_lo -= arg_lo;
                                res += i;
                        }
                } while(n_hi);
        }

        res += (n_lo / inArg1);
        #ifndef FIXMATH_NO_ROUNDING
        if((n_lo % inArg1) >= (inArg1 >> 1))
                res++;
        #endif
        res += (r_hi << 16);

        return (neg ? -res : res);
}

fix16_t fix16_mul(fix16_t inArg0, fix16_t inArg1) {
    __int16_t hi[2] = { (inArg0 >> 16),    (inArg1 >> 16)    };
    __uint16_t lo[2] = { (inArg0 & 0xFFFF), (inArg1 & 0xFFFF) };

     __int32_t r_hi = hi[0] * hi[1];
     __int32_t r_md = (hi[0] * lo[1]) + (hi[1] * lo[0]);
    __uint32_t r_lo = lo[0] * lo[1];
    #ifndef FIXMATH_NO_ROUNDING
    r_lo += 0xFFFF;
    #endif

    r_md += (r_hi & 0xFFFF) << 16;
    r_md += (r_lo >> 16);

    return r_md;

}

Original comment by birdla...@gmail.com on 3 Mar 2011 at 5:04

[GoogleCodeExporter]

Running the code I get -49552.  This is using gcc 4.3.1.  I don't have time 
right now to look into the excecution speed, but ill take a look later today

Original comment by joe.scha...@gmail.com on 3 Mar 2011 at 5:35

[GoogleCodeExporter]

Yes I get the same value but the accurate one is -52352. cycle count seems very 
high for this code. 

let see what you get for cycle count

Original comment by birdla...@gmail.com on 3 Mar 2011 at 6:29

[GoogleCodeExporter]

How are you measuring cycle count

Original comment by joe.scha...@gmail.com on 3 Mar 2011 at 6:36

[GoogleCodeExporter]

Birdland, are you targetting ARM, if so then we could make a much faster 
assembler optimized version, especially if it has a clz instruction.
This implementation is very generic and aimed at accuracy, there's a similar 
divide function to divide a 64-bit number by a 32-bit number in the newly 
submitted int64.h header which may work better, though it is largely the same.

Original comment by Flatmush@googlemail.com on 3 Mar 2011 at 7:52

[GoogleCodeExporter]

Ok so I have never run a profiler before so I am not entirely sure what I am 
doing here, but here it goes.

so as a test I took birdlands code and copied it into fix16.c appending an _b 
to the function name.  
The standard fix16_div is using 64 bit math.
I threw in mul for good measure
I modified the main code to be:

    int x = 5063265; //77.2593;
    int y = -6338360;//-96.7157;
    d = fix16_div_b(x, y);
    k = fix16_div(x,y);
    t = fix16_mul(x,y);
    printf("%d \n",d);
    printf("%d \n",k);
    printf("%d \n",t);
    return 0;

compiled it with gcc 4.2.1 with the -g option. Then used valgrind to look at 
the number of instructions per function. Here is what I came up with.

#############################################################################
Results for Core Duo @ 2.0Ghz(32 bit)


valgrind --tool=callgrind --toggle-collect=fix16_div_b ./main
No optimization =  267
O3                        = 198
Os                         = 209

valgrind --tool=callgrind --toggle-collect=fix16_div ./main
No optimization = 2341
O3                        = 2334
Os                         = 2334

valgrind --tool=callgrind --toggle-collect=fix16_mul ./main
No optimization = 56
O3                         = 10
Os                         = 10

#############################################################################
Results for Core 2 Duo @ 2Ghz(64 bit)

valgrind --tool=callgrind --toggle-collect=fix16_div_b ./main
No optimization = 268
O3                         = 188
Os                         = 200

valgrind --tool=callgrind --toggle-collect=fix16_div ./main
No optimization = 23
O3                        = 14
Os                         = 14

valgrind --tool=callgrind --toggle-collect=fix16_mul ./main
No optimization = 18
O3                         = 10
Os                         = 10

#############################################################################
Results for Arm Cortex M0 @ 48 Mhz (32 bit)
(this uses a timer, which is running at the system clock speed) since it is 
single threaded 
it always runs in the same amount of time.
Also I only compiled with Os, because the code grows too large without it.

also keep in mind that this processor has NO division instruction.  So some 
magic is happening
which is probably why the 32bit and 64 bit functions are similar

fix16_div_b
Os                     = 880

fix16_div
Os                     = 1104

fix16_mul 
Os                     = 145

#############################################################################

Wow.  So the fix16_div routine is very slow for 32 bit systems even if the 
compiler can handle 64 bit numbers
The ARM doesnt take as large of a hit, but that is because the div_b is slow as 
it is, mostly due to the lack of a division routine

Hope all of this makes sense...

Original comment by joe.scha...@gmail.com on 4 Mar 2011 at 2:26

[GoogleCodeExporter]

It's expected that 64-bit division would be slow (though concrete results are 
very nice to have) because there is no trivial way to perform division unlike 
multiplication.
I think that the algorithm I wrote could do with some work, but I'm not sure 
there is going to be a very fast way of doing it.

Original comment by Flatmush@googlemail.com on 4 Mar 2011 at 9:27

[GoogleCodeExporter]

I've found a good way to do fixed point division on arm: 
http://me.henri.net/fp-div.html
It'll be very useful for you birdland, I don't want to steal it to put it in 
the SDK though so I'm going to have to contact him and check if he wants to 
join in or if he's happy for his work to go into our project.

Original comment by Flatmush@googlemail.com on 4 Mar 2011 at 12:27

[GoogleCodeExporter]

Well we have Henri's permission (and encouragement) to use his code as part of 
our library, so I'll add it when I get time (unless one of you want to).

Henri:
"By all means please use the code. It took a lot of work and research to figure 
out, and at the time, I didn't want it all to go to waste - so I wrote it all 
down in that article - but no-one was interested in publishing it at the time. 
For this reason I would be very happy to see it part of an open source library. 
That being said, I don't do ARM work any more - I don't even have the build 
tools installed, so if you can use or adapt the code 'as-is' please let that be 
my contribution."

Original comment by Flatmush@googlemail.com on 4 Mar 2011 at 2:09