Chaste  Release::2018.1
MathsCustomFunctions.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "MathsCustomFunctions.hpp"
37 
38 #include <cmath>
39 #include <iostream>
40 
41 double SmallPow(double x, unsigned exponent)
42 {
43  switch (exponent)
44  {
45  case 0:
46  {
47  return 1.0;
48  }
49  case 1:
50  {
51  return x;
52  }
53  case 2:
54  {
55  return x*x;
56  }
57  case 3:
58  {
59  return x*x*x;
60  }
61  default:
62  {
63  if (exponent % 2 == 0)
64  {
65  // Even power
66  double partial_answer = SmallPow(x, exponent/2);
67  return partial_answer*partial_answer;
68  }
69  else
70  {
71  // Odd power
72  return SmallPow(x, exponent-1)*x;
73  }
74  }
75  }
76 }
77 unsigned SmallPow(unsigned x, unsigned exponent)
78 {
79  switch (exponent)
80  {
81  case 0:
82  {
83  return 1u;
84  }
85  case 1:
86  {
87  return x;
88  }
89  case 2:
90  {
91  return x*x;
92  }
93  case 3:
94  {
95  return x*x*x;
96  }
97  default:
98  {
99  if (exponent % 2 == 0)
100  {
101  // Even power
102  unsigned partial_answer = SmallPow(x, exponent/2);
103  return partial_answer*partial_answer;
104  }
105  else
106  {
107  // Odd power
108  return SmallPow(x, exponent-1)*x;
109  }
110  }
111  }
112 }
113 
114 bool Divides(double smallerNumber, double largerNumber)
115 {
116  double remainder = fmod(largerNumber, smallerNumber);
117  /*
118  * Is the remainder close to zero? Note that the comparison is scaled
119  * with respect to the larger of the numbers.
120  */
121  if (remainder < DBL_EPSILON*largerNumber)
122  {
123  return true;
124  }
125  /*
126  * Is the remainder close to smallerNumber? Note that the comparison
127  * is scaled with respect to the larger of the numbers.
128  */
129  if (fabs(remainder-smallerNumber) < DBL_EPSILON*largerNumber)
130  {
131  return true;
132  }
133 
134  return false;
135 }
136 
137 unsigned CeilDivide(unsigned numerator, unsigned denominator)
138 {
139  if (numerator == 0u)
140  {
141  return 0u;
142  }
143  else
144  {
145  // Overflow-safe for large numbers, but not valid for numerator==0.
146  return ((numerator - 1u) / denominator) + 1u;
147  }
148 }
149 
150 double Signum(double value)
151 {
152  return (0.0 < value) - (value < 0.0);
153 }
154 
155 bool CompareDoubles::IsNearZero(double number, double tolerance)
156 {
157  return fabs(number) <= fabs(tolerance);
158 }
159 
160 double SafeDivide(double numerator, double divisor)
161 {
162  // Avoid overflow
163  if (divisor < 1.0 && numerator > divisor*DBL_MAX)
164  {
165  return DBL_MAX;
166  }
167 
168  // Avoid underflow
169  if (numerator == 0.0 || (divisor > 1.0 && numerator < divisor*DBL_MIN))
170  {
171  return 0.0;
172  }
173 
174  return numerator/divisor;
175 }
176 
177 bool CompareDoubles::WithinRelativeTolerance(double number1, double number2, double tolerance)
178 {
179  double difference = fabs(number1 - number2);
180  double d1 = SafeDivide(difference, fabs(number1));
181  double d2 = SafeDivide(difference, fabs(number2));
182 
183  return d1 <= tolerance && d2 <= tolerance;
184 }
185 
186 bool CompareDoubles::WithinAbsoluteTolerance(double number1, double number2, double tolerance)
187 {
188  return fabs(number1 - number2) <= tolerance;
189 }
190 
191 bool CompareDoubles::WithinAnyTolerance(double number1, double number2, double relTol, double absTol, bool printError)
192 {
193  bool ok = WithinAbsoluteTolerance(number1, number2, absTol) || WithinRelativeTolerance(number1, number2, relTol);
194  if (printError && !ok)
195  {
196  std::cout << "CompareDoubles::WithinAnyTolerance: " << number1 << " and " << number2
197  << " differ by more than relative tolerance of " << relTol
198  << " and absolute tolerance of " << absTol << std::endl;
199  }
200  return ok;
201 }
202 
203 bool CompareDoubles::WithinTolerance(double number1, double number2, double tolerance, bool toleranceIsAbsolute)
204 {
205  bool ok;
206  if (toleranceIsAbsolute)
207  {
208  ok = WithinAbsoluteTolerance(number1, number2, tolerance);
209  }
210  else
211  {
212  ok = WithinRelativeTolerance(number1, number2, tolerance);
213  }
214  if (!ok)
215  {
216  std::cout << "CompareDoubles::WithinTolerance: " << number1 << " and " << number2
217  << " differ by more than " << (toleranceIsAbsolute ? "absolute" : "relative")
218  << " tolerance of " << tolerance << std::endl;
219  }
220  return ok;
221 }
222 
223 double CompareDoubles::Difference(double number1, double number2, bool toleranceIsAbsolute)
224 {
225  if (toleranceIsAbsolute)
226  {
227  return fabs(number1 - number2);
228  }
229  else
230  {
231  double difference = fabs(number1 - number2);
232  double d1 = SafeDivide(difference, fabs(number1));
233  double d2 = SafeDivide(difference, fabs(number2));
234  return d1 > d2 ? d1 : d2;
235  }
236 }
static bool WithinAnyTolerance(double number1, double number2, double relTol=DBL_EPSILON, double absTol=DBL_EPSILON, bool printError=false)
double SmallPow(double x, unsigned exponent)
unsigned CeilDivide(unsigned numerator, unsigned denominator)
double SafeDivide(double numerator, double divisor)
bool Divides(double smallerNumber, double largerNumber)
double Signum(double value)
static double Difference(double number1, double number2, bool toleranceIsAbsolute)
static bool WithinRelativeTolerance(double number1, double number2, double tolerance)
static bool WithinAbsoluteTolerance(double number1, double number2, double tolerance)
static bool IsNearZero(double number, double tolerance)
static bool WithinTolerance(double number1, double number2, double tolerance, bool toleranceIsAbsolute)