23

I am trying to implement the range reduction operation for trigonometry. But instead I think it might be better to just perform a modulo pi/2 operation on incoming data. I was wondering what algorithms exist and are efficient for this operation for 32-bit IEEE 754 floating-point?

I have to implement this in assembly, so fmod, division, multiplication, etc. aren't available to me with just one instruction. My processor uses 16-bit words and I have implemented 32-bit floating point addition, subtraction, multiplication, division, square root, cosine, and sine. I just need range reduction (modulus) for inputting values to cosine and sine.

Greg
  • 9,068
  • 6
  • 49
  • 91
Veridian
  • 3,531
  • 12
  • 46
  • 80
  • 5
    Actually there are a lot of clever algorithms for instance google for "payne hanek range reduction", but I think thats not what you want – Gunther Piez Feb 29 '12 at 20:34
  • 1
    The paper by Ng you linked to in a previous related question of yours actually explains the Payne-Hanek algorithm, which AFAIK is still the state of the art for accurate range reduction. You just have to adapt it to single precision. – janneb Feb 29 '12 at 21:09
  • 1
    @Everyone, please delete/edit your answer so that it applies to my actual question. I am looking for the algorithm within a floating-point modulus. I need to implement what fmod does and minimize the number of divides I perform. – Veridian Mar 01 '12 at 16:41
  • Thanks - fmod was actually exactly what I was looking for (on another project). – Danny Staple Nov 22 '14 at 18:00
  • `remainder()` is like `fmod()` but with IEEE-standard round-to-nearest. – Peter Cordes Apr 09 '16 at 04:54
  • 1
    Please note: any modulo technique involving a floating point approximation will be worthless for larger numbers. If you have an approximation of pi to 16 digits, then accurately dividing a 17-digit number by your approximation can have an error greater than 1, meaning the remainder could be absolutely anywhere in the range 0..pi, revealing nothing about the remainder you were really looking for. – mwfearnley Oct 17 '16 at 23:04

4 Answers4

18

I think standard library's fmod() will be the best choice in most cases. Here's a link to a discussion of several simple algorithms.

On my machine, fmod() uses optimized inline assembly code (/usr/include/bits/mathinline.h):

#if defined __FAST_MATH__ && !__GNUC_PREREQ (3, 5)
__inline_mathcodeNP2 (fmod, __x, __y, \
  register long double __value;                           \
  __asm __volatile__                                  \
    ("1:    fprem\n\t"                            \
     "fnstsw    %%ax\n\t"                             \
     "sahf\n\t"                                   \
     "jp    1b"                               \
     : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc");           \
  return __value)
#endif

So it actually uses a dedicated CPU instruction (fprem) for the calculation.

Helicase
  • 33
  • 3
Michał Kosmulski
  • 9,855
  • 1
  • 32
  • 51
  • Oh, I'm actually trying to implement what fmod does. That is the issue, I am looking for the modulus algorithm for floating point. – Veridian Feb 29 '12 at 20:11
  • The most direct form is probably (code taken from link in my post, but it's kind of the definition of modulus for floating point and thus the obvious way of doing it): template< typename T > T fmod( T x, T y ) { T a = (T)(long long)( x / y ); return x - a * y; } – Michał Kosmulski Feb 29 '12 at 20:17
  • I'm a little worried about rounding on that a*y product, but I'm not sure how to mitigate it. – zmccord Feb 29 '12 at 21:00
  • 2
    The obvious way is unfortunately very inaccurate for large x. fprem is better, but doesn't provide "last bit" accuracy either, for that the Payne-Hanek algorithm is the tool of choice. – janneb Feb 29 '12 at 21:06
  • For future readers: [`fmod()`](http://man7.org/linux/man-pages/man3/fmod.3.html) rounds towards zero, while [`remainder()` rounds to nearest](http://man7.org/linux/man-pages/man3/remainder.3.html). For `remainder`, use `fprem1` instead of `fprem` if using x87 at all. – Peter Cordes Apr 09 '16 at 04:53
  • I found my compiler wouldn't vectorize a loop calling fmod – Old Badman Grey Sep 11 '17 at 11:56
16

Maybe I'm missing the point here, but do you have anything against simply using fmod?

double theta = 10.4;
const double HALF_PI = 2 * atan(1);
double result = fmod(theta, HALF_PI);
Prashant Kumar
  • 20,069
  • 14
  • 47
  • 63
  • 2
    Oh, I'm actually trying to implement what fmod does. That is the issue, I am looking for the modulus algorithm for floating point. – Veridian Feb 29 '12 at 20:11
  • 2
    `fmod` is ok as long as you don't care about precision with large arguments. – Gunther Piez Feb 29 '12 at 20:31
  • 2
    Unless the OP is talking about an environment where fmod isn't available. – Prashant Kumar Feb 29 '12 at 20:32
  • 1
    Unless your math library is correctly rounded (that is, error less than 0.5 ulp, and no, many math libraries aren't), it's better to just use a literal for pi/2. – janneb Feb 29 '12 at 21:01
  • I have to implement this in assembly, so fmod, division, multiplication, etc. aren't available to me with just one instruction. My processor uses 16-bit words and I have implemented 32-bit floating point addition, subtraction, multiplication, division, square root, cosine, and sine. I just need range reduction (modulus) for inputting values to cosine and sine. – Veridian Feb 29 '12 at 21:30
10

The algorithm you want, to limit a floating point value between 0 and some modulus n:

Double fmod(Double value, Double modulus)
{
    return value - Trunc(value/modulus)*modulus;
}

for example pi mod e (3.14159265358979 mod 2.718281828459045)

3.14159265358979 / 2.718281828459045 
   = 1.1557273497909217179

Trunc(1.1557273497909217179)
   = 1

1.1557273497909217179 - 1
   = 0.1557273497909217179

0.1557273497909217179 * e
   = 0.1557273497909217179 * 2.718281828459045
   = 0.42331082513074800

pi mod e = 0.42331082513074800

Ian Boyd
  • 246,734
  • 253
  • 869
  • 1,219
  • 4
    This was particularly helpful for me because - although the initial question was asked in the context of C/C++ programming, I came to this particular question needing a general formula for doing this in a fixed-point number system I'm working in. I'm glad you posted this, because fmod() wasn't a solution to my needs, even though it may have been for the OP. There are quite a few people who need this particular formula in other contexts. – Richard Kettering Aug 18 '16 at 06:25
  • 1
    This can be [very inaccurate for large `value`s](https://stackoverflow.com/questions/9505513/floating-point-modulo-operation/37977740#comment12038890_9505761). There are more sophisticated algorithms. But this will usually be fairly fast, so it's a good choice if it's numerically ok for your use case. – Peter Cordes Oct 08 '17 at 02:10
0

Exact fmod is implemented with long division. The exact remainder is always representable as the dividend and the divisor share the same format. You can look into open-source implementations like glibc and musl. I have also made one in metallic. (shameless plug)

Payne–Hanek range reduction is for constant divisors like π, whose reciprocal we store in advance. Hence, it is not applicable here.

jdh8
  • 3,030
  • 1
  • 21
  • 18