Subversion Repositories programming

/* Ira Snyder

 * Project #2

 * Due 2005-10-28

 */

/* Copyright: Ira W. Snyder

 * Start Date: 2005-10-23

 * End Date: 2005-10-27

 * License: Public Domain

 */

#include <iostream>

#include <iomanip>

#include <cmath>

#include <cstdio>

using namespace std;

// Global Variables

enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20}; // Constants' Indices

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 };

// Global Defines

#define ARRAY_SIZE 21

// Function Prototypes

float *Bessel_BottomUp (int n, float x, float *constants);

float *Bessel_TopDown (int n, float x, float *constants);

void print_table_bottomup (float x);

void print_table_topdown (float x);

void print_equals (int num);

/* 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.

 *

 * Return values:

 * this returns an array of values, in the form:

 * array[0] = calculated value

 * array[1] = part before the subtraction: I[n-2](x)

 * array[2] = part after  the subtraction: (2(n-1)/x) * I[n-1](x)

 */

float *Bessel_BottomUp (int n, float x, float *constants)

{

    int i = 0;

    float *answer = new float[3];

    float *storage = new float[ARRAY_SIZE];

    // Prime the array

    storage[0] = constants[IDX0];

    storage[1] = constants[IDX1];

    // Calculate each value in turn

    for (i=2; i<=n; i++)

    {

        storage[i] = storage[i-2] - (2.0 * (i-1) / x) * storage[i-1];

#ifdef DEBUG

        // This will print out some useful information to find out why this

        // is going right or wrong in it's calculation.

        printf ("storage[%d] = %e - %e * %e = %e - %e = %e\n", i, storage[i-2],

                (2.0 * (i-1)  / x), storage[i-1],

                storage[i-2], (2.0 * (i-1) / x) * storage[i-1], storage[i]);

#endif

    }

    // Store the answer, delete the memory

    answer[0] = storage[n];

    answer[1] = storage[n-2];

    answer[2] = (2.0 * (n-1) / x) * storage[n-1];

    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.

 *

 * Return values:

 * this returns an array of values in the form:

 * array[0] = calculated value

 * array[1] = part before the addition: I[n+2](x)

 * array[2] = part after  the addition: (2(n+1)/x) * I[n+1](x)

 */

float *Bessel_TopDown (int n, float x, float *constants)

{

    int i = 18;

    float *answer = new float[3];

    float *storage = new float[ARRAY_SIZE];

    // Prime the array

    storage[20] = constants[IDX20];

    storage[19] = constants[IDX19];

    // Calculate each value in turn

    for (i=18; i>=n; i--)

    {

        storage[i] = (2.0 * (i+1) / x) * storage[i+1] + storage[i+2];

#ifdef DEBUG

        // This will print out some useful information to find out why this

        // is going right or wrong in it's calculation.

        printf ("storage[%d] = %e * %e + %e = %e + %e = %e\n", i,

                (2.0 * (i+1) / x), storage[i+1], storage[i+2],

                (2.0 * (i+1) / x) * storage[i+1], storage[i+2], storage[i]);

#endif

    }

    // Store the answer, free the memory

    answer[0] = storage[n];

    answer[1] = storage[n+2];

    answer[2] = (2.0 * (n+1) / x) * storage[n+1];

    delete[] storage;

    // Return the result

    return answer;

}

/* Print a table of the Bessel function (for PART #2)

 * for the given value of x */

void print_table_bottomup (float x)

{

    int n, tblwidth = 6+16+16+16+16+16+20;

    float true_val, calc_val, abs_err, rel_err, part1, part2;

    float *constants, *answer;

    // Find the right constants to use for the x value we're using.

    if (x == 1.0) constants = x1_const;

    if (x == 2.0) constants = x2_const;

    if (x == 5.0) constants = x5_const;

    // Reset the output flags

    cout << resetiosflags (ios_base::scientific | ios_base::showpoint);

    // A line of equal signs

    print_equals (tblwidth);

    // Print the x value we're using

    cout << "Using x = " << x

         << " and working Bottom-up (PART #1)" << endl;

    // A line of equal signs

    print_equals (tblwidth);

    // Set the output flags we're using

    cout << setiosflags (ios_base::scientific | ios_base::showpoint);

    // Print the table header

    cout << setw(6) << "n Val."

         << setw(16) << "TRUE VAL"

         << setw(16) << "COMPUTED"

         << setw(16) << "ABS ERROR"

         << setw(16) << "REL ERROR"

         << setw(16) << "I[n-2](x)"

         << setw(20) << "(2(n-1)/x)I[n-1](x)"

         << endl;

    // A line of equal signs

    print_equals (tblwidth);

    // Run once for n = 5,10,15,20

    for (n=5; n<=20; n+=5)

    {

        // Choose the true values

        switch (n)

        {

            case 5:  true_val = constants[IDX5];  break;

            case 10: true_val = constants[IDX10]; break;

            case 15: true_val = constants[IDX15]; break;

            case 20: true_val = constants[IDX20]; break;

        }

        // Get all of the pieces to print

        answer = Bessel_BottomUp(n, x, constants);

        calc_val = answer[0];

        abs_err  = abs (true_val - calc_val);

        rel_err  = abs_err / abs (true_val);

        part1    = answer[1];

        part2    = answer[2];

        delete[] answer;

        // Print a table line

        cout << setw(6) << n

             << setw(16) << true_val

             << setw(16) << calc_val

             << setw(16) << abs_err

             << setw(16) << rel_err

             << setw(16) << part1

             << setw(20) << part2

             << endl;

    }

    cout << endl;

}

/* Print a table of the Bessel function (for PART #2)

 * for the given value of x */

void print_table_topdown (float x)

{

    int n, tblwidth = 6+16+16+16+16+16+20;

    float true_val, calc_val, abs_err, rel_err, part1, part2;

    float *constants, *answer;

    bool done = false;

    // Choose the correct constants for the value of x we're using

    if (x == 1.0) constants = x1_const;

    if (x == 2.0) constants = x2_const;

    if (x == 5.0) constants = x5_const;

    // Reset the output flags

    cout << resetiosflags (ios_base::scientific | ios_base::showpoint);

    // A line of equal signs

    print_equals (tblwidth);

    // Print the x value we're using

    cout << "Using x = " << x

         << " and working Top-down (PART #2)" << endl;

    // A line of equal signs

    print_equals (tblwidth);

    // Set the output flags we're using

    cout << setiosflags (ios_base::scientific | ios_base::showpoint);

    // Print the table header

    cout << setw(6) << "n Val."

         << setw(16) << "TRUE VAL"

         << setw(16) << "COMPUTED"

         << setw(16) << "ABS ERROR"

         << setw(16) << "REL ERROR"

         << setw(16) << "I[n+2](x)"

         << setw(20) << "(2(n+1)/x)I[n+1](x)"

         << endl;

    // A line of equal signs

    print_equals (tblwidth);

    n = 0;

    // Run once for n = 0,1,5,10,15

    while (!done)

    {

        // Choose the correct true values of the function

        switch (n)

        {

            case 0:  true_val = constants[IDX0];  break;

            case 1:  true_val = constants[IDX1];  break;

            case 5:  true_val = constants[IDX5];  break;

            case 10: true_val = constants[IDX10]; break;

            case 15: true_val = constants[IDX15]; break;

        }

        // Get all of the values to be printed

        answer = Bessel_TopDown(n, x, constants);

        calc_val = answer[0];

        abs_err  = abs (true_val - calc_val);

        rel_err  = abs_err / abs (true_val);

        part1    = answer[1];

        part2    = answer[2];

        delete[] answer;

        // Print a line in the table

        cout << setw(6) << n

             << setw(16) << true_val

             << setw(16) << calc_val

             << setw(16) << abs_err

             << setw(16) << rel_err

             << setw(16) << part1

             << setw(20) << part2

             << endl;

        switch (n) // Set n for the next loop

        {

            case 0:  n = 1;  break;

            case 1:  n = 5;  break;

            case 5:  n = 10; break;

            case 10: n = 15; break;

            case 15: done = true; break;

        }

    }

    cout << endl;

}

/* Print a line of equal signs on the screen */

void print_equals (int num)

{

    int i;

    for (i=0; i<num; i++)

        cout << "=";

    cout << endl;

}

int main (void)

{

#ifdef DEBUG

    // Print out lots of good information to find out why these functions are

    // going right or wrong.

    float *answer;

    cout << "Bottom-up" << endl << endl;

    answer = Bessel_BottomUp(20, 1.0, x1_const);

    cout << endl << endl << "Top-down" << endl << endl;

    answer = Bessel_TopDown (0,  1.0, x1_const);

#else

    // Print the first line of output as my name

    cout << "Name: Ira Snyder" << endl << endl;

    // Print all of the tables we need

    print_table_bottomup (1.0);

    print_table_bottomup (2.0);

    print_table_bottomup (5.0);

    print_table_topdown  (1.0);

    print_table_topdown  (2.0);

    print_table_topdown  (5.0);

#endif

    return 0;

}