Algorithms for

Scientic Computation
(An Introduction Using C)
Division of Applied Mathematics
Korea Advanced Institute of Science and Technology
2004
1
Contents
1 Introduction 4
2 Basic Data Types 8
2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Int . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Float, Double, and Long Double . . . . . . . . . . . . . . . . . . 11
2.4 Char . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Arithmetic operators 15
4 Input and Output Functions 18
4.1 printf() and scanf() statements . . . . . . . . . . . . . . . . . . . 18
4.2 File Input and Output . . . . . . . . . . . . . . . . . . . . . . . . 20
5 Control Flow Statements 26
5.1 Conditional Expression . . . . . . . . . . . . . . . . . . . . . . . . 26
5.2 if Statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.3 switch Statements . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.4 while Statements . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.5 for Statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.6 Jump Statements . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6 Function statements 34
6.1 Local and Global Variables . . . . . . . . . . . . . . . . . . . . . 35
6.2 Storage Classes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.3 Recursive Functions . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.4 Macros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
7 Arrays 41
7.1 Arrays and Pointers . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.2 Dynamic Memory Allocation Function . . . . . . . . . . . . . . . 43
7.3 Multi-Dimensional Arrays and Pointers. . . . . . . . . . . . . . . 44
8 Strings 46
8.1 Input and Output functions of Strings . . . . . . . . . . . . . . . 46
8.2 String Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
8.3 Command Line Arguments . . . . . . . . . . . . . . . . . . . . . 50
9 Structures 51
10 Linear Linked Lists 55
10.1 basic list operations . . . . . . . . . . . . . . . . . . . . . . . . . 56
10.2 Stacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
10.3 Queues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
2
11 Trees 64
12 Sorting Methods 69
12.1 Insertion Sort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
12.2 Quicksort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
12.3 Heapsort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
13 Searching Methods 81
13.1 linear and binary search methods . . . . . . . . . . . . . . . . . . 81
13.2 binary search tree . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
14 Hashing 91
14.1 Linear Probing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
14.2 Double Hashing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
14.3 Chaining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
15 Root Finding Methods 99
15.1 Root Finding Based on Search Methods . . . . . . . . . . . . . . 99
15.2 Root Finding Based on Target Functions . . . . . . . . . . . . . . 101
16 Numerical Integration 106
16.1 Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
16.2 Simpsons Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
17 Numerical Solutions of
Ordinary Dierential Equations 112
17.1 Eulers Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
17.2 Runge-Kutta method . . . . . . . . . . . . . . . . . . . . . . . . . 113
18 Numerical Solution of
Simultaneous Equations 117
18.1 simple Gauss-Jordan Elimination . . . . . . . . . . . . . . . . . . 117
18.2 LU Decomposition Methods . . . . . . . . . . . . . . . . . . . . . 121
18.3 Solutions of m equations with n unknowns . . . . . . . . . . . . . 125
19 Random Number Generation 129
19.1 Creating Uniform Random Numbers . . . . . . . . . . . . . . . . 129
19.2 Creating Non-uniform Random Numbers . . . . . . . . . . . . . . 134
3
1 Introduction
Computer Systems : information processing
Structure
why C language?
- Dennis Ritchie, 1973 (AT&T Bell lab) for Unix O.S.
- American National Standards Institute(ANSI) : 1989, ANSI C
- Pros :
portability
general purpose
compact code
structured programming
exible control structures. mid-level language
4
C programming(in Unix environment)
Structure of C program :
#header - preprocessor
main()
{
statements;
.
.
.
}
function a()
{
statements;
.
.
.
}
function b()
{
statements;
.
.
.
}
5
5 types of C statements:
declaration statements.
assignment statements.
function statements.
control statements.
null statements.
Problem solving methodology
1. problem statement
give a clear and concise problem statement.
(to avoid any misunderstanding)
- eg. compute the straight-line distance between two points in a plane.
2. input/ output description
the I/O diagram.
- eg.
3. hand example
simple set of data for testing the program.
4. algorithm development
algorithm : a step-by-step outline of the problem solution
- eg.
(a) give values to the two points.
(b) compute the lengths of the two sides.
(c) compute the distance between the two points using the Pythagorean
theorem.
(d) print the distance between the two points.
Pseudo-code
5. testing
Compare with the hand solution.
6
Sample Program
- Pseudo-code
1. read (x1, y1) and (x2, y2)
2. side_1 = x1 - x2, side_2 = y1 - y2
3. distance = sqrt (side_1*side_1 + side_2*side_2)
4. print distance
- Sample Program
/*---------------------------------------------------------------*/
/* This program computes the distance between two points. */
/* */
#include <stdio.h>
#include <stdlib.h>
#incldue<math.h>
main()
{
/* Declare and initialize variables. */
float x1, y1, x2, y2, side_1, side_2, distance;
/* Read two points. */
printf("two points: x1 y1 x2 y2 = ");
scanf("%f %f %f %f", &x1, &y1, &x2, &y2);
/* Compute two sides. */
side_1 = x1 - x2;
side_2 = y1 - y2;
/* Compute the distance */
distance = sqrt(side_1*side_1 + side_2*side_2);
/* Print the distance */
printf("distance = %6.2f\n", distance);
}
/*---------------------------------------------------------------*/
7
2 Basic Data Types
2.1 Data
Data
data
_

_
characters : char
numbers
_

_
integers : int(short, int, long)
real numbers : oat, double, long double
Qualiers - short, long, unsigned
- int : int, short(int), long(int), unsigned(int),
unsigned short, unsigned long
- double : long double
Number of bytes for data types : sizeof operator
printf("short has %d bytes\n", sizeof(short));
short has 2 bytes
Example program
/* ...*/ - comments
#include<stdio.h> - input/output functions
(< > : standard C library)
<stdlib.h> - constants for exiting the program
<math.h> - mathematical functions
(.h : header files)
main()
{...} - block
statements : end with ;
- declaration statements : data_type list of variables;
- assignment statements : identifier = expression;
- control flow statements
- function statements
- null statements
8
The declaration statements dene memory locations for variables.
eg.
short x, y=1; - short : data_type.
- x, y : variables.
- y=1 : initialization.
...
?
!
... ...
address leld
conlenl leld
memory ...
x
y
! byle8bls)
2 byles lor shorl
Selection of the name of variables.
- small letter, capital letter, underline(

), number.
- number cannot come rst.
- maximum number of characters(831).
- no reserved words.
eg.
valid not valid
wiggly $z**
cat1 1cat
Hot tub Hot-tub
KcaB dont
9
2.2 Int
int constants
- decimal number : 123 {0,1,2,...,9}
- octal number : O123 {0,1,2,...,7}
- hexadecimal number : Ox123 {0,...,9,a,b,...,f}
range of values
short (2) 2
15
2
15
1 (32768 32767)
unsigned short (2) 0 2
16
1 (0 65535)
int (4) 2
31
2
31
1
unsigned int (4) 0 2
32
1
eg. unsigned (range : 0 15)
binary number decimal number
0000 0
0001 1
0010 2
.
.
.
.
.
.
1111 15
eg. int (range : 8 7)
2s complement method for negative integers.
x = x + 1
where, x is a binary number.
binary number decimal number
0111 7
.
.
.
.
.
.
0001 1
0000 0
1111 -1
.
.
.
.
.
.
1000 -8
eg. range of short(2bytes)(2
15
2
15
1)
binary number 1000 0000 0000 0000 0111 1111 1111 1111
hexa-decimal number 8 0 0 0 7 f f f
10
2.3 Float, Double, and Long Double
constants :
1. 2.75 3.14 1.6e 19 (= 1.6 10
19
)
- Always be stored double form.
binary number (IEEE 754 standard) - oat(32bits)
!) 8) 23)
snqle bl exponenl parl lraclon parl
eg. the precision digits and the range of exponents for oats (32bits)
i. valid numbers
- fraction part (23bits)
- range of values : 2
0
2
1
2
22

1
2
22
=
1
(2
10
)
2.2

=
1
(10
3
)
2.2
= 10
6.6
= precision digits = 6 (decimal digits)
ii. range of exponents
- exponent part (8bits)
- range of values : 0 2
8
1(255)
- bias : 2
7
1(127) = 127 128 i.e. 2
127
2
128
- transform to the exponent of 10 : 2
128
= (2
10
)
12.8
= (10
3
)
12.8
= 10
38.4
= the range of the exponent of 10 : 38 38
11
The limitations of data type (in SUNWS)
/*---------------------------------------------------------------*/
/* This program prints the system limitations. */
/* */
#include <stdio.h>
#include <limits.h>
#include <float.h>
main() {
/* byte size for each data type */
printf("short has %d bytes\n", sizeof(short));
printf("int has %d bytes\n", sizeof(int));
printf("long has %d bytes\n\n", sizeof(long));
printf("float has %d bytes\n", sizeof(float));
printf("double has %d bytes\n", sizeof(doluble));
printf("long double has %d bytes\n\n", sizeof(long double));
/* Print integer type maximums. */
printf("short maximum: %d\n", SHRT_MAX);
printf("int maximum: %d\n", INT_MAX);
printf("long maximum: %ld\n\n", LONG_MAX);
/* Print float precision, range, and maximum. */
printf("float precision digits: %d\n", FLT_DIG);
printf("float maximum exponent: %d\n", FLT_MAX_10_EXP);
printf("float maximum: %e\n\n", FLT_MAX);
/* Print double precision, range, and maximum */
printf("double precision digits: %d\n", DBL_DIG);
printf("double maximum exponent: %d\n", DBL_MAX_10_EXP);
printf("double maximum: %e\n\n", DBL_MAX);
/* Print long precision, range, and maximum */
printf("long double precision digits: %d\n", LDBL_DIG);
printf("long double maximum exponent: %d\n", LDBL_MAX_10_EXP);
printf("long double maximum: %e\n\n", LDBL_MAX);
}
/*---------------------------------------------------------------*/
12
short has 2 bytes
int has 4 bytes
long has 4 bytes
float has 4 bytes
double has 8 bytes
long double has 16 bytes
short maximum : 32767
int maximum : 2147483647
long maximum :2147483647
float precision digits : 6
float maximum exponent : 38
float maximum : 3.402823e+38
double precision digits : 15
double maximum exponent : 308
double maximum : 1.797693e+308
long double precision digits : 33
long double maximum exponent : 4932
long double maximum : 1.189731e+4932
13
2.4 Char
Express by ASCII code
(7 or 8 bits)
4bls 4bls
zone
bls
dql
bls
bl bnary dql
eg. 7 bit ASCII
character binary number decimal number
A 1000001 65
B 1000010 66
.
.
.
.
.
.
.
.
.
P 1010000 80
.
.
.
.
.
.
.
.
.
Z 1011010 90
- constants :
A B C a - single quotation.
cf. 1 : 0000001
1 : 0110001
- Non-existing character in the keyboard.
Use \
eg.
char beep = \007
\n \012 (newline)
\t \011 (tab)
\b \010 (backspace)
\r \015 (carriage return)
\ ~~ \\
~~ \
" ~~ \"
14
3 Arithmetic operators
arithmetic operators
_

_
unary operators (one operand) : =, (data_type), -, +
binary operators (two operands)
_

_
additive operators : +, -
multiplicative operators : *, /, %
= (assignment)
fahr = 0; (correct)
0 = fahr; (wrong)
The direction of assignment is from left to right.
- (sign change)
x = -1;
One or more space between = and -1.
+, -, *
x = a+b;
/
- type of integer numbers : 5/2 2 (quotient)
- type of real numbers : 5./2 2.5
% (remainder, modulus)
3%2 1, 8%2 0
(data type) cast operator
eg.
int i, j,;
i = 1.6+1.7; (i=3)
j = (int)1.6+(int)1.7; (j=2)
15
increment and decrement operators : ++, --
The priority is same as unary operators.
- prefix operators
x = ++n : n = n+1; x = n;
x = --n : n = n-1; x = n;
- postfix operators
x = n++ : x = n; n = n+1;
x = n-- : x = n; n = n-1;
eg.
int x, n=1;
x = n++; (x=1, n=2)
x = ++n; (x=3, n=3)
abbreviated assignment operators : op=,
The priority is same as = operators.
- op : +, -, *, /, % and bitwise operators
eg.
x += 3 (x = x+3)
x *= y+1 (x = x*(y+1))
the general form of assignment statement: identifier = expression;
- expression: a constant, another variable, or the results of an operation.
eg.
5
4+20
c=3+8
6+(c=3+8)
16
Priority of operators
- basic rule:
unary operators except =

