0

I want to find the intersection point of two lines and then check if that point is equal to a predefined value (here I take it the exact value).

Problem schematic

This is my minimum working example,

PROGRAM foo
    IMPLICIT NONE

    ! Real kind number
    INTEGER, PARAMETER :: RKN = 8

    ! Parameters of line formulae
    REAL(RKN) :: a1, b1, c1, a2, b2, c2

    ! Exact values
    REAL(RKN), DIMENSION(2) :: ans

    ! Points of first line
    REAL(RKN), DIMENSION(2) :: pnt_a, pnt_b

    ! Slope of first line
    REAL(RKN) :: m

    ! Intersection
    REAL(RKN), DIMENSION(2) :: x

    ! Output format
    CHARACTER(*), PARAMETER :: fmt = '(X, A, X, GO)'

    pnt_a = [- 4._RKN, - 1._RKN]
    pnt_b = [- 1._RKN ,  0._RKN]

    m = (pnt_b(2) - pnt_a(2)) / (pnt_b(1) - pnt_a(1))

    a1 = m
    b1 = - 1._RKN
    c1 = m * pnt_a(1) - pnt_a(2)

    ! Second line is horizontal and has no vertical intercept
    a2 =   0._RKN
    b2 =   1._RKN
    c2 =   0._RKN

    ! Cramer's rule
    x = [(c1 * b2 - b1 * c2) / (a1 * b2 - b1 * a2), (a1 * c2 - c1 * a2) / (a1 * b2 - b1 * a2)]

    ans = [- 1._RKN, 0._RKN]

    WRITE(*, fmt) 'x(1)   = ', x(1)
    WRITE(*, fmt) 'ans(1) = ', ans(1)
    PRINT *, x(1) == ans(1)
    PRINT *, ABS(x(1) - ans(1)) < EPSILON(1._RKN)

END PROGRAM

After compilation, output is,

x(1)   =     -.9999999999999998
ans(1) =     -1.000000000000000
F
F

Though reals are double precision, the output is not close enough to exact answer, even using EPSILON()intrinsic function. If RKN is set to 16, it produces,

x(1)   =    -1.00000000000000000000000000000000
ans(1) =    -1.00000000000000000000000000000000
F
F

How can I handle this type of problems?

Shaqpad
  • 107
  • 2
  • 13
  • 1
    What is the problem? It really isn't clear. What output were you expecting? – John Coleman Mar 26 '17 at 14:22
  • I expect an answer close enough to exact value. – Shaqpad Mar 26 '17 at 14:28
  • 2
    If your problem is that carrying out a comparison on two `REAL` numbers is not finding them equal when you want them to, define your own "equality" condition. I.e. set a difference that you consider to be sufficiently small that any differences are just down to floating point rounding errors in the calculation. – SteveES Mar 26 '17 at 14:41
  • @SteveES Do you mean using `EPSILON()` function? – Shaqpad Mar 26 '17 at 14:43
  • Using `G0` for output needn't give you a representation of the entire bit pattern of the output item. So, just because you're seeing something output as `-1` it really doesn't mean that it is. – francescalus Mar 26 '17 at 14:44
  • 2
    You are doing floating point division. That introduces round-off errors that can't be wished away. If it really is an issue, you could use exact rational arithmetic (at the cost of many more cpu cycles). – John Coleman Mar 26 '17 at 14:46
  • @francescalus Yes of course. I am aware of that. – Shaqpad Mar 26 '17 at 15:04
  • 1
    you might use quad precision for the calculation, convert the result to double, then do your comparison. Then i think at least your ` – agentp Mar 26 '17 at 15:33
  • @agentp I added a few lines to top of the question. Actually, I am trying to check if the intersection point is close enough to a predefined value (say, the exact value). The problem is simple, but the difference between numerical and exact values is greater than `EPSILON()`. – Shaqpad Mar 26 '17 at 15:45
  • Your question seems to be more about **precision in formatting or printing the results**, rather than in the actual calculation itself. See existing duplicates. – smci Mar 26 '17 at 16:14
  • 1
    @smci My take is that the question isn't about formatting, but is questioning why the equality conditions are evaluating to `False`. – SteveES Mar 26 '17 at 16:33
  • 1
    **[Is the use of machine epsilon appropriate for floating-point equality tests?](http://stackoverflow.com/questions/3281237/is-the-use-of-machine-epsilon-appropriate-for-floating-point-equality-tests)** I still say the OP wants to round their result, rather than worry about getting exact float/double comparisons to work. – smci Mar 26 '17 at 16:36
  • @smci, I never suggested it wasn't a duplicate ;) – SteveES Mar 26 '17 at 16:40
  • 1
    @smci, that was essentially what I was getting at - using the `EPSILON()`-style comparison, but setting a more appropriate value (that you control) in place of `EPSILON()`. Rounding is the other way of achieving this. – SteveES Mar 26 '17 at 16:43
  • @SteveES: ok cool. Well we've triaged the OP's confusion. Other people can clean up the wording, find a canonical question and close-as-duplicate. agentp's suggestion to internally use quad-precision just to make sure that a double-precision result passes a comparison using epsilon, is insane overkill! Round the results when displaying them, is all. – smci Mar 26 '17 at 16:47
  • if you need the result to be good to machineepsilon at some precision then you need to do the calculation at higher precision. It is of course over kill if you are just being retentive about using mechineepsilon as your tolerance parameter. – agentp Mar 26 '17 at 18:37
  • This question never needed the result to be accurate to machineepsilon and that's been obvious all along... and you can always round the result before displaying it, or use a format spec. Every language has tons of duplicates on these, which I referenced. – smci Mar 26 '17 at 22:14
  • @smci Thanks for your comments. That machine epsilon accuracy is not needed here may be correct, but my question is not about displaying or formatting results. – Shaqpad Mar 26 '17 at 22:21

0 Answers0