1

I'm using gfortran and I'm starting to write a new program that I know will have billions and billions of float equality checks (think "Monte-Carlo simulation"). And I need to know if two real numbers are "close enough".

It could be a simple number (a == b) or 2D/3D points, and even complex numbers. (I'll mostly do simple numbers, however.)

Are you aware of a gfortran intrinsic that could help me here, or a perhaps a function in BLAS/LAPACK that could do the trick. Or maybe some neat (possibly non portable) trick for blazing fast comparison? An epsilon of 0.0001 should be enough for me, but maybe there is a magic number that gfortran on x86_64 could optimize with bitwise operation.

I don't really know and google couldn't find a solution for me. Can you ?

If you're not aware of a neat trick, what do you think would be the fastest implementation of the traditional abs(a-b) <= epsilon? Or is it simply the fastest and there is nothing I can do about it? (luckily it's a "pure" function so it will probably be inlined)

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
ker2x
  • 440
  • 3
  • 11
  • 1
    dropping `=` should already be faster.. by how much? I'd recommend using a [profiler](https://stackoverflow.com/questions/18191635/good-profiler-for-fortran-and-mpi) but I would not expect big numbers – norok2 Sep 27 '19 at 17:47
  • interesting, i'll give it a try. i don't care if it's < or <=, so it's worth to give it a try. i'll go check the AMD developer manual (they're both asm instruction) and profile it. thx! – ker2x Sep 27 '19 at 17:59
  • 1
    actually, now that I think it better, it will not make any difference whatsoever. And probably there is no *magic blazing fast comparison trick*, if you think about what could be happening at gating level https://stackoverflow.com/questions/12135518/is-faster-than – norok2 Sep 27 '19 at 19:18
  • According to godbolt (in C with gcc9) one call seta (Set byte if above) and the other setnb (Set byte if not below). It's doing the exact opposite but it doesn't matter (it just have to test if the bit is not set instead of checling if the byte is set, and there you have it) – ker2x Sep 27 '19 at 19:53
  • the compiler is really smart btw. if you disassemble "return (fabsf(a - b) < 0.0000);" the asm code is "xor eax, eax; ret;" :D – ker2x Sep 27 '19 at 19:55
  • The compiler isn't necessarily really smart. The code is actually rather questionable. The absolute value of a number can never be less than zero, so for non-exceptional values the result is always false. Why the xor? Because NaN. – evets Sep 28 '19 at 06:11
  • No, the xor is there to zero the register. Nothing to do with NaN. Or well, a comparison with NaN is defined to always return false, and fabs(NaN) = NaN, so whether a or b is NaN doesn't have any effect on the result. – janneb Sep 28 '19 at 20:05

0 Answers0