> binary operators > . . .
operators Associativity
(when same priority)
( )
-, (data_type), ++, --
*, /, %
+, -
=
Priority of operands
- Basic rule : Change to the many number of bytes operands.
(double > oat > long > int > short)
eg.
int i, j;
float x, y;
i = 5./2*5; (i=12)
x = 5./2*5; (x=12.5)
j = 5/2*5; (j=10)
y = 5/2*5; (y=10)
eg.
f =
x
3
2x
2
+x 6.3
x
2
+ 0.05005x 3.14
f = (x*x*x-2*x*x+x-6.3)/(x*x+0.05005*x-3.14);
17
4 Input and Output Functions
4.1 printf() and scanf() statements
printf() functions
printf("control string", arguments);
cf. Arguments can be omitted.
eg.
printf("output");
output
printf("%c %d\n", A, 10);
A 10
%c, %d - numeric conversion specifier
variable-type specier meaning
character %c character
%s string
integer %d decimal no.
%o octal no.
%x hexadecimal no.
%u unsigned no.
oating-point %f oating point type
values %e exponent type
%g shorter between %f and %e
cf. qualiers - short : h, long : l, double : l, long double : L
- The number of conversion speciers should be the number of outputs.
eg.
printf("%c %d\n", 65, 65);
A 65
printf("%d", A-8+5);
4
printf("%d %o %x", 10, 10, 10);
10, 12, a
printf("%f %e\n", 2.4, 2.4);
2.400000 2.400000e+00
printf("%s", string);
string
18
- user dened output forms
integer real number string
%nd or %n.mf or %n.ms or
%-nd %-n.mf %-n.ms
n : no. of spaces
- : starting from left
m(in real number) : no of floating points
m(in string) : first m characters
eg.
printf("%d\n", 126); : 126
printf("%10d\n", 126); : #######126
printf("%-10d\n", 126); : 126#######
eg.
printf("%f\n", 1234.56); : 1234.56####
printf("%e\n", 1234.56); : 1.23456#e+03
printf("%3.1f\n", 1234.56); : 1234.6
printf("%10.3f\n", 1234.56); : ##1234.560
printf("%10.3e\n", 1234.56); : #1.235e+03
eg.
printf("%2s\n", "string"); : string
printf("%10s\n", "string"); : ####string
printf("%-10.5s\n", "string"); : strin#####
19
scanf() functions
scanf("control string", arguments);
- Arguments should be memory addresses of variables.
- Conversion specifiers are same as printf().
eg.
char ch;
scanf("%c %d", &ch, &i);
cf. & operator: returns the memory address of assigned variable.
4.2 File Input and Output
unix command :
a.out <in.file> out.file
declaration of data type of le :
FILE *fp; (*fp : file pointer - the variable that indicate
the memory address of file)
fopen() and fclose() functions
fopen("filename", "r");
r : read
w : write
a : append
- fopen() returns the memory address of le.
20
eg.
FILE *in;
in = fopen("test". r");
fclose(in) : fclose(file-pointer)
- notifying the end of a task to the file
Input, Output function related to les.
- character
char ch; FILE *fp;
ch = getchar();
ch = getc(fp); : read one character in the FILE
which is indicated by fp.
ch = getc(stdin) : stdin - keyboard
putchar(ch);
putc(ch, fp); : write ch to the fp
putc(ch, stdout); : stdout - screen
- string
char str[100]; FILE *fp;
gets(str);
fgets(str, max, fp); : read the first max characters in
a string of the FILE indicated by fp
and assign to str.
puts(str);
fputs(str, fp); : write the string to the FILE
indicated by fp
- printf & scanf
scanf("control-string", argument_list);
fscanf(fp, "control-string", argument_list); : read from the FILE
indicated by fp
printf("control-string", argument_list);
fprintf(fp, "control-string", argument_list) : write to the FILE
indicated by fp
21
eg.
main()
{
FILE *fp;
int age;
fp=fopen("sam", "r");
fscanf(fp, "%d", &age);
fclose(fp);
fp=fopen("sam", "w");
fprintf(fp, "sam is %d years old\n", age);
fclose(fp)
}
eg.
if(fp==NULL)
printf("error");
else{
...
Programming Linear Modeling (Regression)
the process that determines the linear equation that is the best t to a
set of data points in terms of minimizing the sum of the squared distances
between the line and the data points.
For the given sample :
(x
1
, y
1
), (x
2
, y
2
), . . . , (x
n
, y
n
)
linear model(estimator) :
y(x) = mx +b
problem : for the given sample and estimator, nd m and b which mini-
mizes
E =
1
2
n

i=1
(y
i
y(x
i
))
2
.
22
If y is linear, E has a quadratic form.
_
E
m
_
m=m

= 0,
_
E
b
_
b=b

= 0
m

n
k=1
x
k

n
k=1
y
k
n

n
k=1
x
k
y
k
(

n
k=1
x
k
)
2
n

n
k=1
x
2
k
b

n
k=1
x
k

n
k=1
x
k
y
k

n
k=1
x
2
k

n
k=1
y
k
(

n
k=1
x
k
)
2
n

n
k=1
x
2
k
23
- ozone measurements progrmamming
1. problem statement: for a given data (le),
determine a linear model for estimating the ozone mixing ratio at a
specied altitude and
range of altitudes
2. I/O diagram

ranqe ol allludes
lnear model lor ozone mxnq ralo
dala
'zone!.dal'

ranqe ol allludes
lnear model lor ozone mxnq ralo
dala
'zone!.dal'
3. hand example
- data: (ppmv: parts per million value)
altitude (x
k
) ppmv (y
k
)
20 3
24 4
26 5
28 6
4

k=1
x
k
= 98,
4

k=1
y
k
= 18,
4

k=1
x
2
k
= 2436,
4

k=1
x
k
y
k
= 464
denominator = (
4

k=1
x
k
)
2
4
4

k=1
x
2
k
= 140
m

= (
4

k=1
x
k
4

k=1
y
k
4
4

k=1
x
k
y
k
)/denominator = 0.37
b

= (
4

k=1
x
k
4

k=1
x
k
y
k

4

k=1
x
2
k
4

k=1
y
k
)/denominator = 4.6
4. algorithm:
(a) read data le
(b) compute range of altitudes, and parameters(m & b) of linear model
(c) print range of altitudes and linear model
24
- variables
FILE fp
x, y : data file
float sumx, sumy, sumx2, sumxy, denom, m, b
: linear model
first, last : range
integer i, j, k
- C program
#include<stdio.h>
main()
{
int i=0;
float x, y, sumx=0, sumy=0, sumx2=0, sumxy=0,
denom, m, b, first, last;
FILE *fp;
/* read data file */
fp=fopen("zone1.dat","r");
while((fscanf(fp,"%f %f", &x, &y))==2)
{
i++;
if(i==1) first=x;
sumx += x;
sumy += y;
sumx2 += x*x;
sumxy += x*y;
}
last=x;
fclose(fp);
/* compute m & b */
denom=sumx*sumx-i*sumx2;
m=(sumx*sumy-i*sumxy)/denom;
b=(sumx*sumxy-sumx2*sumy)/denom;
/* print range, m & b */
printf("range of altitudes in km: %.2f to %.2f\n", first, last);
printf("linear model: m=%.2f, b=%.2f\n", m, b);
}
25
5 Control Flow Statements
1. selection statements: if, if-else, if-else if, switch
2. repetition statements: while, do-while, for
3. jump statements: goto, continue, break, return
Selection and repetition statements are conditionally controlled.
Jump statements are unconditionally controlled.
5.1 Conditional Expression
evaluated to be true (1) or false (0)
relational operators
> : greater than
>= : greater than or equal to
< : less than
<= : less than or equal to
equality operators
== : equal to
!= : not equal to
logical operators
&& : and
|| : or
! : not
precedence and associativity
arith. op. > rel. op. > eq. op. > logical op.
26
operators associativity
( )
++, --, -, !, (data_type)
*, /, %
+, -
<, <=, >, >=
==, !=
&&
||
=, op=
eg.
x>y+2 : (x>(y+2))
ex!=xye==zee : ((ex!=xye)==zee)
a>c==d!=e : (((a>c)==d)!=e)
5>2&&4>7 : ((5>2)&&(4>7))
!0&&2>3||5>4&&-1 : (((!0)&&(2>3))||((5>4)&&(-1)))
5.2 if Statements
if statements
1) if(expr)
statement; : binary decision.
2) if(expr)
statement;
else
statement; : binary decision
3) if(expr)
statement;
else if(expr)
statement;
...
else
statement; : multiway decision.
eg1.
a=1; sum=0;
if(a<10)
{
++a; : a=2
sum += a; : sum=2
}
27
eg2.
j=k=m=1;
if(x>y)
if(y<z)
k++; x y z j k m
else 1 2 3 2 1 1
m++; 3 1 2 1 2 1
else 3 2 1 1 1 2
j++;
j=k=m=1;
if(x>y)
{
if(y<z)
k++; x y z j k m
} 1 2 3 2 1 1
else 3 1 2 1 2 1
j++; 3 2 1 1 1 1
eg3.
X
X TX
TX

y =
_
_
_
1 if x > 1
x if 1 x 1
1 if x < 1
if(x>1) if(x<-1)
y = 1; y = -1;
else if(x>=-1) else if(x<=1)
y = x; y = x;
else else
y = -1; y = 1;
28
conditional operator (ternary)
E1?E2:E3 operator
if E1
E2;
else
E3;
eg1.
if(y<0) : x = (y<0)?-y:y;
x = -y;
else
x = y;
eg2.
x = max(a,b)
x = (a>b)?a:b;
29
5.3 switch Statements
switch statements : multiway decision.
switch(expr){ : (expr) - integral type (char or int)
case const_expr: statement; - const_expr (char or int)
...
case const_expr: statement;
...
default : statement;
...
}
eg.
int i;
scanf("%d", &i); i
switch(i){ 1 grater than 4
case 2: printf("%d\n", 2); 2 2
break; 3 3
case 3: printf("%d\n", 3); 4 4
break;
case 4: printf("%d\n", 4); 5 greater than 4
break; 6 greater than 4
default: printf("greater than 4\n"); ... ...
}
30
5.4 while Statements
while statements
1) while(expr) 2) do
statement; statement;
while(expr);
eg.
i=1; sum=0; i=1; sum=0;
while(i<=10){ do{
sum += i; sum += i;
i++; i++;
} }while(i<=10);
i 11 i 11
sum
10

i=1
i = 55 sum
10

