Quadruple Precision in C++ (Gcc)

Quadruple Precision in C++ (GCC)

Apparently, this seems to have been an installation error on my part.

While the core C/C++ portion of the GCC includes libquadmath.so, the Fortran version supplies libquadmath.a and quadmath.h, which can be included to access the functions.

#include <quadmath.h>
#include <iostream>
int main()
{
  char* y = new char[1000];
  quadmath_snprintf(y, 1000, "%Qf", 1.0q);
  std::cout << y << std::endl;
  return 0;
}

Epsilon in quadruple precision (gcc)

long double is not guaranteed to be implemented as IEEE-745 quadruple precision. C++ reference reads:

long double - extended precision floating point type. Does not necessarily map to types mandated by IEEE-754. Usually 80-bit x87 floating point type on x86 and x86-64 architectures.

If long double is implemented as 80-bits x86 extended precision, then epsilon is 2-63 = 1.0842e-19. This is the value you get as the output.

Some compilers support __float128 type that has quadruple precision. In GCC long double becomes an alias for __float128 if -mlong-double-128 command line option is used, and on x86_64 targets __float128 is guaranteed to be IEEE quadruple precision type (implemented in software).

std::numeric_limits is not specialized for __float128. To get the value of epsilon the following trick can be used (assuming a little-endian machine):

__float128 f1 = 1, f2 = 1;      // 1.q       -> ...00000000
std::uint8_t u = 1;
std::memcpy(&f2, &u, 1);        // 1.q + eps -> ...00000001
std::cout << double(f2 - f1);   // Output: 1.9259e-34

With GCC you can use libquadmath:

#include <quadmath.h>
...

std::cout << (double)FLT128_EPSILON;

to get the same output.

Converting a working code from double-precision to quadruple-precision: How to read quadruple-precision numbers in FORTRAN from an input file

Almost all of this is already in comments by @IanH and @Vladimi.

I suggest mixing in a little Fortran 90 into your FORTRAN 77 code.

Write all of your numbers with "E". Change your other program to write the data this way. Don't bother with "D" and don't try to use the infrequently supported "Q". (Using "Q" in constants in source code is an extension of gfortran -- see 6.1.8 in manual.)

Since you want the same source code to support two precisions, at the top of the program, have:

use ISO_FORTRAN_ENV
WP = real128

or

use ISO_FORTRAN_ENV
WP = real64

as the variation that changes whether your code is using double or quadruple precision. This is using the ISO Fortran Environment to select the types by their number of bits. (use needs to between program and implicit none; the assignment statement after implicit none.)

Then declare your real variables via:

real (WP) :: MyVar

In source code, write real constants as 1.23456789012345E+12_WP. The _type is the Fortran 90 way of specifying the type of a constant. This way you can go back and forth between double and quadruple precision by only changing the single line defining WP

WP == Working Precision.

Just use "E" in input files. Fortran will read according to the type of the variable.

Why not write a tiny test program to try it out?

---

Related Topics