ulp.c

Go to the documentation of this file.
00001 /*                           U L P . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 2010-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @file libbn/ulp.c
00021  *
00022  * Routines useful for performing comparisons and dynamically
00023  * calculating floating point limits including the Unit in the Last
00024  * Place (ULP).
00025  *
00026  * In this context, ULP is the distance to the next normalized
00027  * floating point value larger that a given input value.
00028  *
00029  * TODO: handle NaN, +-Inf, underflow, overflow, non-IEEE, float.h
00030  *
00031  * This file is completely in flux, incomplete, limited, and subject
00032  * to drastic changes.  Do NOT use it for anything.
00033  *
00034  * It also assumes an IEEE 754 compliant floating point
00035  * representation.
00036  */
00037 
00038 #include "common.h"
00039 
00040 #include <float.h>
00041 
00042 
00043 double
00044 bn_epsilon()
00045 {
00046 #if defined(DBL_EPSILON)
00047     return DBL_EPSILON;
00048 #elif defined(HAVE_IEEE754)
00049     static const double val = 1.0;
00050     register long long next = *(long long*)&val + 1;
00051     return val - *(double *)&next;
00052 #else
00053     /* must be volatile to avoid long registers */
00054     volatile double tol = 1.0;
00055     while (1.0 + (tol * 0.5) != 1.0) {
00056         tol *= 0.5;
00057     }
00058 #endif
00059 }
00060 
00061 
00062 float
00063 bn_epsilonf()
00064 {
00065 #if defined(FLT_EPSILON)
00066     return FLT_EPSILON;
00067 #elif defined(HAVE_IEEE754)
00068     static const float val = 1.0;
00069     register long next = *(long*)&val + 1;
00070     return val - *(float *)&next;
00071 #else
00072     /* must be volatile to avoid long registers */
00073     volatile float tol = 1.0f;
00074     while (1.0f + (tol * 0.5f) != 1.0f) {
00075         tol *= 0.5f;
00076     }
00077 #endif
00078 }
00079 
00080 
00081 double
00082 bn_dbl_min()
00083 {
00084     register long long val = (1LL<<52);
00085     return *(double *)&val;
00086 }
00087 
00088 
00089 double
00090 bn_dbl_max()
00091 {
00092     static const double val = INFINITY;
00093     register long long next = *(long long*)&val - 1;
00094     return *(double *)&next;
00095 }
00096 
00097 
00098 double
00099 bn_flt_min()
00100 {
00101     register long val = (1LL<<23);
00102     return *(float *)&val;
00103 }
00104 
00105 
00106 double
00107 bn_flt_max()
00108 {
00109     static const float val = INFINITY;
00110     register long next = *(long*)&val - 1;
00111     return *(float *)&next;
00112 }
00113 
00114 
00115 double
00116 bn_ulp(double val)
00117 {
00118     double next;
00119     register long long up, dn, idx;
00120 
00121     if (isnan(val) || !isfinite(val))
00122         return val;
00123 
00124     if (val >=0) {
00125         up = *(long long*)&val + 1;
00126         return *(double *)&up - val;
00127     }
00128     dn = *(long long*)&val - 1;
00129     return *(double *)&dn - val;
00130 }
00131 
00132 
00133 /*
00134  * Local Variables:
00135  * tab-width: 8
00136  * mode: C
00137  * indent-tabs-mode: t
00138  * c-file-style: "stroustrup"
00139  * End:
00140  * ex: shiftwidth=4 tabstop=8
00141  */
Generated on Tue Dec 11 13:14:28 2012 for LIBBN by  doxygen 1.6.3