i=1
i = 55
i = 11; sum = 55; i = 11; sum = 0;
.
.
.
.
.
.
i 11 i 12
sum 0 sum 11
31
5.5 for Statements
for statements
for(expr1; expr2; expr3)
statement;
* expr1 : initialization
* expr2 : condition for repetition
* expr3 : modification to the loop-control values
eg.
i=1; sum=0; : expr1
while(i<=10){ : (i<=10) - expr2
sum+=i;
i++ : (i++) - expr3
}
for(i=1, sum=0; i<=10; i++) : , comma operator
sum+=i; (the lowest priority)
cf.
int i, j, k; distinguish the variable
func(a, b, c); the direction of argument assignments:
eg. func(a, a++, c);
variety of for statements
for(ch=a; ch<=q; ch++) : ch=a, ..., q
for(n=10; n>0; n--) : n=10, 9, ..., 1
for(n=2; n<60; n*=2) : n=2, 4, 8, 16, 32
for(n=1; n*n<100; n+=2) : n=1, 3, 5, 7, 9
nested for loop :
for(i=5; i>=1; i--){ 5 4 3 2 1
for(j=i; j>=1; j--) 4 3 2 1
printf("%d", j); 3 2 1
printf("\n"); 2 1
} 1
32
5.6 Jump Statements
- break statements: an early exit from
"for, while, do_while, switch"
i=1; sum=0; i=1; sum=0;
while(1){ while(i<=10){
if(i>10) --> sum += i;
break; i++;
sum += i; }
i++;
}
- continue statements: next iteration of "for, while, do-while"
for(i=1, sum=0; i<=10; i++){ for(i=1, sum=0; i<=10; i++)
if(i%2==0) continue; --> if(i%2!=0) sum += i;
sum += i;
}
- goto statements: unconditional jump to the specified location
goto label; label can be character or number.
if(x==0) goto overflow
else
z = y/x;
...
overflow : printf("Infinity\n");
33
6 Function statements
function statements
function: generates an output for the given input
procedure: a part of whole program function in C
denition :
return_data_type function_name(arguments)
data_type arguments;
{
statements;
:
return(expr);
}
return_data_type:
- It does not need to be declared when there is no return statement.
- If it is declared, return_data_type should be declared in main().
return statement: one value can be returned.
eg.
return(1);
return(x + y);
return(x, y); not valid
eg.
main() | int sum(a)
{ | int a;
int i, j; | {
int sum( ); | int i=1, j=0;
scanf("%d", &i); | while (i <= a)
j = sum(i); | {
printf("%d\n", j); | j += i;
} | i++;
| }
| return(j);
| }
34
eg.
int j;
main() | void sum(a)
{ | int a;
int i; | {
void sum( ); | int i=1;
scanf("%d", &i); | while (i <= a)
sum(i); | {
printf("%d\n", j); | j += i;
} | i++;
| }
| }
6.1 Local and Global Variables
Variables are eective in the declared block ( ).
If the names of local and global variables are same, the local variables are
eective, that is, global variables are not valid.
eg.
int global = 1;
main( )
{
int local = 2;
...
{
int very_local;
very_local = global +local;
}
...
}
...
eg.
{
int a = 5;
{
int a = 2;
printf("%d\n", a);
}
printf("%d\n", a)
}
35
eg.
{
int a = 5, b = -3, c = -7
{
int b = 8;
float c = 9.9;
a = b;
{
int c;
c = b;
printf("%d %d %d \n", a, b, c); -> 8 8 8
}
printf("%d %d %f \n", a, b, c); -> 8 8 9.9
}
printf("%d %d %d \n", a, b, c); -> 8 -3 -7
}
6.2 Storage Classes
Storage classes are specied in the declaration statement.
denition: storage class data type variables;
auto - dynamic storage in the main memory
register - temporary storage in the CPU
static - static storage in the main memory
extern - referred variables from the outside of function
cpu
reqsler
MEMCPY
proces
s!
process
2
lexl
proqram)
dynomc
sloraqe
slalc
sloraqe
slalc
aulo
auto: no declaration of storage class auto
eg. int i, j; auto int i, j;
register:
- Variables are stored in the register located at CPU
fast access time.
- If not possible, variables are stored in the dynamic storage, that is, same
as auto.
eg. register int i, j;
36
static:
- Variables are maintained until the program terminates.
- Initialization of variables is performed only once.
- If variables are not initialized, variables are set to 0.
eg.
int i;
main()
{
int j;
for(j=1; j<=3; j++) {
int k=0;
static int l = 0;
i++; k++; l++;
printf("%d %d %d %d \n", i, j, k, l);
}
}
extern : declaration to use variables outside of function
eg.
int global;
main()
{
extern int global;
{
extern int global;
...
- extern can be used for variables in dierent le.
- If referred variables are declared in the same le, extern is not necessary.
37
6.3 Recursive Functions
A function that invokes itself.
eg. factorial computation
n! = 1 2 n =

n
i=1
i iterative form
n! = n (n 1)! if n 1 recursive form
1 if n = 0
iterative form : recursive form long :
long factorial(int k) | long factorial(int k)
{ | {
int j; long term; | if (k == 0)
if(k ==0) return(1); | return(1);
else{ | else
term = 1; | return(k*factorial(k-1));
for(j=2; j<= k; j++) | }
term *= j; |
return(term); |
} |
} |
P
U
S
H
l!)
l2)
l4)
l5)
l3)
compulalon usnq slack
l0)!
!l0) !
2l!) 2
3l2) 6
4l3) 24
5l4) !20
P
C
P
STACK
38
eg.
main()
{
void up_down( ); a1
up_down(1); a2
} a3
void up_down(int n) a4
{ b4
printf("a%d\n", n); b3
if (n < 4) up_down(n+1); b2
printf("b%d\n", n); b1
}
eg. Fibonacci sequence
a
n
= a
n1
+a
n2
if n 2
1 if n = 1, n = 0.
lim
n
a
n
a
n1
=

