Applied Numerical Analysis, YDL and computers

Just want to have a nice chat with your beer? Go ahead ;-)

Applied Numerical Analysis, YDL and computers

Postby aguilarojo » 28 May 2009, 05:19

Not too long ago consumers considered as essential information regarding a computer's processing speed, it's RAM and storage capacity. Certain principles of advanced mathematics (and even fundamental arithmetic) can be used as a means of extracting useful information about computers not commonly discussed or published. One of the concepts I've found very interesting, which I believe nearly everyone can benefit from understanding is an algorithm originally discussed and explained in the book: "Applied Numerical Methods in C" by Shoichiro Nakamura

It is understood that any integer imagined is a number bearing infinite digits which signify the value of that integer. Any computer's memory can only maintain a finite representation of any number. The discrepancy between how we humans think and how computers process numbers can be considered by referring to the basic number line we all learned during our childhood and later years. Most remember how the number line ranging at least from 0 to 10 was evenly marked; some remember something of the quadrants of the entire Cartesian Plane and more. Regardless how one remembers the number line it's key mnemonic feature was it's consistent and regular markings, and subdivisions.

Dr. Nakamura explains that computers have gaps in the range of numbers they work with and these gaps are irregular and vary according to that computer's capacity to represent a number within a finite amount of bits. Unlike human number line constructs there exists within any computer between the gaps of one and another digital (and therefore approximate) representation of a number - no numerical representation of any kind. Dr. Nakamura gives this gap a name, calling it "machine epsilon".

The computer industry was aware of the problem representing numbers digitally and this led to two primarily different means of CPU design, namely Intel and PowerPC. Here are some references which discuss some details for those who would appreciate some greater depth:


The prevalence of human error has developed the nearly universal common practice multiple reviews of nearly any activity where measurement and accounting are important. Within the human world, if there is an error the measurement or assessment is done again from a different perspective or procedure. This correction process is not so straightforward with computers as the errors involved include errors built into the hardware and often bundled within the software. The odd reality is that although a computer can process millions of calculations a second only a human can determine where the error actually exists and take the necessary steps or procedures to correct or at least minimize them.

However, identifying and focusing upon machine epsilon is one of the more difficult and usually hidden causes of error. It's correction also, is not obvious.

Dr. Nakamura, explains how machine epsilon not only affects round off errors, but can also skew computations of linear and other equations which can in turn become propagated into all kinds of business and scientific reports. I recommend reading Dr. Nakamura's work to learn more details however he did produce a program in C (called eps.c, it appears on p. 12) which theoretically can be run on any computer to determine it's machine epsilon. In his book he included in a table the machine epsilon of several computers.

I modified his program only in so far that, my version:


  • is clearly labeled.
  • includes the above mentioned table.
  • includes computer's on which I ran his program.
  • outputs it's results to a file.

