From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugzilla@gcc.gnu.org>
Received: by sourceware.org (Postfix, from userid 48)
	id D0C08384D0F0; Sat, 19 Nov 2022 20:14:01 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org D0C08384D0F0
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org;
	s=default; t=1668888841;
	bh=uVfw0mQhNKjdyizpu3SZBd+xvyfwT0joFbyXu0LChus=;
	h=From:To:Subject:Date:In-Reply-To:References:From;
	b=TfHdlvntPg9+/HDqw6U6yPC3W1de8WBcxtvq7ldkomQWYnt8lsFmqIvePWPr3BZak
	 YOJ9w7ITZVbdC0HnSw2ldLqp2/rND3n1P4mb1VeDO3A5C6Zu4SVcuA7emkQ0YPRShw
	 RLDYl0X3nDYUGAowk+0PR7AfmffDShXQX937CeWo=
From: "anlauf at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>
To: gcc-bugs@gcc.gnu.org
Subject: [Bug fortran/107753] gfortran returns NaN in complex divisions
 (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)
Date: Sat, 19 Nov 2022 20:14:01 +0000
X-Bugzilla-Reason: CC
X-Bugzilla-Type: changed
X-Bugzilla-Watch-Reason: None
X-Bugzilla-Product: gcc
X-Bugzilla-Component: fortran
X-Bugzilla-Version: 12.2.0
X-Bugzilla-Keywords: 
X-Bugzilla-Severity: normal
X-Bugzilla-Who: anlauf at gcc dot gnu.org
X-Bugzilla-Status: WAITING
X-Bugzilla-Resolution: 
X-Bugzilla-Priority: P4
X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org
X-Bugzilla-Target-Milestone: ---
X-Bugzilla-Flags: 
X-Bugzilla-Changed-Fields: 
Message-ID: <bug-107753-4-Lvw8witxJs@http.gcc.gnu.org/bugzilla/>
In-Reply-To: <bug-107753-4@http.gcc.gnu.org/bugzilla/>
References: <bug-107753-4@http.gcc.gnu.org/bugzilla/>
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/
Auto-Submitted: auto-generated
MIME-Version: 1.0
List-Id: <gcc-bugs.sourceware.org>

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D107753

--- Comment #11 from anlauf at gcc dot gnu.org ---
(In reply to Weslley da Silva Pereira from comment #7)
> More data for the discussion:
> 1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization
> flags `-O` but still reproduce the wrong result for complex divisions with
> huge numbers. See

It is possible that gfortran's dependence on optimization level depends
on the version.  If one wants to test run-time behavior and avoid
compile-time simplification, it may be helpful to add:

  volatile :: x, y, z

I then get consistent results for -O0 / -O1.

> 4. My Ubuntu 20.04.5 LTS with compiler ifort 2021.7.1 computes the complex
> division `x/x` accurately even for the case of huge numbers. Scenarios
> tested:
>    - I tested the program in
> https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/
> test_zcomplexdiv.f and the one in https://godbolt.org/z/b3WKWodvn.
>    - I tested ifort with flags -fp-model precise and -fp-model fast. The
> latter enables more aggressive optimizations on floating-point data.
>    - I tested compilation with optimization flags -O0, -O, -O1, -O2, -O3.=
=20

Intel might be fine, but at least some current llvm-based compilers (Nvidia,
AMD flang) show more or less similar behavior to gfortran.

E.g. nvfortran 22.11:

 (1.7976931348623157E+308,1.7976931348623157E+308)
 (8.9884656743115795E+307,8.9884656743115795E+307)
 (4.4942328371557898E+307,4.4942328371557898E+307)
 (NaN,0.000000000000000)
 (NaN,0.000000000000000)
 (1.000000000000000,0.000000000000000)

As a sidenote: we are really discussing borderline cases here, valid
but only rarely occuring in normal code execution.  If I replace

  x =3D cmplx( huge(0.0d0), huge(0.0d0), dp )
  y =3D cmplx( b**(E-1), b**(E-1), dp )

by

  x =3D cmplx( nearest(huge(0.0d0),-1.d0), nearest(huge(0.0d0),-1.d0), dp )
  y =3D cmplx( nearest(b**(E-1),   -1.d0), nearest(b**(E-1),   -1.d0), dp )

then I get

               (1.0000000000000000,0.0000000000000000)
               (1.0000000000000000,0.0000000000000000)

instead of NaN.=