Files

main.cpp

#include <iostream>
#include <cmath>
#include <stdexcept>
#include <limits>

using namespace std;


double integrateSin( double a, double b, int N ) {
    double h = (b-a)/N;
    double total = 0.0;
    for (int i=0; i<N; i++) {
        double x = (i+0.5)*h + a;
        double f = sin(x);
        total +=f;
    }
    return total/N;
}


double infiniteIntegral( double x ) {
    // we perform the substitution
    // x + 1 - 1/s;
    // to change the infinite integral to an integral between 0 and 1
    double a = 0;
    double b = 1;
    int N = 1000;
    double h = (b-a)/N;
    double total = 0.0;
    for (int i=0; i<N; i++) {
        double s = (i+0.5)*h + a;
        double t = x + 1 - 1/s;
        double f = pow(s,-2) * exp( - 0.5*t*t );
        total +=f;
    }
    return total/N;
}


double hornerFunction( double x, double a0, double a1) {
    return a0 + x*a1;
}

double hornerFunction( double x, double a0, double a1, double a2) {
    return a0 + x*hornerFunction( x, a1, a2);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3) {
    return a0 + x*hornerFunction( x, a1, a2, a3);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6, double a7) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6, a7);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6, double a7, double a8) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6, a7, a8);
}

const double a0 = 2.50662823884;
const double a1 = -18.61500062529;
const double a2 = 41.39119773534;
const double a3 = -25.44106049637;
const double b1 = -8.47351093090;
const double b2 = 23.08336743743;
const double b3 = -21.06224101826;
const double b4 = 3.13082909833;
const double c0 = 0.3374754822726147;
const double c1 = 0.9761690190917186;
const double c2 = 0.1607979714918209;
const double c3 = 0.0276438810333863;
const double c4 = 0.0038405729373609;
const double c5 = 0.0003951896511919;
const double c6 = 0.0000321767881768;
const double c7 = 0.0000002888167364;
const double c8 = 0.0000003960315187;

double norminvNoCheck( double x ) {
    // We use Moro's algorithm
    double y = x - 0.5;
    if (y<0.42 && y>-0.42) {
        double r = y*y;
        return y*hornerFunction(r,a0,a1,a2,a3)/hornerFunction(r,1.0,b1,b2,b3,b4);
    } else {
        double r;
        if (y<0.0) {
            r = x;
        } else {
            r = 1.0 - x;
        }
        double s = log( -log( r ));
        double t = hornerFunction(s,c0,c1,c2,c3,c4,c5,c6,c7,c8);
        if (x>0.5) {
            return t;
        } else {
            return -t;
        }
    }
}

double norminv( double x, bool checkRange=1) {
    // Note the #include <stdexcept> at the top of the file
    if (checkRange && (x<0 || x>1.0)) {
        throw logic_error("Parameter x is out of range for norminv. It should be between 0 and 1");
    } else {
        return norminvNoCheck(x);
    }
}




int main() {
	
	cout << "Integrate sin using rectangle rule\n";
	cout << integrateSin(1,3,1000) << "\n";
	cout << "Correct answer is\n";
	cout << (cos(3)-cos(1)) << "\n";

    // question 2 effectively asks us to calculate normcdf except
    // for the factor of sqrt( 2 pi )
    double expected = infiniteIntegral(1.96);
    double actual = 0.975 * sqrt( 2 * 3.141 );
    cout << "Expected = "<<expected<<", actual="<<actual<<"\n";

}