Numerical Recipes

Numerical Recipes The Art of Scientific Computing (with applications in computational physics)

Computer Architecture CPU Memory von Neumann architecture: both data and instructions are stored in memory. External Storage

Program Organization int main(void) { … } double func(double x) { … }

First Example #include <stdio.h> int main(void) { printf(”hello, world\n”); return 0; } gcc hello.c (to get a.out) [Or other way depending on your OS]

Data Types #include <stdio.h> int main(void) { int i,j,k; double a,b,f; char c, str[100]; j = 3; a = 1.05; c = ’a’; str[0] = ’p’; str[1] = ’c’; str[2] = ’\0’; printf(”j=%d, a=%10.6f, c=%c, str=%s\n”, j, a, c, str); }

Equal is not equal • x = 10 is not 10 = x • x = x + 1 made no sense if it is math • x = a + b OK, but a+b = x is not C. • In general, left side of ‘=’ refers to memory location, right side expression can be evaluated to numerical values

Expressions • Expressions can be formed with +, -, *, / with the usual meaning • Use parentheses ( …) if meaning is not clear, e.g., (a+b)*c • Be careful 2/3 is 0, not 0.666… • Large number of operators exists in C, such as ++i, --j, a+=b, &, !, &&, ||, ^, ?a:b, etc

Use a Clear Style (A) k=(2-j)*(1+3*j)/2; k=j+1; if(k == 3) k=0; switch(j) { case 0: k=1; break; case 1: k=2; break; case 2: k=0; break; default: { fprintf(stderr, ”unexpected value for j”); exit(1); } } (D) k=(j+1)%3; (B) (C)

Control Structures in C - loop for (j=0; j < 10; ++j) { a[j] = j; } while (n < 1000) { n *= 2; } do { n *= 2; } while (n < 1000);

Control Structure - conditional if (b > 3) { a = 1; } if (n < 1000) { n *= 2; } else { n = 0; }

Control Structure - break for( ; ; ) { ... if(. . .) break; }

Pointers • Pointer is a variable in C that stores address of certain type • int *p; double *ap; char *str; • You make it pointing to something by (1) address operator &, e.g., p = &j, (2) malloc() function, (3) or assignment, str = ”abcd” • Get the value the pointer is pointing to by dereferencing, *p

1D Array in C • int a[4]; defines elements a[0],a[1],a[2], and a[3] • a[j] is same as *(a+j),a has a pointer value • double b[4], *bb; bb=b-1; then valid range of index for b is from 0 to 3, but bb is 1 to 4

1D Array Argument Passing void routine(double bb[], int n) // bb[1..n] (range is 1 to n) • We can use as double a[4]; routine(a-1, 4);

2D Array in C double m[3][5]; defines fixed size array. Example below defines dynamic 2D array. double **m; m = (double **) malloc(3*sizeof(double *)); for(i=0; i<3; ++i) { m[i] = (double *)malloc(5*sizeof(double)); }

Representation of 2D Array, fixed size vs dynamic

Special Treatment of Array in NR • float *vector(long nl, long nh) allocates a float vector with index [nl..nh] • float **matrix(long nrl, long nrh, long ncl, long nch) allocates a 2D matrix with range [nrl..nrh] by [ncl..nch]

Call by Value vs Call by Reference • In a function call, a copy of the value is passed to a function argument. E.g.: void func(int a, int b[]){ a = a + 1; b[0] = b[0] + 1; } int main(void) { int i, A[5]; i = 1; A[0] = 1; func(i, A); i = ?; A[0] = ?; }

Header Files in NR #include ”nr.h” #include ”nrutil.h” Distinction of file names in <….> and in ”… ”

Precedence and Association

Pre/post Increment, Address of • Consider f(++i) vs f(i++), what is the difference? • &a vs *a • Conditional expression x = (a < b) ? c : d;

Macros in C #define DEBUG #define PI 3.141592653 #define SQR(x) ((x)*(x))

Computer Representation of Numbers • Unsigned or two’s complement integers (e.g., char) 1000 0000 = 128 or -128 1000 0001 = 129 or -127 1000 0010 = 130 or -126 1000 0011 = 131 or -125 . . . 1111 1100 = 252 or -4 1111 1101 = 253 or -3 1111 1110 = 254 or -2 1111 1111 = 255 or -1 0000 0000 = 0 0000 0001 = 1 0000 0010 = 2 0000 0011 = 3 0000 0100 = 4 0000 0101 = 5 0000 0110 = 6 . . . 0111 1111 = 127

Real Numbers on Computer Example for β=2, p=3, emin= -1, emax=2 ε 0 0.5 1 2 3 4 5 6 7 ε is called machine epsilon.

Floating Point, sMBe-E, not IEEE

IEEE 754 Standard (32-bit) • The bit pattern represents If e = 0: (-1)sf 2-126 If 0<e<255: (-1)s (1+f) 2e-127 If e=255 and f = 0: +∞ or -∞ and f ≠ 0: NaN … … s b-1 b-2 b-23 e f = b-12-1 + b-22-2 + … + b-232-23

Error in Numerical Computation • Integer overflow • Round off error • E.g., adding a big number with a small number, subtracting two nearby numbers, etc • How does round off error accumulate? • Truncation error (i.e. discretization error) • The field of numerical analysis is to control truncation error

(Machine) Accuracy • Limited range for integers (char, int, long int, and long long int) • Limited precision in floating point. We define machine ε as such that the next representable floating point number is (1 + ε) after 1. ε 10-7 for float (32-bit) and  10-15 for double (64-bit)

Stability • An example of computing Φnwhere • We can compute either by Φn+1 = ΦnΦ or Φn+1 = Φn-1 – Φn Results are shown in stability.c program

Reading Materials • “Numerical Recipes”, Chap 1. • “What every computer scientist should know about floating-point arithmetic”. Can be downloaded from http://perso.ens-lyon.fr/jean-michel.muller/goldberg.pdf • “The C Programming Language”, Kernighan & Ritchie

Problems for Lecture 1 (C programming, representation of numbers in computer, error, accuracy and stability, assignment to be handed in next week thursday) 1. (a) An array is declared as char *s[4] = {”this”, ”that”, ”we”, ”!”}; What is the value of s[0][0], s[0][4], and s[2][1]? (b) If the above array is to be passed to a function, how should it be used, i.e., the declaration of the function, and use of the function in calling program? If the array is declared as fixed size: char t[4][5]; instead, then how should it be passed to a function? 2. (a) Study the IEEE 754 standard floating point representation for 32-bit single precision numbers (float in C) and write out the bit-pattern for the float numbers 0.0, 1.0, 0.1, and 1/3. (b) For the single precision floating point representation (32-bit number), what is the exact value of machine epsilon? What is the largest positive number? 3. For the recursion relation: F n+1 =Fn-1 – Fn with F0 and F1 arbitrary, find the general solution Fn. Based on its solution, discuss why is it unstable for computing the power of golden mean Φ? (Hint: consider trial solution of the form Fn = Arn).

Numerical Recipes