Strange behavior of Fmod operation in Float32Full

tosjduenfs · 2013-10-12 22:49

I'm trying to do a modulo operation in the following:

'Calculate the Local Hour Angle:  LHA = GMST + Longitude - RA  where 0 <= LHA <=360

 LHA := FLT.Fmod(FLT.Fsub(FLT.Fadd(GMST,LongDecimal),RaDecimal),360.0)

parallax serial terminal says
LHA = -147.0883
GMST = 24.375
LongDecimal = -87.83828
RaDecimal = 83.625

The correct result should be about 212.9 but it gives me -147.0883.

Has anyone else had an issue with using the modulo operation on a negative number? Am I using it incorrectly?

Thanks,
Mike

kuroneko · 2013-10-13 05:12

tosjduenfs wrote: »

The correct result should be about 212.9 but it gives me -147.0883.

They are both valid remainders for this operation (see below), it depends on the implementation which one is returned. IOW if you need it to be between 0.0 and 360.0 check for <0 and add 360.0.

24.375 + (-87.83828) - 83.625 = -147.08828

[COLOR="blue"]-147.08828 mod 360.0[/COLOR]

-147.08828 = -147.08828 - 360.0 * (+0)
+212.91172 = -147.08828 - 360.0 * (-1)

Mark_T · 2013-10-13 06:36

From what I recall IEEE float operations are fully specified, no ambiguity exists.

Heater. · 2013-10-13 10:05

Mark_T,

Standards are great. Let's see what the different languages actually do:

Python:

# python
Python 2.7.3 (default, Jan  2 2013, 13:56:14) 
[GCC 4.7.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> -147.08828 % 360.0
212.91172

JavaScript:

# node
> -147.08828 % 360.0
-147.08828
>

C:

# cat  main.c 
#include <stdio.h>
#include <math.h>

int main(int argc, char* argv[])
{
        printf ("%f\n", fmodf(-147.08828, 360.0));
}
# ./main 
-147.088287

C#: (Using IEEERemainer function)

                      IEEERemainder              Modulus
-147.08828 / 360 =          -147.08828           -147.08828
147.08828 / 360 =           147.08828            147.08828

bc: Arbitrary precision calculator on Linux

# bc
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
-147.08828 % 360.0
-147.08828

Turns out only Python agrees with you here.

tosjduenfs · 2013-10-13 10:09

I guess I thought a modulo operation should always give a positive result. Regardless, I implemented my own.

  LHA := FLT.Fsub(FLT.Fadd(GMST,LongDecimal),RaDecimal)


   IF LHA < 0.0
     repeat until (LHA > 0.0) OR (LHA == 0.0)
       LHA := FLT.Fadd(LHA,360.0)
       
   IF LHA > 360.0
     repeat until (LHA < 360.0) OR (LHA == 360.0)
       LHA := FLT.Fsub(LHA,360.0)

Mark_T · 2013-10-13 12:09

Why mention language semantics? I'm talking about IEEE 754 floating point specification, where the behaviour
is specified as far as I remember to the bit. The only thing that might be causing variation is the rounding mode
which affects various operations, but in a precisely determined way. I don't have a copy of the spec to hand so
can't lookup Fmod.

Duane Degn · 2013-10-13 12:32

Mike, Do you know about F32? It's an improved floating point object. I think it has all the methods Float32Full has but only uses one cog. I also think Jonathan managed to increase the speed of many of the methods.

I just check and the OBEX version hasn't been updated with the latest bug fix. The latest version can be found in post #98 of Jonathan's thread.

Heater. · 2013-10-13 12:38

Mark_T,

Why mention language semantics? I'm talking about IEEE 754 floating point specification

I agree with you. I would have naively assumed that on modern computers such operations were passed on to a floating point unit and the language would accept whatever came back as it is, which should be IEEE 754 compliant. Or at least the languages would perform IEEE 754 in software. Seems this is not true.

So what does IEE 754 say about Floating point remainder?

So far I have found from wikipedia: "This is not like a normal modulo operation, it can be negative for two positive numbers. It returns the exact value of x(round(x/y)·y)". That sounds even worse than I would expect!

From here: http://msdn.microsoft.com/en-us/library/system.math.ieeeremainder.aspx I find that the IEEE floating point remainder is defined as:

IEEERemainder = dividend - (divisor * Math.Round(dividend / divisor))

where as modulus is defined as:

Modulus = (Math.Abs(dividend) - (Math.Abs(divisor) * 
          (Math.Floor(Math.Abs(dividend) / Math.Abs(divisor))))) * 
          Math.Sign(dividend)

I can't even find that IEEE has a modulus operator. And the remander operator it does define may not do what you would like, for example:

Math.IEEERemainder(3, 2) = -1

tosjduenfs · 2013-10-13 20:47

Duane Degn wrote: »

Mike, Do you know about F32? It's an improved floating point object. I think it has all the methods Float32Full has but only uses one cog. I also think Jonathan managed to increase the speed of many of the methods.

I just check and the OBEX version hasn't been updated with the latest bug fix. The latest version can be found in post #98 of Jonathan's thread.

Mark, I've been using the F32 object I found on the obex, I'll update to the new version to see if it give me the result I need rather than using the loops. Thanks

Heater. · 2013-10-14 01:52

re: modulus in IEE 754

I just found the following note about a certain _frem function in some ARM libraries:

Functions that perform the IEEE 754 remainder operation. This is defined to take two numbers, x and y, and return a number z so that z = x - n * y, where n is an integer. To return an exactly correct result, n is chosen so that z is no bigger than half of x (so that z might be negative even if both x and y are positive). The IEEE 754 remainder function is not the same as the operation performed by the C library function fmod, where z always has the same sign as x. Where the IEEE 754 specification gives two acceptable choices of n, the even one is chosen. This behavior occurs independently of the current rounding mode.

So there we have it.

Should Float32 do what most languages modulus operator does or it should it do what IEEE does? Should it provide both?

tosjduenfs · 2013-10-14 17:48

It would be nice if you could specify the behavior of the modulo operation.

Heater. · 2013-10-15 01:16

The net is awash with such specifications of modulo, modulus and remainder.

Have a google for fmod (as used in C) and IEEremainder. I'm sure you can even find the actual 754 specification somewhere.

Strange behavior of Fmod operation in Float32Full

Comments