/* Learn you some FPU limits
* by Pegasus Epsilon <pegasus@pimpninjas.org>
*
* build: gcc epsilon.c -lquadmath -o epsilon
* run: ./epsilon
* check hardware math or software math:
* $ gcc -S epsilon.c -lquadmath
* $ grep div epsilon.s
* hardware results:
* fdivrp %st, %st(1)
* fdivrp %st, %st(1)
* software results:
* .globl __divtf3
* call __divtf3@PLT
* call __divtf3@PLT
* Different architectures will yield different results.
*/
#include <quadmath.h>
#include <stdio.h>
/* You should get different results with float, double, long double, __float80,
* and __float128. __float128 is obviously meant to be the most precise.
* __float128 is also probably done in software, so you know... All this
* confusion... The subnormal limit and
* epsilon are *slightly* better using __float128 instead of long double.
* They're both quite a lot better than double. I can see no difference
* whatsoever between __float80 and __float128. Interesting is that __float80
* and __float128 give the same subnormal limit, where __float80 is implemented
* in hardware, and __float128 is implemented in software, on this machine.
* This means that I can get better than long double precision, for no extra
* execution time. I may have to render some super-deep zooms. YMMV, though,
* you probably have a different system.
*/
#define PRECISION 5
#if PRECISION == 1
/* GCC 64bit, storage size seems to be 64 bits */
# define FLOAT float /* ISO */
# define FMT "g"
#elif PRECISION == 2
/* ISO C98, "at least as good as float" */
/* GCC 64bit, storage size seems to be 128 bits */
# define FLOAT double
# define FMT "g"
#elif PRECISION == 3
/* GCC extension, meant to be 80 bits precise */
/* GCC 64bit, storage size seems to be 256 bits */
/* done in hardware here */
# define FLOAT __float80 /* GCC */
# define FMT "Lg"
#elif PRECISION == 4
/* ISO C98, "at least as good as double" */
/* MSVC, exactly the same as double */
/* GCC 32bit, 80 bits precise, storage size depends on flags */
/* GCC 64bit, storage size seems to be 256 bits */
/* done in hardware here */
# define FLOAT long double
# define FMT "Lg"
#else
/* GCC extension, meant to be 128 bits precise */
/* GCC 64bit, storage size seems to be 256 bits */
/* done in software here */
# define FLOAT __float128 /* GCC */
# define FMT "Qg"
#endif
int main (void) {
#if PRECISION == 5
char buf[256] = { 0 };
#endif
FLOAT limit, epsilon, zero;
printf("size of type: %lu bits\n", sizeof(FLOAT) * 16);
printf("format: %%"FMT"\n");
for (zero = (FLOAT)1; (FLOAT)1 != ((FLOAT)1 + zero); zero /= (FLOAT)2)
epsilon = zero;
for (zero = (FLOAT)1; zero; zero /= (FLOAT)2)
limit = zero;
#define EPSILON_MSG "machine epsilon (ϵ): "
#define ROUNDOFF_MSG "unit roundoff (u): "
#define LIMIT_MSG "subnormal limit: "
#if PRECISION == 5
quadmath_snprintf(buf, 256, "%"FMT, epsilon);
printf(EPSILON_MSG"%s\n", buf);
quadmath_snprintf(buf, 256, "%"FMT, epsilon / 2);
printf(ROUNDOFF_MSG"%s\n", buf);
quadmath_snprintf(buf, 256, "%"FMT, limit);
printf(LIMIT_MSG"%s\n", buf);
#else
printf(EPSILON_MSG"%"FMT"\n", epsilon);
printf(ROUNDOFF_MSG"%"FMT"\n", epsilon / 2);
printf(LIMIT_MSG"%"FMT"\n", limit);
#endif
puts("and just in case machine epsilon is displaying as zero above:");
printf("epsilon < limit (should be false): %s\n", epsilon < limit ? "true" : "false");
printf("epsilon == limit (should be false): %s\n", epsilon == limit ? "true" : "false");
printf("epsilon > limit (should be true): %s\n", epsilon > limit ? "true" : "false");
return 0;
}