// Arup Guha
// 1/13/2012
// COT 4500 - Numerical Calculus

// Illustration of textbook example of rounding and chopping on page 27.

// Known bugs: The Horner's method results do NOT agree with the text.
//             I am not sure whether the problem is in the chop and rnd
//             functions or the horner's functions.

// Body of the code still needs to be commented.

#include <stdio.h>
#include <math.h>

double chop(double input, int precision);
double rnd(double input, int precision);
double eval_horners_chop(double poly[], int length, double x, int digits);
double eval_horners_rnd(double poly[], int length, double x, int digits);
double reg_eval_chop(double poly[], int length, double x, int digits);
double my_pow_rnd(double base, int exp, int digits);
double reg_eval_rnd(double poly[], int length, double x, int digits);

int main() {

    double poly[] = {1.5, 3.2, -6.1, 1};
    printf("horner chop = %lf\n", eval_horners_chop(poly, 4, 4.71, 3));
    printf("horner rnd = %lf\n", eval_horners_rnd(poly, 4, 4.71, 3));
    printf("reg chop = %lf\n", reg_eval_chop(poly, 4, 4.71, 3));
    printf("reg rnd = %lf\n", reg_eval_rnd(poly, 4, 4.71, 3));
    return 0;
}

double chop(double input, int precision) {

    int exponent = 0;

    if (fabs(input) >= 1) {
        while (fabs(input) >= 1) {
            input /= 10;
            exponent++;
        }
    }
    else if (fabs(input) < .1 && input != 0) {
        while (fabs(input) < .1) {
            input *= 10;
            exponent--;
        }
    }

    int digits = (int)(input*pow(10, precision));

    return (digits)/pow(10, precision - exponent);
}

double rnd(double input, int precision) {

    int exponent = 0;

    if (fabs(input) >= 1) {
        while (fabs(input) >= 1) {
            input /= 10;
            exponent++;
        }
    }
    else if (fabs(input) < .1 && input != 0) {
        while (fabs(input) < .1) {
            input *= 10;
            exponent--;
        }
    }

    int digits = (int)(input*pow(10, precision));
    double leftover = input*pow(10, precision) - digits;

    if (leftover < .5)
        return (digits)/pow(10, precision - exponent);
    else
        return (digits+1)/pow(10, precision - exponent);
}

double eval_horners_chop(double poly[], int length, double x, int digits) {

    int i;
    double ans = 0;
    for (i=length-1; i>=0; i--) {

        ans = chop(chop(x, digits)*chop(ans, digits),digits) + chop(poly[i], digits);

    }

    return chop(ans, digits);
}

double eval_horners_rnd(double poly[], int length, double x, int digits) {

    int i;
    double ans = 0;
    for (i=length-1; i>=0; i--) {
        ans = rnd(rnd(x,digits)*rnd(ans,digits),digits) + rnd(poly[i],digits);
    }

    return rnd(ans, digits);
}

double my_pow_chop(double base, int exp, int digits) {

    double ans = 1;
    int i;
    for (i=0; i<exp; i++) {
        ans = chop(chop(ans, digits)*chop(base, digits), digits);
    }

    return chop(ans, digits);
}

double my_pow_rnd(double base, int exp, int digits) {

    double ans = 1;
    int i;
    for (i=0; i<exp; i++) {
        ans = rnd(rnd(ans, digits)*rnd(base, digits), digits);
    }

    return rnd(ans, digits);
}

double reg_eval_chop(double poly[], int length, double x, int digits) {

    double sum = 0;
    int i;
    for (i=0; i<length; i++) {
        sum = sum + chop(chop(poly[i], digits)*my_pow_chop(x, i, digits), digits);
    }

    return chop(sum, digits);
}

double reg_eval_rnd(double poly[], int length, double x, int digits) {

    double sum = 0;
    int i;
    for (i=0; i<length; i++) {
        sum = sum + rnd(rnd(poly[i], digits)*my_pow_rnd(x, i, digits), digits);
    }

    return rnd(sum, digits);
}