Feature/float128 a

[wuzhiqiao-ninecube]

• Please target the develop branch not main.
• Review our Contributing Guide, and ensure that you sign your last commit (at minimum) as per the guide.
• Provide a concise description of what your pull request does, and why it is needed/benefical.
• Add a note about your change to the CHANGELOG.md and docs/shared/RecentChanges.rst files. Notice that the former is a markdown file and the latter is reStructuredText, so the formatting is slightly different.
• After your PR is opened, ensure that all of the tests are passing (a SUNDIALS developer will have to allow the testing to run).

Add comprehensive __float128 support with quadmath dependency and ostream integration

This commit implements full support for the __float128 data type across SUNDIALS:

1. Core Implementation

   • Added dependency on libquadmath for quad-precision math operations
   • New CMake configuration options:
     • -D SUNDIALS_PRECISION=float128
     • -D SUNDIALS_MATH_LIBRARY=SUNDIALS_MATH_LIBRARY=/usr/lib/x86_64-linux-gnu/libm.so;/usr/lib/gcc/x86_64-linux-gnu/13/libquadmath.so
   • Automatic definition of SUNDIALS_FLOAT128_PRECISION macro when enabled

2. Precision System Enhancements

   • Added SUN_DIGITS10 constant to sundials_type.h (avoiding std::numeric_limits::digits10)
   • Implemented complete set of quadmath math functions:
     • SUNRsqrt, SUNRabs, SUNRisnan, SUNRexp, SUNRlog
     • SUNRceil, SUNRcopysign, SUNRpowerR, SUNRround
     • SUNRsin, SUNRcos, SUNRsinh, SUNRcosh, SUNRasin, SUNRacos, SUNRatan
   • three macros — SUNRcosh, SUNRsinh, and SUNRlog — were moved from the src/arkode/arkode_lsrkstep_impl.h file to include/sundials/sundials_math.h
   • Increased precision of PI constant to 40+ digits:

     #define PI SUN_RCONST(3.141592653589793238462643383279502884)

     The output files (*.out) of the related examples have been changed as a result, and a corresponding commit was made.

3. Code Unification

   • Replaced all local math macros (ABS, SIN, etc.) with unified SUNR* macros
   • Updated C/C++ examples to use precision-portable math operations:

     // Before
     #define ATAN(x) (atan((x)))
     
     // After
     y = SUNRatan(x);

   • Modified C++ examples to use SUNR* macros instead of std functions

4. I/O System

   • Implemented printf formatting for __float128:

     • Defined GSYM="Qg", ESYM="Qe", FSYM="Qf" for float128
     • Added precision-specific format specifiers:

       #if defined(SUNDIALS_FLOAT128_PRECISION)
         printf("Value = %Qg", float128_val);
       #endif

       However, the compilation warning "warning: unknown conversion type character ‘Q’ in format [-Wformat=]" could not be eliminated.

   • Added std::ostream support for __float128:

     • Full stream state capture (width/precision/flags)
     • Dynamic quadmath_snprintf formatting
     • Scientific/fixed/general format detection
     • Example usage:

       __float128 q = 1.234q;
       std::cout << std::scientific << std::setprecision(15) << q;

5. Environment Support Matrix

   Environment	Support Status
   C	✓ Implemented
   C++	✓ Implemented
   OpenMP	✓ Implemented
   Pthread	✓ Implemented
   MPI	✓ Implemented
   Fortran	✗ Not supported
   OpenMP Device	✗ Not supported
   CUDA	✗ Not supported
   HIP	✗ Not supported

Key Technical Decisions:

• Avoided std::numeric_limits in favor of SUN_DIGITS10
• Unified math operations through SUNR* macros
• Maintained printf compatibility with quadmath specifiers
• Implemented native C++ stream integration
• Ensured precision constants meet float128 requirements

This implementation enables quad-precision computations while maintaining API consistency across all precision modes. The solution provides seamless integration with both C-style printf and C++ ostream interfaces, with special attention to cross-environment compatibility and precision portability.

Signed-off-by: Wu Zhiqiao wuzhiqiao79@gmail.com

[drreynolds]

I could be mistaken, but I believe that __float128 is only supported by GCC, and even then, only on x86-64 and i386 targets. In other words, __float128 is not supported by other compilers (e.g., Clang, MSVC) or on other architectures (AArch64, most GPUs). One of our top goals with SUNDIALS is to ensure that it is portable, and not specifically tied to a small set of compilers/architectures.

Furthermore, the existing heuristic constants, and many of the method coefficients, throughout SUNDIALS have only been tested for at most long double, which generally has less precision (80-bit on many compilers/architectures), so it is unclear which of these would need to be adjusted for true 128-bit floating-point computations. If this is nontrivial, then our team will struggle to take on the level of "new work" that would be required to truly support __float128 throughout the SUNDIALS.

@wuzhiqiao-ninecube, can you provide a compelling reason why you think that SUNDIALS should add comprehensive __float128 support, including dependency on the [GCC-specific] libquadmath library? I realize that this implementation likely took considerable effort, and in perusing the changes it looks like you fixed a number of minor issues throughout the repository, so I hate to disregard it entirely, but without a very clear motivation we cannot accept a contribution that rules out so many modern computing systems and compilers, and that would significantly add so much unfunded effort onto our team.

