/* Simple Integrator; version 1.0 */
/* (c) Ian Hickson 2000, distributed under the GNU GPL */
#include /* for input/output */
#include /* for maths */
#define MaxSignificantFigures 16
#define MinSignificantFigures 3
#define Function exp(-pow(x,2))
#define KnownAnswer 0.5*sqrt(M_PI)
#define LowLimit 0 /* lower limit -- set by brief */
#undef HighLimit /* upper limit -- use findHighLimit() */
/* cannot use "long double" mathematics because exp() and pow() don't
support long doubles!!! */
/* Wrapper around the function f(x) which we will integrate */
double f(double x) {
return Function;
}
/* Returns the number to use as the high limit, typically, the number
to use as infinity. This is very roughly the point at which the
value of the function to be integrated becomes negligible (f(xN)=0)
*/
double findHighLimit() {
#ifndef HighLimit /* otherwise we can simply use given value */
/* Sweeps values of xN to find at what point f(x) becomes negligible */
/* WARNING! This assumes that the function is asymptotic about y=0
and does not cross the x axis. More importantly, it assumes that
the high limit is infinity! #define HighLimit if these
assumptions do not apply. */
/* Note: Granularity is 1.0 */
double xN;
for (xN=LowLimit; f(xN); xN++);
return xN;
#else
return HighLimit;
#endif
}
/* Returns the number to use as the low limit, typically 0 */
double findLowLimit() {
#ifndef LowLimit
#error no low limit algorithm defined!
#else
return LowLimit;
#endif
}
/* equivalent of abs() but for doubles (abs() is for ints only) */
double absd(double x) {
return x < 0.0 ? -x : x;
}
/* returns the number of common signifcant figures between doubles a
and b */
int getCommonSignificantFigures(double a, double b) {
/* Method:
+ multiply or divide a and b until one is 0.1 <= x < 1
and the other is y < 1
+ multiply both by ten, comparing their integer portions until
they differ, or until the number of comparisons exceeds the
number of significant digits available in the type used (as
specified by MaxSignificantFigures).
Possible bugs XXX: This will overflow if a or b are close to the
maximum or minimum exponents already (won't occur with the
equation of assessment version 54).
*/
int MatchingSignificantFigures;
while (absd(a) < 0.1 || absd(b) < 0.1) {
a *= 10;
b *= 10;
}
while (absd(a) >= 1 || absd(b) >= 1) {
a /= 10;
b /= 10;
}
MatchingSignificantFigures = -1;
do {
MatchingSignificantFigures++;
a *= 10;
b *= 10;
} while (((long long int)a == (long long int)b) &&
(MatchingSignificantFigures < MaxSignificantFigures));
return MatchingSignificantFigures;
}
/* Integrate from x1 to xN in lastN steps */
/* Uses the Extended Trapezoidal Numerical method. */
double Integrate(double x1, double xN, int lastN) {
double I = 0.5 * (f(x1)+f(xN)); /* First and last terms */
double h = (xN - x1) / (lastN - 1); /* Find the trapezium width */
double stepDelta;
int i;
for (i=1; i MinSignificantFigures) {
Result = EstimateIntegral(x1, xN, SignificantFiguresRequired);
printf(" It would seem the integral is: %.*g to %d significant figures\n",
SignificantFiguresRequired, Result, SignificantFiguresRequired);
printf(" My full estimate is..........: %.*g\n", MaxSignificantFigures, Result);
printf(" The actual answer is.........: %.*g to %d significant figures [KNOWN]\n",
SignificantFiguresRequired, KnownAnswer, SignificantFiguresRequired);
printf(" The full actual answer is....: %.*g [KNOWN]\n",
MaxSignificantFigures, KnownAnswer);
printf(" The first error is in significant digit number %d.\n",
getCommonSignificantFigures(KnownAnswer, Result)+1);
printf(" The error is +/-%g or %g%%.\n\n", absd(Result-KnownAnswer),
absd(100.0*((Result-KnownAnswer)/KnownAnswer)));
} else {
printf("Sorry, I cannot be confident about numbers with such a low\
accuracy (less than %d).\n", MinSignificantFigures+1);
}
} else {
printf("Sorry, I cannot be confident about that kind of high accuracy\
(more than %d).\n", MaxSignificantFigures-1);
}
}