| #include <limits> |
| |
| extern "C" { |
| float slamch_(char* cmach); |
| double dlamch_(char* cmach); |
| float slamc3_(char* cmach); |
| double dlamc3_(char* cmach); |
| } |
| |
| // Determines machine parameters. |
| template <typename T> |
| T xlamch(const char c) { |
| if (c == 'E' || c == 'e') { |
| // eps = relative machine precision |
| return std::numeric_limits<T>::epsilon()/std::numeric_limits<T>::radix; |
| } |
| |
| if (c == 'S' || c == 's') { |
| // sfmin = safe minimum, such that 1 / sfmin does not overflow |
| T sfmin = std::numeric_limits<T>::min(); |
| const T small = T(1) / std::numeric_limits<T>::max(); |
| if (small >= sfmin) { |
| sfmin = small * (T(1) + std::numeric_limits<T>::epsilon() / |
| std::numeric_limits<T>::radix); |
| } |
| return sfmin; |
| } |
| |
| if (c == 'B' || c == 'b') { |
| // base = base of the machine. |
| return std::numeric_limits<T>::radix; |
| } |
| |
| if (c == 'P' || c == 'p') { |
| // prec = eps * base |
| return std::numeric_limits<T>::epsilon(); |
| } |
| |
| if (c == 'N' || c == 'n') { |
| // t = number of (base) digits in the mantissa |
| return std::numeric_limits<T>::digits; |
| } |
| |
| if (c == 'R' || c == 'r') { |
| // rnd = 1.0 when rounding occurs in addition. |
| return 1; |
| } |
| |
| if (c == 'M' || c == 'm') { |
| // emin = minimum exponent before (gradual) underflow |
| return std::numeric_limits<T>::min_exponent; |
| } |
| |
| if (c == 'U' || c == 'u') { |
| // rmin = underflow threshold - base**(emin-1) |
| return std::numeric_limits<T>::min(); |
| } |
| |
| if (c == 'L' || c == 'l') { |
| // emax = largest exponent before overflow |
| return std::numeric_limits<T>::max_exponent; |
| } |
| if (c == 'O' || c == 'o') { |
| // rmax = overflow threshold - (base**emax)*(1-eps) |
| return std::numeric_limits<T>::max(); |
| } |
| |
| return 0; |
| } |
| |
| float slamch_(char* cmach) { return xlamch<float>(*cmach); } |
| double dlamch_(char* cmach) { return xlamch<double>(*cmach); } |
| |
| // xLAMC3 is intended to force A and B to be stored prior to doing |
| // the addition of A and B , for use in situations where optimizers |
| // might hold one of these in a register. |
| float slamc3_(float *a, float *b) { |
| volatile float retval = *a + *b; |
| return retval; |
| } |
| |
| double dlamc3_(double *a, double *b) { |
| volatile double retval = *a + *b; |
| return retval; |
| } |
| |