modpi function

[StefanKarpinski]

The naive definition modpi(x) = mod(x,pi) is wildly off for large values of x. We should have version of modpi that computes this correctly. cc: @alanedelman.

[andrioni]

Is this solvable without resorting to BigFloat?

[andrioni]

This SO discussion seems interesting.

[stevengj]

My understanding is that it is solvable without resorting to arbitrary-precision arithmetic. The term to google for is "range reduction", because this problem most commonly shows up in range reduction of trigonometric functions to reduce their argument to be in [0,π/4]. See e.g. this paper.

[vtjnash]

openlibm includes __ieee754_rem_pio2 which I think does this. it might just need to be exported?

[pao]

I already mentioned this where Stefan first brought this up, but modtau would also be useful for a heck of a lot of dynamics problems, too.

[fabianlischka]

See also JuliaMath/openlibm#30

[fabianlischka]

I'm working on exposing ieee754_rem_pio2, modpi, mod2pi, and modpio2 = mod(_,pi/2). Just finishing up the testing.

[fabianlischka]

Incidentally, does anyone know the difference between k_rem_pio2.c and e_rem_pio2.c?
(see
https://github.com/JuliaLang/openlibm/blob/master/src/k_rem_pio2.c?source=cc and https://github.com/JuliaLang/openlibm/blob/master/src/e_rem_pio2.c?source=c )

Thanks!

[JeffBezanson]

k_rem_pio2 seems to implement a lower-level kernel function called by e_rem_pio2. e_rem_pio2 I suspect has the more useful entry point.

[fabianlischka]

Thanks. That's the one I've been using. I'll try the whole clone/push/submit pull request thing soon-ish (first time, yay), though I must first get Julia made on Mavericks (doesn't work for me yet :-/ )

[pao]

Do see the comments on the OS 10.9 issue; there's a few things worth trying in there: #4614. Also this more recent issue; I think you need to be using bash right now: #4776.

[fabianlischka]

Thanks pao, it's working now. I needed to a) run "xcode-select --install" to update the Xcode command line tools, and b) put the julia directory in the PATH, rather than just put a soft link to the julia executable from /usr/local/bin. I'm sort of done with testing as well, and will clone now.

One problem: my modpi will currently take values (such as 0.7853881633974482) that should remain un-modified and change the last bit (because it internally devolves to a mod pi/2 routine and then adds the correct multiple of pi/2). I'll put in in a check to catch these cases and just return the number.

However, this highlights the fact that the implementers of ieee754_rem_pio2 really maximised precision (retaining one more bit) by returning the result in the range [-pi/4, pi/4] rather than [0, pi/2]... :-)

[fabianlischka]

see here for a summary: http://nbviewer.ipython.org/7443293

Some questions:

questions:

• are the pi-related constants defined elsewhere already?
• can the pi-related constants defintions be wrapped up in a macro? Should they?
• should mod(x,pi) redirect to modpi(x), or issue a warning, or do nothing special?

After
julia> mod(x::FloatingPoint,y::MathConst{:π}) = modpi(x)

we have:

julia> mod(5706674932067741.0,pi) # correct, (new) modpi called
4.237546464512562e-16

julia> mod(5706674932067741,pi) # first arg int: original "mod" called
0.2224559947753093

julia> mod(5706674932067741,pi*1) # second arg Float64: original "mod" called
0.2224559947753093

julia> mod(5706674932067741.0,pi*1) # second arg Float64: original "mod" called
0.2224559947753093

• should all the tests be included in the test dir, or just a selection, or none?
• currently, mod2pi doesn't work for integer arguments:
  julia> mod2pi(123)
  ERROR: no method mod2pi(Int64,)

Could introduce a method that converts ints to floats. The problem being that ints > 2^53 might be truncated, and thus the modpi result be absolutely wrong - should one issue a warning, or throw an error, or silently ignore it?

[JeffBezanson]

I think all we should do is add the modpi functions for float arguments,
since that is all that is actually supported in a meaningful way.
On Nov 13, 2013 12:44 PM, "fabianlischka" notifications@github.com wrote:

see here for a summary: http://nbviewer.ipython.org/7443293

Some questions:

questions:

• are the pi-related constants defined elsewhere already?
• can the pi-related constants defintions be wrapped up in a macro?
  Should they?
• should mod(x,pi) redirect to modpi(x), or issue a warning, or do
  nothing special?

After
julia> mod(x::FloatingPoint,y::MathConst{:ð}) = modpi(x)

we have:

julia> mod(5706674932067741.0,pi) # correct, (new) modpi called
4.237546464512562e-16

julia> mod(5706674932067741,pi) # first arg int: original "mod" called
0.2224559947753093

julia> mod(5706674932067741,pi*1) # second arg Float64: original "mod"
called
0.2224559947753093

julia> mod(5706674932067741.0,pi*1) # second arg Float64: original "mod"
called
0.2224559947753093

should all the tests be included in the test dir, or just a selection,
or none?

currently, mod2pi doesn't work for integer arguments:
julia> mod2pi(123)
ERROR: no method mod2pi(Int64,)

