Subversion Repositories programming

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20};

#define ARRAY_SIZE 21
float x1_const[7] = { 1.266065978,     5.651591040E-01, 2.714631560E-04,
                      2.752948040E-10, 2.370463051E-17, 1.587678369E-23,
                      3.966835986E-25 };

float x2_const[7] = { 2.279585302,     1.590636855,     9.825679323E-03,
                      3.016963879E-07, 8.139432531E-13, 8.641603385E-18,
                      4.310560576E-19 };

float x5_const[7] = { 2.723987182E+01, 2.433564214E+01, 2.157974547,
                      4.580044419E-03, 1.047977675E-6 , 4.078415017E-10,
                      5.024239358E-11 };

/* Copied from the Project #2 Handout.
 *
 * FOR TESTING ONLY
 */
float Bessel_Recursive_BU (int n, float x, float base0, float base1)
{
    if (n == 0)
        return base0;
    else if (n == 1)
        return base1;

    return Bessel_Recursive_BU (n-2, x, base0, base1)
           + 2.0F * (float)(n-1) * Bessel_Recursive_BU (n-1, x, base0, base1) / x;
}

float Bessel_Recursive_TD (int n, float x, float base19, float base20)
{
    if (n == 19)
        return base19;
    else if (n == 20)
        return base20;

    return (2.0F * (float)(n+1) / x) * Bessel_Recursive_TD (n+1, x, base19, base20)
           + Bessel_Recursive_TD (n+2, x, base19, base20);
}


/* Find the Bessel Function.
 * Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
 * by simple substitution.
 *
 * This is the recurrence relation for PART #1.
 */
float Bessel_BottomUp (int n, float x, float base0, float base1)
{
    int i = 0;
    float  answer;
    float *storage = new float[ARRAY_SIZE];

    // Prime the array
    storage[0] = base0;
    storage[1] = base1;

    // Calculate each value in turn
    for (i=2; i<=20; i++)
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];

    // Store the answer, delete the memory
    answer = storage[n];
    delete[] storage;

    // Return the result
    return answer;
}

/* Find the Bessel Function, from the Top -> Down.
 * Algorithm: I[n-1](x) = (2*n/x)*I[n](x) + I[n+1](x)
 * Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)
 * by simple substitution.
 *
 * This is the recurrence relation for PART #2.
 */
float Bessel_TopDown (int n, float x, float base19, float base20)
{
    int i = 18;
    float answer;
    float *storage = new float[ARRAY_SIZE];

    // Prime the array
    storage[20] = base20;
    storage[19] = base19;

    // Calculate each value in turn
    for (i=18; i>=0; i--)
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];

    // Store the answer, free the memory
    answer = storage[n];
    delete[] storage;

    // Return the result
    return answer;
}

/* Print out a table for a given x value, and function.
 * Call like this: print_table (1.0, &BesselFunc);
 *
 * It chooses the correct constants appropriately.
 */
void print_table (float x, float (*bessel_func)(int, float, float, float))
{
    int n;
    float base0, base1, base19, base20;
    float true_val, calc_val, abs_err, rel_err;
    float *constants;

    if (x == 1.0)
        constants = x1_const;

    if (x == 2.0)
        constants = x2_const;

    if (x == 5.0)
        constants = x5_const;

    // Set the output flags
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);

    // Print the table header
    cout << setw(10) << "n Val."
         << setw(20) << "TRUE"
         << setw(20) << "COMPUTED"
         << setw(20) << "ABS ERROR"
         << setw(20) << "REL ERROR"
         << endl;

    // A line of equal signs
    for (n=0; n<(20+20+20+20+10); n++)
        cout << "=";

    cout << endl;

    // Print out the line for n = 5
    n = 5;
    true_val = constants[IDX5];
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
    abs_err = abs (true_val - calc_val);
    rel_err = abs_err / abs (true_val);

    cout << setw(10) << n
         << setw(20) << true_val
         << setw(20) << calc_val
         << setw(20) << abs_err
         << setw(20) << rel_err
         << endl;

    // Print out the line for n = 10
    n = 10;
    true_val = constants[IDX10];
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
    abs_err = abs (true_val - calc_val);
    rel_err = abs_err / abs (true_val);

    cout << setw(10) << n
         << setw(20) << true_val
         << setw(20) << calc_val
         << setw(20) << abs_err
         << setw(20) << rel_err
         << endl;

    // Print out the line for n = 15
    n = 15;
    true_val = constants[IDX15];
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
    abs_err = abs (true_val - calc_val);
    rel_err = abs_err / abs (true_val);

    cout << setw(10) << n
         << setw(20) << true_val
         << setw(20) << calc_val
         << setw(20) << abs_err
         << setw(20) << rel_err
         << endl;

    // Print out the line for n = 20
    n = 20;
    true_val = constants[IDX20];
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
    abs_err = abs (true_val - calc_val);
    rel_err = abs_err / abs (true_val);

    cout << setw(10) << n
         << setw(20) << true_val
         << setw(20) << calc_val
         << setw(20) << abs_err
         << setw(20) << rel_err
         << endl;
}

int main (void)
{
    print_table (1.0, &Bessel_BottomUp);
    cout << endl;
    print_table (2.0, &Bessel_BottomUp);
    cout << endl;
    print_table (5.0, &Bessel_BottomUp);
    cout << endl;


    return 0;
}