92

I recently ran into an issue that could easily be solved using modulus division, but the input was a float:

Given a periodic function (e.g. sin) and a computer function that can only compute it within the period range (e.g. [-π, π]), make a function that can handle any input.

The "obvious" solution is something like:

#include <cmath>

float sin(float x){
    return limited_sin((x + M_PI) % (2 *M_PI) - M_PI);
}

Why doesn't this work? I get this error:

error: invalid operands of types double and double to binary operator %

Interestingly, it does work in Python:

def sin(x):
    return limited_sin((x + math.pi) % (2 * math.pi) - math.pi)
Community
  • 1
  • 1
Brendan Long
  • 53,280
  • 21
  • 146
  • 188
  • 20
    π is not equal to 3.14, and in fact its not representable as any floating point type. Computing `sin(x)` for large values of `x` actually requires a very difficult transcendental argument reduction process that cannot get by with any finite approximation of pi. – R.. GitHub STOP HELPING ICE May 23 '11 at 21:03
  • While all of the answers are without doubt helpful nobody seems to know the answer to this question. – Nikolai Ruhe May 23 '11 at 21:04
  • 3
    This is almost certainly a homework assignment, so floating point errors are either outside the scope of the assignment, or this is meant to lead to a discussion of more rigorous numerical analysis. Either way, `fmod` is likely what the instructor is looking for. – Dennis Zickefoose May 23 '11 at 21:16
  • 2
    It's not homework, it's just something that came up while reading another SO question (http://stackoverflow.com/questions/6091837/sin-and-cos-are-slow-is-there-an-alternatve/6091846#6091846) – Brendan Long May 23 '11 at 22:10
  • In that case, the comments are quite on topic. – Dennis Zickefoose May 24 '11 at 01:51
  • @R.. It seems to me that a finite, 1000-decimal approximation of pi should be enough for the naive implementation of double precision argument reduction. We are at worst talking about reducing a number near DBL_MAX (in the 10^308s) into a denormal (in the 10^-324s). What could go wrong? (not that `fmod` will be very useful in implementing the argument reduction) – Pascal Cuoq Jul 14 '12 at 20:50
  • 2
    OK, I should have been more precise in my statement. My point was that if the argument can grow unboundedly large (not just double precision exponent size), no finite approximation of pi will suffice. For double, yes, a very very long approximation of pi will suffice. – R.. GitHub STOP HELPING ICE Jul 15 '12 at 00:21
  • @R.: But if `x` is so large in magnitude that `x + M_PI` has a noticeable precision loss, then you probably need a more precise representation for `x` anyway if you want any reasonable estimate of its sine. – aschepler Jul 10 '13 at 01:12
  • 1
    @aschepler: I don't think you understood the issue. – R.. GitHub STOP HELPING ICE Jul 10 '13 at 01:30

9 Answers9

85

Because the normal mathematical notion of "remainder" is only applicable to integer division. i.e. division that is required to generate integer quotient.

In order to extend the concept of "remainder" to real numbers you have to introduce a new kind of "hybrid" operation that would generate integer quotient for real operands. Core C language does not support such operation, but it is provided as a standard library fmod function, as well as remainder function in C99. (Note that these functions are not the same and have some peculiarities. In particular, they do not follow the rounding rules of integer division.)

AnT stands with Russia
  • 312,472
  • 42
  • 525
  • 765
  • 9
    Of worth, from the definition of % in the 98 standard: "(a/b)*b + a%b is equal to a." For floating point types, `(a/b)*b` already equals `a` [insofar as such a statement can be made for floating point types], so `a%b` would never be particularly useful. – Dennis Zickefoose May 23 '11 at 21:10
  • 1
    @Dennis: Indeed, algebraically, in a field the remainder is always 0. The most appropriate definition of the `%` operator for floating point, I suppose, would be `a-(a/b)*b`, which would be either 0 or a very small value. – R.. GitHub STOP HELPING ICE May 23 '11 at 21:17
  • 7
    @Dennis: You can easily fix this formula by requiring that "floor(a/b)*b + a%b = a". Note that for integers, floor(a/b) = a/b. – vog May 23 '11 at 21:22
  • C-style integer division uses trunc, not floor, but the point remains. – dan04 Jul 10 '13 at 02:23
  • 29
    **-1** Re "the normal mathematical notion of "remainder" is only applicable to integer division", the mathematical notion of modulo arithmetic works as well for floating point values, and this is one of the first issues that Donald Knuth discusses in his classic *The Art of Computer Programming* (volume I). I.e. it was once basic knowledge. Today, students don't receive the education they pay for, IMHO. – Cheers and hth. - Alf Aug 23 '14 at 10:51
60

You're looking for fmod().

I guess to more specifically answer your question, in older languages the % operator was just defined as integer modular division and in newer languages they decided to expand the definition of the operator.

EDIT: If I were to wager a guess why, I would say it's because the idea of modular arithmetic originates in number theory and deals specifically with integers.

Doug Stephen
  • 7,181
  • 1
  • 38
  • 46
  • 11
    "older languages" - APL dates back to the 1960's, and its modulo operator "|" works with both integers and floating point data (also with scalar, vector, matrix, tensor, ...) . There's no good reason that C's "%" modulo operator couldn't have performed the same function as fmod if used with floating point numbers. – rcgldr Aug 23 '14 at 01:16
  • @rcgldr The design goals didn't require floating point modulo. C was implemented to compile Unix and to limit the amount of assembly language needed for the OS. "C is an imperative procedural language. It was designed to be compiled using a relatively straightforward compiler, to provide low-level access to memory, to provide language constructs that map efficiently to machine instructions, and to require minimal run-time support. " https://en.wikipedia.org/wiki/C_(programming_language) – harper Feb 08 '19 at 08:52
  • 1
    @harper - Since C includes floating point arithmetic such as add, subtract, multiply, and divide, using the same syntax that it does for integers, I don't see why it could not also have included modulo using the same syntax (%). The choice to include it or not seems arbitrary. – rcgldr Feb 08 '19 at 09:28
18

I can't really say for sure, but I'd guess it's mostly historical. Quite a few early C compilers didn't support floating point at all. It was added on later, and even then not as completely -- mostly the data type was added, and the most primitive operations supported in the language, but everything else left to the standard library.

Jerry Coffin
  • 476,176
  • 80
  • 629
  • 1,111
  • 1
    +1 for being the first reasonable answer I see as I'm reading down the list. Actually, having now read them all, this is the *only* reasonable answer. – Cheers and hth. - Alf Aug 23 '14 at 10:54
  • A belated +1 from me too. I used to write in C for 6809 and Z80 embedded systems. No way could I afford the space to include the c runtime library. Even had to write my own startup code. Floating point was a luxury I could not afford :) – Richard Hodges Jan 02 '16 at 15:44