Let me also point out half of the first paragraph from our Contributing Guide, where we say: "At this time, the SUNDIALS team does not have the resources to review and take in large additions to the code or significant new features. Larger additions can be contributed as a SUNDIALS "addon" which is a component that may be optionally downloaded by users and then compiled and installed with SUNDIALS."

[wuzhiqiao-ninecube]

I could be mistaken, but I believe that __float128 is only supported by GCC, and even then, only on x86-64 and i386 targets. In other words, __float128 is not supported by other compilers (e.g., Clang, MSVC) or on other architectures (AArch64, most GPUs). One of our top goals with SUNDIALS is to ensure that it is portable, and not specifically tied to a small set of compilers/architectures.

Furthermore, the existing heuristic constants, and many of the method coefficients, throughout SUNDIALS have only been tested for at most long double, which generally has less precision (80-bit on many compilers/architectures), so it is unclear which of these would need to be adjusted for true 128-bit floating-point computations. If this is nontrivial, then our team will struggle to take on the level of "new work" that would be required to truly support __float128 throughout the SUNDIALS.

@wuzhiqiao-ninecube, can you provide a compelling reason why you think that SUNDIALS should add comprehensive __float128 support, including dependency on the [GCC-specific] libquadmath library? I realize that this implementation likely took considerable effort, and in perusing the changes it looks like you fixed a number of minor issues throughout the repository, so I hate to disregard it entirely, but without a very clear motivation we cannot accept a contribution that rules out so many modern computing systems and compilers, and that would significantly add so much unfunded effort onto our team.

Let me also point out half of the first paragraph from our Contributing Guide, where we say: "At this time, the SUNDIALS team does not have the resources to review and take in large additions to the code or significant new features. Larger additions can be contributed as a SUNDIALS "addon" which is a component that may be optionally downloaded by users and then compiled and installed with SUNDIALS."

Thank you very much for your detailed and thoughtful feedback on my pull request. I truly appreciate the time you took to review the changes and to explain the portability and maintainability concerns.

I understand that __float128 and libquadmath are traditionally associated with GCC on x86-64/i386, and that portability is a top priority for SUNDIALS. However, I’d like to share a few additional observations and motivations that might help clarify the potential value of this feature.

1. Compiler and Platform Support for __float128
   I have done some preliminary testing on Ubuntu with GCC, Clang, and Intel compilers, and all three successfully compiled the code using __float128 and libquadmath. While I fully acknowledge that this is not a guarantee for all platforms (e.g., MSVC, AArch64, GPUs), it does suggest that support may be broader than initially assumed. Still, I agree that any integration must be done in a way that does not break the build on unsupported platforms.

2. Motivation for Adding __float128 Support
   My interest in this feature stems from several scientific computing needs:

High-precision reference solutions: In ODE/DAE problems, having a very accurate reference solution (computed with __float128) can help validate the correctness of lower-precision simulations and numerical methods.

Semiconductor device simulation: In my previous work on semiconductor physics and device simulation, extremely small quantities (e.g., carrier concentrations, tunneling probabilities) often require precision beyond double or even long double. __float128 provides a valuable tool for such problems.

Mixed-precision iterative solvers: In the context of linear system solvers, mixed-precision strategies (e.g., using __float128 for residuals or selective refinement) are an active research area. Having native __float128 support in SUNDIALS would enable experimentation and potential advances in this direction.

3. Nature of the Modifications
   I should also note that the majority of changes in the pull request are concentrated in the examples directory. The core modifications are limited to two files: sundials_math.h and sundials_types.hpp. In the examples, I replaced locally-defined macros for mathematical functions (e.g., SUNRsqrt, SUNRabs, SUNRisnan, etc.) with the ones already provided by sundials_math.h, and similarly unified the printf format specifier macros (GSYM, ESYM, FSYM). These changes are not strictly tied to __float128—they improve code consistency and could be submitted as a separate, independent pull request. If that would make review easier, I’m happy to split the work and first contribute these cleanups.

4. Exploring the Addon Approach
   Regarding the suggestion to implement this feature as an addon (as mentioned in the Contributing Guide), I agree that this could be a viable path. However, I’m not entirely sure how to structure an addon that would allow __float128 to be used as a basic sunrealtype option. Typically, addons are more self-contained components (e.g., new solvers or vector implementations). I’m open to discussing how we might design such an addon, or if there are alternative ways to integrate high-precision support without compromising the core’s portability. If the team has any guidance or examples, I would be glad to follow them.

In the meantime, I would be happy to:

Split off the non‑__float128 improvements into a separate PR for immediate inclusion.

Offer to help maintain any 128‑bit related code in the future, to reduce the burden on your team.

Please let me know what you think would be best. Again, thank you for your honest and constructive feedback—I really appreciate it, and I’m committed to aligning this contribution with SUNDIALS’ goals.

Best regards,
wuzhiqiao