Fortran round-off errors

Question

I have simple code, which flags nodes with in region enclosed by cylinder. On implementing the code, the result is mild tilt of the cylinder observed case with 90 degrees

The actual issue: The above algorithm is implemented in Fortran. The code checks for points in Cartesian grid if inside the cylinder. Following being the test case: The cylinder makes an angle 90 degrees in the yz-plane with respect to y-axis. Therefore, the orientation vector $\vec{o}$ is (0, 1, 0).

Case 1: Orientation vector is assigned directly with $\vec{o}=(0.0,1.0,0.0)$. This results in perfect cylinder with $\theta=90.$

Case 2: Orientation vector is specified with intrinsic Fortran functions with double precision accuracy dsin and dcos with $\vec{o}=(0.0, \sin(\pi/2.0), \cos(\pi/2.0))$ with $\pi$ value assigned with more than 20 significant decimal points. The resulting cylinder results in mild tilt.

The highlighted region indicates the extra material due to tilt of cylinder with respect to Cartesian axes. I also tried architecture specific maximum precision "pi" value. This also results in same problem.

This shows like the actual angle made by cylinder is not 90 degrees. Can anyone suggest valid solution for this problem. I need to use the inbuilt trigonometric functions for arbitrary angles and looking for accurate cell flagging method.

Note: All operations are performed with double precision accuracy.

The actual function is below. rk is defined parameter with value 8

  pure logical function in_particle(p,px,x)
  type(md_particle_type),intent(in) :: p
  real(kind=rk),intent(in) :: px(3),x(3)
  real(kind=rk) :: r(3),rho(3),rop(2),ro2,rdiff,u
  rop = particle_radii(p)  ! (/R_orth,R_para/)
  ro2 = rop(1)**2
  rdiff = rop(2) - rop(1)

  r = x-px

! Case 1:
!  u = dot_product((/0.0_rk,-1.0_rk,0.0_rk/),r)
!  rho = r-u*(/0.0_rk,-1.0_rk,0.0_rk/)

! Case 2:
  u = dot_product((/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/),r)
  rho = r-u*(/0.0_rk,-dsin(pi/2.0_rk),dcos(pi/2.0_rk)/)

  if((u.le.rdiff).and.(u.ge.-rdiff)) then
    in_particle = dot_product(rho,rho) < ro2
  else
    in_particle = .false.
  end if
end function in_particle

Note: The trigonometric operations are done inside the code to explain the problem better. However the original code reads the orientation in vector form from user. Then converts this information to quaternions for particle-particle collision operations. On converting quaternions back to orientation vector, this error is even more amplified. Even before the start of collision, the orientation of cylinder tends to be disoriented by 2 lattice cells.

Do you mean that if you comment the two lines below "Case 2" and uncomment the two lines below "Case 1", you get the expected result, while in the above code you get an unexpected result? What do you get if you insert a line like `print *, dsin(pi/2.0_rk); stop` (to verify that the value of pi is OK)? — roygvib, Jun 18 '16 at 23:46
Just a style point - dsin and dcos is very Fortran 66. Simply sin and cos is perfectly OK, and has been for almost 40 years. — Ian Bush, Jun 19 '16 at 07:06
Cross-posting -- http://scicomp.stackexchange.com/questions/24247/fortran-round-off-error-with-floating-point-operations -- is rather frowned upon on Stack Exchange. — High Performance Mark, Jun 19 '16 at 09:59
@HighPerformanceMark: I posted the question there first as its more relevant there. But later realized that this community could have better reach and more solutions. If required, I can delete the question there. — SKPS, Jun 19 '16 at 14:45
@roygvib: Yes. If I uncomment lines below "Case 1" and comment "Case 2" lines, I get the expected result. I printed the values of the angle, only the numbers after 20th significant digit is little different. — SKPS, Jun 20 '16 at 08:06
One more confirmation: Do you get cos(pi/2.0_rk) = 0.00000000000000006123... ? (the number of zero is 17.) If so, I have no clue at the moment for those two lines... As for other lines, I'm wondering if rdiff is guaranteed to be positive (because the IF conditional seems to assume it). — roygvib, Jun 20 '16 at 08:43
@roygvib: rdiff is always positive. And I am also getting cos(pi/2.0_rk) = 0.00000000000000006123. Can such a small deviation from true value cause such strong disorientation? — SKPS, Jun 20 '16 at 09:00
No, I don't think so (just wanted to confirm that pi retains double precision). The projection itself looks stable. There might be some corner cases where in_particle changes true/false because of the slight numerical errors, but I don't think it would lead to systematic tilt. There might be some other problems but no idea at the moment...(sorry!). — roygvib, Jun 20 '16 at 09:43
aside from the value of `pi`, it seems your error should be uniform over the `x` coordinate because you are taking the dot product against `[0,-sin,cos]`. A precision error in `pi` should cause a uniform distortion of the cross section, not a "tilt". — agentp, Jun 20 '16 at 18:05
Though a wild guess, giving some slight buffer region to the IF conditional, say, `- ( rdiff + eps ) < u .and. u < ( rdiff + eps )` with eps = 1.0d-3 etc might help...? (to avoid possible cases of particles exactly located on the boundary). — roygvib, Jun 23 '16 at 04:51

Sneftel · Answer 1 · 2016-06-27T19:17:56.793

cos(pi/2) is not necessarily going to give you exactly 0, no matter how exact you make the cos calculation, and no matter how many digits of pi you have, because:

pi, as an irrational number, will contain up to 1/2 ulp of error when represented as an FP number; and
sin and cos are not guaranteed by the IEEE-754 standard to be correctly rounded (or even implemented).

Now, sin(pi/2) is extremely likely to come out as 1 regardless of precision and FP architecture, simply because sin has such a low derivative around 1; with single-precision floats, it should come out to 1 if you're anywhere within about 3e-4 of the exact value of pi/2. The problematic call is the cos, which has lots of precision to play with around 0 and a derivative of about -1 in the neighborhood.

Still, we're talking about extremely small values here. I think what's really potentiating the problem here is the in/out test you're doing, combined with ordinary FP rounding rules. I would guess, in fact, that if you were to bias your test points by, say, a quarter of the grid quantum, you'd see all straight verticals in your voxelization (though it might not be symmetrical around the minor axes).

Another option would be to actually discard some precision from your sin/cos calculation before doing the dot product, effectively quantizing your axes.

score 3 · Answer 2 · answered Jun 23 '16 at 21:02

Short answer: Create a table of sin and cos of common angles (0, pi/6, pi/4, pi/3, pi/2, pi and their multiples) and compute only for uncommon angles. The reason being that errors with uncommon angles will be tolerated by most people while errors with common angles will likely not be tolerated.

Explanation: Because floating point computation is not exact (that is its nature), you sometime need a little bit of compromise between the accuracy and the readability of the code.

One way of doing that is to avoid to compute something that is known exactly. To do that, you can check the value of the angle and do the actual computation only if it is not an obvious angle. For example angle 0, 90, 180 and 270 degrees have obvious values of sin and cos. More generally, the cos and sin of common angles (0, pi/6, pi/4, pi/3, pi/2, pi and their multiples) are known exactly (even if they are irrational numbers).

Rather than just using a table, I'd suggest using a function that computes trig functions for arguments that the number of radians in a circle is the nearest floating-point value to 2pi, to accommodate the common usage cases where the argument is formed by multiplying a number by the fp representation of pi. — supercat, Mar 30 '18 at 00:50

Fortran round-off errors

2 Answers2