12

The modulo operator % in C and C++ is defined for two integers, however, there is an fmod() function available for usage with doubles.

Mark Elliot
  • 75,278
  • 22
  • 140
  • 160
  • 5
    This is the answer to OP's question, but ignores the fundamental problem in what OP is trying to do: `sin(fmod(x,3.14))` or even `sin(fmod(x,M_PI))` is not equal to `sin(x)` for large values of `x`. In fact the values may differ by as much as 2.0. – R.. GitHub STOP HELPING ICE May 23 '11 at 21:05
  • 2
    @R..: Right, but that's a different question and I'm not entirely sure there's an accepted answer, though there's plenty of [research on the subject](http://scholar.google.com/scholar?hl=en&q=The+evaluation+of+periodic+functions+with+large+input+arguments&btnG=Search&as_sdt=1%2C21&as_ylo=&as_vis=0) – Mark Elliot May 23 '11 at 21:13
  • @R - I fixed the equation to do it correctly. The actual equation wasn't the point (it was fairly easy to figure out once I had the function to test it with). – Brendan Long May 23 '11 at 23:37
  • Is not `%` the remainder operator and not a modulo operator? – chux - Reinstate Monica Mar 03 '14 at 19:51
7

The constraints are in the standards:

C11(ISO/IEC 9899:201x) §6.5.5 Multiplicative operators

Each of the operands shall have arithmetic type. The operands of the % operator shall have integer type.

C++11(ISO/IEC 14882:2011) §5.6 Multiplicative operators

The operands of * and / shall have arithmetic or enumeration type; the operands of % shall have integral or enumeration type. The usual arithmetic conversions are performed on the operands and determine the type of the result.

The solution is to use fmod, which is exactly why the operands of % are limited to integer type in the first place, according to C99 Rationale §6.5.5 Multiplicative operators:

The C89 Committee rejected extending the % operator to work on floating types as such usage would duplicate the facility provided by fmod

Yu Hao
  • 119,891
  • 44
  • 235
  • 294
2

The % operator does not work in C++, when you are trying to find the remainder of two numbers which are both of the type Float or Double.

Hence you could try using the fmod function from math.h / cmath.h or you could use these lines of code to avoid using that header file:

float sin(float x) {
    float temp;
    temp = (x + M_PI) / ((2 *M_PI) - M_PI);
    return limited_sin((x + M_PI) - ((2 *M_PI) - M_PI) * temp));
}
AustinWBryan
  • 3,249
  • 3
  • 24
  • 42
Shashank
  • 36
  • 4
2

"The mathematical notion of modulo arithmetic works for floating point values as well, and this is one of the first issues that Donald Knuth discusses in his classic The Art of Computer Programming (volume I). I.e. it was once basic knowledge."

The floating point modulus operator is defined as follows:

m = num - iquot*den ; where iquot = int( num/den )

As indicated, the no-op of the % operator on floating point numbers appears to be standards related. The CRTL provides 'fmod', and usually 'remainder' as well, to perform % on fp numbers. The difference between these two lies in how they handle the intermediate 'iquot' rounding.

'remainder' uses round-to-nearest, and 'fmod' uses simple truncate.

If you write your own C++ numerical classes, nothing prevents you from amending the no-op legacy, by including an overloaded operator %.

Best Regards

