Add new rational number library

7 files changed, 886 insertions, 3 deletions

diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index 247b925e..bd4f44ca 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -62,6 +62,7 @@ set(LIBCRYSTFEL_SOURCES
     src/peakfinder8.c
     src/taketwo.c
     src/xgandalf.c
+    src/rational.c
 )
 
 if (HAVE_FFTW)
@@ -101,6 +102,7 @@ set(LIBCRYSTFEL_HEADERS
     src/peakfinder8.h
     src/taketwo.h
     src/xgandalf.h
+    src/rational.h
 )
 
 add_library(${PROJECT_NAME} SHARED
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index 1e616e43..6e2ade6a 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -35,7 +35,9 @@
 #include <string.h>
 #include <assert.h>
 
+#include "rational.h"
 #include "integer_matrix.h"
+#include "utils.h"
 
 
 /**
@@ -95,7 +97,7 @@ IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols)
  *
  * Returns: a newly allocated copy of @m, or NULL on error/
  **/
-IntegerMatrix *intmat_copy(IntegerMatrix *m)
+IntegerMatrix *intmat_copy(const IntegerMatrix *m)
 {
 	IntegerMatrix *p;
 	int i, j;
@@ -127,6 +129,19 @@ void intmat_free(IntegerMatrix *m)
 }
 
 
+void intmat_size(const IntegerMatrix *m, unsigned int *rows, unsigned int *cols)
+{
+	if ( m == NULL ) {
+		*rows = 0;
+		*cols = 0;
+		return;
+	}
+
+	*rows = m->rows;
+	*cols = m->cols;
+}
+
+
 /**
  * intmat_set:
  * @m: An %IntegerMatrix
@@ -162,6 +177,96 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j)
 }
 
 
+/* intmat_intvec_mult:
+ * @m: An %IntegerMatrix
+ * @vec: An array of floats
+ * @ans: An array of floats in which to store the result
+ *
+ * Multiplies the matrix @m by the vector @vec.  The size of @vec must equal the
+ * number of columns in @m, and the size of the result equals the number of rows
+ * in @m.
+ *
+ * Returns: non-zero on error
+ **/
+int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, float *ans)
+{
+	unsigned int i;
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		unsigned int j;
+
+		ans[i] = 0.0;
+		for ( j=0; j<m->cols; j++ ) {
+			ans[i] += intmat_get(m, i, j) * vec[j];
+		}
+
+	}
+
+	return 0;
+}
+
+
+/* intmat_rationalvec_mult:
+ * @m: An %IntegerMatrix
+ * @vec: An array of %Rational
+ * @ans: An array of %Rational in which to store the result
+ *
+ * Multiplies the matrix @m by the vector @vec.  The size of @vec must equal the
+ * number of columns in @m, and the size of the result equals the number of rows
+ * in @m.
+ *
+ * Returns: non-zero on error
+ **/
+int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec,
+                            Rational *ans)
+{
+	unsigned int i;
+
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		unsigned int j;
+
+		ans[i] = rtnl_zero();
+		for ( j=0; j<m->cols; j++ ) {
+			Rational t;
+			t = rtnl_mul(vec[j], rtnl(intmat_get(m, i, j), 1));
+			ans[i] = rtnl_add(ans[i], t);
+		}
+
+	}
+
+	return 0;
+}
+
+
+/* intmat_solve_rational:
+ * @m: An %IntegerMatrix
+ * @vec: An array of %Rational
+ * @ans: An array of %Rational in which to store the result
+ *
+ * Solves the matrix equation m*ans = vec, where @ans and @vec are
+ * vectors of %Rational.
+ *
+ * This is just a convenience function which creates a %RationalMatrix out of
+ * @m and then calls rtnl_mtx_solve().
+ *
+ * Returns: non-zero on error
+ **/
+int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec,
+                          Rational *ans)
+{
+	RationalMatrix *rm;
+	int r;
+
+	rm = rtnl_mtx_from_intmat(m);
+	r = rtnl_mtx_solve(rm, vec, ans);
+	rtnl_mtx_free(rm);
+	return r;
+}
+
+
 /**
  * intmat_intvec_mult:
  * @m: An %IntegerMatrix
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index c4983aef..aeea38dd 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -40,17 +40,22 @@
  **/
 typedef struct _integermatrix IntegerMatrix;
 
+#include "rational.h"
+
 #ifdef __cplusplus
 extern "C" {
 #endif
 
 /* Alloc/dealloc */
 extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols);
-extern IntegerMatrix *intmat_copy(IntegerMatrix *m);
+extern IntegerMatrix *intmat_copy(const IntegerMatrix *m);
 extern IntegerMatrix *intmat_identity(int size);
 extern void intmat_free(IntegerMatrix *m);
 
 /* Get/set */
+extern void intmat_size(const IntegerMatrix *m, unsigned int *rows,
+                        unsigned int *cols);
+
 extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j,
                        signed int v);
 extern signed int intmat_get(const IntegerMatrix *m,
@@ -65,14 +70,21 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i
                                         signed int m21, signed int m22, signed int m23,
                                         signed int m31, signed int m32, signed int m33);
 
-/* Matrix-(int)vector multiplication */
+/* Matrix-vector multiplication */
+extern int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec,
+                                float *ans);
 extern signed int *intmat_intvec_mult(const IntegerMatrix *m,
                                       const signed int *vec);
+extern int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec,
+                                   Rational *ans);
 
 /* Matrix-matrix multiplication */
 extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
                                          const IntegerMatrix *b);
 
+extern int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec,
+                                 Rational *ans);
+
 /* Inverse */
 extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m);
 
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
new file mode 100644
index 00000000..bc783711
--- /dev/null
+++ b/libcrystfel/src/rational.c
@@ -0,0 +1,472 @@
+/*
+ * rational.c
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <locale.h>
+
+#include "rational.h"
+#include "integer_matrix.h"
+#include "utils.h"
+
+
+/**
+ * SECTION:rational
+ * @short_description: Rational numbers
+ * @title: Rational numbers
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational number library
+ */
+
+
+/* Eucliden algorithm for finding greatest common divisor */
+static signed int gcd(signed int a, signed int b)
+{
+	while ( b != 0 ) {
+		signed int t = b;
+		b = a % b;
+		a = t;
+	}
+	return a;
+}
+
+
+static void squish(Rational *rt)
+{
+	signed int g;
+
+	if ( rt->num == 0 ) {
+		rt->den = 1;
+		return;
+	}
+
+	g = gcd(rt->num, rt->den);
+	assert(g != 0);
+	rt->num /= g;
+	rt->den /= g;
+
+	if ( rt->den < 0 ) {
+		rt->num = -rt->num;
+		rt->den = -rt->den;
+	}
+}
+
+
+Rational rtnl_zero()
+{
+	Rational r;
+	r.num = 0;
+	r.den = 1;
+	return r;
+}
+
+
+Rational rtnl(signed long long int num, signed long long int den)
+{
+	Rational r;
+	r.num = num;
+	r.den = den;
+	squish(&r);
+	return r;
+}
+
+
+double rtnl_as_double(Rational r)
+{
+	return (double)r.num/r.den;
+}
+
+
+static void overflow(long long int c, long long int a, long long int b)
+{
+	setlocale(LC_ALL, "");
+	ERROR("Overflow detected in rational number library.\n");
+	ERROR("%'lli < %'lli * %'lli\n", c, a, b);
+	abort();
+}
+
+
+static void check_overflow(long long int c, long long int a, long long int b)
+{
+	if ( (a==0) || (b==0) ) {
+		if ( c != 0 ) overflow(c,a,b);
+	} else if ( (llabs(c) < llabs(a)) || (llabs(c) < llabs(b)) ) {
+		overflow(c,a,b);
+	}
+}
+
+
+Rational rtnl_mul(Rational a, Rational b)
+{
+	Rational r;
+	r.num = a.num * b.num;
+	r.den = a.den * b.den;
+	check_overflow(r.num, a.num, b.num);
+	check_overflow(r.den, a.den, b.den);
+	squish(&r);
+	return r;
+}
+
+
+Rational rtnl_div(Rational a, Rational b)
+{
+	signed int t = b.num;
+	b.num = b.den;
+	b.den = t;
+	return rtnl_mul(a, b);
+}
+
+
+Rational rtnl_add(Rational a, Rational b)
+{
+	Rational r, trt1, trt2;
+
+	trt1.num = a.num * b.den;
+	trt2.num = b.num * a.den;
+	check_overflow(trt1.num, a.num, b.den);
+	check_overflow(trt2.num, b.num, a.den);
+
+	trt1.den = a.den * b.den;
+	trt2.den = trt1.den;
+	check_overflow(trt1.den, a.den, b.den);
+
+	r.num = trt1.num + trt2.num;
+	r.den = trt1.den;
+	squish(&r);
+	return r;
+}
+
+
+Rational rtnl_sub(Rational a, Rational b)
+{
+	b.num = -b.num;
+	return rtnl_add(a, b);
+}
+
+
+/* -1, 0 +1 respectively for a<b, a==b, a>b */
+signed int rtnl_cmp(Rational a, Rational b)
+{
+	Rational trt1, trt2;
+
+	trt1.num = a.num * b.den;
+	trt2.num = b.num * a.den;
+
+	trt1.den = a.den * b.den;
+	trt2.den = a.den * b.den;
+
+	if ( trt1.num > trt2.num ) return +1;
+	if ( trt1.num < trt2.num ) return -1;
+	return 0;
+}
+
+
+Rational rtnl_abs(Rational a)
+{
+	Rational r = a;
+	squish(&r);
+	if ( r.num < 0 ) r.num = -r.num;
+	return r;
+}
+
+
+/**
+ * rtnl_format
+ * @rt: A %Rational
+ *
+ * Formats @rt as a string
+ *
+ * Returns: a string which should be freed by the caller
+ */
+char *rtnl_format(Rational rt)
+{
+	char *v = malloc(32);
+	if ( v == NULL ) return NULL;
+	if ( rt.den == 1 ) {
+		snprintf(v, 31, "%lli", rt.num);
+	} else {
+		snprintf(v, 31, "%lli/%lli", rt.num, rt.den);
+	}
+	return v;
+}
+
+
+/**
+ * SECTION:rational_matrix
+ * @short_description: Rational matrices
+ * @title: Rational matrices
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational matrix library
+ */
+
+
+struct _rationalmatrix
+{
+	unsigned int rows;
+	unsigned int cols;
+	Rational *v;
+};
+
+
+/**
+ * rtnl_mtx_new:
+ * @rows: Number of rows that the new matrix is to have
+ * @cols: Number of columns that the new matrix is to have
+ *
+ * Allocates a new %RationalMatrix with all elements set to zero.
+ *
+ * Returns: a new %RationalMatrix, or NULL on error.
+ **/
+RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols)
+{
+	RationalMatrix *m;
+
+	m = malloc(sizeof(RationalMatrix));
+	if ( m == NULL ) return NULL;
+
+	m->v = calloc(rows*cols, sizeof(Rational));
+	if ( m->v == NULL ) {
+		free(m);
+		return NULL;
+	}
+
+	m->rows = rows;
+	m->cols = cols;
+
+	return m;
+}
+
+
+RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m)
+{
+	RationalMatrix *n;
+	int i;
+
+	n = rtnl_mtx_new(m->rows, m->cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<m->rows*m->cols; i++ ) {
+		n->v[i] = m->v[i];
+	}
+
+	return n;
+}
+
+
+Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j)
+{
+	assert(m != NULL);
+	return m->v[j+m->cols*i];
+}
+
+
+void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v)
+{
+	assert(m != NULL);
+	m->v[j+m->cols*i] = v;
+}
+
+
+RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m)
+{
+	RationalMatrix *n;
+	unsigned int rows, cols;
+	int i, j;
+
+	intmat_size(m, &rows, &cols);
+	n = rtnl_mtx_new(rows, cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<rows; i++ ) {
+		for ( j=0; j<cols; j++ ) {
+			n->v[j+cols*i] = rtnl(intmat_get(m, i, j), 1);
+		}
+	}
+
+	return n;
+}
+
+
+void rtnl_mtx_free(RationalMatrix *mtx)
+{
+	if ( mtx == NULL ) return;
+	free(mtx->v);
+	free(mtx);
+}
+
+
+/* rtnl_mtx_solve:
+ * @m: A %RationalMatrix
+ * @vec: An array of %Rational
+ * @ans: An array of %Rational in which to store the result
+ *
+ * Solves the matrix equation m*ans = vec, where @ans and @vec are
+ * vectors of %Rational.
+ *
+ * In this version, @m must be square.
+ *
+ * The number of columns in @m must equal the length of @ans, and the number
+ * of rows in @m must equal the length of @vec, but note that there is no way
+ * for this function to check that this is the case.
+ *
+ * Returns: non-zero on error
+ **/
+int rtnl_mtx_solve(const RationalMatrix *m, const Rational *ivec, Rational *ans)
+{
+	RationalMatrix *cm;
+	Rational *vec;
+	int h, k;
+	int i;
+
+	if ( m->rows != m->cols ) return 1;
+
+	/* Copy the matrix and vector because the calculation will
+	 * be done in-place */
+	cm = rtnl_mtx_copy(m);
+	if ( cm == NULL ) return 1;
+
+	vec = malloc(m->rows*sizeof(Rational));
+	if ( vec == NULL ) return 1;
+	for ( h=0; h<m->rows; h++ ) vec[h] = ivec[h];
+
+	/* Gaussian elimination with partial pivoting */
+	h = 0;
+	k = 0;
+	while ( h<=m->rows && k<=m->cols ) {
+
+		int prow = 0;
+		Rational pval = rtnl_zero();
+		Rational t;
+
+		/* Find the row with the largest value in column k */
+		for ( i=h; i<m->rows; i++ ) {
+			if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) {
+				pval = rtnl_abs(rtnl_mtx_get(cm, i, k));
+				prow = i;
+			}
+		}
+
+		if ( rtnl_cmp(pval, rtnl_zero()) == 0 ) {
+			k++;
+			continue;
+		}
+
+		/* Swap 'prow' with row h */
+		for ( i=0; i<m->cols; i++ ) {
+			t = rtnl_mtx_get(cm, h, i);
+			rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i));
+			rtnl_mtx_set(cm, prow, i, t);
+		}
+		t = vec[prow];
+		vec[prow] = vec[h];
+		vec[h] = t;
+
+		/* Divide and subtract rows below */
+		for ( i=h+1; i<m->rows; i++ ) {
+
+			int j;
+			Rational dval;
+
+			dval = rtnl_div(rtnl_mtx_get(cm, i, k),
+			                rtnl_mtx_get(cm, h, k));
+
+			for ( j=0; j<m->cols; j++ ) {
+				Rational t = rtnl_mtx_get(cm, i, j);
+				Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j));
+				t = rtnl_sub(t, p);
+				rtnl_mtx_set(cm, i, j, t);
+			}
+
+			/* Divide the right hand side as well */
+			Rational t = vec[i];
+			Rational p = rtnl_mul(dval, vec[h]);
+			vec[i] = rtnl_sub(t, p);
+		}
+
+		h++;
+		k++;
+
+
+	}
+
+	/* Back-substitution */
+	for ( i=m->rows-1; i>=0; i-- ) {
+		int j;
+		Rational sum = rtnl_zero();
+		for ( j=i+1; j<m->cols; j++ ) {
+			Rational av;
+			av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]);
+			sum = rtnl_add(sum, av);
+		}
+		sum = rtnl_sub(vec[i], sum);
+		ans[i] = rtnl_div(sum, rtnl_mtx_get(cm, i, i));
+	}
+
+	free(vec);
+	rtnl_mtx_free(cm);
+
+	return 0;
+}
+
+
+/**
+ * rtnl_mtx_print
+ * @m: A %RationalMatrix
+ *
+ * Prints @m to stderr.
+ *
+ */
+void rtnl_mtx_print(const RationalMatrix *m)
+{
+	unsigned int i, j;
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		fprintf(stderr, "[ ");
+		for ( j=0; j<m->cols; j++ ) {
+			char *v = rtnl_format(m->v[j+m->cols*i]);
+			fprintf(stderr, "%4s ", v);
+			free(v);
+		}
+		fprintf(stderr, "]\n");
+	}
+}
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
new file mode 100644
index 00000000..8b182163
--- /dev/null
+++ b/libcrystfel/src/rational.h
@@ -0,0 +1,94 @@
+/*
+ * rational.h
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef RATIONAL_H
+#define RATIONAL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * Rational
+ *
+ * The Rational is an opaque-ish data structure representing a rational number.
+ *
+ * "Opaque-ish" means that the structure isn't technically opaque, allowing you
+ * to assign and allocate them easily.  But you shouldn't look at or set its
+ * contents, except by using the accessor functions.
+ **/
+typedef struct {
+	/* Private, don't modify */
+	signed long long int num;
+	signed long long int den;
+} Rational;
+
+
+/**
+ * RationalMatrix
+ *
+ * The RationalMatrix is an opaque data structure representing a matrix of
+ * rational numbers.
+ **/
+typedef struct _rationalmatrix RationalMatrix;
+
+#include "integer_matrix.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern Rational rtnl_zero(void);
+extern Rational rtnl(signed long long int num, signed long long int den);
+extern double rtnl_as_double(Rational r);
+
+extern Rational rtnl_mul(Rational a, Rational b);
+extern Rational rtnl_div(Rational a, Rational b);
+extern Rational rtnl_add(Rational a, Rational b);
+extern Rational rtnl_sub(Rational a, Rational b);
+
+extern signed int rtnl_cmp(Rational a, Rational b);
+extern Rational rtnl_abs(Rational a);
+
+extern char *rtnl_format(Rational rt);
+
+extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols);
+extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m);
+extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j);
+extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v);
+extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m);
+extern void rtnl_mtx_free(RationalMatrix *mtx);
+extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec,
+                          Rational *ans);
+extern void rtnl_mtx_print(const RationalMatrix *m);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif	/* RATIONAL_H */
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 02f3ca02..beb141c3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -83,3 +83,7 @@ if (HAVE_OPENCL)
   add_test(gpu_sim_check gpu_sim_check)
 endif (HAVE_OPENCL)
 
