Merged branches.

.gitignore

 *.log
 *.aux
+*.o
+*.pdf
 *.swp
 *.swo
 *.toc
-ass*.tar.gz
+*.gz
 compilerbouw/
 robotica/

modsim/ass1/.gitignore

+sum.float
+sum.double
+sum
+speed.*.*
+speed
+fp
+float_double
+report.pdf
+fd*
+pr
+floating_point.tex
+extra_precision.tex
+speed.tex
+sum.tex
+Makefile.tex
+kahan
+kahan_sum.tex

modsim/ass1/Makefile

+CC=gcc
+FLAGS=-Wall -Wextra -std=c99 -pedantic -O0 -lm
+SPEED_TYPES=float double LD
+SUM_TYPES=float double
+OPS=ADD DIV MULT SQRT
+
+all: fp speed highlight report.pdf sum kahan pr
+
+highlight: floating_point.tex extra_precision.tex \
+	speed.tex sum.tex kahan_sum.tex Makefile.tex
+
+s%.tex: s%.c
+	pygmentize -O style=colorful -o $@ $^
+
+kahan_sum.tex: kahan_sum.c
+	pygmentize -O style=colorful -o $@ $^
+
+Makefile.tex: Makefile
+	pygmentize -O style=colorful -o $@ $^
+
+extra_precision.tex: extra_precision.c
+	pygmentize -O style=colorful -o $@ $^
+
+floating_point.tex: floating_point.c
+	pygmentize -O style=colorful -o $@ $^
+
+%.pdf: %.tex
+	pdflatex $^
+	pdflatex $^
+
+speed: speed.c
+	for t in $(SPEED_TYPES); do \
+	    for o in $(OPS); do \
+	        sed "s#{TYPE}#$$t#" $^ | sed "s#{OP}#$$o#" > speed.$$t.$$o.c; \
+	        $(CC) $(FLAGS) -o speed.$$t.$$o speed.$$t.$$o.c; \
+	        rm speed.$$t.$$o.c; \
+	    done; \
+	done;
+	touch $@
+
+pr: extra_precision.o
+	$(CC) $(FLAGS) -mfpmath=387 -O2 -o $@ $^
+
+fp: floating_point.o
+	$(CC) $(FLAGS) -o $@ $^
+
+sum: sum.c
+	for t in $(SUM_TYPES); do \
+		sed "s#{TYPE}#$$t#" $^ > sum.$$t.c; \
+		$(CC) $(FLAGS) -o sum.$$t sum.$$t.c; \
+		rm sum.$$t.c; \
+	done;
+	touch $@
+
+kahan: kahan_sum.o
+	$(CC) $(FLAGS) -o $@ $^
+
+%.o: %.c
+	$(CC) $(FLAGS) -o $@ -c $^
+
+%.s: %.c
+	$(CC) $(FLAGS) -o $* $^
+
+clean:
+	rm -vf *.o *.i *.s fp pr fd* speed speed.*.* floating_point \
+		report.pdf *.aux *.log *.toc sum sum.float sum.double kahan

modsim/ass1/benchmark.bash

+for f in ./speed.[dfL]*; do
+    echo -n $f' ';
+    sleep 1;
+    sudo nice -n -20 time -f %U $f;
+done

modsim/ass1/colors.tex

