aboutsummaryrefslogtreecommitdiff
path: root/tests/rational_check.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-10-29 10:05:24 +0100
committerThomas White <taw@physics.org>2019-03-11 16:49:36 +0100
commit9a3ef1de0b661085a26d9be60bcff9e261783acd (patch)
treeb1861fe26f2bbc0e9642ba8921c5e1f826a92512 /tests/rational_check.c
parent4337cafe052c4ad238c969dfa4cb7c7ac52f5e07 (diff)
Add new rational number library
Diffstat (limited to 'tests/rational_check.c')
-rw-r--r--tests/rational_check.c194
1 files changed, 194 insertions, 0 deletions
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;
+}