BRL-CAD
ulp.c
Go to the documentation of this file.
1 /* U L P . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2010-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file libbn/ulp.c
21  *
22  * NOTE: This is a work-in-progress and not yet published API to be
23  * used anywhere. DO NOT USE.
24  *
25  * Routines useful for performing comparisons and dynamically
26  * calculating floating point limits including the Unit in the Last
27  * Place (ULP).
28  *
29  * In this context, ULP is the distance to the next normalized
30  * floating point value larger than a given input value.
31  *
32  * TODO: handle NaN, +-Inf, underflow, overflow, non-IEEE, float.h
33  *
34  * This file is completely in flux, incomplete, limited, and subject
35  * to drastic changes. Do NOT use it for anything.
36  *
37  * It also assumes an IEEE 754 compliant floating point
38  * representation.
39  */
40 
41 #include "common.h"
42 
43 #include <math.h>
44 #include <limits.h>
45 #ifdef HAVE_FLOAT_H
46 # include <float.h>
47 #endif
48 
49 /* #define HAVE_IEEE754 1 */
50 
51 double
53 {
54 #if defined(DBL_EPSILON)
55  return DBL_EPSILON;
56 #elif defined(HAVE_IEEE754)
57  static const double val = 1.0;
58  long long next = *(long long*)&val + 1;
59  return val - *(double *)&next;
60 #else
61  /* must be volatile to avoid long registers */
62  volatile double tol = 1.0;
63  while (1.0 + (tol * 0.5) > 1.0) {
64  tol *= 0.5;
65  }
66  return tol;
67 #endif
68 }
69 
70 
71 float
73 {
74 #if defined(FLT_EPSILON)
75  return FLT_EPSILON;
76 #elif defined(HAVE_IEEE754)
77  static const float val = 1.0;
78  long next = *(long*)&val + 1;
79  return val - *(float *)&next;
80 #else
81  /* must be volatile to avoid long registers */
82  volatile float tol = 1.0f;
83  while (1.0f + (tol * 0.5f) > 1.0f) {
84  tol *= 0.5f;
85  }
86  return tol;
87 #endif
88 }
89 
90 
91 double
93 {
94  long long val = (1LL<<52);
95  return *(double *)&val;
96 }
97 
98 
99 double
101 {
102 #if defined(DBL_MAX)
103  return DBL_MAX;
104 #elif defined(INFINITY)
105  static const double val = INFINITY;
106  long long next = *(long long*)&val - 1;
107  return *(double *)&next;
108 #else
109  return 1.0/bn_dbl_min();
110 #endif
111 }
112 
113 
114 double
116 {
117  long val = (1LL<<23);
118  return *(float *)&val;
119 }
120 
121 
122 double
124 {
125 #if defined(FLT_MAX)
126  return FLT_MAX;
127 #elif defined(INFINITY)
128  static const float val = INFINITY;
129  long next = *(long*)&val - 1;
130  return *(float *)&next;
131 #else
132  return 1.0/bn_flt_min();
133 #endif
134 }
135 
136 
137 double
139 {
140  return sqrt(bn_flt_min());
141 }
142 
143 
144 double
146 {
147  return sqrt(bn_flt_max());
148 }
149 
150 
151 double
153 {
154  return sqrt(bn_dbl_min());
155 }
156 
157 
158 double
160 {
161  return sqrt(bn_dbl_max());
162 }
163 
164 
165 double
166 bn_ulp(double val)
167 {
168  long long up, dn;
169 
170  if (isnan(val) || isinf(val))
171  return val;
172 
173  if (val >=0) {
174  up = *(long long*)&val + 1;
175  return *(double *)&up - val;
176  }
177  dn = *(long long*)&val - 1;
178  return *(double *)&dn - val;
179 }
180 
181 
182 /*
183  * Local Variables:
184  * tab-width: 8
185  * mode: C
186  * indent-tabs-mode: t
187  * c-file-style: "stroustrup"
188  * End:
189  * ex: shiftwidth=4 tabstop=8
190  */
Definition: db_flip.c:35
double bn_dbl_min(void)
Definition: ulp.c:92
double bn_flt_min_sqrt(void)
Definition: ulp.c:138
float bn_epsilonf(void)
Definition: ulp.c:72
double bn_flt_max(void)
Definition: ulp.c:123
Header file for the BRL-CAD common definitions.
double bn_dbl_max_sqrt(void)
Definition: ulp.c:159
float f
Definition: db_flip.c:36
double bn_dbl_max(void)
Definition: ulp.c:100
double bn_dbl_min_sqrt(void)
Definition: ulp.c:152
double bn_epsilon(void)
Definition: ulp.c:52
double bn_flt_min(void)
Definition: ulp.c:115
double bn_flt_max_sqrt(void)
Definition: ulp.c:145
double bn_ulp(double val)
Definition: ulp.c:166