Love
  • 21
  • 3
  • A somewhat better definition would use the mathematical value floor(num/den+0.5). Such a definition would always yield a result that was precisely representable, while Knuth's definition would not do so for e.g. (1E-37, 1.0), since the mathematically-precise result for Knuth's definition would be (1E-37 - 1.0). – supercat Jan 27 '22 at 23:39
  • From a practical perspective, I think Mr. Knuth's definition would probably be essentially as useful. Defining the operator in such a way that it never requires rounding may be aesthetically desirable, but it's not nearly as important as saying that two numbers which differ by some exact multiple of the denominator should belong to the same modular equivalence class. Indeed, while lack of rounding might be aesthetically desirable, in cases where (x+den)-den and (x-den)+den would belong to the same equivalence class as each other, but not to x, it may be better to put x into that class. – supercat Feb 11 '22 at 16:42
  • BTW, what do you think of the concept of having "rounding modes", rather than having separate "add, taking ceiling of result" and "add, taking floor of the result" operations? How often is it useful to have a piece of code that is sometimes executed in one mode and sometimes the other? If e.g. a program does `x = 16.0; y=0.1; z=x+y;`, replacing the latter statement with `z=1.1;` would be an improvement if the execution-time rounding semantics were known to match those of the compiler, but an implementation that supports rounding modes wouldn't be able to do that. – supercat Feb 11 '22 at 18:17
  • @supercat I find it preferrable, in most cases, to let the caller decide how to run things. That's why I opted to let them control the FPU rounding instead of handing them a precooked "fix". As You probably know, the FPU has _four_ rounding modes. The behavior of the modulus operator depends on the FPU rounding mode. If the RM is RM_NEAREST (the default), `%=` behaves like `remainder`, and if it's RM_TRUNC,` %=` behaves like `fmod`. The behavior is not specified for RM_FLOOR and RM_CEIL, so feel free to explore. – Love Feb 18 '22 at 04:59
  • Addendum.. Definition of rounding mode flags: `#define RM_NEAREST (0<<10) #define RM_FLOOR (1<<10) #define RM_CEIL (2<<10) #define RM_TRUNC (3<<10)` And here's how to set the rounding mode (in x64 asm): `_setrm87 PROC ;; extern "C" void __stdcall _setrm87( ushort rm ); fnstcw [rsp-4] mov ax, [rsp-4] and ax, not CW_RM_MASK ; Mask away the RM bits and cx, CW_RM_MASK ; Just the RM bits (rcx=rm) or ax, cx mov [rsp-4], ax fldcw [rsp-4] ret _setrm87 ENDP` Sorry about the format. Cut/paste it to your editor and reformat. Best Regards – Love Feb 18 '22 at 05:24
  • If a function performs multiple FPU operations (which would be the case with most functions that perform any), how often is it useful to globally set an FPU mode that will be applied to all of them? I would think that if one wants to compute the minimum possible value for `a*b-c*d`, one would want to have `a*b` rounded toward -inf, and `c*d` rounded toward +inf. Rounding everything the same direction would yield a meaningless result. – supercat Feb 18 '22 at 08:22
  • @supercat In case you need different rounding for different terms of an expression you must of course use the likes of `floor` and `ceil` for the terms. Otherwise the active rounding mode applies to everything. The default is round-to-nearest. However, you can change the active rounding mode to fit your purpose, using the `FLDCW` instruction. But even then you'll have to use RTL rounding functions explicitly if you need to do what you said. Best Regards – Love Feb 19 '22 at 18:05
  • @supercat Oh, BTW... `int( 1E-37 / 1.0 )` is 0, not 1, so Knuth's definition would yield 1E-37. Best Regards – Love Feb 19 '22 at 18:36
  • Sorry--I meant `(-1E-37,1.0)`. The IEEE definition yields -1E-37, which would be fine if for any value y there were only one value of x such that `fmod(x,1.0)` would yield `y`. If one were to specify that `fmod(-0.5,1.0)` and `fmod(0.5,1.0)` both yield `0.5`, and all other cases use round-to-nearest, that could work nicely. – supercat Feb 19 '22 at 18:56
2

try fmod

Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680
justin
  • 104,054
  • 14
  • 179
  • 226
1

For C/C++, this is only defined for integer operations.

Python is a little broader and allows you to get the remainder of a floating point number for the remainder of how many times number can be divided into it:

>>> 4 % math.pi
0.85840734641020688
>>> 4 - math.pi
0.85840734641020688
>>> 
AustinWBryan
  • 3,249
  • 3
  • 24
  • 42
Andrew
  • 2,530
  • 16
  • 9
  • 3
    Remainder is *not* 'another name for modulus' !! See: http://stackoverflow.com/questions/13683563/whats-the-difference-between-mod-and-remainder or from a math-pov: http://math.stackexchange.com/questions/801962/difference-between-modulus-and-remainder Put simply: modulo and remainder are *only* the same for positive numbers and another example is that remainder doesn't let you go around the compass (counter-clockwise). Please correct that as I am to frugal to downvote `:P` – GitaarLAB Mar 21 '16 at 21:31