aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/CMakeLists.txt2
-rw-r--r--libcrystfel/src/integer_matrix.c107
-rw-r--r--libcrystfel/src/integer_matrix.h16
-rw-r--r--libcrystfel/src/rational.c472
-rw-r--r--libcrystfel/src/rational.h94
-rw-r--r--tests/CMakeLists.txt4
-rw-r--r--tests/rational_check.c194
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;
+}