Code: Select all
 /* A simple form of this program appears on p. 12 of Nakamura's text on "Applied Numerical Methods in C".  Any number stored in computer memory will suffer from error due to attempting to represent a real number (with infinite decimal places) in a data space with a finite number of bits.  This program measures such error in single and double-precision values.  The numbers generated as machine epsilon is that computer's error.

    The modifications made:

        1. This version outputs to a file, res3.txt
   
        2. This version uses void functions in manipulating output to a file.
 
        3. Dr. Nakamura's table reports results for the IBM PC, IBM 370, the VAX 11
           and Cray XMP.
 
        4. This version permits the user to easily compare the current system against
           others listed in the tables mentioned above.
 
        5. Res3.c was run utilizing YDL on the Mac IIci and the G3 Minitower; res3.c
           was also run on an R/S6000 utilizing cc within AIX.
 
    Program Name: res3.c
    Programmer: Derick Centeno
    Completion Date: 2.26.2003
    Last modification: 5.57.2009
  */
 
 /*Required libraries*/
  #include<stdio.h>
  #include<stdlib.h>
 
 /*Functions used*/
  void prt_table_one(FILE *f);
  void prt_hdr_table_two(FILE *f);
  void prt_end_table_two(FILE *f);
  void prt_new_hdr(FILE *f);
 
  int main(void) {
 
  float z, g = 1, s_eps;
  double q, f = 1, d_eps;
 
  /*Declare p_outfile as a pointer so that FILE streams output to it. */
  FILE *p_outfile;
 
  /*Proper declaration ensuring opening of the file -- res3.txt*/
  p_outfile = fopen("res3.txt", "w");
 
  /*The following tests that res3.txt does open AND prints out a warning to the screen if it doesn't.*/
  if (p_outfile == NULL) {
 
     /*This message prints to screen alerting user immediately of problem.*/
      printf("ERROR: Could not open the file: \"res3.txt\"\n");
      exit(EXIT_FAILURE);
  }
 
  fprintf(p_outfile, "\nCalculation of single precision values: \n");
  prt_new_hdr(p_outfile);
 
      do {
          g = g/2;
          z = g * 0.98 + 1;
          z = z - 1;
          fprintf(p_outfile, "g = %15.8E  z = %15.8E\n", g, z);
          if (z > 0)
                    s_eps =  z;
        } while (z > 0);
 
  fprintf(p_outfile, "\nCalculation of double precision values: \n");
  prt_new_hdr(p_outfile);
 
        do {
            f = f/2;
            q = f * 0.98 + 1;
            q = q - 1;
            fprintf(p_outfile, "f = %15.8E  q = %15.8E\n", f, q);
            if (q > 0)
                     d_eps = q;
        } while (q > 0);
 
    prt_table_one(p_outfile);
    prt_hdr_table_two(p_outfile);
    fprintf(p_outfile, "Single    |       1.19E-07                 1.19E-07           %1.2E\n", s_eps);
    fprintf(p_outfile, "Double   |       2.22E-16                 2.22E-16           %1.2E\n", d_eps);
    prt_end_table_two(p_outfile);
 
    /*Close the file*/
    fclose(p_outfile);
 
    return (EXIT_SUCCESS);
    }
 
   void prt_table_one(FILE *p_outfile) {
       fprintf(p_outfile, "\n                               Machine epsilon                            \n");
       fprintf(p_outfile, "---------------------------------------------------------------------------\n");
       fprintf(p_outfile, "Precision |     IBM-PC*    IBM-370     VAX-11    Cray-XMP    Mac IIci**\n");
       fprintf(p_outfile, "===========================================================================\n");
       fprintf(p_outfile, "Single    |     1.19E-7    9.53E-7     1.19E-7   3.55E-15    1.19E-07\n");
       fprintf(p_outfile, "Double    |     2.77E-17   2.22E-16    2.77E-17  1.26E-29    1.08E-19\n");
       fprintf(p_outfile, "---------------------------------------------------------------------------\n");
       fprintf(p_outfile, "   *Quick C and Basic                **Symantec's Think C\n");
  }
 
  void prt_hdr_table_two(FILE *p_outfile) {
      fprintf(p_outfile, "\n                               Machine epsilon (continued)\n");
      fprintf(p_outfile, "---------------------------------------------------------------------------\n");
      fprintf(p_outfile, "Precision |        Mac PowerPC            IBM RS/6000       Current System\n");
      fprintf(p_outfile, "          |    233 MHz G3 Minitower        Model 59H\n");
      fprintf(p_outfile, "===========================================================================\n");
  }
 
 void prt_end_table_two(FILE *p_outfile) {
     fprintf(p_outfile, "---------------------------------------------------------------------------\n");
  }
 
 void prt_new_hdr(FILE *p_outfile) {
      fprintf(p_outfile, "     Maximum  Error       Minimum Error\n");
      fprintf(p_outfile, "========================================\n");
  }


The output of this code res3.txt follows:

Code: Select all
   
   Calculation of single precision values:
    Maximum  Error         Minimum Error
 ========================================
 g =  5.00000000E-01  z =  4.90000010E-01
 g =  2.50000000E-01  z =  2.45000005E-01
 g =  1.25000000E-01  z =  1.22499943E-01
 g =  6.25000000E-02  z =  6.12499714E-02
 g =  3.12500000E-02  z =  3.06249857E-02
 g =  1.56250000E-02  z =  1.53125525E-02
 g =  7.81250000E-03  z =  7.65621662E-03
 g =  3.90625000E-03  z =  3.82816792E-03
 g =  1.95312500E-03  z =  1.91402435E-03
 g =  9.76562500E-04  z =  9.57012177E-04
 g =  4.88281250E-04  z =  4.78506088E-04
 g =  2.44140625E-04  z =  2.39253044E-04
 g =  1.22070312E-04  z =  1.19686127E-04
 g =  6.10351562E-05  z =  5.98430634E-05
 g =  3.05175781E-05  z =  2.99215317E-05
 g =  1.52587891E-05  z =  1.49011612E-05
 g =  7.62939453E-06  z =  7.51018524E-06
 g =  3.81469727E-06  z =  3.69548798E-06
 g =  1.90734863E-06  z =  1.90734863E-06
 g =  9.53674316E-07  z =  9.53674316E-07
 g =  4.76837158E-07  z =  4.76837158E-07
 g =  2.38418579E-07  z =  2.38418579E-07
 g =  1.19209290E-07  z =  1.19209290E-07
 g =  5.96046448E-08  z =  0.00000000E+00
 
  Calculation of double precision values:
   Maximum  Error         Minimum Error
 ========================================
 f =  5.00000000E-01  q =  4.90000000E-01
 f =  2.50000000E-01  q =  2.45000000E-01
 f =  1.25000000E-01  q =  1.22500000E-01
 f =  6.25000000E-02  q =  6.12500000E-02
 f =  3.12500000E-02  q =  3.06250000E-02
 f =  1.56250000E-02  q =  1.53125000E-02
 f =  7.81250000E-03  q =  7.65625000E-03
 f =  3.90625000E-03  q =  3.82812500E-03
 f =  1.95312500E-03  q =  1.91406250E-03
 f =  9.76562500E-04  q =  9.57031250E-04
 f =  4.88281250E-04  q =  4.78515625E-04
 f =  2.44140625E-04  q =  2.39257812E-04
 f =  1.22070312E-04  q =  1.19628906E-04
 f =  6.10351562E-05  q =  5.98144531E-05
 f =  3.05175781E-05  q =  2.99072266E-05
 f =  1.52587891E-05  q =  1.49536133E-05
 f =  7.62939453E-06  q =  7.47680664E-06
 f =  3.81469727E-06  q =  3.73840332E-06
 f =  1.90734863E-06  q =  1.86920166E-06
 f =  9.53674316E-07  q =  9.34600830E-07
 f =  4.76837158E-07  q =  4.67300415E-07
 f =  2.38418579E-07  q =  2.33650208E-07
 f =  1.19209290E-07  q =  1.16825104E-07
 f =  5.96046448E-08  q =  5.84125519E-08
 f =  2.98023224E-08  q =  2.92062758E-08
 f =  1.49011612E-08  q =  1.46031380E-08
 f =  7.45058060E-09  q =  7.30156891E-09
 f =  3.72529030E-09  q =  3.65078456E-09
 f =  1.86264515E-09  q =  1.82539228E-09
 f =  9.31322575E-10  q =  9.12696141E-10
 f =  4.65661287E-10  q =  4.56348070E-10
 f =  2.32830644E-10  q =  2.28173924E-10
 f =  1.16415322E-10  q =  1.14086962E-10
 f =  5.82076609E-11  q =  5.70434810E-11
 f =  2.91038305E-11  q =  2.85218515E-11
 f =  1.45519152E-11  q =  1.42608148E-11
 f =  7.27595761E-12  q =  7.13051840E-12
 f =  3.63797881E-12  q =  3.56514818E-12
 f =  1.81898940E-12  q =  1.78257409E-12
 f =  9.09494702E-13  q =  8.91287044E-13
 f =  4.54747351E-13  q =  4.45643522E-13
 f =  2.27373675E-13  q =  2.22932783E-13
 f =  1.13686838E-13  q =  1.11466392E-13
 f =  5.68434189E-14  q =  5.57331958E-14
 f =  2.84217094E-14  q =  2.77555756E-14
 f =  1.42108547E-14  q =  1.39888101E-14
 f =  7.10542736E-15  q =  6.88338275E-15
 f =  3.55271368E-15  q =  3.55271368E-15
 f =  1.77635684E-15  q =  1.77635684E-15
 f =  8.88178420E-16  q =  8.88178420E-16
 f =  4.44089210E-16  q =  4.44089210E-16
 f =  2.22044605E-16  q =  2.22044605E-16
 f =  1.11022302E-16  q =  0.00000000E+00
 
                                 Machine epsilon
  ---------------------------------------------------------------------------
Precision |   IBM-PC*    IBM-370   VAX-11       Cray-XMP       Mac IIci**
  ===========================================================================
Single    |  1.19E-7     9.53E-7   1.19E-7        3.55E-15        1.19E-07
Double   |  2.77E-17     2.22E-16 2.77E-17        1.26E-29       1.08E-19
  -----------------------------------------------------------------------------------------------------------------------------------------------------
     *Quick C and Basic                **Symantec's Think C
 
                                 Machine epsilon (continued)
  ---------------------------------------------------------------------------------------------------------------------------------------------------
Precision  |      Mac PowerPC           IBM RS/6000     Current System
            |    233 MHz G3 Minitower     Model 59H
==========================================================================
Single    |       1.19E-07                 1.19E-07              1.19E-07
Double    |       2.22E-16                 2.22E-16              2.22E-16
 ---------------------------------------------------------------------------
                                                                                                                       


I invite you to run res3.c on your system, be it a PS3 or something else. Discover your machine epsilon. When I ran res3.c last it was on my current system running YDL 6.1 on a Titanium Powerbook G4. The machine epsilon for your system will appear under the Current System column.

I have a proposal which I believe would be interesting to the YDL community. Let interested persons run res3.c on their systems and determine the machine epsilon of their systems posting what they discover here; we can expand upon the tables already mentioned as a community allowing this data to be available to all. Going beyond running res3.c however, is a different variation of this proposal running res3.c on Cell systems under standard mode and varying arrangements of the SPEs into different configurations in order to discover whether there is any variation (an increase or decrease) of the machine epsilon. If there is, let it be documented here as to what was done to the Cell's configuration which caused the machine epsilon to change. I'm not sure this proposal is challenging enough to become a graduate thesis in itself, but as a community project aimed at exploring principles of advanced computing the concept is a way to expand each persons understanding and/or interest.

Let the fun begin...
Last edited by aguilarojo on 05 Nov 2010, 01:46, edited 7 times in total.

Everything on the Earth has a purpose.
Every disease an herb to cure it.
And every person has a mission.
This is the Indian Theory of Existence.
-- Morning Dove, Salish (1888-1936)
User avatar
aguilarojo
ydl guru
ydl guru
 
Posts: 227
Joined: 06 May 2009, 14:50
Location: New York City

Re: Applied Math, YDL and computers

Postby ppietro » 28 May 2009, 05:29

Thanks for taking the time to write that up and for your efforts to post fun and interesting articles. :)

I do have a tiny little suggestion. If you could paste your source code without the line numbers, it would be easier for us to cut and paste into editors it so that we can try it ourselves.

e.g.:

Code: Select all
/* A simple form of this program appears on p. 12 of Nakamura's text on "Applied Numerical Methods in C".  Any number stored in computer
memory will suffer from error due to attempting to represent a real number (with infinite decimal places) in a data space with a finite
number of bits.  This program measures such error in single and double-precision values.  The numbers generated as machine epsilon is
that computer's error.

The modifications made:

1. This version outputs to a file, res3.txt

2. This version uses void functions in manipulating output to a file.

3. Dr. Nakamura's table reports results for the IBM PC, IBM 370, the VAX 11 and Cray XMP.

4. This version permits the user to easily compare the current system against others listed in the tables mentioned above.

5. The res3.c was run utilizing YDL on the Mac IIci and the G3 Minitower; res3.c was also run on an R/S6000 utilizing cc within AIX.

Program Name: res3.c
Programmer: Derick Centeno
Completion Date: 2.26.2003
Last modification: 5.57.2009
*/

/*Required libraries*/
#include<stdio.h>
#include<stdlib.h>

/*Functions used*/
void prt_table_one(FILE *f);
void prt_hdr_table_two(FILE *f);
void prt_end_table_two(FILE *f);
void prt_new_hdr(FILE *f);

