00001 # ifndef CPPAD_NEAR_EQUAL_INCLUDED
00002 # define CPPAD_NEAR_EQUAL_INCLUDED
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 # include <complex>
00152 # include <cppad/local/cppad_assert.hpp>
00153 # include <cppad/check_numeric_type.hpp>
00154
00155 namespace CppAD {
00156
00157
00158 template <class Type>
00159 bool isfinite(const Type &z1, const Type &z2, const Type &x , const Type &y)
00160 { Type infinity = Type(1) / z1;
00161 Type nan = z1 / z2;
00162
00163
00164 bool xNan = ( x != x || x == nan );
00165 bool yNan = ( y != y || y == nan );
00166
00167
00168 bool xInf = (x == infinity || x == - infinity);
00169 bool yInf = (x == infinity || x == - infinity);
00170
00171 return ! (xNan | yNan | xInf | yInf);
00172 }
00173
00174 template <class Type>
00175 bool NearEqual(const Type &x, const Type &y, const Type &r, const Type &a)
00176 {
00177 CheckNumericType<Type>();
00178 Type zero(0);
00179
00180 CPPAD_ASSERT_KNOWN(
00181 zero <= r,
00182 "Error in NearEqual: relative error is less than zero"
00183 );
00184 CPPAD_ASSERT_KNOWN(
00185 zero <= a,
00186 "Error in NearEqual: absolute error is less than zero"
00187 );
00188
00189
00190 if( ! isfinite(zero, zero, x, y) )
00191 return false;
00192
00193 Type ax = x;
00194 if( ax <= zero )
00195 ax = - ax;
00196
00197 Type ay = y;
00198 if( ay <= zero )
00199 ay = - ay;
00200
00201 Type ad = x - y;
00202 if( ad <= zero )
00203 ad = - ad;
00204
00205 if( ad <= a )
00206 return true;
00207
00208 if( ad <= r * (ax + ay) )
00209 return true;
00210
00211 return false;
00212 }
00213
00214 template <class Type>
00215 bool NearEqual(
00216 const std::complex<Type> &x ,
00217 const std::complex<Type> &y ,
00218 const Type &r ,
00219 const Type & a )
00220 {
00221 CheckNumericType<Type>();
00222 Type zero(0);
00223
00224 CPPAD_ASSERT_KNOWN(
00225 zero <= r,
00226 "Error in NearEqual: relative error is less than zero"
00227 );
00228 CPPAD_ASSERT_KNOWN(
00229 zero <= a,
00230 "Error in NearEqual: absolute error is less than zero"
00231 );
00232
00233
00234 if( ! isfinite(zero, zero, x.real(), x.imag()) )
00235 return false;
00236 if( ! isfinite(zero, zero, y.real(), y.imag()) )
00237 return false;
00238
00239 std::complex<Type> d = x - y;
00240
00241 Type ad = std::abs(d);
00242 if( ad <= a )
00243 return true;
00244
00245 Type ax = std::abs(x);
00246 Type ay = std::abs(y);
00247 if( ad <= r * (ax + ay) )
00248 return true;
00249
00250 return false;
00251 }
00252
00253 template <class Type>
00254 bool NearEqual(
00255 const std::complex<Type> &x ,
00256 const Type &y ,
00257 const Type &r ,
00258 const Type & a )
00259 {
00260 return NearEqual(x, std::complex<Type>(y, Type(0)), r, a);
00261 }
00262
00263 template <class Type>
00264 bool NearEqual(
00265 const Type &x ,
00266 const std::complex<Type> &y ,
00267 const Type &r ,
00268 const Type & a )
00269 {
00270 return NearEqual(std::complex<Type>(x, Type(0)), y, r, a);
00271 }
00272
00273 }
00274
00275 # endif