Fortran provides quad-precision floating point numbers in certain implementations. Declaring a number in the following manner, will create a variable with IEEE 128-bit quad precision (the 16 implies 16 bytes):
real (kind=16) :: quad_x real (kind=8) :: double_x real (kind=4) :: single_x
In order to see its affect, consider a piece of Fortran code which calculates π using the equation of Leibniz. Run with double precision provides the following output:
How many terms to calculate pi to? 10000 The value of pi is: 3.1416925847879611
Run with quad precision on the other hand:
How many terms to calculate pi to? 10000 The value of pi is: 3.14169258478796109557151794433593750
In C, 128 bits can be accessed via the type __float128. This may not be availabe in all compilers, but GCC (4.6 and newer) does seem to support it – however, according to the GCC specs, it is said to be “an order of magnitude or two slower” than double precision. To make it even *easier* to use, they provide a bunch of functions, like quadmath_snprintf which converts the number into a string so it can be output. But why make it so difficult? Fortran deals with every variable output in the same way (i.e. transparently).
Here’s what to do – call it a double double, and make the conversion specifier for I/O %TH.