int main(void) {

float z, g = 1, s_eps;
double q, f = 1, d_eps;

/*Declare p_outfile as a pointer so that FILE streams output to it. */
FILE *p_outfile;

/*Proper declaration ensuring opening of the file -- res3.txt*/
p_outfile = fopen("res3.txt", "w");

/*The following tests that res3.txt does open AND prints out a warning to the screen if it doesn't.*/
if (p_outfile == NULL) {
    /*This message prints to screen alerting user immediately of problem.*/
    printf("ERROR: Could not open the file: \"res3.txt\"\n");
    exit(EXIT_FAILURE);
}

fprintf(p_outfile, "\nCalculation of single precision values: \n");
prt_new_hdr(p_outfile);

    do {
        g = g/2;
        z = g * 0.98 + 1;
        z = z - 1;
        fprintf(p_outfile, "g = %15.8E  z = %15.8E\n", g, z);
        if (z > 0)
            s_eps =  z;
    } while (z > 0);

fprintf(p_outfile, "\nCalculation of double precision values: \n");
prt_new_hdr(p_outfile);

    do {
        f = f/2;
        q = f * 0.98 + 1;
        q = q - 1;
        fprintf(p_outfile, "f = %15.8E  q = %15.8E\n", f, q);
        if (q > 0)
            d_eps = q;
    } while (q > 0);

prt_table_one(p_outfile);
prt_hdr_table_two(p_outfile);
fprintf(p_outfile, "Single    |       1.19E-07                 1.19E-07           %1.2E\n", s_eps);
fprintf(p_outfile, "Double    |       2.22E-16                 2.22E-16           %1.2E\n", d_eps);
prt_end_table_two(p_outfile);

/*Close the file*/
fclose(p_outfile);

return (EXIT_SUCCESS);
}

void prt_table_one(FILE *p_outfile) {
    fprintf(p_outfile, "\n                               Machine epsilon                            \n");
    fprintf(p_outfile, "---------------------------------------------------------------------------\n");
    fprintf(p_outfile, "Precision |     IBM-PC*    IBM-370     VAX-11    Cray-XMP    Mac IIci**\n");
    fprintf(p_outfile, "===========================================================================\n");
    fprintf(p_outfile, "Single    |     1.19E-7    9.53E-7     1.19E-7   3.55E-15    1.19E-07\n");
    fprintf(p_outfile, "Double    |     2.77E-17   2.22E-16    2.77E-17  1.26E-29    1.08E-19\n");
    fprintf(p_outfile, "---------------------------------------------------------------------------\n");
    fprintf(p_outfile, "   *Quick C and Basic                **Symantec's Think C\n");
}

void prt_hdr_table_two(FILE *p_outfile) {
    fprintf(p_outfile, "\n                               Machine epsilon (continued)\n");
    fprintf(p_outfile, "---------------------------------------------------------------------------\n");
    fprintf(p_outfile, "Precision |        Mac PowerPC            IBM RS/6000       Current System\n");
    fprintf(p_outfile, "          |    233 MHz G3 Minitower        Model 59H\n");
    fprintf(p_outfile, "===========================================================================\n");
}

void prt_end_table_two(FILE *p_outfile) {
    fprintf(p_outfile, "---------------------------------------------------------------------------\n");
}

void prt_new_hdr(FILE *p_outfile) {
    fprintf(p_outfile, "     Maximum  Error       Minimum Error\n");
    fprintf(p_outfile, "========================================\n");
}


Cheers,
Paul
User avatar
ppietro
Site Admin
Site Admin
 
Posts: 4965
Joined: 13 Sep 2007, 22:18

Re: Applied Math, YDL and computers

Postby aguilarojo » 28 May 2009, 17:03

Thanks for your kind comments, Paul.

By the way, I implemented your request and hopefully added an interesting twist for those who have access to a Cell or node of Cells.
Hopefully others find what I've written as interesting as you have.

All the best...

Everything on the Earth has a purpose.
Every disease an herb to cure it.
And every person has a mission.
This is the Indian Theory of Existence.
-- Morning Dove, Salish (1888-1936)
User avatar
aguilarojo
ydl guru
ydl guru
 
Posts: 227
Joined: 06 May 2009, 14:50
Location: New York City


Return to Speaker's Corner

Who is online

Users browsing this forum: No registered users and 4 guests