Listing 3 Computes Machine Floating-point Parameters
#include <stdio.h> #include <math.h> #include <float.h> main() { int beta, p; double a, b, eps, epsp1, sigma, nums; /* Discover radix */ a = 1.0; do { a = 2.0 * a; b = a + 1.0; } while ((b - a) == 1.0); b = 1.0; do b = 2.0 * b; while ((a + b) == a); beta = (int) ((a + b) - a); printf("radix:\n"); printf("\talgorithm: %d\n",beta); printf("\tprovided: %d\n",FLT_RADIX); /* Compute precision in bits */ p = 0; a = 1.0; do { ++p; a *= (double) beta; b = a + 1.0; } while ((b - a) == 1.0); printf("precision:\n"); printf("\talgorithm: %d\n",p); printf("\tprovided: %d\n",DBL_MANT_DIG); /* Compute machine epsilon */ eps = 1.0; do { eps = 0.5 * eps; epsp1 = eps + 1.0; } while (epsp1 > 1.0); epsp1 = 2.0 * eps + 1.0; eps = epsp1 - 1.0; printf("machine epsilon:\n"); printf("\talgorithm: %g\n",eps); printf("\tformula: %g\n",pow(FLT_RADIX,1.0-DBL_MANT_DIG)); printf("\tprovided: %g\n",DBL_EPSILON); /* Compute smallest normalized magnitude */ printf("smallest non-zero magnitude:\n"); printf("\tformula: %g\n",pow(FLT_RADIX,DBL_MIN_EXP-1)); printf("\tprovided: %g\n",DBL_MIN); /* Compute larget normalized magnitude */ printf("largest non-zero magnitude:\n"); printf("\tformula: %g\n", pow(FLT_RADIX,DBL_MAX_EXP-1) * FLT_RADIX * (1.0 - pow(FLT_RADIX,-DBL_MANT_DIG))); printf("\tprovided: %g\n",DBL_MAX); printf("smallest exponent: %d\n",DBL_MIN_EXP); printf("largest exponent: %d\n",DBL_MAX_EXP); nums = 2 * (FLT_RADIX - 1) * pow(FLT_RADIX,DBL_MANT_DIG-1) * (DBL_MAX_EXP - DBL_MIN_EXP + 1); printf("This system has %g numbers\n",nums); return 0; } /* Output (Borland C++): radix: algorithm: 2 provided: 2 precision: algorithm: 53 provided: 53 machine epsilon: algorithm: 2.22045e-16 formula: 2.22045e-16 provided: 2.22045e-16 smallest non-zero magnitude: formula: 2.22507e-308 provided: 2.22507e-308 largest non-zero magnitude: formula: 1.79769e+308 provided: 1.79769e+308 smallest exponent: -1021 largest exponent: 1024 This system has 1.84287e+19 numbers */ /* End of File */