https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96711
--- Comment #11 from Steve Kargl <sgk at troutmask dot apl.washington.edu> --- On Thu, Aug 20, 2020 at 01:47:58PM +0000, bre08 at eggen dot co.uk wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96711 > > --- Comment #10 from B Eggen <bre08 at eggen dot co.uk> --- > I've been experimenting with the suggested work-around > > m = anint(y) > > This works for larger numbers, even in quad precision, however, it breaks down > a long way before the integer*16 range is exhausted, consider the code below, > which starts with 2^113 and tries to double it, minus 1. The minus 1 is not > taking effect: > Of course, that is going to fail if you want an exact result. REAL(16) is a floating point format, which has 113 digits of precision. This means that integers less than 2**113 are exactly representable as REAL(16) floating point numbers. When you double (2._16)**113 to (2._16)**114, that is a prefectly fine REAL(16) floating point number. Subtracting 1._16 from (2._16)**114 is a valid floating point operation. Your problem lies in that result is rounded, and 1._16 is less than epsilon(1._16) in comparison to (2._16)**114. > I guess at some point NINT() will be fixed, can anyone suggest a robust > workaround that is valid until 2^127-1 ? If you're trying to use floating point arithmetic in computations and you need exact integral values up to 2**127-1, then there isn't a REAL type that works. If you're looking for arbitrary precision arithmetic and you want to use Fortran, I suggest that you google "David Bailey arithmetic".