+add_executable(rational_check rational_check.c)
+target_include_directories(rational_check PRIVATE ${COMMON_INCLUDES})
+target_link_libraries(rational_check ${COMMON_LIBRARIES})
+add_test(rational_check rational_check)
diff --git a/tests/rational_check.c b/tests/rational_check.c
new file mode 100644
index 00000000..8957cfdd
--- /dev/null
+++ b/tests/rational_check.c
@@ -0,0 +1,194 @@
+/*
+ * rational_check.c
+ *
+ * Check that rational numbers work
+ *
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2019 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <gsl/gsl_rng.h>
+
+#include <rational.h>
+#include <utils.h>
+
+
+static Rational gen_rtnl(gsl_rng *rng, double *dptr, int sz)
+{
+	Rational r;
+	signed int num, den;
+	do {
+		num = gsl_rng_uniform_int(rng, 2*sz) - sz;
+		den = gsl_rng_uniform_int(rng, 2*sz) - sz;
+	} while ( den == 0 );
+	r = rtnl(num, den);
+	*dptr = (double)num/den;
+	return r;
+}
+
+
+static int test_rational(gsl_rng *rng)
+{
+	Rational r1, r2, r3;
+	double rd1, rd2, rd3;
+	char op;
+
+	r1 = gen_rtnl(rng, &rd1, 100);
+	r2 = gen_rtnl(rng, &rd2, 100);
+
+	switch ( gsl_rng_uniform_int(rng, 3) ) {
+
+		case 0:
+		r3 = rtnl_mul(r1, r2);
+		rd3 = rd1*rd2;
+		op = '*';
+		break;
+
+		case 1:
+		r3 = rtnl_div(r1, r2);
+		rd3 = rd1/rd2;
+		op = '/';
+		break;
+
+		case 2:
+		r3 = rtnl_add(r1, r2);
+		rd3 = rd1+rd2;
+		op = '+';
+		break;
+
+		case 3:
+		r3 = rtnl_sub(r1, r2);
+		rd3 = rd1-rd2;
+		op = '-';
+		break;
+
+		default:
+		abort();
+
+	}
+
+	if ( isinf(rd3) && isinf(rtnl_as_double(r3)) ) return 0;
+
+	if ( fabs(rtnl_as_double(r3) - rd3) > 0.001 ) {
+		ERROR("%10s %c %10s = %10s", rtnl_format(r1), op,
+		      rtnl_format(r2), rtnl_format(r3));
+		ERROR(" --->  %10f %c %10f = %10f\n", rd1, op, rd2, rd3);
+		return 1;
+	}
+
+	return 0;
+}
+
+
+static int test_rational_matrix(gsl_rng *rng)
+{
+	const int size = 3;
+	RationalMatrix *rm;
+	Rational *rvec;
+	Rational *rans;
+	gsl_matrix *gm;
+	gsl_vector *gvec;
+	gsl_vector *gans;
+	int i, j;
+	int filt;
+	int fail = 0;
+
+	rm = rtnl_mtx_new(size, size);
+	gm = gsl_matrix_alloc(size, size);
+	rvec = malloc(size*sizeof(Rational));
+	gvec = gsl_vector_alloc(size);
+	rans = malloc(size*sizeof(Rational));
+	if ( (rm==NULL) || (gm==NULL) || (rvec==NULL)
+	  || (gvec==NULL) || (rans==NULL) ) return 1;
+
+	/* Find a matrix equation which actually works */
+	do {
+		for ( i=0; i<size; i++ ) {
+			double d;
+			Rational r;
+			for ( j=0; j<size; j++ ) {
+				r = gen_rtnl(rng, &d, 5);
+				rtnl_mtx_set(rm, i, j, r);
+				gsl_matrix_set(gm ,i, j, d);
+			}
+			r = gen_rtnl(rng, &d, 5);
+			rvec[i] = r;
+			gsl_vector_set(gvec ,i, d);
+		}
+
+		gans = solve_svd(gvec, gm, &filt, 0);
+	} while ( (gans==NULL) || (filt!=0) || (isnan(gsl_vector_get(gans,0))) );
+
+	rtnl_mtx_solve(rm, rvec, rans);
+
+	for ( i=0; i<size; i++ ) {
+		if ( fabs(rtnl_as_double(rans[i]) - gsl_vector_get(gans, i) > 0.001) )
+		{
+			STATUS("%25s %10f %10f\n", rtnl_format(rans[i]),
+			                           rtnl_as_double(rans[i]),
+			                           gsl_vector_get(gans, i));
+			fail = 1;
+		}
+	}
+
+	gsl_matrix_free(gm);
+	rtnl_mtx_free(rm);
+	free(rans);
+	free(rvec);
+	gsl_vector_free(gvec);
+	gsl_vector_free(gans);
+	return fail;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int fail = 0;
+	gsl_rng *rng;
+	int i;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+	gsl_set_error_handler_off();
+
+	STATUS("Arithmetic test...\n");
+	for ( i=0; i<1000; i++ ) {
+		fail += test_rational(rng);
+	}
+
+	STATUS("Matrix test...\n");
+	for ( i=0; i<1000; i++ ) {
+		fail += test_rational_matrix(rng);
+	}
+
+	gsl_rng_free(rng);
+
+	if ( fail != 0 ) return 1;
+	return 0;
+}