5+1
2
: golden ratio
int fibonacci(int k)
{
if (k == 0 || k ==1) return(1);
else (fibonacci(k-1) + fibonacci(k-2));
}
fibonacci(5) = ?
OXP
O[P
OYP
OZP
OYP OZP
OXP
OWP
O\P
X X
R
Y
Z
\
Z
R
Y
R
R
R
R
X
_
39
6.4 Macros
Macros are specied by a preprocessing directive.
form: #define macro_name(parameters) macro_text
usage:
The macro_name is replaced by the macro_text.
It should be described in one line.
If we need to describe the macro in the next line,

is attached at
the end of line.
; should not be described.
no space between macto_name and parameters
eg.
#define PI 3.141592
#define degree_C(x) ((x)-32)*(5.0/9.0))
#define MAX(x,y) ((x)>y?(x):(y))
characteristics
faster execution than function statement
If the program segment is large size and frequently used, macros are
inecient.
Sometimes, we need to check the meaning of macros.
cc -E test.c
eg.
#define SQ(X) X*X -> ((X)*(X))
#define PR(X) printf("X is %d\n", X)
main()
{
int x=4;
PR(SQ(X+2)); X+2 * X+2 -> 14
PR(100/SQ(2)); 100/2*2 -> 100
PR(SQ(++X)); ++X * ++X -> 36
}
#include<stdio.h>: the le is includes from the specied path.
#include "myfile.h": include le in the current directory.
#include "/user/kil/am320/myfile.h": the path of le is specied.
40
7 Arrays
Declaration
int a[10]; 1 dimensional array, a[0], a[1], ..., a[9]
int b[10][10]; 2 dimensional array
b[0][0], b[0][1], ..., b[0][9]
b[1][0], b[1][1], ..., b[1][9]
.
.
.
b[9][0], b[9][1], ..., b[9][9]
Initialization
Only external array and static array can be initialized.
If arrays are not initialized, the default value of each element is zero.
eg.
static int a[10]={1,2,3,4,5,6,8,9,10};
static int b[2][10]={{1,2,3,4,5,6,7,8,9,10},{2,4,6,8,10}};
7.1 Arrays and Pointers
The use of pointer.
- name of array = memory address of the rst element,
that is, pointer constant.
eg.
int data[10]; data == &data[0]
pointer variable and pointer constant.
eg.
int *pti, dates[10];
dates == &dates[0]
dates+2 == &dates[2]
*(dates+2) == dates[2]
cf. * operator fetches the values of the assigned memory address.
41
pti = dates; memory address of dates is assigned to the pointer pti
pti+2 == &dates[2]
*(pti+2) == dates[2]
*pti+2 == dates[0]+2
R W
X
a
`
cf.
- pti+1: the memory address of the next data to the data
indicated by pti
- dates+1: the memory address of the next data to
dates[0], that is, &dates[1]
- in general,
*(pti+i) == dates[i]
*(dates+i) == dates[i]
passing the array to the argument of function
memory address of an array is passed
eg.
main() {
int ages[50]
convert(ages);
}
convert(a) | convert(a)
int a[]; | int *a;
{ | {
: | :
a[i] | *(a+i) = a[i]
} | }
a[i] == *(a+i) == ages[i]
&a[i] == a+i == &ages[i]
42
memory address is returned
eg.
main()
{
char *gets();
...
}
char *gets();
{
...
return( ); : return the memory address
}
7.2 Dynamic Memory Allocation Function
Memory allocation functions
#include<stdlib.h> is required.
- malloc(size): allocation of storage for an object of size bytes
returns a pointer.
- calloc(n, size): allocation of storage for an array,
n objects of size bytes
returns a pointer.
Storage is initialized to all zeros.
- free(ptr): free preciously allocated storage for an object
indicated by ptr.
eg. programming an array using pointers.
int x[10]; the index value 10 is xed.
int *x;
x = (int *) calloc (10, sizeof(*x));
*(x+i) == x[i]
43
eg.
int n, *x;
printf("number of elements in an array = ");
scanf("%d", &n)
x = (int *)calloc(n, sizeof(*x));
printf("scanning element values\n");
for(i=0; i<n; i++){
printf("element %d = ", i);
scanf("%d", x+i);
}
...
printf("printing element values\n");
for(i=0; i<n; i++)
printf("x[%d]: %d ", i, *(x+i));
...
}
7.3 Multi-Dimensional Arrays and Pointers.
int zippo[4][2]; /* initialization of two-dimensional arrays */
int *pti;
zppo0|0| zppo0|!| zppo!|0| zppo!|!| . . . .
the sequence of memory storage.
memory address of two dimensional array:
zippo == zippo[0] == &zippo[0][0]
zippo+1 == zippo[1] == &zippo[1][0]
zippo+2 == zippo[2] == &zippo[2][0]
44
pointers of pointers :
the indirection operator * is used more than one time.
eg.
main()
{
static short x[2][2]={{1,2},{3,4}};
short *px[2], **pp;
px[0] = x[0];
px[1] = x[1];
pp = px;
}
eg.
**pp == X[0][0] == 1
**(pp+1) == X[1][0] == 3
*(*pp+1) == X[0][1] == 2
*(*(pp+1)+1) == X[1][1] == 4
in general,
*(*(pp+i)+j) == x[i][j]
cf. 3 dimensional array
*(*(*(pp+i)+j)+k) = x[i][j][k]
eg. programming two dimensional array using pointers
int y[10][10];
...
int **y;
y = (int **)calloc(10, sizeof(*y));
for(i=0; i<10; i++)
*(y+i) = (int *)calloc(10, sizeof(**y));
...
*(*(y+i)+j) == y[i][j]
45
8 Strings
denition of strings:
- one dimensional arrays of char
- terminated by a null character \0
- String constants are written between double quotes.
eg. string
- If " is described, \ is used, that is, \"
declaration of strings:
1. array size is specied:
eg.
char name[10]="string"; or
char name[10]={s,t,r,i,n,g,\0 };
2. array size is not specied:
eg.
char name[]="string";
3. using pointers:
eg.
char *p="string"; or
char *p; p="string";
array of strings: multi-dimensional array of char
eg.
static char fruit[3][7]={"Apple","Pear","Orange"};
static char *fruit[3]={"Apple","Pear","Orange"};
8.1 Input and Output functions of Strings
character input and output functions:
putchar(): print one character to the screen
getchar(): read one character from the keyboard
eg.
char ch;
ch = getchar();
putchar(ch);
putchar(a);
putchar(65);
46
string input and output functions:
gets(s) functions:
- read a string into s from stdin.
- characters are placed in s until a new line character is read.
- the value of s is returned.
puts(s) functions:
- the string s is copied to stdout.
- a new line character is placed at the end of string.
- the value of s is returned.
eg.
char name[80];
char *ptr;
ptr = gets(name); Hong Kil Dong
puts(ptr); Hong Kil Dong
puts(ptr+5); Kil Dong
puts(name); Hong Kil Dong
cf.
scanf("%s", name); (!= gets(name);)
printf("%s\n", name); (== puts(name);)
47
8.2 String Functions
string functions:
- arguments of string functions are memory addresses of strings.
- #include<string.h> is required.
strlen(s): a count of the number of characters before \0 is returned.
strcmp(s1,s2): an integer is returned which is less than, equal to, or
greater than zero depending on whether s1 is lexicographically less than,
equal to, or greater then s2.
strcat(s1,s2): a copy of s2, including the null character, is appended
to the end of s1 and the value of s1 is returned.
strcpy(s1,s2): s2 is copied into s1 until \0 is moved; whatever exits in
s1 is overwritten and the value of s1 is returned.
eg.
char *p, *s1, *s2, *s3;
s1="abcde"; s2="fghij"; s3="hij";
printf("%d\n", strlen(s1+2)); 3
printf("%d\n", strcmp(s2+2, s3)); 0
strcat(s2, s3);
puts(s2); fghijhij
strcpy(s1+1, s2+2);
puts(s1); ahijhij
48
eg.
static char *p[2][3]= {"abc","defg", "hi","jklmo","pgstuvw","xyz"};
printf("%c\n", ***p);
***p == p[0][0][0] == a
**p[1] == p[1][0][0] == j
**(p[1]+2) == p[1][2][0] == x
(*(*(p+1)+1))[7] == p[1][1][6] == w
*(p[1][2]+2) == p[1][2][2] == z
**(p[1]+1) == **(*(p+1)+1) == p[1][1][0] == p
*(*(*(p+1)+1)+1) == p[1][1][1] == q
cf.
p == memory address of pointer array
*p == value of p[0] = memory address of the 0th row
**p == value of p[0][0] == memory address of "abc"
***p == a
49
8.3 Command Line Arguments
Two arguments, conventionally called argc and argv, can be used with main()
to communicate with the operating system.
The variable argc provides a count of the number of command line argu-
ments.
The array argv is an array of pointers to char, and can be thought of an
array of strings.
eg.
main(argc, argv)
int argc;
char *argv[];
{
int i;
printf("argc = %d\n", argc);
for (i=0; i<argc; i++)
printf("argv[%d] = %s\n", i, argv[i]);
}
If we compile the program and then give the following command:
a.out in.data out.data
The following is printed on the screen:
argc = 3
argv[0] = a.out
argv[1] = in.data
argv[2] = out.data
The command line arguments are convenient to pass the le name or optional
values as arguments to main().
50
9 Structures
hierarchy of a le:
file
record
field
. . .
. . .
declaration of structures:
struct structure_name{
field_type field_name;
...
} variable_name;
eg.
#define MAX_TITLE 80
#define MAX_AUTHOR 40
struct book{
char title[MAX_TITLE];
char author[MAX_AUTHOR];
float value;
};
struct book libry;
initialization of structure variables:
static or global variables are initialized.
eg.
static struct book libry = {
"The pirate at the damsel",
"Tenes Vivottee",
1.95
};
member operators:
To access the values of members, we use the structure member operator
".".
eg.
libry.title = "KAIST Annual Report";
51
array of structures:
eg.
#define MAX_BOOK 1000
static struct book libry[MAX_BOOK];
libry[0].value = 2.0;
libry[4].title = "LaTex";
libry[2].title[4]
nested structures:
another structure is dened within the structure.
eg.
#define MAX_LENGTH 20
struct names {
char first[MAX_LENGTH];
char last[MAX_LENGTH];
};
struct guy {
struct names handle; /* nested structure */
char favfood[MAX_LENGTH];
char job[MAX_LENGTH];
float income;
}
main()
{
/* initialization of nested structure */
static struct guy fellow = {
{"Franco", "Wathall"},
"eggplant",
"doormat customizer",
15435.00
};
...
}
52
pointers of structures:
pointers of structures are useful in the following points:
1) member of structure can be easily accessed by -> operator,
2) structure variables can be easily passed to arguments of function, and
3) data structures (eg. list, tree, ... ) can be described.
eg.
main()
{
static struct guy fellow[2] = {
{{"Franco","Wathall"},
"eggplant",
"doormat customizer",
15435.00
},
{{"Rodney","Swillbelly"},
"salmon mousse",
"interior decorator",
35000.00
}
};
struct guy *him;
him = fellow;
printf("him->income is $%.2f\n\n", him->income);
printf("(*him).income is $%.2f\n\n", (*him).income);
him++;
printf("him->favfood is %s\n\n", him->favfood);
printf("him->names.last is %s\n\n", him->names.last);
}
him->income is $15435.00
(*him).income is $15435.00
him->favfood is salmon mousse
him->names.last is Swillbelly
cf.
him->income == (*him).income
Here, . operator has higher priority than * operator.
53
unions: unions follow the same syntax as structures but have members
that share storage.
eg.
union int_or_float {
int n;
float x;
}temp;
temp.n = 4444;
printf("%10d %10.4e\n", temp.n, temp.x); 0000004444 4.3387e-29
temp.x = 4444.0;
printf("%10d %10.4e\n", temp.n, temp.x); -0536852854 4.4440e+03
...
the use of typedef: the typedef declaration allows a data type to be
explicitly associated with user dened identier.
eg.
typedef float REAL;
REAL x, y[25], *pr; (== float x, y[25], *pr;)
variations of typedef:
eg.
typedef char *string;
string name, sign; (== char *name, *sign;)
typedef double scalar;
typedef scalar vector[10];
vector a, b, c; (== a[10], b[10], c[10];)
54
10 Linear Linked Lists
a linear linked list consists of a series of items (or nodes), where each node
contains a pointer to the next node.
struct list {
char d;
struct list *next;
};
typedef struct list element;
typedef element *link;
link head; (= struct list *head;)
head = (link) malloc(sizeof(element));
head -> d = a;
head -> next = NULL;

head -> next = (link) malloc(sizeof(element));

55
10.1 basic list operations
creating a list
counting the elements
looking up an element
concatenation two lists
inserting an element
deleting an element
creating a list
cf.
- stack: nodes are accessed in a last in, rst out (LIFO) manner.
- queue: nodes are accessed in a rst in, rst out (FIFO) manner.
- creating a list
a function that will produce a list from a string
(1) iterative function
link s_to_l(s)
char s[ ];
{
link head = NULL, tail;
int i;
if (s[0]! = \0) {
head = (link)malloc(sizeof(element));
head -> d = s[0];
tail = head;
for (i=1; s[i]! = \0; i++){
tail -> next = (link) malloc(sizeof(element));
tail = tail -> next;
tail -> d = s[i];
}
tail -> next = NULL;
}
return(head);
}
56
(2) recursive function
link s_to_l(s)
char s[ ];
{
link head;
if (s[0] == \0) return(NULL);
else {
head = (link) malloc(sizeof(element));
head -> d = s[0];
head -> next = s_to_l(s+1);
return(head);
}
}
- counting the elements
int count(head)
link head;
{
if (head == NULL) return (0);
else return(1+count(head -> next));
}
- printing the elements
void print_list(head)
link head;
{
if (head == NULL) prinf ("NULL");
else{
printf("%c ->", head ->d);
print_list(head -> next);
}
}
- concatenating list a and list b
void cat_list[a, b)
link a, b;
{
if (a -> next == NULL) a -> next = b;
else cat_list (a -> next, b);
}
57
- inserting an element
void insert(p1, p2, q)
link p1, p2, q;
{
p1 -> next = q;
q -> next = p2;
}
- deleting an element
void delete_from_list(head, ch)
link head; char ch;
{
list temp, prev; int found = 0;
temp = head;
while (temp != NULL && !found) {
if (temp -> d == ch) found = 1;
else {
prev = temp;
temp = temp -> next;
}
if found{
if(temp == head) head = temp -> next;
else {prev -> next = temp -> next}; free(temp);
}
}
- recursive deletion of a list
void destroy_list(head)
link head;
{
if (head != NULL){
destroy_list (head -> next);
free(head);
}
}
58
10.2 Stacks
Stacks are data structures to which entries can be added or removed in a
last in rst out (LIFO) manner.

a
a
standard stack operations:
isempty(t): return 1 if stack is empty
otherwise 0
vtop(t): return the value of the top element
pop(t): return the value of the top element and remove it
push(t, d): place the value d onto the top of the stack
declaration of stack data type:
struct stack{
char d;
struct stack *next;
};
typedef struct stack element;
typedef element *top;
59
- int isempty(t)
top t;
{
return(t == NULL);
}
- char vtop(t)
top t;
{
return(t -> d);
}
- char pop(t)
top *t;
{
top t1 = *t; char ch;
if (! isempty(t1)) {
ch = t1 -> d;
*t = t1 -> next;
free(t1);
return(ch);
}
else
return(\0);
}
- void push(t, x)
top *t ; char x;
{
top temp;
temp = (top) malloc (sizeof (element));
temp -> d= x;
temp -> next = *t;
*t = temp;
}
60
eg. reversing a string using a stack
void reverse(s)
char s[ ];
{
int i;
top t = NULL;
for(i=0; s[i]! = \0, i++) push(&t, s[i]);
for(i=0; s[i]! = \0, i++) putchar(pop(&t));
}
61
10.3 Queues
Queues are data structures that can be used to store information that
needs to be accessed in a rst in, rst out (F1, F0) manner.
quene
lronl rear
. . .
add a node
declaration of queue data type:
struct list{
char d;
struct list *next;
};
typedef struct list element;
typedef element *link;
struct queue_struct{
link front;
link rear;
};
typedef queue_struct *queue;
queue operations:
make queue: set up the queue pointers.
add to queue: add a new item to the rear of queue.
remove from queue: remove an item from the front of a queue.
62
queue make_queue( )
{
queue q;
q = (queue)malloc(sizeof(struct queue_struct));
if (q != NULL){
q -> front = NULL;
q -> rear = NULL;
}
return(q);
}
void add_to_quene(q, ch)
queue q; char ch;
{
link node;
node = (link)malloc(sizeof(element));
node -> d = ch;
node -> next = NULL;
if (q -> rear == NULL)
q -> front = q -> rear = node;
else{
q -> rear -> next = node;
q -> rear = node;
}
}
char remove_from_quene(q)
queue q;
{
char ch = \0; link temp;
if (q -> front != NULL) {
ch = q -> front -> d;
temp = q -> front;
q -> front = q -> front -> next;
free(temp);
}
if (q -> front == NULL)
q -> rear = NULL;
return(ch);
}
63
11 Trees
p
r
rool
branch
q o
level 0
level !
level 2
level 3
rrool node
leal nodes
o. p. q chldren ol r sblnqs
nodes lhal share a common parenl)
k-ary trees: every node has at most k children.
balanced tree: the leaf nodes are all either at the lowest
or next-to-lowest levels.
declaration of binary tree:
struct node{
int d
struct node *left;
struct node *rignt;
};
typedef struct node element;
typedef element *btree;
d
lell rqhl
64
create a binary tree from an array
int a[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
btree create_tree( );
create_tree(a, 0, 10)
a|
a2!| a22|
0
! 2
3 4 5 6
7 8 9
btree init_node(d, p1, p2)
int d; btree p1, p2;
{
btree t;
t = (btree)malloc(sizeof(element));
t -> data = d
t -> left = p1;
t -> right = p2;
return(t);
}
btree creat_tree(a, i, size)
int a[ ], i, size;
{
if (i >= size) return(NULL);
else return(int_node(a[i], create_tree(a, 2*i+1, size),
creat_tree(a, 2*i+2, size)));
}
65
create_tree(a, 0, 10);

wX wY
W
X Y
Z [ \ ]
^ _ ` u|ss
u|ss
UGU UGU UGU UGU
UGU UGU UGU UGU UGU
66
binary tree traversal
- inorder : left, root, right
- preorder : root, left, right
- postorder : left, right, root
inorder(root)
btree root;
{
if (root != NULL){
inorder(root -> left);
printf("%d ", root ->d);
inorder(root -> right);
}
}
preorder(root)
btree root;
{
if(root != NULL){
printf("%d ", root -> d);
preorder(root -> left);
preorder(root -> right);
}
}
postorder(root)
btree rooot;
{
if (root != NULL){
postorder(root -> left);
postorder(root -> right);
printf("%d ", root -> d);
}
}
67
eg.
W
X
Y
Z [ \ ]
^ _ `

inorder(root); 7 3 8 1 9 4 0 5 2 6
preorder(root); 0 1 3 7 8 4 9 2 5 6
postorder(root); 7 8 3 9 4 1 5 6 2 0
68
12 Sorting Methods
records: data structures which include several elds of information.
keys: the data attributes which are used for sorting
elementary methods: require O(n
2
) steps on the average to sort n ran-
domly arrayed records
bubble sort
insertion sort
selection sort
shell sort
advanced methods: require O(nlog n) steps on the average
merge sort
quick sort
heap sort
12.1 Insertion Sort
insertion sort(A) algorithm: an array A[1, . . ., n] is sorted in ascending
order.
1. for j 2 to length[A]
2. do key A[ j ]
3. i j-1
4. while i>0 and A[ i ] > key
5. do A[i+1] A[ i ]
6. j i-1
7. A[i+1] key
Steps 3 through 7 indicate that insert A[ j ] into
the sorted sequence A[1, . . ., j-1]
69
eg. A = 5, 2, 4, 6, 1, 3 )
Y
[
]
X
Z
running time analysis cost time
1. for j 2 to length[A] c
1
n
2. do key A[ j ] c
2
n 1
3. i j 1 c
3
n 1
4. while i > 0 and A[ i ] > key c
4

n
j=2
t
j
5. do A[i + 1] A[ i ] c
5

n
j=2
(t
j
1)
6. i i 1 c
6

n
j=2
(t
j
1)
7. A[i + 1] key c
7
n 1
T(n) = c
1
n +c
2
(n 1) +c
3
(n 1) +c
4

n
j=2
t
j
+c
5

n
j=2
(t
j
1)
+c
6

n
j=2
(t
j
1) +c
7
(n 1)
t
j
= j for j = 1, . . . , n 1

n
j=2
t
j
=

n
j=2
j =
n(n+1)
2
1

n
j=2
(t
j
1 =
n(n1)
2
T(n) = (
c
4
2
+
c
5
2
+
c
6
2
)n
2
+(c
1
+c
2
+c
3
+
c
4
2

c
5
2

c
6
2
+c
7
)n(c
2
+c
3
+c
4
+c
7
)
(worst-case) running tim of o(n
2
)
average-case running time of o(n
2
)
70
C program of insertion sort:
void insertion_sort(a, n)
double a[ ]; int n;
{
int i, j;
double key;
for(j=1; j<n; j++){
key = a[j];
i = j-i;
while(i >= 0 && a[i] > key){
a[i+1] = a[i];
i = i-1;
}
a[i+1] = key;
}
}
71
12.2 Quicksort
Sorting Performance
worst-case running time: O(n
2
)
average running time: O(nlog n)
Algorithm:
1. divide procedure:
the array A[p,r] is partitioned into two nonempty subarrays A[p,q]
and A[q+1,r] such that each element of A[p,q] is less than or equal
to each element of A[q+1,r].
2. conquer procedure:
A[p,q] and A[q+1,r] are sorted by recursive calls to quicksort.
void quicksort(a, n)
int a[ ], n;
{
int k, pivot;
if(findpivot(a, n, &pivot)!= 0){
k = partition(a, n, pivot);
quicksort(a, k);
quicksort(a+k, n-k);
}
}
int findpivot(a, n, pivot_ptr)
int a[ ], n, *pivot_ptr;
{
int i;
for(i=1; i<n; i++)
if(a[0] != a[i]){
*pivot_ptr = (a[0]>a[i])?a[0]:a[i];
return(1);
}
return(0);
}
72
int partition(a, n, pivot)
int a[ ], n, pivot;
{
int i=0, j=n-1, temp;
while(i<=j){
while(a[i] < pivot) i++;
while(a[j] >=pivot) j--;
if(i<j){
temp = a[j];
a[i] = a[j];
a[j] = temp;
i++; j--;
}
}
return(i);
}
eg.
int a[12]={4,7,5,2,5,9,2,1,9,-5,-3};
i: 0 1 2 3 4 5 6 7 8 9 10 11
a[i]: 4 7 9 5 2 5 9 2 1 9 -5 -3
4 -3 -5 5 2 5 1 2 9 9 9 7
2 -3 -5 1 2 5 5 4 7 9 9 9
1 -3 -5 2 2 4 5 5
-5 -3 1
-5 -3
Figure 1: Execution of quicksort
73
Performance of quicksort
1. Worst-case partitioning:
T(n) = T(n 1) +O(n)
where O(n) represents partitioning cost.
Since T(1) = O(1),
T(2) = O(1) +O(2), ,
T(n) =

O(k) = O(

k) = O(n
2
).
In general,

O(f(k)) = O(

f(k)).
) (
2
n O
n
n -1
n -2
n -3
1
1
1
1
1 1
2
2
n
n
n -1
no. of operations
Figure 2: Worst-Case Partitioning
74
2. Best-case partitioning:
T(n) = 2T(n/2) +O(n) = O(nlog
2
n)
In general, T(n) = f(n) +aT(n/b)
T(n) = O(n
log
b
a
log
b
n).
(proof)
T(n) = f(n) +aT(n/b)
= f(n) +af(n/b) +a
2
T(n/b
2
)
= f(n)+af(n/b)+a
2
f(n/b
2
)+...+a
log
b
n1
f(n/b
log
b
n1
)+a
log
b
n
T(1).
Since a
log
b
n
T(1) = O(n
log
b
a
),
T(n) = O(n
log
b
a
) +

a
j
f(n/b
j
). ... (1)
f(n) = O(n
log
b
a
).

a
j
(n/b
j
)
log
b
a
= n
log
b
a

(a/b
log
b
a)
j
= n
log
b
a

1 = n
log
b
a
log
b
n.
... (2)
From (1) and (2),
T(n) = O(n
log
b
a
log
b
n).
If a = b = 2, T(n) = O(nlog
2
n).
4
n
4
n
4
n
4
n n
2
log
) log (
2
n n O
2
n
2
n
n n
n
n
n
1 1 1 1
Figure 3: Best-case partitioning
75
3. Average-case partitioning:
The average-case running time is much closer to the best case.
eg.
T(n) = T(9n/10) +T(n/10) +O(n)
log
10
n = O(log
2
n) (log
10
n < log
2
n)
log
10/9
n = O(log
2
n) (log
10/9
n > log
2
n)
the total cost = O(nlog
2
n)
n
10
1
n
10
9
100
1
100
9
100
9
n n n n
100
81
1
1
n n
n
n
n
n
10
log
n
9
10
log
<= n
<= n
Figure 4: Average-case partitioning
76
12.3 Heapsort
heap:
1. ascending heap
a balanced binary tree in which the key value for every node is greater
than or equal to the key value for its parent node.
2. descending heap
heapsort algorithm: making a heap from an array a[1..n]
heapify(a, n)
{
for i <- int(n/2) downto 1
bubbledown(a, i, n)
}
bubbledown(a, i, n)
{
child <- 2i
if (child < n) and (a[child+1] > [child])
child <- child+1
if (child <= n) and (a[i] < a[child]){
swap(a[i], a[child])
bubbledown(a, child, n)
}
5
3 7
9 1 7 9
18 16 18
18
18 9
16 5 7 3
9 7 1
heapify(a,n)
1 2 3 4 5 6 7 8 9 10
5 7 3 9 1 7 9 18 16 18
18 18 9 16 5 7 3 9 7 1
eg.
77
heapsort(a, n)
{
heapify(a, n)
i <- n
while (i > 1) {
swap(a[1], a[i])
i <- i-1
bubbledown(a, 1, i)
}
}
eg.
1 2 3 4 5 6 7 8 9 10
18 18 9 16 5 7 3 9 7 1
18 16 9 9 5 7 3 1 7 18
16 9 9 5 5 7 3 1 18 18
18
18 9
16 5 7 3
9 7 1
18
16 9
9 5 7 3
1 7 18
16
9 9
7 5 7 3
1 18 18
9
7 9
1 5 7 3
16 18 18
1
3 5
7 7 9 9
16 18 18
1 3 5 7 7 9 9 16 18 18
78
Performance of heapsort
running time of heapsort :
no. of calls to bubbledown
. . .
<= 3 <= 2 0
8
n
4
n
2
n
n 1
2 3 1 1 0
<= 1
1 2 3 4 5 6 7 8
2
1
8
7
6
5 4
3
eg. n =8
79
heapsort(a, n)
{
heapify(a, n) -> O(n)
i <- n
while (i > 1) {
swap(a[1], a[i])
i <- i-1
bubbledown(a, 1, i) -> O(log n)
}
}
total running time of heapsort: O(nlog
2
n)
80
13 Searching Methods
13.1 linear and binary search methods
linear search: the sequential process of scanning through the records,
starting at the rst record
binary search: the procedure for locating an element in an already sorted
array
eg.
int bin_search(a, n, target)
int a[ ], n, target;
{
int low, high, middle;
low = 0; high = n-1;
while(low <= high){
middle = (low+high)/2;
if (target < a[middle])
high = middle-1;
else if (target > a[middle])
low = middle+1;
else
return(middle);
}
return(not_found);
}
it requires O(log
2
n) iterations in the worst case.
eg.
int a[8]={1,4,8,10,11,12,16,17};
bin_search(a, 8, 8);
low=0, high=7
----------------
middle=3, high=2
----------------
middle=1, low=2
----------------
middle=2, found
----------------
81
13.2 binary search tree
every node in the left subtree has a key value no larger than the root node.
every node in the right subtree has a key value no smaller than the root
node.
this property applies recursively to each subtree of the complete tree.
eg.
5
4 9
1 8 10
10
9
8
5
4
1
tree traversal: a process of visiting all nodes in a tree.
pre-order : root,left,right
in-order : left,root,right
post-order : left,right,root
searching in a binary search tree
searching procedure :
1. check to see if the root node is empty
2. compare the search key with the root key
3. if it matches, end the search
4. if it is smaller, search the left subtree
5. otherwise, search the right subtree
82
struct node{
int data;
struct node *left;
struct node *right;
};
typedef struct node *bin_tree;
bin_tree find_node(t, value)
bin_tree t; int value;
{
if(empty_tree(t)) return(t);
else{
if(t->data == value) return(t);
else if(t->data > value)
return(find_node(t->left, value));
else return(find_node(t->right, value));
}
}
int empty_tree(t)
bin_tree t;
{
if(t == NULL) return(1);
else return(0);
}
83
Adding a node to a binary search tree
Only one location at which the new node can be added.
Finding the location and adding the new node is required.
void add_node(t, new_node)
bin_tree t, new_node;
{
bin_tree prev;
while(t != NULL){
if(new_node->data < t->data){
prev = t;
t = t->left;
}
else{
prev = t;
t = t->right;
}
}
if(new_node->data < prev->data)
prev->left = new_node;
else
prev->right = new_node;
}
eg.
7
add
5
4 9
1 8 10
7
84
Construction of a binary search tree
add_node()
the structure of the resulting tree will depend heavily on the sequence
in which the nodes are added.
make a binary search tree from a sorted array
bin_tree creat_bst(a, n)
bin_tree a; int n;
{
int mid;
if(n == 0) return(NULL);
mid = n/2;
(a+mid)->left = creat_bst(a, mid);
(a+mid)->right = creat_bst(a+mid+1, n-mid-1);
return(a+mid);
}
0 1 11
0 1 2 3 4 5 6 7 8 9 10 11
a
data
0 1
left right
11
l m r
l m m r r l
l l l l m m m m
6
3 9
1 5 8 11
0 2 4 7 10
(balanced) banary search tree
. . .
85
Node deletion
1. locate the node to be deleted
2. node deletion
(a) a leaf node: just remove the node.
(b) a node with only one child:
the node is replaced with its child node.
(c) a node with two children:
- two options:
i. replace it with the largest-key-valued node from
the left subtree.
ii. replace it with the smallest-key-valued node from
the right subtree.
55
91 77 70
86 74 37
80 51 40
63 48
(a)
(b)
(c)
functions to be implemented:
delete_node( ) includes the following functions:
find_parent( )
delete_leaf( )
delete_one_child( )
delete_two_children( )
86
#define true 1
#define false 0
bin_tree delete_node(t, value, flag)
bin_tree t; int value; int *flag;
{
bin_tree parent, node;
node = find_parent(t, &parent, value);
if(node == NULL) {
*flag = false;
return(t);
}
if((node->left == NULL) && (node->right == NULL))
t = delete_leaf(t, parent, node);
else if((node->left == NULL) || (node->right == NULL))
t = delete_one_child(t, parent, node);
else delete_two_children(t, node);
*flag = true;
return(t);
}
bin_tree find_parent(t, p, value)
bin_tree t; bin_tree *p; int value;
{
*p = NULL;
while(t != NULL && value != t->data){
*p = t;
if(value < t->data)
t = t->left;
else
t = t->right;
}
return(t);
}
87
bin_tree delete_leaf(t, p, node)
bin_tree t, p, node;
{
if(p == NULL) t = NULL;
else if(p->data < node->data)
p->right = NULL;
else
p->left = NULL;
free(node);
return(t);
}
88
bin_tree delete_one_child(t, p, node)
bin_tree t, p, node;
{
if(p == NULL){
if(node->left == NULL) node = t->right;
else node = t->left;
free(t);
return(node);
}
if(node->left == NULL)
if(p->data < node->data)
p->right = node->right;
else
p->left = node->right;
else
if(p->data < node->data)
p->right = node->left;
else
p->left = node->left;
free(node);
return(t);
}
89
void delete_two_children(t, node)
bin_tree t, node;
{
int temp;
bin_tree small, small_parent;
/* find the node which has the smallest
key value in the right subtree */
small = node->right;
small_parent = node;
while(small->left != NULL){
small_parent = small;
small = small->left;
}
temp = small->data;
if(small->left == NULL && small->right == NULL)
delete_leaf(t, small_parent, small);
else
delete_one_child(t, small_parent, small);
node->data = temp;
}
90
14 Hashing
Hashing involves some transformation of the search key which is used to
access an element of the array
eg. telephone directories, on-line dictionaries
Hash Table : an array of active records
Hash Function : mapping the record key into the appropriate array index
spreading the records fairly evenly over the hash table, and minimizing
the number of collisions (duplicate hashes)
collision resolution process
linear probing
double hashing
chaining
student ID number % 8
960324
970344
980285
(in am320)
eg. h(key_value)=key_value % M
where M is the size of the hash table.
h
0
1
7
hash table
** hash entry can be a pointer, or a
data structure containing all info.
91
14.1 Linear Probing
1. set up a hash table
key_value
h 0
M-1
if null,
set that entry in the hash table
else
scan the hash table entries
to find the first null entry.
(no null entry -> hash table is full)
2. search the hash table
key_value
h 0
M-1
if the hash table entry points
to the item, search is successful.
else if null, search is unsuccessful.
else scan subsequent entries
until a null entry is found or searching
item is found.
init_hash(): allocates an initial, empty hash table
make_hash(): 1(hash table -> linear linked list)
search_hash(): 2
#define table_full -1
struct node{
int data;
struct node *next;
};
typedef struct node element;
typedef element *list
92
list *init_hash(n); int n;
{
list *hash_table;
int i;
hash_table = calloc(n, sizeof(list));
if(hash_table == NULL) exit(-1);
for(i=0; i<n; i++)
*(hash_table+i) = NULL;
return(hash_table);
}
int search_hash(table, n, key)
list table[ ]; int n, key;
{
int entry = key%n;
int count = 0;
while((count < n) && (table[entry] != NULL) &&
(table[entry]->data != key)){
entry = (++entry)%n;
count++;
}
if(count == n) return(table_full);
return(entry);
}
cf. exit(status)
this function terminates a program.
all buffers are finished and files are closed
the function returns the value of status to the calling
process.
status: 0 the program ran properly
(integer) non-zero otherwise
93
list *make_hash(l, n)
list l; int n;
{
list *hash_table;
int j;
hash_table = init_hash(n);
while(l != NULL){
j = search_hash(hash_table, n, l->data);
if(j == table_full) exit(2);
if(*(hash_table+j) == NULL) *(hash_table+j) = l;
else printf("Record duplicate exits\n");
l=l->next;
}
return(hash_table);
}
h
l
ht
ht+j
94
14.2 Double Hashing
Repeating the hash process itself with the second hash function.
eg. settling up a hash table
entry number = h(k)
if null, set that entry in the hash table
else apply the second hash function
h(k)=M-2-k%(M-2) [1,...,M-2]
entry number = entry no.+h(k)
repeat the above procedure until a null entry is found
or table_full is found.
int search_hash(table, n, key)
list table[ ]; int n, key;
{
int entry, count;
entry = key%n;
for((count = 0; (count < n) && (table[entry] != NULL)&&
(table[entry]->data != key); count++)
entry = (entry+(n-2)-key%(n-2))%n;
if(count == n) return(table_full);
return(entry);
}
95
14.3 Chaining
Set up a linked list whenever a duplicate hash is encountered.
If there is a lot of uncertainty about the problem size, chaining may be
the preferred approach.
key_value
h 0
n-1
1
int search_hash_chain(table, n, key, pre1)
list table[ ]; int n, key; list *pre1;
{
int entry;
list temp;
entry = key%n;
temp = table[entry];
while((temp != NULL) && (temp->data != key)){
*pre1 = temp;
temp = temp->next;
}
return(entry);
}
96
list *make_hash_chain(l, n)
list l; int n;
{
list *hash_table; int i, j;
list temp, pre2;
hash_table = init_hash(n);
while(l != NULL){
temp = l->next;
j = search_hash_chain(hash_table, n, l->data, &pre2);
if(*(hash_table+j) == NULL) ==>(i)
{*(hash_table+j) = l; l->next = NULL;}
else if(pre2->next == NULL) ==>(ii)
{pre2->next = l; l->next = NULL;}
else
{printf("Record duplicate exits\n"); free(l);}
l = temp;
}
return(hash_table);
}
temp l
l
a b c
temp l
pre2
pre1
pre2
pre1
a
0
(i)
(ii)
b c
n-1
97
Deleting a record
linear probing & double hashing: a dummy node can be allocated to
the delete node to permit searches for duplicated hash record lying
beyond the removed record.
chaining: remove a node from a linked list
98
15 Root Finding Methods
Finding the real root of univariate function
assumptions:
1. known interval [x
1
, x
2
]
initial starting value, x
0
[x
1
, x
2
]
2. at most one root
3. approximate solution
- oating point arithmetic
- iterative methods (asymptotic convergence)
basic idea: nd the interval [x
1
, x
2
] where f(x
1
) f(x
2
) 0
the interval [x
1
, x
2
] must contain at least one root.
(odd number of roots)
15.1 Root Finding Based on Search Methods
1. linear search:
find_interval( ) returns
1 if a root is found
0 otherwise
f: target function; xstart, xend: given interval;
xleft, xright: interval containing a root;
int find_interval(f, xstart, xend, stepsize, xleft, xright)
double (*f)(), xstart, xend, stepsize, *xleft, *xright;
{
int found_root = TRUE;
*xleft = xstart;
*xright = *xleft + stepsize;
while(*xleft < xend && f(*xleft)*f(*xright) > 0.0){
*xleft = *xright;
*xright += stepsize;
if(*xright > xend) *xright = xend;
}
if(f(*xleft)*f(*xright) > 0.0) found_root = FALSE;
return(found_root);
}
99
nds only the rst interval containing the root.
fails to nd an interval if an ever number of roots are suciently close
to one another
2. binary search (bisection method):
divide and conquer method -
the initial interval [x
1
, x
2
] is successively halved
the size of interval at the nth iteration = (x
2
x
1
)/2
n
f: target function; x1, x2: given interval;
epsilon: tolerance (accuracy) of solution
double bisect(f, x1, x2, epsilon)
double (*f)(), x1, x2, epsilon;
{
double y;
for(y=(x1+x2)/2.0; fabs(x1-y)>epsilon; y=(x1+x2)/2.0)
if(f(x1)*f(y) <= 0.0) x2 = y;
else x1 = y;
return(y);
}
100
15.2 Root Finding Based on Target Functions
1. secant method
algorithm
(a) tting a straight line
between the points, (x1, f(x1)) and (x2, f(x2)).
(b) nd the root, x

of the straight line.

is close to x2, stop.

, and go to step1.
pros and cons:
root can be the value outside the initial interval.
there is no way to detect that it will fail to nd a root
(x
2
, y
2
)
y = f(x)
x
2
x
3
x
1
(x
1
, y
1
)
x
4
101
C program
f: target function; x1, x2: given interval;
epsilon: tolerance; max_tries: maximum number of iterations;
found_flag: flag of finding a solution
double secant(f, x1, x2, epsilon, max_tries, found_flag)
double (*f)(), x1, x2, epsilon;
int max_tries, *found_flag
{
int count = 0;
double root = x1;
*found_flag = 1;
while(fabs(x2-x1) > epsilon && count < max_tries){
root = x1 - f(x1)*(x2-x1)/(f(x2)-f(x1));
x1 = x2;
x2 = root;
count++;
}
if(fabs(x2-x1) > epsilon) *found_flag = 0;
return(root);
}
102
2. Newtons method
algorithm:
(a) initial guess x
n
(n = 0).
(b) nd the root of the straight line that ts a tangent to f(x):
x
n+1
= x
n

f(x
n
)
f

(x
n
)
_
f

(x
n
) =
f(x
n
) f(x
n+1
)
x
n
x
n+1
and f(x
n+1
) 0 as n
_
(c) if [x
n+1
x
n
[ < , stop
otherwise, go to 2.
cf.
f(x
n
)
f

(x
n
)
could result in a oating-point overow error message
dierent initial guess, start again
x
0
x
1
x
2
103
C program
f: target function; fprime: f(x_n);
dmin: minimum value of f(x_n); x0: starting point;
epsilon: tolerance; error: flag of finding a solution
double newton(f, fprime, dmin, x0, epsilon, error)
double (*f)(), (*fprime)(), dmin, x0, epsilon; int *error;
{
double deltax;
deltax = 2.0*epsilon;
*error = 0;
while(!(*error) && fabs(deltax) > epsilon)
if(fabs(fprime(x0) > dmin){
deltax = f(x0)/fprime(x0);
x0 = x0 - deltax;
}
else
*error = 1;
return(x0);
}
104
problems in applying root nding methods
(a) false convergence in secant and Newtons methods -
very small change in x
n
and x
n+1
but far from a root
check whether the returned value is a root
x
i1
x
i
x
i+1
(b) cycling (in Newtons method)
f(x) =
1
1+e
x

1
2
selecting a dierent starting value.
xn+1
xn
x
y
105
16 Numerical Integration
Integration of f(x) in the range of x [a, b] is approximately computed by
the summation of weighted function values.
16.1 Trapezoidal Rule
Rectangular Rule: Integration of f(x) in the range of x [a, b] is approxi-
mated by Riemannian integration formula. rarely used in practice.
y = f(x)
x
y
y = f(x)
y
x
106
Trapezoidal Rule: Integration of f(x) in the range of x [a, b] is approx-
imated by the summation of trapezoidal segments more accurate ap-
proximation than the rectangular rule.
y = f(x)
y
x
double trapezoidal(f, a, b, n)
double (*f)(), a, b; int n;
{
double answer, h;
int i;
answer = f(a)/2;
h = (b-a)/n;
for(i=1; i<=n; i++)
answer = answer + f(a+i*h);
answer = answer - f(b)/2;
return(h*answer);
}
x
f(x)
f(x + h)
sum=
h
2
(f(x) +f(x + h))
x + h
107
The error bound of trapezoidal rule is given by O(h
2
).
Here, the error is given by
E =
_
b
a
f(x)dx
b a
2
(f(a) +f(b))
.
The Taylor series expansion of f(x) around x is
f(x) = f( x) +zf

( x) +
z
2
2
f

( x) +
where z = x x.
That is,
_
b
a
f(x)dx =
_
h/2
h/2
(f( x) +zf

( x) +
z
2
2
f

( x) + )dz
= hf( x) +
1
24
h
3
f

( x) + (1)
On the other hand,
b a
2
(f(a) +f(b)) =
h
2
_
f( x)
h
2
f

( x) +
1
2
_
h
2
_
2
f

( x) +f( x) +
h
2
f

( x) +
_
= hf( x) +
1
8
h
3
f

( x) + (2)
From (1) and (2),
E =
_
b
a
f(x)dx
b a
2
(f(a) f(b))

=
1
12
h
3
f

( x)
Therefore, total error is given by
E
T

=
1
12
(B A)
3
N
3
N

i=1
f

( x)
where h =
BA
N
and x
i
=
1
2
(x
i
+x
i+1
).
Let

f

=
1
N

N
i=1
f

( x). Then,
E
T

=
1
12
(B A)h
2

f

O(h
2
).
108
16.2 Simpsons Rule
Simpsons Rule: approximating the function within each of n intervals by
some polynomials
the second-order polynomials (Simpsons rule):
Sum

=
h
6
[f(x) + 4f(x +
h
2
) +f(x +h)]
error bound: O(h
4
)
f(x)
f(x + h)
x + h x
f(x +
h
2
)
the third-order polynomials (Simpsons 3/8 rule):
Sum

=
h
8
[f(x) + 3f(x +
h
3
) + 3f(x +
2
3
h) +f(x +h)]
error bound : O(h
5
)
f(x)
f(x + h)
x + h x
f(x +
2h
3
)
f(x +
h
3
)
109
double simpsion(f, a, b, n)
double (*f)(double), a, b; int n;
{
double answer, h, x;
int i;
answer = f(a);
h = (b-a)/n;
for(i=1; i<=n; i++){
x = a+i*h;
answer = answer + 4*f(x-h/2) + 2*f(x);
}
answer = answer - f(b);
return(h*answer/6);
}
How to decide the number of intervals n?
large n more accurate integration, but requires more computation time.
1. increase n until the change in approximated integral is small
2. wide intervals on portions of function that are nearly linear,
smaller intervals on portions of function that are more complicated.
110
- The First Method:
double new_simpson(f, a, b, no, tolerance)
double (*f)(double), a, b; int no; double tolerance;
{
double check = tolerance + 1.0;
double lowval, val;
lowval = simpson(f, a, b, no);
while(check > tolerance){
no = 2*no;
val = simpson(f, a, b, no);
check = fabs((val-lowval)/val);
lowval = val;
}
return(val);
}
- The Second Method
double adapt(f, a, b, n, tolerance)
double(*f)(double), a, b; int n; double tolerance;
{
double val, check;
val = simpson(f, a, b, 2*n);
check = fabs((simpson(f,a,b,n)-val)/val);
if(check > tolerance)
val = adapt(f, a, (a+b)/2.0, n, tolerance)+
adapt(f, (a+b)/2.0, b, n, tolerance);
return(val);
}
Other Methods
Quadrature method (Refer Numerical Recipes in C)
Monte Carlo integration
111
17 Numerical Solutions of
Ordinary Dierential Equations
the solution of the rst-order ordinary dierential equations (O.D.E.s)
dy
dx
= g(x, y)
conversion of an math-order O.D.E. to the rst-order O.D.E.
eg. y

= g(x, y, y

, y

)
let y
0
= y, y
1
= y

, y
2
= y

, then
y

2
= g(x, y
0
, y
1
, y
2
)
where y

1
= y
2
, y

0
= y
1
y

=
_
_
y

2
y

1
y

0
_
_
=
_
_
g(x, y
0
, y
1
, y
2
)
y
2
y
1
_
_
17.1 Eulers Method
Eulers method is the most straightforward method, but not accurate.
Problem: y

= g(x, y)
algorithm:
Step 1. set up the variables :
x = x
0
, y = y
0
(initial value of x and y)
xlast (the largest value of x)
h (an increment of x)
Step 2. output x and y
Step 3. set x = x +h
Step 4. compute y = y +h g(x, y)
Step 5. if x > last, stop
otherwise, go to Step 2.
112
C program:
void euler(x0, y0, g, h, xlast)
double x0, y0, (*g)(), h, xlast;
{
for(; x0<=xlast; x0=x0+h){
printf("%e %e\n", x0, y0);
y0 = y0 + h*g(x, y);
}
}
sources of error in Eulers method: (x
0
, y
0
)
(x
0
+h, y
0
+hg(x
0
, y
0
) y
1
= [y(x
0
+h) y
1
[ = O(h
2
)
rst-order term in Taylor expansion is considered.
the local error per step: O(h
2
)
Here, errors are accumulated.
17.2 Runge-Kutta method
better estimation of the derivative by evaluating the derivative at more
than one point for each step and combining the various computed deriva-
tives to obtain a better approximation of the solution
the derivative is computed at the point x +
h
2
(the 2nd-order Runge-Kutta method)
algorithm:
Step 1. set up the variables:
x = x
0
, y = y
0
, xlast, h
temp (a temporary variable)
Step 2. output x and y
Step 3. compute temp = h g(x, y)
Step 4. compute y = y +h g(x +
h
2
, y +temp/2)
Step 5. set x = x +h
Step 6. if x > xlast, stop
otherwise, go to step 2
the local error per step: O(h
3
).
113
Second-order Runge-Kutta method
dy
dx
= g(x, y), y(0) = y
0
y
n+1
= y
n
+
_
x
n+1
x
n
g(x, y)dx
x
n+1
= x
n
+h
applying the trapezoidal rule to the integration term:
_
x
n+1
x
n
g(x, y)dx

=
1
2
h[g(x
n
, y
n
) +g(x
n+1
, y
n+1
)]
_
y
n+1
= y
n
+hg(x
n
, y
n
), (estimate of y
n+1
using Eulers method)
y
n+1
= y
n
+
h
2
[g(x
n
, y
n
) +g(x
n+1
, y
n+1
)] (1)
the local error par step: O(h
3
)
Error Analysis:
dy
dx
= g(x, y)
Taylor series expansion around x
n
and y
n
:
y
n+1
= y
n
+ hg +
h
2
2
[g
x
+ g
y
g] +
h
3
6
[g
xx
+ 2g
xy
g + g
yy
g
2
+ g
x
g
y
+ g
2
y
g] +
o(h
4
) (2)
where g
x
=
g
x
[
x=x
n
,y=y
n
, g
y
=
g
y
[
x=x
n
,y=y
n
applying Taylor series expansion of g(x
n+1
, y
n+1
) around (x
n
, y
n
) in (1)
y
n+1
= y
n
+hg +
h
2
2
[g
x
+g
y
g] +
h
3
4
[g
xx
+ 2g
xy
g +g
yy
g
2
] +o(h
4
) (3)
by comparing (2) and (3)
error in the order of h
3
O(h
3
)
114
Fourth-order Runge-Kutta method (most widely used method)
dy
dx
= g(x, y), y(0) = y
0
y
n+1
= y
n
+
_
x
n+1
x
n
g(x, y)dx
x
n+1
= x
n
+h
1. applying the simpsons 1/3 rule in the integral term:
_
x
n+1
x
n
g(x, y)dx

=
h
6
[g(x
n
, y
n
) + 4g(x
n
+
h
2
, y
n
+ ) +g(x
n+1
, y
n+1
)]
k
1
= hg(x
n
, y
n
)
k
2
= hg(x
n
+
h
2
, y
n
+
k
1
2
)
k
3
= hg(x
n
+
h
2
, y
n
+
k
2
2
)
k
4
= hg(x
n
+h, y
n
+k
3
)
y
n+1
= y
n
+
1
6
[k
1
+ 2k
2
+ 2k
3
+k
4
]
the local error per step: O(h
5
)
2. applying the Simpsons 3/8 rule in the integral term:
_
x
n+1
x
n
g(x, y)dx

=
h
8
[g(x
n
, y
n
)+3g(x
n
+
h
3
, y
n
+)+3g(x
n
+
2h
3
, y
n
)+
g(x
n
, y
n
)]
k
1
= hg(x
n
, y
n
)
k
2
= hg(x
n
+
h
3
, y
n
+
k
1
3
)
k
3
= hg(x
n
+
2h
3
, y
n
+
k
1
3
+
k
2
3
)
k
4
= hg(x
n
+h, y
n
+k
1
k
2
+k
3
)
y
n+1
= y
n
+
1
8
[k
1
+ 3k
2
+ 3k
3
+k
4
]
the local error per step: O(h
5
)
115
algorithm:
Step 1. set up the variables
x = x0, y = y0, xlast, h, ta, tb, tc, td
Step 2. output x and y
Step 3. compute
ta = h g(x, y) estimation of y evaluated at x
tb = h g(x +h/2, y +ta/2),
tc = h g (x + h/2, y + tb/2) estimation of y evaluated at
x +h/2
td = h g(x +h, y +tc) estimation of y evaluated at x +h
Step 4. compute y = y +
1
6
(ta + 2tb + 2tc +td)
Step 5. set x = x +h
Step 6. if x > xlast, terminate
otherwise, go to step 2
C program:
void runge_kutta(x0, y0, g, h, xlast)
double x0, y0, (*g)(), h, xlast;
{
double ta, tb, tc, td;
for(; x0<=xlast; x0=x0+h){
printf("%e %e\n", x0, y0);
ta = h*g(x0, y0);
tb = h*g(x0+h/2.0, y0+ta/2.0);
tc = h*g(x0+h/2.0, y0+tb/2.0);
td = h*g(x0+h, y0+tc);
y0 = y0+(ta+2.0*tb+2.0*tc+td)/6.0;
}
}
cf. error bound: O(h
5
)
cf. 2nd version based on the Simpsons 3/8 rule:
ta = h g(x, y)
tb = h g(x +
h
3
, y +
ta
3
)
tc = h g(x +
2h
3
, y +
ta
3
+
tb
3
)
td = h g(x +h, y +ta tb +td)
y = y +
1
8
(ta + 3tb + 3tc +td)
116
18 Numerical Solution of
Simultaneous Equations
Consider only the case where the number of linear equations is same as the
number of unknowns: n equations and n unknowns
_
_
_
_
_
a
00
a
01
a
0,n1
a
10
a
11
a
1,n1
.
.
.
.
.
.
.
.
.
.
.
.
a
n1,0
a
n1,1
a
n1,n1
_
_
_
_
_
_
_
_
_
_
x
0
x
1
.
.
.
x
n1
_
_
_
_
_
=
_
_
_
_
_
b
0
b
1
.
.
.
b
n1
_
_
_
_
_
Ax = b
18.1 simple Gauss-Jordan Elimination
Let us consider the following simultaneous equations:
_
_
_
2x
0
+ 3x
1
+ 5x
2
= 10
x
0
+ 2x
1
+ 4x
2
= 5
4x
0
+ 7x
1
+ 3x
2
= 15

_
_
_
2x
0
+3x
1
+ 5x
2
= 10
0.5x
1
+ 1.5x
2
= 0
1x
1
7x
2
= 5

_
_
_
2x
0
4x
2
= 10
0.5x
1
+1.5x
2
= 0
10x
2
= 5

_
_
_
2x
0
= 12
0.5x
1
= 0.75
10x
2
= 5
x
0
= 6, x
1
= 1.5, and x
2
= 0.5
G-J elimination by the ith diagonal element (a
ii
):
row
j
= row
j

a
ji
a
ii
row
i
(j ,= i)
117
To program G-J elimination method, st make a working matrix by
augmenting A matrix with b (n rows, n + 1 columns).
e.g.
w =
_
_
A b
_
_
=
_
_
2 3 5 10
1 2 4 5
4 7 3 15
_
_
C program:
void simple_gauss(a, b, x, n)
double **a, b[], x[]; int n;
{
double **work;
double m; int i, j, k;
/* working matrix */
work = (double **)calloc(n, sizeof(*work));
for(i=0; i<n; i++)
*(work+i) = (double *)calloc(n+1, sizeof(**work));
for(i=0; i<n; i++){
for(j=0; j<n; j++)
*(*(work+i)+j) = *(*(a+i)+j);
*(*(work+i)+n) = b[i];
}
/* G-J elimination */
for(i=0; i<n; i++)
for(j=0; j<n; j++){
if(j != i){
m = *(*(work+j)+i)/*(*(work+i)+i);
for(k=i; k<=n; k++)
*(*(work+j)+k) -= m* *(*(work+i)+k);
}
}
/* compute x */
for(i=0; i<n; i++)
x[i] = *(*(work+i)+n)/*(*(work+i)+i);
}
118
problem in G-J elimination:
this method fails if any of the pivot elements is zero.
Select a nonzero pivot element at each step (partial pivoting).
Select the pivot element with the greatest absolute magnitude
(for the numerical stability).
If all of candidate pivot elements are zero,
no solution or an innite number of solutions exists.
other pivoting methods (in order to provide numerical stability)
implicit pivoting:
rescale all equations so that they are all same order of magnitude
use partial pivoting
the nal solution is rescaled back to the original units
full pivoting :
exchanges of both rows and columns are considered in order to
select the pivot element.
(exchanging columns permutes the subscripts in vector x.)
the subscripts in the nal solution should be reordered.
119
matrix inversion:
G-J elimination can be used to nd the matrix inverse
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
_
x
0
x
1
x
n1
_
_
_
_
_
_
=
_
_
_
_
_
_
_
1 0 0
0 1 0
0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 1
_
_
_
_
_
_
_
A A
1
= I
Ax
0
=
_
_
_
_
_
1
0
.
.
.
0
_
_
_
_
_
, Ax
1
=
_
_
_
_
_
0
1
.
.
.
0
_
_
_
_
_
, , Ax
n1
=
_
_
_
_
_
0
.
.
.
0
1
_
_
_
_
_
make a working matrix such as
_
_
[
A [ I
[
_
_
120
18.2 LU Decomposition Methods
Any non-singular matrix A can be decomposed to the LU form.
A = LU
where L is a lower triangular matrix and U is a upper triangular matrix.
eg.
_
_
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
_
_
=
_
_
1 0 0
l
21
1 0
l
31
l
32
1
_
_
_
_
u
11
u
12
u
13
0 u
22
u
23
0 0 u
33
_
_
where l
ii
= 1.
algorithm
1. the rst row of U:
u
1j
= a
1j
, j = 1, . . . , N.
the rst column of L:
l
i1
= a
i1
/u
11
, i = 2, . . . , N.
2. the nth row of U (n 2):
u
nj
= a
nj

n1

k=1
l
nk
u
kj
, j = n, . . . , N.
the nth column of L (n 2):
l
in
=
_
a
in

n1

k=1
l
ik
u
kn
_
/u
nn
, i = n + 1, . . . , N.
121
eg.
A =
_
_
2 1 3
1 3 2
3 1 3
_
_
u
11
= 2, u
12
= 1, u
13
= 3,
l
11
= 1, l
21
= 0.5, l
31
= 1.5,
u
22
= 3 (1/2) 1 = 3.5, u
23
= 2 (1/2)(3) = 0.5,
l
22
= 1, l
32
= [1 (1.5) 1]/3.5 = 0.14,
u
33
= 3 (1.5)(3) (0.14)(0.5) = 1.57.
Therefore,
L =
_
_
1 0 0
0.5 1 0
1.5 0.14 1
_
_
and
U =
_
_
2 1 3
0 3.5 0.5
0 0 1.57
_
_
.
122
Solving simultaneous equations
Change simultaneous equations of Ax = y with LUx = y.
Let Ux = z, then Lz = y.
Solve z Solve x.
eg.
- forward elimination:
_
_
1 0 0
l
21
1 0
l
31
l
32
1
_
_
_
_
z
1
z
2
z
3
_
_
=
_
_
y
1
y
2
y
3
_
_

z
1
= y
1
z
2
= y
2
z
1
l
21
z
3
= y
3
z
1
l
31
z
2
l
32
- backward substitution:
_
_
u
11
u
12
u
13
0 u
22
u
23
0 0 u
33
_
_
_
_
x
1
x
2
x
3
_
_
=
_
_
z
1
z
2
z
3
_
_

x
3
= z
3
/u
33
x
2
= (z
2
u
23
x
3
)/u
22
x
1
= (z
1
u
12
x
2
u
13
x
3
)/u
11
algorithm
1. Forward elimination step:
z
1
= y
1
,
z
i
= y
i

_

i1
j=1
l
ij
z
j
_
, i = 2, . . . , N.
2. Backward substitution step:
x
N
= z
N
/u
NN
,
x
i
= (z
i

N
j=i+1
u
ij
x
j
)/u
ii
, i = N 1, . . . , 1.
123
Computing the determinant
- Determinant of matrix A:
det(A) =

a
i1
a
j2
a
k3
a
rN
where the summation is extended over
all permutations of the rst subscripts of a.
eg.
det(A) = det
_
_
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
_
_
= a
11
a
22
a
33
+a
21
a
32
a
13
+a
31
a
12
a
23
a
11
a
32
a
23
a
21
a
12
a
33
a
31
a
22
a
13
.
5 5 matrix: 5! = 120 terms (each term need 4 multiplications)
10 10 matrix: 10! 2 10
8
terms (each term need 9 multiplications)
- Two important rules of determinants:
1. det(BC) = det(B)det(C)
2. det(M) = the product of all diagonal elements of M
if M is an upper or lower triangular matrix:
det(A) = det(LU) = det(L)det(U) = det(U)
i.e. det(A) =

N
i=1
u
ii
- Checking ill-conditioned problems
1. A number of problems which are solvable, yet their solutions become
very inaccurate because of severe round-o errors.
2. Compare det(A)det(A
1
) with 1.
3. Compare A
1
(A
1
)
1
with AA
1
.
If there is signicant deviation,
increase the precision of computation.
124
18.3 Solutions of m equations with n unknowns
Compute the solution of Ax = y:
Case 1. if m = n and rank(A) = n or det(A) ,= 0,
unique solution of x exists.
Case 2. if m < n, rank(A) < n and det(A) = 0,
A: singular matrix.
no unique solution, innite number of solutions.
eg. x +y = 1, m = 1, n = 2
Solution:
x = 1 y or y = 1 x
if we have m linearly independent equations for n unknowns and m < n,
we have m basic variable n m free variables.
eg. G-J elimination for m < n:
2u + 3v + 1w + 4x +y = 6
2u + 3v + 1w + 1x y = 1
4u + 6v 1w + 1x + 2y = 5
m = 3, n = 5
u v w x y RHS
2 3 1 4 1 6
2 3 1 1 1 1
4 6 1 1 2 5
G-J elimination:
u v w x y RHS
1 1.5 0 0 0.05 0.4
1 0 1.5 1.5
1 0.6 1.6
125
basic variables: u, w, x, free variables: v, y.
_
u = 0.4 1.5v + 0.05y
w = 1.5 + 1.5y
x = 1.6 0.6y
Case 3. if m > n, no solution. We need to solve the approximate solution.
eg. linear regression
N samples are given:
(x
1
, y
1
), (x
2
, y
2
), . . . , (x
N
, y
N
).
The samples are generated by
y
N
=
L

i=0
w
i
x
ik
+
where
y
k
: the output for the kth sample,
w
i
: the ith weight, and
x
ik
: the ith input for the kth sample.
: the random noise, usually Gaussian noise having
the distribution of N(0,
2
).
if L = 2, the estimation of output for the kth sample:
y
k
= w
0
+w
1
x
1k
126
y
k
= w
T
x
k
: the estimation of output for the kth sample
where x
k
= [x
0k
, x
1k
, , x
Lk
]
T
and w = [w
0
, w
1
, , w
L
]
T
.
the error for the kth pattern:
e
k
= y
k
y
k
= y
k
w
T
x
k
We want to estimate w such that minimizing the mean square error (MSE):
E[e
2
k
] =
_
e
2
k
p(y[x)dy.
In practice, we use
1
N
N

k=1
e
2
k
instead of E[e
2
k
].
MSE = E[e
2
k
] = E[(y
k
y
k
)
2
]
= E[(y
k
w
T
x
k
)
2
]
= E[y
2
k
] +w
T
E[x
k
x
T
k
]w 2E[y
k
x
T
k
]w
Let R = E[x
k
x
T
k
] (input correlation matrix) and P = E[y
k
x
k
].
Then,
MSE = E[e
2
k
] = E[y
2
k
] +w
T
Rw 2P
T
w.
E
w
[
w=w
= 0
2Rw

2P = 0
w

= R
1
P
minimum mean square error (MMSE):
MMSE = E[y
2
k
] +w
T
Rw

2P
T
w

= E[y
2
k
] P
T
w

127
regression algorithm:
1. Read (x
1
, y
1
), , (x
N
, y
N
).
2. Compute R and P:
R =
1
N

N
k=1
x
k
x
T
k
P =
1
N

N
k=1
y
k
x
k
3. Compute R
1
by G-J elimination method
(or LU decomposition method).
4. Compute w

:
w

= R
1
P
128
19 Random Number Generation
19.1 Creating Uniform Random Numbers
linear congruential method
- The i + 1th random number (integer) is generated by
r
i+1
= (ar
i
+b)%m (0, , m1)
where a, b and m are xed constants.
- r
0
: seed, m: range.
eg.
a = 13, b = 1, and m = 17
r
0
= 1
r
1
= (13 1 + 3)%17 = 16
r
2
= (13 16 + 3)%17 = 7
r
3
= (13 7 + 3)%17 = 9( cycle length=4)
r
4
= (13 9 + 3)%17 = 1
.
.
.
- Select values of a, b and m so that the cycle length is
as long as possible.
- Independent random number and uniform distribution.
- No absolute set of rules for selecting the values of a, b and m.
a typical form:
a = 8x + 5
where x is a positive integer and b is an odd number.
129
ANSI C Standard(1988): rand() generate the random number
between 0 and RAND MAX (32767)
static unsigned long next=1;
int rand(void)
{
next = next*1103515245 + 12345;
return((unsigned long) (next/65536)%32768);
}
eg. generating random numbers between 0 and INT MAX:
#include <limits.h>
static int next=1;
int rand(void)
{
next = next*84314861 + 453816693;
if (next < 0)
next += INT_MAX + 1;
return (next);
}
eg. setting the seed:
default: srand(1)
void srand(int seed)
{
next = seed;
}
130
additive congruential method
1. Make a table of random integer.
2. Select the two table entry except the given entry.
3. Update the random number of the given table entry by adding
two random numbers of the selected table entries.
better property of uniform random numbers.
#include <stdlib.h>
#define TSIZE 55
static int table_index;
static unsigned int table[ISIZE];
void rinit (void)
{
for (table_index=0; table_index<TSIZE; table_index++)
table[table_index] = rand();
}
#include <limits.h>
#define B 23
#define D 54
int rand_ac(void)
{
table_index = (table_index+1)%TSIZE;
table[table_index] = (table[(table_index+B)%TSIZE]+
table[(table_index+D)%TSIZE])%INT_MAX;
return(table[table_index]);
}
131
Monte-Carlo integration
Step 1. Generate the random number between [a, b] and between [0, A].
Step 2. Count the number under f(x) : C.
Step 3. Compute the integrated value as
Area =
C
N
A(b a)
where N is total random numbers.
- The basic theorem of Monte Carlo integration:
_
fdV < f > V V
_
< f
2
> < f >
2
N
where
< f >=
1
N
N

i=1
f(x
i
)
< f
2
>=
1
N
N

i=1
f
2
(x
i
)
- error bound is given by O(N
1/2
).
cf. If we are use quasi-random numbers instead of pseudo-random numbers,
the error bound is redundant to O(N
1
). (refer Numerical Recipes in C).
132
eg. Monte-Carlo integration to approximate .
main()
{
int i, count=0, maxdraws;
double x, y, simple_uniform();
printf("Enter number of points to be drawn: \n");
scanf("%d", &maxdraws);
for (i=1; i<=maxdraws; i++) {
x = simple_uniform();
y = simple_uniform();
if (x*x + y*y < 1.0) count++;
}
printf("Estimate of PI is %lf\n", 4*count/maxdraws);
}
double simple_uniform(void)
{
return ((double) rand()/RAND_MAX);
}
133
19.2 Creating Non-uniform Random Numbers
inverse cumulative method
- Let f
X
(x) be probability density function.
- The cumulative distribution function is dened by
F
X
(x) = Prob[X x] =
_
x
0
f(x)dx
We want to generate random values with F
X
(x):
Step 1. Draw a random value u between 0 and 1 (uniform distribution).
Step 2. Compute x = F
1
X
(u).
F
1
X
(u) always exists since F
X
(u) is monotonic increasing function.
This method is possible when the inverse of F
X
(x) is known.
double rand_icm(f_inverse)
double (*f_inverse)();
{
double u;
u = (double) rand()/RAND_MAX;
return (f_inverse(u));
}
134
construction method
- Generating random variables with dierent density or mass functions using
known mathematical relationship.
eg.
random variable X: uniformly distributed between 0 and 1.
Y = g(x) = a + (b a)X
Y will be uniformly distributed between a and b.
double uniform(a, b)
double a, b;
{
return ((double) rand()/RAND_MAX*(b-a)+a);
}
135
eg. Bernoulli drawing
Let Prob[X = 1] = p and Prob[X = 0] = 1 p.
Toss a fair coin: p =
1
2
r heads in a sample of n independent coin tosses
a binomial distribution:
P(r) = nC
r
p
r
(1 p)
nr
Let Y = X
1
+X
2
+ +X
N
(no. of heads in n trials). Then,
E[Y ] = np
V ar(Y ) = np(1 p)
If n is large, binomial distribution is closely approximated by
normal distribution (eg. np(1 p) 5).
int bernoulli(double p)
{
if ((double) rand()/RAND_MAX < p) return(1);
else return(0);
}
int binomial(double p, int n)
{
int i, count=0;
for (i=0; i<n; i++)
count += bernoulli(p);
return(count);
}
136
eg. Gaussian distribution
Let X
1
, X
2
, , X
N
be independent random variables that are
uniformly distributed between
1
2
and
1
2
and
Y =
N

i=1
X
i
.
Then, Y is approximated as a random variable with Gaussian
distribution having
E[Y ] = 0 and
V ar(Y ) =
N
12
cf. the central limit theorem:
Let random variables X be of continuous type and independent with
E[X
i
] =
i
and
2
X
i
=
2
i
.
Then,
Y = X
1
+ +X
N
has mean
Y
and variance
2
Y
given by

Y
= E[Y ] =
1
+
2
+ +
N
,

2
Y
= V ar(Y ) =
2
1
+
2
2
+ +
2
N
,
and its density by
f(Y ) = f
1
(x) f
N
(x)
where is convolution operator, that is, i.e.,
f
1
(x) f
2
(x) =
_

f
1
(x

)f
2
(x x

)dx

.
137
double standard_gaussian(void)
{
int i;
double sum = 0.0;
for (i=1; i<=12; i++)
sum += (double) rand()/RAND_MAX;
return (sum - 6.0);
}
double general_gaussian (double m, double sigma)
{
return (standard_gaussian()*sigma + m);
}
rejection method
- same idea as Monte-Carlo simulation.
Step 1. Identify the largest value of f
X
(x), f
max
and the interval [x
min
, x
max
].
Step 2. Draw a uniform distributed random value x

in [x
min
, x
max
].
Step 3. Draw a uniform distributed random value y

in [0, f
max
].
Step 4. If y

f(x

), return x

. Otherwise, go to Step 2.
This method is good when the cumulative distribution function has
no closed form or when it can not be inverted.
double rejection(f, fmax, xmin, xmax)
double (*f)(), fmax, xmin, xmax;
{
double xstat = uniform(xmin, xmax);
while (bernoulli(f(xstat)/fmax) == 0)
xstat = uniform(xmin, xmax);
return (xstat);
}
138