Hi!
On Wed, May 19, 2021 at 08:35:26PM +0200, Tobias Burnus wrote:
> On 19.05.21 17:15, Segher Boessenkool wrote:
> >> real(16) :: y ! 128bit REAL
> >> integer(16), parameter :: k2 = nint (2 / epsilon (y), kind(k2))
> >> integer(16), parameter :: m2 = 10384593717069655257060992658440192_16
> >> !2**113
> >> if (k2 /= m2) stop 3
> >>
> >>On x86_64-linux-gnu, k2 == m2 — but on powerpc64le-linux-gnu,
> >>k2 == 2**106 instead of 2**113.
> >>
> >>My solution is to permit also 2**106 besides 2**113.
> >I do not understand Fortran well enough, could you explain what the code
> >is supposed to do?
>
> First, 2_16 means the integer '2' of the integer kind '16', i.e. int128_t
> type.
Ah ok.
> And I think this testcase tries to ensure that the result of 'nint'
> both at compile time and at runtime matches what should be the result.
That is the main thing I missed :-)
> ('a**b' is the Fortran syntax for: 'a' raised to the power of 'b'.)
I know, I use it all the time :-) Many other languages have this
nowadays btw. It is such a succinct notation :-)
> nint(2/epsilon(y)). Here, 'epsilon' is the
> "Model number that is small compared to 1."
> Namely: b**(p-1) = '(radix)**(1-digits)'
> alias 'real_format *fmt = REAL_MODE_FORMAT (mode)'
> with radix = fmt->b and digits = fmt->p;
>
> [b**(p-1) is from the Fortran standard but 'b' and 'p' also match the
> ME/target names, while radix/digits matches the FE names and also the
> Fortran intrinsic inquiry function names.]
>
> This is for radix = 2 equivalent to:
>
> 2/2**(1-digits) = 2*2**(digits-1) = 2**(digits)
>
> On x86-64, digits == fmt->p == 113.
And the same on powerpc64le with -mabi=ieeelongdouble, presumably.
> Our powerpc64le gives digits == 106.
It stil defaults to -mabi=ibmlongdouble.
> Having written all this, I wonder why we don't just
> rely on the assumption that '2**digit(x)' works – and use this
> to generate the valid.
If that works (and I have no reason to believe it won't), that is as
simple a solution as you will find :-)
> As I like that patch and believe it is obvious, I intent to
> commit it as such – unless there are further comments.
> It passes on both x86-64-gnu-linux and powerpc64le-none-linux-gnu.
> I think the radix == 2 is a good bet, but if we ever run into issues,
> it can also be changed to use radix(...) as well ...
I don't think any GCC target does 10 or 16 by default, and less chance
in Fortran even :-)
Thanks!
Segher
