Inconsistent value of Pi in the source

[serguei-patchkovskii]

Hi,

I am trying to track down a numerical inconsistency in a real16 (quadruple-precision) build of LAPACK, and I noticed that 2Pi and Pi/2 is hard-coded in multiple locations throughout the source. The two values I see are:

PIOVER2, which is set to 1.57079632679489662 (either D0 or E0, depending on the precision) and

TWOPI, which is set to 6.2831853071795864769252867663 (again either D0 or E0). Parenthetically, this value is off in the last decimal digit - it should have been ...7666 if properly rounded.

These values are sufficiently accurate and consistent for real4 and real8 math. The value of PIOVER2 is already not quite accurate enough for real10 math, and neither value is sufficiently accurate for real16, which requires about 34 decimal digits.

The fix is to set

PIOVER2 = 1.570796326794896619231321691639751442099

TWOPI = 6.283185307179586476925286766559005768394

These values are given to 40 significant digits, which should be more than adequate for real16 and complex32 use.

The PIOVER2 value is found in:

./SRC/cbbcsd.f: $ ONE = 1.0E0, PIOVER2 = 1.57079632679489662E0,
./SRC/dbbcsd.f: $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
./SRC/sbbcsd.f: $ ONE = 1.0E0, PIOVER2 = 1.57079632679489662E0,
./SRC/zbbcsd.f: $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
./TESTING/EIG/cckcsd.f: $ PIOVER2 = 1.57079632679489662E0,
./TESTING/EIG/ccsdts.f: PARAMETER ( PIOVER2 = 1.57079632679489662E0,
./TESTING/EIG/dckcsd.f: $ PIOVER2 = 1.57079632679489662D0,
./TESTING/EIG/dcsdts.f: PARAMETER ( PIOVER2 = 1.57079632679489662D0,
./TESTING/EIG/sckcsd.f: $ PIOVER2 = 1.57079632679489662E0,
./TESTING/EIG/scsdts.f: PARAMETER ( PIOVER2 = 1.57079632679489662E0,
./TESTING/EIG/zckcsd.f: $ PIOVER2 = 1.57079632679489662D0,
./TESTING/EIG/zcsdts.f: PARAMETER ( PIOVER2 = 1.57079632679489662D0,

The TWOPI value is found in:

./SRC/clarnv.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./SRC/dlarnv.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./SRC/slarnv.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./SRC/zlarnv.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/clarnd.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/clatms.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/clatmt.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/dlarnd.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/dlatms.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/dlatmt.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/slarnd.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/slatms.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/slatmt.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
./TESTING/MATGEN/zlarnd.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/zlatms.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )
./TESTING/MATGEN/zlatmt.f: PARAMETER ( TWOPI = 6.2831853071795864769252867663D+0 )

[mu578]

Hello. you have C double exported constants here: https://github.com/moe123/macadam/blob/master/macadam/mcconsts.h and in your system <math.h> you may relay on them. You can also export real16 or new ones using https://www.wolframalpha.com. However ...7666 vs ...7663 won't change anything it will be rounded as normalized by rt.

[mu578]

Hello, yes I understand, meanwhile, I am not confident with 128-bit defs and its support. After checking in and out those concerns, the most straightforward step I see, would be to use a compiler implementing or emulating such binary storage. Then computing estimate of each values needed and dumping them to their respective hexadecimal true representation. I am not aware of what fortran today support for const declaration. Anyhow it departs from f77 compatibility forcibly.

[mu578]

"hexadecimal true representation" will give you a true normalized representation. There is nothing conflated here. AIX is an OS what's important is 128-bit floating-point hardware support or emulation. You know, I am a pretty regular guy; I like synthesized chitchat.

[zerothi]

I think there is no reason not to have the precision increased in case users ask for it. @serguei-patchkovskii could you make a PR for this?

[langou]

Thanks Nick. Thanks Sergei. This is merged.