Could introduce a method that converts ints to floats. The problem being
that ints > 2^53 might be truncated, and thus the modpi result be
absolutely wrong - should one issue a warning, or throw an error, or
silently ignore it?

Reply to this email directly or view it on GitHubhttps://github.com//issues/2524#issuecomment-28416433
.

[fabianlischka]

So, if we do NOT add the mod method for Pi:
mod(x::FloatingPoint,y::MathConst{:π}) = modpi(x)
we would have:
mod(100,pi) == mod(100.0,pi) != modpi(100.0) # modpi(100) not defined
If we add it, we would have:
mod(100,pi) != mod(100.0,pi) == modpi(100.0)

The value of mod(100.0, pi) would move in the "correct" direction. But it might be confusing that it matters whether the first arg is a float or an int....

[ivarne]

Why can't you raise InexactError() if the int can't be converted lossless to a float and give the improved results for ints also?

[StefanKarpinski]

Thinking about this again, InexactError isn't bad. Maybe that's the best approach. I looked into if MPFR has a modpi-like function and they don't seem to.

[fabianlischka]

What about this:

modpi(x::Int32)  = modpi(float64(x)) 
modpi(x::Int64)  = if int(float(x))==x modpi(float64(x)) else error("Integer argument to modpi() is too large.") end

?

[StefanKarpinski]

That seems reasonable. I would write it like this:

modpi(x::Int32) = float32(modpi(float64(x)))

function modpi(x::Int64)
  fx = float64(x)
  fx == x || error("argument to modpi is too large: $x")
  modpi(fx)
end

[ivarne]

@StefanKarpinski Why would you reduce the precision on 32 bit systems? If someone wants a float32, he can convert it himself?

[vtjnash]

Stefan, Why not use the two argument version of assert (@asset x "msg") instead of the x || error("msg") form?

[stevengj]

@vtjnash, the usual definition of assertion is restricted to predicates that can never fail unless there is a bug in the software; it seems like bad style to use it to throw ordinary exceptions.

[ivarne]

Sorry for the negative feedback but error() does not raise a InexactError

[kmsquire]

Why would you reduce the precision on 32 bit systems? If someone wants a float32, he can convert it himself?

@invarne, two things:

1. there's no reduction in precision on 32-bit systems, only with 32-bit floats--all floats are Float64 by default, even on 32-bit systems
2. for type stability, if someone is already calculating with 32-bit floats, suddenly switching to 64-bit in the middle of a calculation is not really a good idea

[ivarne]

@kmsquire I was talking about the conversion to float32 when modpi is used with int32 argument. I think that the most probable reason someone use a int32 is because they are on a 32 bit system. Because 32 bit systems use float64, by default i think that modpi should return that. Sorry for the confusion.

I totally agree that modpi should return float32 when the argument is float32

[kmsquire]

@ivarne, my bad, I hadn't noticed the int32 in the signature. I agree with you--it probably should be Float64 by default.

[StefanKarpinski]

Sure, I'm not really particularly stuck on that. I believe that we generally convert Int32 values to Float32s but this case is certainly debatable.

[fabianlischka]

yes, so currently there is no method modpi() for Float32, only for Float64.

The int cases are handled differently because we know that every Int32 can be represented as a Float64 exactly, while with large Int64 we might be off by eps=1 or 2 or more, and then modpi becomes quite meaningless (though, come to think of it, one could just add the difference between x::Int64 and fx = int(float(x)) to the result of modpi(float(x)) to get the correct result (might need to modpi it once more).... Might be a bit overkill though).

@ivarne, no worries, my first pull request, happy to get feedback. How would I raise InexactError?

[StefanKarpinski]

The easiest thing to do for Float32 is just modpi(x::Float32) = float32(modpi(float64(x)).

[StefanKarpinski]

I guess we return Float64 for sin(::Int32) and such, so clearly my behavior above of returning Float32 for modpi(::Int32) is not consistent. Please ignore that bit.

[ivarne]

throw(InexactError())

Unfortunatly InexactError does not have a field for message, so if you want to help those who hit the error, you need to do the https://github.com/JuliaLang/julia/blob/master/base/repl.jl#L75 trick, or convince Stefan/Jeff to add messages to all exceptions.

[fabianlischka]

sure, will do.
BTW, this whole "x==y || error()", or "some_condition && return 5" is really preferred to if statements? I'm surprised, it strikes me as less legible, at least initially. (On the other hand, given that there is no "then" or ":" separating the condition from the "then" part, it might actually be easier to read)

[StefanKarpinski]

I write it that way, others don't. It's a matter of taste. I think it's silly for an error check to take up multiple lines. I certainly think that writing x == y || error("x not equal to y") is clearer than cramming an if-else onto a single line.

[quinnj]

+1 for || error syntax. One of my favorite tricks I've picked up from Stefan.

[fabianlischka]

@ivarne - for now, I've left the "normal" error, with the informative error msg. Not sure how to do the InexactError with informative error msg. Open for suggestions, obviously.
For integers, we return Float64 (and Int64 inputs that can't be represented as a Float64 throw an error). For Float32, we return Float32.
I've added tests, see the pull request #4799
Anything else?