From 9a3ef1de0b661085a26d9be60bcff9e261783acd Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 29 Oct 2018 10:05:24 +0100 Subject: Add new rational number library --- libcrystfel/CMakeLists.txt | 2 + libcrystfel/src/integer_matrix.c | 107 ++++++++- libcrystfel/src/integer_matrix.h | 16 +- libcrystfel/src/rational.c | 472 +++++++++++++++++++++++++++++++++++++++ libcrystfel/src/rational.h | 94 ++++++++ 5 files changed, 688 insertions(+), 3 deletions(-) create mode 100644 libcrystfel/src/rational.c create mode 100644 libcrystfel/src/rational.h (limited to 'libcrystfel') 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 #include +#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; irows; i++ ) { + + unsigned int j; + + ans[i] = 0.0; + for ( j=0; jcols; 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; irows; i++ ) { + + unsigned int j; + + ans[i] = rtnl_zero(); + for ( j=0; jcols; 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 + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#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 ab */ +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; irows*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; iv[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; hrows; 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; irows; 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; icols; 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; irows; i++ ) { + + int j; + Rational dval; + + dval = rtnl_div(rtnl_mtx_get(cm, i, k), + rtnl_mtx_get(cm, h, k)); + + for ( j=0; jcols; 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; jcols; 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; irows; i++ ) { + + fprintf(stderr, "[ "); + for ( j=0; jcols; 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 + * + * 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 . + * + */ + +#ifndef RATIONAL_H +#define RATIONAL_H + +#ifdef HAVE_CONFIG_H +#include +#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 */ -- cgit v1.2.3