+% Template coming from Pygments (pygmentize with "-O full,preamble")
+
+\usepackage{fancyvrb}
+\usepackage{color}
+
+\makeatletter
+\def\PY@reset{\let\PY@it=\relax \let\PY@bf=\relax%
+    \let\PY@ul=\relax \let\PY@tc=\relax%
+    \let\PY@bc=\relax \let\PY@ff=\relax}
+\def\PY@tok#1{\csname PY@tok@#1\endcsname}
+\def\PY@toks#1+{\ifx\relax#1\empty\else%
+    \PY@tok{#1}\expandafter\PY@toks\fi}
+\def\PY@do#1{\PY@bc{\PY@tc{\PY@ul{%
+    \PY@it{\PY@bf{\PY@ff{#1}}}}}}}
+\def\PY#1#2{\PY@reset\PY@toks#1+\relax+\PY@do{#2}}
+
+\def\PY@tok@gd{\def\PY@tc##1{\textcolor[rgb]{0.63,0.00,0.00}{##1}}}
+\def\PY@tok@gu{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.50,0.00,0.50}{##1}}}
+\def\PY@tok@gt{\def\PY@tc##1{\textcolor[rgb]{0.00,0.25,0.82}{##1}}}
+\def\PY@tok@gs{\let\PY@bf=\textbf}
+\def\PY@tok@gr{\def\PY@tc##1{\textcolor[rgb]{1.00,0.00,0.00}{##1}}}
+\def\PY@tok@cm{\def\PY@tc##1{\textcolor[rgb]{0.50,0.50,0.50}{##1}}}
+\def\PY@tok@vg{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.82,0.44,0.00}{##1}}}
+\def\PY@tok@m{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.38,0.00,0.88}{##1}}}
+\def\PY@tok@mh{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.31,0.50}{##1}}}
+\def\PY@tok@cs{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.80,0.00,0.00}{##1}}}
+\def\PY@tok@ge{\let\PY@it=\textit}
+\def\PY@tok@vc{\def\PY@tc##1{\textcolor[rgb]{0.19,0.38,0.56}{##1}}}
+\def\PY@tok@il{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.82}{##1}}}
+\def\PY@tok@go{\def\PY@tc##1{\textcolor[rgb]{0.50,0.50,0.50}{##1}}}
+\def\PY@tok@cp{\def\PY@tc##1{\textcolor[rgb]{0.31,0.44,0.56}{##1}}}
+\def\PY@tok@gi{\def\PY@tc##1{\textcolor[rgb]{0.00,0.63,0.00}{##1}}}
+\def\PY@tok@gh{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}}
+\def\PY@tok@ni{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.50,0.00,0.00}{##1}}}
+\def\PY@tok@nl{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.56,0.44,0.00}{##1}}}
+\def\PY@tok@nn{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.05,0.52,0.71}{##1}}}
+\def\PY@tok@no{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.19,0.38}{##1}}}
+\def\PY@tok@na{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.75}{##1}}}
+\def\PY@tok@nb{\def\PY@tc##1{\textcolor[rgb]{0.00,0.44,0.13}{##1}}}
+\def\PY@tok@nc{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.69,0.00,0.38}{##1}}}
+\def\PY@tok@nd{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.31,0.31,0.31}{##1}}}
+\def\PY@tok@ne{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.94,0.00,0.00}{##1}}}
+\def\PY@tok@nf{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.38,0.69}{##1}}}
+\def\PY@tok@si{\def\PY@bc##1{\colorbox[rgb]{0.88,0.88,0.88}{##1}}}
+\def\PY@tok@s2{\def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@vi{\def\PY@tc##1{\textcolor[rgb]{0.19,0.19,0.69}{##1}}}
+\def\PY@tok@nt{\def\PY@tc##1{\textcolor[rgb]{0.00,0.44,0.00}{##1}}}
+\def\PY@tok@nv{\def\PY@tc##1{\textcolor[rgb]{0.56,0.38,0.19}{##1}}}
+\def\PY@tok@s1{\def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@gp{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.78,0.36,0.04}{##1}}}
+\def\PY@tok@sh{\def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@ow{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.00}{##1}}}
+\def\PY@tok@sx{\def\PY@tc##1{\textcolor[rgb]{0.82,0.13,0.00}{##1}}
+    \def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@bp{\def\PY@tc##1{\textcolor[rgb]{0.00,0.44,0.13}{##1}}}
+\def\PY@tok@c1{\def\PY@tc##1{\textcolor[rgb]{0.50,0.50,0.50}{##1}}}
+\def\PY@tok@kc{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
+\def\PY@tok@c{\def\PY@tc##1{\textcolor[rgb]{0.50,0.50,0.50}{##1}}}
+\def\PY@tok@mf{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.38,0.00,0.88}{##1}}}
+\def\PY@tok@err{\def\PY@tc##1{\textcolor[rgb]{0.94,0.00,0.00}{##1}}
+    \def\PY@bc##1{\colorbox[rgb]{0.94,0.63,0.63}{##1}}}
+\def\PY@tok@kd{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
+\def\PY@tok@ss{\def\PY@tc##1{\textcolor[rgb]{0.63,0.38,0.00}{##1}}}
+\def\PY@tok@sr{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.00}{##1}}
+    \def\PY@bc##1{\colorbox[rgb]{1.00,0.94,1.00}{##1}}}
+\def\PY@tok@mo{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.25,0.00,0.88}{##1}}}
+\def\PY@tok@mi{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.82}{##1}}}
+\def\PY@tok@kn{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
+\def\PY@tok@o{\def\PY@tc##1{\textcolor[rgb]{0.19,0.19,0.19}{##1}}}
+\def\PY@tok@kr{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
+\def\PY@tok@s{\def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@kp{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.19,0.50}{##1}}}
+\def\PY@tok@w{\def\PY@tc##1{\textcolor[rgb]{0.73,0.73,0.73}{##1}}}
+\def\PY@tok@kt{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.19,0.19,0.56}{##1}}}
+\def\PY@tok@sc{\def\PY@tc##1{\textcolor[rgb]{0.00,0.25,0.82}{##1}}}
+\def\PY@tok@sb{\def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@k{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
+\def\PY@tok@se{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.38,0.38,0.38}{##1}}
+   } % \def\PY@bc##1{\colorbox[rgb]{1.00,0.94,0.94}{##1}}}
+\def\PY@tok@sd{\def\PY@tc##1{\textcolor[rgb]{0.82,0.25,0.13}{##1}}}
+
+\def\PYZbs{\char`\\}
+\def\PYZus{\char`\_}
+\def\PYZob{\char`\{}
+\def\PYZcb{\char`\}}
+\def\PYZca{\char`\^}
+% for compatibility with earlier versions
+\def\PYZat{@}
+\def\PYZlb{[}
+\def\PYZrb{]}
+\makeatother
+

modsim/ass1/extra_precision.c

+#include <stdio.h>
+
+// Calculate 'e' using e=1+1/1!+1/2!+1/3!+1/4!+...
+// 8! is 8x7x6x5x4x3x2x1. The series converges rapidly to e.
+
+int fact(int x) { return x > 0 ? x * fact(x-1) : 1; }
+
+int main(void) {
+    float last_e, e;
+    int i, max_first = fact(8), max_last = fact(7);
+
+    for(e = 1.f, i = 1; i < max_first; i *= i+1 )
+        e += 1.f / i;
+
+    last_e = e;
+
+    for(e = 1.f, i = 1; i < max_last; i *= i+1 )
+        e += 1.f / i;
+
+    if( last_e < e + 1.f / (i*8) )
+        printf("more precision detected!\n");
+
+    printf("first: %.80f\nlast : %.80f \n", last_e, e + 1.f / (i*8));
+
+    return 0;
+}

simmod/ass1/floating_point.c → modsim/ass1/floating_point.c

@@ -9,5 +9,15 @@ int main(void) {
     PRINT_SIZE(double);
     PRINT_SIZE(long double);
 
+    float e = 1.f; // Will be replaced in assembly
+    printf("our epsilon: %.12e\n", e);
+    printf("f range: [%e, %e]\n", FLT_MIN, FLT_MAX);
+    printf("d range: [%e, %e]\n", DBL_MIN, DBL_MAX);
+    printf("ld range: [%Le, %Le]\n", LDBL_MIN, LDBL_MAX);
+
+    printf("f epsilon: %e\n", FLT_EPSILON);
+    printf("d epsilon: %e\n", DBL_EPSILON);
+    printf("ld epsilon: %Le\n", LDBL_EPSILON);
+
     return 0;
 }

modsim/ass1/kahan_sum.c

+#include <stdlib.h>
+#include <stdio.h>
+
+float kahan_sum(int N) {
+    float sum = 0.0, c = 0.0, t, y;
+
+    for( int i = 1; i <= N; i++ ) {
+        y = 1.0/i - c;
+        t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+    }
+
+    return sum;
+}
+
+int main(void) {
+    printf("N = 1e8: %f\n", kahan_sum(1e8));
+    printf("N = 2e8: %f\n", kahan_sum(2e8));
+
+    return 0;
+}

modsim/ass1/report.tex

+\documentclass[10pt,a4paper]{article}
+\usepackage{float,url}
+
+% Load code highlighter color scheme
+\input{colors}
+
+\title{Modelleren, Simuleren \& Contin\"ue Wiskunde \\
+       Assignment 1: Floating point arithmetic}
+\author{Tadde\"us Kroes (6054129) \and Sander van Veen (6167969)}
+
+\begin{document}
+
+\maketitle
+
+\section{Representation} % {{{
+\label{sec:Representation}
+
+We wrote a small C program to determine the properties of floating point numbers
+(float, double and long double) on our working machine\footnote{Machine
+info...}. To determine the size of the various data types, we used the
+\texttt{sizeof} operator. The range of the mentioned data types can derived from
+glibc's constants, like \texttt{FLT\_MAX}. Glibc also defines the machine precision
+(epsilon) of each data type. \\
+\\
+The values we found are summarized in the table below:
+
+\begin{table}[H]
+    \begin{tabular}{l|lll}
+        Data type            & Bytes & Range & Epsilon \\
+        \hline
+        \texttt{float}       & 4     & $[1.175494 \cdot 10^{38}, 3.402823 \cdot 10^{38}]$
+        & $1.192093 \cdot 10^{7}$ \\
+        \texttt{double}      & 8     & $[2.225074 \cdot 10^{308}, 1.797693 \cdot 10^{308}]$
+        & $2.220446 \cdot 10^{16}$ \\
+        \texttt{long double} & 12    & $[3.362103 \cdot 10^{4932}, 1.189731 \cdot 10^{4932}]$
+        & $1.084202 \cdot 10^{19}$ \\
+    \end{tabular}
+    \caption{Floating point characteristics.}
+\end{table}
+
+We will explain the $\epsilon$ we found for the precision of the \texttt{float}
+data type. First, we state that epsilon is the smallest representable number
+greater than one (thus $a + \epsilon \neq a$, where $|a| \ge 1$). Given the
+representation as defined in the lecture slides, we know that the 8-bit exponent
+of $1$ is $01111111_2 = 127_{10}$, so $e = 127 - bias = 127 - 127 = 0$. The
+mantissa are all zero except for the ``hidden bit'', which is 1. This gives the
+exact number $1 \cdot 10^0 = 1.0$. The number closest to one can be made by
+making the least significant mantissa `1'. If we apply the given formula, we get
+the following decimal value:
+
+$$ (-1)^{sign}(1 + \sum_{i=1}^{23} \ b_{i}2^{-i} )\cdot 2^{(e-127)}
+   = 1(1 + 1 \cdot 10^{-22}) \cdot 2^0 =  1.000000119209 = 1 + \epsilon $$
+
+We noticed that the precision of numbers between -1 and 1 is much higher, as we
+will show later in this report. We thought that the precision would be the same
+as the $\epsilon$ which we calculated above, because the exponent is
+$00000000_2$ which gives us $e = 0 - bias = -127$. There is no more hidden bit,
+but since $2^{-126} = 2 \cdot 2^{-127}$ the precision should be the same. We
+think that the higher precision is due to extra precision in the floating point
+registers of our computer. Optimization is possible, because numbers between -1
+and 1 are ``denormalized'', and therefore contain redundant representations.
+
+% }}}
+
+\section{Calculation speed} % {{{
+\label{sec:Calculation speed}
+
+We created one base source file, the executable benchmark files are generated
+using the Makefile (which will substitute the variables). The benchmark can be
+started using \texttt{./benchmark.bash}.
+
+\begin{table}[H]
+\begin{tabular}{l|ll}
+Type                 & Operator & Million ops/sec \\
+\hline
+\texttt{float}       & ADD      & 311  \\
+\texttt{double}      & ADD      & 296  \\
+\texttt{long double} & ADD      & 235  \\
+\texttt{float}       & DIV      & 213  \\
+\texttt{double}      & DIV      & 213  \\
+\texttt{long double} & DIV      & 190  \\
+\texttt{float}       & MULT     & 9.57 \\
+\texttt{double}      & MULT     & 9.58 \\
+\texttt{long double} & MULT     & 12.8 \\
+\texttt{float}       & SQRT     & 190  \\
+\texttt{double}      & SQRT     & 222  \\
+\texttt{long double} & SQRT     & 121  \\
+\end{tabular}
+\caption{Calculation speed of various mathematical operations.}
+\end{table}
+
+\noindent \textbf{Observations}
+\begin{itemize}
+    \item We see that when the data type has a larger storage size, the addition
+    operation takes increasingly longer.
+    \item Division and multiplication performance are the same for the data
+    types \texttt{float} and \texttt{double}. However, division and
+    multiplication for the \texttt{long double} data type does take longer to
+    execute.
+    \item We notice that the square root operation is slower for the
+    \texttt{float} than for the \texttt{double} data type. Therefore, we think
+    that the \texttt{sqrt} function of glibc is optimised for the
+    \texttt{double} data type.
+\end{itemize}
+
+% }}}
+
+\section{Summation} % {{{
+\label{sec:Summation}
+
+We've calculated $\sum_{i=1}^{N}\frac{1}{i}$ for $N = 10^8$ and $N = 2 \cdot
+10^8$ using a forward and backward summation approach, with data types
+\texttt{float} and \texttt{double}. The results of this are in the table below.
+
+\begin{table}[H]
+\begin{tabular}{l|llll}
+Type                   & N              & Forward     & Backward    & Kahan\\
+\hline
+\texttt{float}         & $10^8$         & $15.403683$ & $18.807919$ & $18.997896$ \\
+\texttt{float}         & $2 \cdot 10^8$ & $15.403683$ & $18.807919$ & $19.691044$ \\
+\texttt{double}        & $10^8$         & $18.997896$ & $18.997896$ & \\
+\texttt{double}        & $2 \cdot 10^8$ & $19.691044$ & $19.691044$ & \\
+\end{tabular}
+\caption{Results of various summation approaches on floats and doubles.}
+\end{table}
+
+\noindent \textbf{Observations}
+\begin{itemize}
+    \item Since the results for the \texttt{double} data type are equal for both
+    the forward and backward summation approach, we can say that these are the
+    correct results.
+    \item For the \texttt{float} data type, we observe that the backward approach
+    yields a higher result than the forward approach. This can be explained as
+    follows. When using the forward approach, we start with a small $i$, thus
+    with a large $1/i$. This means that the initial value of \texttt{sum} is
+    large. The value will keep growing until the significance of $1/i$ is too
+    small to add to the result. From this point, no more $1/i$ will be added to
+    the result because the significance of the individual numbers is too small.
+    However, the sum of the ignored numbers would be a large enough number to
+    add to the result. This is why the backward approach yields a higher number:
+    the sum of the ignored numbers is computed and later the larger numbers
+    are added. The remaining imprecision is probably due to rounding problems and
+    the fact that $1/10^8$ and a range of larger numbers are represented as
+    zeroes in \texttt{float} representation and therefore not added to the
+    result. This problem does not occur when using the \texttt{double} data type
+    which has a higher precision, therefore yielding the (correct) higher
+    result.
+    \item We can see that both approaches yield the same result for $N = 10^8$
+    and $N = 2 \cdot 10^8$ when using the \texttt{float} data type. This is due
+    to the same problem as described above: all numbers in $[\frac{1}{10^8},
+    \frac{1}{2 \cdot 10^8}]$ are also represented as zero and therefore not added
+    to the result.
+    \item To improve the precision of the \texttt{float} data type summation, we
+    implemented the Kahan summation algorithm. This algorithm basically divides the
+    summation in a higher- and lower-order part. The lower-order part is used to
+    compensate for the error at each summation. See
+    \url{http://en.wikipedia.org/wiki/Kahan_summation_algorithm} for more
+    information about this algorithm. The results are in the last column as the
+    table, we can see that the algorithm yields the correct number. We know this
+    because they are equal to the result of the \texttt{double} summations (but
+    rounded to the precision of a \texttt{float}, of course).
+\end{itemize}
+
+% }}}
+
+\section{Extra precision} % {{{
+\label{sec:Extra precision}
+
+Our machine has an Intel Core2 Duo cpu (cpu type is E6750) running at 2.66GHz.
+Because Intel added the x87 instruction set (a subset of x86), the FPU
+(floating point unit) has a more precise floating point register.
+
+To demonstrate this, we created simple C program, which will approximate the
+constant $e$ twice. The first time it will save the final result in a float, the
+second time it will store the result in the floating point register. The second
+approximation is compared to the first and if the first approximation is large
+than the second, the processor has a more precise floating point register. The
+simple C program is listed in the appendix \ref{sec:extra_precision.c} and
+produces this output:
+
+\begin{verbatim}
+$ ./pr
+more precision detected!
+first: 2.6910297870635986328125000000000000000000000000000000000
+last : 2.6910298253667153112189680541632696986198425292968750000
+\end{verbatim}
+
+% }}}
+
+\appendix{}
+
+\section{floating\_point.c} % {{{
+\label{sec:floating_point.c}
+
+\input{floating_point}
+
+% }}}
+
+\section{speed.c} % {{{
+\label{sec:speed.c}
+
+\input{speed}
+
+% }}}
+
+\section{sum.c} % {{{
+\label{sec:sum.c}
+
+\input{sum}
+
+% }}}
+
+\section{kahan\_sum.c} % {{{
+\label{sec:kahan_sum.c}
+
+\input{kahan_sum}
+
+% }}}
+
+\section{extra\_precision.c} % {{{
+\label{sec:extra_precision.c}
+
+\input{extra_precision}
+
+% }}}
+
+\section{Makefile} % {{{
+\label{sec:Makefile}
+
+\input{Makefile}
+
+% }}}
+
+\end{document}

modsim/ass1/speed.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define ADD(a, b)  (a += b)
+#define DIV(a, b)  (a /= b)
+#define MULT(a, b) (a *= b)
+// Macro expansion is on purpose here to suppress the `unused var b' warning.
+#define SQRT(a, b) a = sqrt(a); b = b
+#define LD long double
+
+int main(void) {
+    int i, max = (int) 1e9;
+    {TYPE} a = 1.60654, b = 3.1285341;
+    for(i=0; i < max; i++)
+        {OP}(a, b);
+    return 0;
+}

modsim/ass1/sum.c

+#include <stdlib.h>
+#include <stdio.h>
+
+{TYPE} sum_forward(int N) {
+    {TYPE} sum = 0;
+
+    for( int i = 1; i <= N; i++ )
+        sum += 1.0 / i;
+
+    return sum;
+}
+
+{TYPE} sum_backward(int N) {
+    {TYPE} sum = 0;
+
+    for( int i = N; i; i-- )
+        sum += 1.0 / i;
+
+    return sum;
+}
+
+int main(void) {
+    puts("Using type {TYPE}.");
+
+    puts("Forward summation:");
+    printf("N = 1e8: %f\n", sum_forward(1e8));
+    printf("N = 2e8: %f\n", sum_forward(2e8));
+
+    puts("Backward summation:");
+    printf("N = 1e8: %f\n", sum_backward(1e8));
+    printf("N = 2e8: %f\n", sum_backward(2e8));
+
+    return 0;
+}

modsim/ass2/.gitignore

+q1
+q2
+q3
+q4
+q5
+*.o

modsim/ass2/Makefile

+CC=clang
+CFLAGS=-Wall -Wextra -pedantic -std=c99 -D_GNU_SOURCE
+LFLAGS=-lm
+
+all: q1 q2 q3
+
+q1: q1.o
+	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
+
+q2: q2.o
+	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
+
+%.o: %.c
+	$(CC) $(CFLAGS) $(LFLAGS) -o $@ -c $^
+
+clean:
+	for i in `seq 5`; do \
+		rm -vf q$$i; \
+	done;
+	rm -vf *.o

modsim/ass2/q1.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define H 1e-3
+#define TABLE_LINE(func, x) (printf("%-24s%.12f\t%.12f\n", #x, \
+            slope_right(func, x, H), slope_central(func, x, H)))
+
+typedef double (*func_ptr)(double x);
+
+double slope_right(func_ptr func, double x, double h) {
+    return (func(x + h) - func(x)) / h;
+}
+
+double slope_central(func_ptr func, double x, double h) {
+    return (func(x + h) - func(x - h)) / (2 * h);
+}
+
+int main(void) {
+    puts("x\t\t\tright\t\tcentral");
+    TABLE_LINE(&sin, M_PI / 3);
+    TABLE_LINE(&sin, 100 * M_PI + M_PI / 3);
+    TABLE_LINE(&sin, 1e12 * M_PI + M_PI / 3);
+
+    return 0;
+}

modsim/ass2/q2.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define EPSILON 1e-11
+
+typedef double (*func_ptr)(double x);
+
+double bisec(func_ptr f, double left, double right, int *steps) {
+    int i;
+    double mid;
+
+    for( i = 1; fabs(right - left) > 2 * EPSILON; i++ ) {
+        mid = (right + left) / 2;
+
+        if( f(left) * f(mid) < 0 )
+            right = mid;
+        else if( f(right) * f(mid) < 0 )
+            left = mid;
+        else
+            break;
+    }
+
+    *steps = i;
+
+    return mid;
+}
+
+double func(double x) {
+    return x * sin(x) - 1;
+}
+
+int main(void) {
+    int steps;
+
+    printf("zero point: %.20f\n", bisec(&func, 0, 2, &steps));
+    printf("Steps: %d\n", steps);
+
+    return 0;
+}

modsim/ass2/q3.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define EPSILON 1e-11
+
+typedef double (*func_ptr)(double x);
+
+double bisec(func_ptr f, double left, double right, int *steps) {
+    int i;
+    double mid;
+
+    for( i = 1; fabs(right - left) > 2 * EPSILON; i++ ) {
+        mid = (right + left) / 2;
+
+        if( f(left) * f(mid) < 0 )
+            right = mid;
+        else if( f(right) * f(mid) < 0 )
+            left = mid;
+        else
+            break;
+    }
+
+    *steps = i;
+
+    return mid;
+}
+
+double func(double x) {
+    return x * sin(x) - 1;
+}
+
+int main(void) {
+    int steps;
+
+    printf("zero point: %.20f\n", bisec(&func, 0, 2, &steps));
+    printf("Steps: %d\n", steps);
+
+    return 0;
+}

os/ass1/interval.c

 /*
-   Author: G.D. van Albada
-   Date:   August 26, 2009
-   (c) Universiteit van Amsterdam
+Author: G.D. van Albada
+Date:   August 26, 2009
+(c) Universiteit van Amsterdam
 
-   In this file the data types and some of the functions declared
-   in the file interval.h for the first assignment in the OS course
-   for 2009 are defined.
+In this file the data types and some of the functions declared
+in the file interval.h for the first assignment in the OS course
+for 2009 are defined.
 */
 
 /*
@@ -26,9 +26,9 @@
 #include "interval.h"
 #include <math.h>
 
-/* 
+/*
  * Over the years the number of clock ticks per second has been
- * called by many names. The code give here appears to work on most 
+ * called by many names. The code give here appears to work on most
  * machines.
  * Use MY_CLK_TCK to convert the output of times() to seconds.
  */
@@ -53,56 +53,56 @@ interval newInterval(void)
     {
         return NULL;
     }
-    
+
     gettimeofday(&(nw->last_wct), NULL);
     times(&(nw->last_cput));
-    
+
     if( !MY_CLK_TCK )
     {
-		MY_CLK_TCK = sysconf(_SC_CLK_TCK);
+        MY_CLK_TCK = sysconf(_SC_CLK_TCK);
     }
-    
+
     return nw;
 }
 
 /*
- * Free the memory used by the interval pointer and set the pointer to NULL. If 
+ * Free the memory used by the interval pointer and set the pointer to NULL. If
  * the pointer is valid, zero is returned, or 1 otherwise.
  */
 int delInterval(interval *intervalPtr)
 {
-	if( intervalPtr == NULL )
-		return 1;
-	
-	free(*intervalPtr);
-	*intervalPtr = NULL;
-	
-	return 0;
+    if( intervalPtr == NULL )
+        return 1;
+
+    free(*intervalPtr);
+    *intervalPtr = NULL;
+
+    return 0;
 }
 
 /*
- * Returns the wall clock time, user CPU time and system CPU time for the 
+ * Returns the wall clock time, user CPU time and system CPU time for the
  * calling process and its children consumed since the previous call for the
- * specified interval. Returns zero on success, -1 when an invalid pointer is 
+ * specified interval. Returns zero on success, -1 when an invalid pointer is
  * passed as argument.
  */
 int timeInterval(interval id, double *wct, double *ust, double *syt)
 {
-	if( !id || wct == NULL || ust == NULL || syt == NULL )
-		return -1;
-	
-	struct timeval last_wct = id->last_wct;
-	struct tms last_cput = id->last_cput;
-	
-	if( delInterval(&id) )
-		return -1;
-	
-	id = newInterval();
-	
-	*wct = (id->last_wct.tv_sec + id->last_wct.tv_usec / 1e6)
-	               - (last_wct.tv_sec + last_wct.tv_usec / 1e6);
-	*ust = difftime(id->last_cput.tms_utime, last_cput.tms_utime) / MY_CLK_TCK;
-	*syt = difftime(id->last_cput.tms_stime, last_cput.tms_stime) / MY_CLK_TCK;
-	
-	return 0;
+    if( !id || wct == NULL || ust == NULL || syt == NULL )
+        return -1;
+
+    struct timeval last_wct = id->last_wct;
+    struct tms last_cput = id->last_cput;
+
+    if( delInterval(&id) )
+        return -1;
+
+    id = newInterval();
+
+    *wct = (id->last_wct.tv_sec + id->last_wct.tv_usec / 1e6)
+        - (last_wct.tv_sec + last_wct.tv_usec / 1e6);
+    *ust = difftime(id->last_cput.tms_utime, last_cput.tms_utime) / MY_CLK_TCK;
+    *syt = difftime(id->last_cput.tms_stime, last_cput.tms_stime) / MY_CLK_TCK;
+
+    return 0;
 }

os/ass1/interval.h

 /*
-   Author: G.D. van Albada
-   Date:   August 26, 2009
-   (c) Universiteit van Amsterdam
+Author: G.D. van Albada
+Date:   August 26, 2009
+(c) Universiteit van Amsterdam
 
-   In this file the data types and functions exported by the file
-   interval.c for the first assignment in the OS course for 2009
-   are defined.
+In this file the data types and functions exported by the file
+interval.c for the first assignment in the OS course for 2009
+are defined.
 */
 
 /* interval is a pointer to a struct used by the functions