3

I want to run a particle simulation with periodic boundary conditions - for simplicity, let's assume a 1D simulation with a region of length 1.0. I could enforce these conditions using the following short snippet:

if (x > 1.0)
    x -= 1.0;
if (x < 0.0)
    x += 1.0;

but it feels "clumsy"1 - especially when generalizing this to higher dimensions. I tried doing something like

x = x % 1.0;

which takes good care of the case x > 1.0 but doesn't do what I want for x < 0.02. A few examples of the output of the "modulus" version and the "manual" version to show the difference:

Value: 1.896440,  mod: 0.896440,  manual: 0.896440
Value: -0.449115, mod: -0.449115, manual: 0.550885
Value: 1.355568,  mod: 0.355568,  manual: 0.355568
Value: -0.421918, mod: -0.421918, manual: 0.578082

1) For my current application, the "clumsy" way is probably good enough. But in my pursuit of becoming a better programmer, I'd still like to learn if there's a "better" (or at least better-looking...) way to do this.

2) Yes, I've read this post, so I know why it doesn't do what I want, and that it's not supposed to. My question is not about this, but rather about what to do instead.

Community
  • 1
  • 1
Tomas Aschan
  • 58,548
  • 56
  • 243
  • 402
  • Is there a reason why just putting the `if/else` in a method, possibly wrapped by a class with a field specifying the boundary if that's a parameter, isn't workable? – chrylis -cautiouslyoptimistic- Oct 30 '13 at 10:19
  • @chrylis: It's definitely a better approach - but it still feels a little clumsy. Also, although I didn't specifically need it this time, it is nice to know a branch-less way of doing this for high-performance applications. – Tomas Aschan Oct 30 '13 at 10:22

2 Answers2

4

You can use % with this slight modification x = (x + 1.0) % 1.0

Ivaylo Strandjev
  • 69,226
  • 18
  • 123
  • 176
0

The best approach is probably to subtract the floor of the value from the value itself. Computing a floating-point remainder in compliance with IEEE standards is rather complicated, and unless one is on a system which can detect and accelerate "easy cases", especially those where the divisor is a power of two, a pair of if statements is apt to be faster.

It might be interesting, though, to consider why fmod was designed the annoying way it was: if fmod were designed to return a value between 0 and the dividend, then the precision of the result when the dividend is a very small positive number would be much better than the precision when the dividend is a very small negative number (the precision would be limited to that of the divisor). The advantages of having fmod's precision be relatively symmetric about zero probably outweighs the advantages of having the results be non-negative, but that doesn't imply the IEEE is the only good way to design a range-limiting function.

An alternative approach which would combine the advantages of the IEEE's approach and the zero-to-divisor approach would be to specify that a mod function must yield a result whose numerical value was (for d positive) less than the numerical value of d/2, but no less than -d/2. Such a definition would always yield a result that was representable in the source operands' type (if D is a very small value such that D/2 is not precisely representable, the range of the modulus would be symmetrical). Unfortunately, I know of no library mod functions that work this way.

supercat
  • 77,689
  • 9
  • 166
  • 211