Commit db6a27c2 authored by Taddeüs Kroes's avatar Taddeüs Kroes

Added Kahan summation algorithm to question 3.

parent d205c1c6
...@@ -9,3 +9,4 @@ report.pdf ...@@ -9,3 +9,4 @@ report.pdf
fd* fd*
pr pr
floating_point.tex floating_point.tex
kahan
...@@ -4,7 +4,7 @@ SPEED_TYPES=float double LD ...@@ -4,7 +4,7 @@ SPEED_TYPES=float double LD
SUM_TYPES=float double SUM_TYPES=float double
OPS=ADD DIV MULT SQRT OPS=ADD DIV MULT SQRT
all: fp speed highlight report.pdf sum pr all: fp speed highlight report.pdf sum kahan pr
highlight: floating_point.tex highlight: floating_point.tex
...@@ -34,11 +34,14 @@ fp: floating_point.o ...@@ -34,11 +34,14 @@ fp: floating_point.o
sum: sum.c sum: sum.c
for t in $(SUM_TYPES); do \ for t in $(SUM_TYPES); do \
sed "s#{TYPE}#$$t#" $^ > sum.$$t.c; \ sed "s#{TYPE}#$$t#" $^ > sum.$$t.c; \
gcc $(FLAGS) -o sum.$$t sum.$$t.c; \ $(CC) $(FLAGS) -o sum.$$t sum.$$t.c; \
rm sum.$$t.c; \ rm sum.$$t.c; \
done; done;
touch $@ touch $@
kahan: kahan_sum.o
$(CC) $(FLAGS) -o $@ $^
%.o: %.c %.o: %.c
$(CC) $(FLAGS) -o $@ -c $^ $(CC) $(FLAGS) -o $@ -c $^
...@@ -47,4 +50,4 @@ sum: sum.c ...@@ -47,4 +50,4 @@ sum: sum.c
clean: clean:
rm -vf *.o *.i *.s fp pr fd* speed speed.*.* floating_point \ rm -vf *.o *.i *.s fp pr fd* speed speed.*.* floating_point \
report.pdf *.aux *.log *.toc sum sum.float sum.double report.pdf *.aux *.log *.toc sum sum.float sum.double kahan
#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;
}
...@@ -112,13 +112,13 @@ We've calculated $\sum_{i=1}^{N}\frac{1}{i}$ for $N = 10^8$ and $N = 2 \cdot ...@@ -112,13 +112,13 @@ We've calculated $\sum_{i=1}^{N}\frac{1}{i}$ for $N = 10^8$ and $N = 2 \cdot
\texttt{float} and \texttt{double}. The results of this are in the table below. \texttt{float} and \texttt{double}. The results of this are in the table below.
\begin{table}[H] \begin{table}[H]
\begin{tabular}{l|lll} \begin{tabular}{l|llll}
Type & N & Forward & Backward \\ Type & N & Forward & Backward & Kahan\\
\hline \hline
\texttt{float} & $10^8$ & $15.40368$ & $18.80792$ \\ \texttt{float} & $10^8$ & $15.40368$ & $18.80792$ & $18.99790$ \\
\texttt{float} & $2 \cdot 10^8$ & $15.40368$ & $18.80792$ \\ \texttt{float} & $2 \cdot 10^8$ & $15.40368$ & $18.80792$ & $19.69104$ \\
\texttt{double} & $10^8$ & $18.99790$ & $18.99790$ \\ \texttt{double} & $10^8$ & $18.99790$ & $18.99790$ & \\
\texttt{double} & $2 \cdot 10^8$ & $19.69104$ & $19.69104$ \\ \texttt{double} & $2 \cdot 10^8$ & $19.69104$ & $19.69104$ & \\
\end{tabular} \end{tabular}
\caption{Results of various summation approaches on floats and doubles.} \caption{Results of various summation approaches on floats and doubles.}
\end{table} \end{table}
...@@ -149,11 +149,13 @@ Type & N & Forward & Backward \\ ...@@ -149,11 +149,13 @@ Type & N & Forward & Backward \\
to the same problem as described above: all numbers in $[\frac{1}{10^8}, 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 \frac{1}{2 \cdot 10^8}]$ are also represented as zero and therefore not added
to the result. to the result.
\item To improve the precision of the \texttt{float} data type summation, we
implemented the Kahan summation algorithm.
\end{itemize} \end{itemize}
% }}} % }}}
\appendix{} \appendix
\section{floating\_point.c} % {{{ \section{floating\_point.c} % {{{
\label{sec:floating_point.c} \label{sec:floating_point.c}
......
...@@ -23,12 +23,12 @@ int main(void) { ...@@ -23,12 +23,12 @@ int main(void) {
puts("Using type {TYPE}."); puts("Using type {TYPE}.");
puts("Forward summation:"); puts("Forward summation:");
printf("N = 1e8: %e\n", sum_forward(1e8)); printf("N = 1e8: %f\n", sum_forward(1e8));
printf("N = 2e8: %e\n", sum_forward(2e8)); printf("N = 2e8: %f\n", sum_forward(2e8));
puts("Backward summation:"); puts("Backward summation:");
printf("N = 1e8: %e\n", sum_backward(1e8)); printf("N = 1e8: %f\n", sum_backward(1e8));
printf("N = 2e8: %e\n", sum_backward(2e8)); printf("N = 2e8: %f\n", sum_backward(2e8));
return 0; return 0;
} }
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment