4

I wrote a basic fixed-point variant of the sine function which utilizes a lookup table (targeting an AVR micro controller w/o FPU). My implementation also accepts negative values and values which exceed 2π, like its floating point pendant in math.h does.

So I need to map the given value to a range between 0 and 2π (i.e. their fixed point counterparts). For positive arguments it's easy to trim them with C's build-in remainder operator %. As this is not an option for negative values, I'm using the following (obvious) approach:

    unsigned int modulus = (a - (INT_FLOOR(a / b) * b) + b) % b;

a and b are integer typed values and INT_FLOOR() is just meant as a hint that the fractional part of (a/b) is truncated. This formula ensures that the calculated modulus (which is used as an index for the table array) is always positive and that also negative arguments are mapped their positive counterparts (maintaining phase shifts in both directions).

My problem with this approach is that it seems to be overly complex as it involves no less than five arithmetic operations. Is there a more efficient approach I am missing?

Timmy O'Mahony
  • 53,000
  • 18
  • 155
  • 177
Developer's Destiny
  • 131
  • 9
  • 7
    The sine is an odd function, so `sin(-x) = -sin(x)`, the cosine is even, `cos(-x) = cos(x)`. That should simplify things. – Daniel Fischer May 03 '12 at 20:34
  • 2
    Can you use [`fmod`](http://www.cplusplus.com/reference/clibrary/cmath/fmod/) to do remainder of floats? It works correctly for negative numbers, as well as positives. – Sergey Kalinichenko May 03 '12 at 20:36
  • 1
    You seem to be wanting essentially the remainder of the absolute values, so why not do that directly: `unsigned in modulus = abs(a) % abs(b);` Alternatively, something like: `modulus = a % b; if (modulus < 0) modulus += b;` – Jerry Coffin May 03 '12 at 20:36
  • @DanielFischer: Thank you very much, I was missing the obvious (which might explain why someone tagged this as homework). Unfortunately, my naïve approach still seems to be the better solution as the resulting code size is shorter by 26 bytes (avr-gcc -O2). dasblinkenlight: fmod is for floating point math only, which is something I want to avoid at all costs on that architecture JerryCoffin: using abs() here and there increases the code size more than using my naïve approach, avr-gcc does strange things sometimes – Developer's Destiny May 03 '12 at 22:42

2 Answers2

4

Unless your integer arguments have been scaled such that they're in terms of multiples of π (e.g. 65536 means 2π), trying to do argument reduction is probably misguided, as 2π is irrational and any reduction mod 2π will introduce error that scales up with the number of periods until the entire result of the reduction becomes error. This is actually a serious problem in many floating point trig implementations.

I would recommend either not doing argument reduction at all, or using an angle scale based on a power of two rather than radians (so, for example, 0x10000 or 0x1000000 corresponds to 2π or 360 degrees). Then argument reduction becomes a single bitwise and operation.

R.. GitHub STOP HELPING ICE
  • 208,859
  • 35
  • 376
  • 711
  • I'm aware of accumulating rounding errors. Fortunately, they are within acceptable ranges regarding my tasks even after several hundred rounds (I don't exceed 1000 rounds and the interim results aren't deeply nested). Since I'm using a lookup tables on architectures with (at most) 65536 bytes of flash memory, I don't see how I could get away without using argument reduction. Your second suggestion is very interesting, though. It involves modifications of all surrounding routines who use that sine function but it might increase precision a lot. Definitely something I need more thinking about. – Developer's Destiny May 03 '12 at 23:58
  • I think it would increase both precision and performance a lot. – R.. GitHub STOP HELPING ICE May 04 '12 at 00:15
0

I see your formula correct, although not so pretty.

If you make b a power of 2 then some operations can be done with bit masks. Something along the lines of:

unsigned int modulus = (a - (a & ~(b-1)) + b) & (b-1);

And if b is a constant or a macro it should optimize quite a bit.

rodrigo
  • 94,151
  • 12
  • 143
  • 190
  • `b` is indeed a macro but it isn't practical to modify it or the surrounding calculations to handle powers of two (which results in both code size increase and precision loss). In my case, `b` is the amount of temporal quantization steps of my approximated sine curve. It has to be a divisor of my fixed-point version of 2π so that I can efficiently map fixed-point angles to their corresponding quantization step (which indeed involves just one bit shift). But thanks for the hint nonetheless. – Developer's Destiny May 03 '12 at 23:23
  • And that's why old school math programs (MS-DOS games for example) divided the cirle in 2^x (512 or 256) _binary degrees_. That old 360 degree custom has no use in computers. – rodrigo May 04 '12 at 09:35