Plan for Implementing Quadruple Precision (REAL(KIND=16)) in LFortran

Overview

LFortran will support full IEEE 754 binary128 quadruple precision (REAL(KIND=16) and COMPLEX(KIND=16)) using LLVM’s native fp128 type for all arithmetic and most operations. This provides exact binary128 semantics (15-bit exponent, 113-bit mantissa) on every target LLVM supports.

Transcendental intrinsics will be provided through an optional, thin runtime wrapper that can link to libquadmath (or the platform’s _Float128 / long-double functions when available). Basic operations, array intrinsics, and most language features will work with zero extra dependencies.

The implementation is designed to be incremental, portable, and consistent with existing LFortran infrastructure (ASR → LLVM IR lowering, runtime library, multiple backends, JIT).

1. Frontend and ASR Level

• Extend the parser and semantics to fully recognize KIND=16:
  • Literals: 1.0_q0, 1.q0, REAL(1.0, KIND=16)
  • SELECTED_REAL_KIND(precision=33, range=4931) → 16
  • ISO_FORTRAN_ENV: expose REAL128 (value 16)
• Update all type inquiry intrinsics (KIND, PRECISION, RANGE, EPSILON, etc.) to return correct values for KIND=16.
• Add RealKind::Quad (or equivalent) in ASR and ensure COMPLEX(KIND=16) is represented as two fp128 values.
• Support SELECT TYPE, TRANSFER, and all generic interfaces for quad precision.

2. Lowering and LLVM Code Generation

• Map ASR real(kind=16) directly to LLVM fp128 type (llvm::Type::getFP128Ty).
• Use ConstantFP::get with APFloat::IEEEquad semantics for all literals and compile-time constants.
• All basic operations lower automatically:
  • Arithmetic: +, −, ×, /, ** (real**real)
  • Comparisons, casts, ABS, SIGN, MOD, MODULO, AINT, ANINT, FLOOR, CEILING, FRACTION, EXPONENT, SET_EXPONENT, NEAREST, SPACING, SCALE, etc.
  • SQRT (via llvm.sqrt.f128)
• Array intrinsics lower to loops or LLVM vectorized operations on fp128:
  • MATMUL, DOT_PRODUCT
  • SUM, PRODUCT, MAXVAL, MINVAL, MAXLOC, MINLOC
  • TRANSPOSE, PACK, UNPACK, reductions, and any elemental array operation
• I/O formatting, READ/WRITE with quad values, and all runtime support use the existing LFortran I/O infrastructure (already partially implemented).

These features require no extra libraries and work immediately in the LLVM backend, JIT, and interactive mode.

3. Math Intrinsics and Runtime Library

Default Operations (No Extra Library Required)

The following are handled purely by LLVM fp128 and existing runtime code:

• All arithmetic, comparisons, and conversions listed above
• Array intrinsics (MATMUL, DOT_PRODUCT, SUM, PRODUCT, etc.)
• Most elemental intrinsics: ABS, SIGN, MOD, FLOOR, CEILING, SQRT, FRACTION, EXPONENT, etc.
• All I/O, type inquiry, and language features

Operations Requiring the Quad Math Runtime Wrapper

These map to transcendental functions and will be provided by a thin C++ wrapper in LFortran’s runtime (one file per intrinsic):

Real functions (suffix q in the library):

• Trigonometric & hyperbolic: acos, asin, atan, atan2, cos, sin, tan, acosh, asinh, atanh, cosh, sinh, tanh, sincos
• Exponential & logarithmic: exp, exp2, expm1, log, log2, log10, log1p, pow
• Roots: sqrt (fallback to LLVM if possible), cbrt, hypot
• Rounding: ceil, floor, round, trunc, rint, nearbyint, lrint, llrint, lround, llround
• Manipulation: frexp, ldexp, modf, scalbn, scalbln, ilogb
• Classification: isnan, isinf, finite, issignaling, signbit, copysign, fmax, fmin, fdim, fabs, fma
• Special: erf, erfc, tgamma, lgamma
• Bessel: bessel_j0, bessel_j1, bessel_jn, bessel_y0, bessel_y1, bessel_yn
• Other: nan, nextafter, fmod, remainder, remquo

Complex functions:

• abs, arg, imag, real, conj, proj
• exp, log, log10, pow, sqrt
• sin, cos, tan, sinh, cosh, tanh
• asin, acos, atan, asinh, acosh, atanh

The runtime wrapper (lfortran_rt.quadmath or similar) will be built only when the CMake flag is enabled and will call the underlying library functions.

Platform Fallback for libc/libm

On platforms where long double is already IEEE 128-bit (LDBL_MANT_DIG == 113), the wrapper can alias directly to standard libm long-double functions (or C23 _Float128 suffixed functions when available in glibc ≥ 2.26). No extra library is needed.

4. Other Backends and Interoperability

• C backend: Emit __float128 (GCC) or _Float128 (C23).
• WASM / experimental backends: Use LLVM’s software fp128 emulation or emit a clear compile-time error if unsupported.
• COMPLEX(KIND=16): Stored as two consecutive fp128 values (matches ABI expectations).
• Maintain full interoperability with C/C++ via iso_c_binding (REAL(KIND=16) ↔ __float128).

5. Build System and Configuration

• Add CMake flag: -DLFORTRAN_F128_MATH_LIB=libquadmath (default: auto-detect or off)
• When enabled, build the thin wrapper library and link it automatically (similar to how other runtime components are handled).
• Document the optional nature of the dependency and any LGPL implications if libquadmath is used.
• Provide a pure-LFortran fallback mode (optional) using MPFR at 113-bit precision only if licensing requires avoiding external libraries.

6. Testing and Validation

• Extend existing test suite:
  • Literals, mixed-precision expressions, SELECTED_REAL_KIND
  • All listed intrinsics (both real and complex)
  • Array operations (MATMUL, DOT_PRODUCT, reductions)
  • I/O round-trip accuracy
  • Edge cases: denormals, NaN, Inf, under/overflow
• Add cross-compilation and platform-specific tests (x86_64 with/without libquadmath, PowerPC with native 128-bit long double, etc.).
• Ensure JIT and interactive mode support quad literals and intrinsics out of the box.

Why This Approach Is Optimal

• Minimal effort and maximal reuse: Leverages LLVM’s mature fp128 support for everything possible → arithmetic and array intrinsics work immediately.
• Portability: Optional dependency; zero-cost on platforms with native 128-bit long double.
• Performance: Uses highly optimized library routines for transcendentals.
• User experience: Matches expectations from other Fortran compilers; no surprises at link time for basic code.
• Future-proof: Works with all LFortran backends and scales to new targets as LLVM improves fp128 support.

This plan gives LFortran production-grade quadruple precision support with a clean, maintainable architecture. Implementation can start with frontend/ASR + LLVM lowering (quick win), then add the runtime wrapper incrementally.