tools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2024 Andreas Buchheit <buchheit@num.uni-sb.de>
// SPDX-FileCopyrightText: 2024-2026 Jonathan Busse <jonathan@jbusse.de>
// SPDX-FileCopyrightText: 2024 Ruben Gutendorf
// <ruben.gutendorf@uni-saarland.de>
//
// SPDX-License-Identifier: AGPL-3.0-only
 
#ifndef EPSTEIN_TOOLS
#define EPSTEIN_TOOLS
#include <complex.h>
#include <stdbool.h>
 
static inline double real_int_pow(double base, unsigned int exp) {
    double b2;
    switch (exp) {
    case 0:
        return 1.0;
    case 1:
        return base;
    case 2:
        return base * base;
    case 3:
        return base * base * base;
    case 4:
        b2 = base * base;
        return b2 * b2;
    case 5:
        b2 = base * base;
        return b2 * b2 * base;
    case 6:
        b2 = base * base;
        return b2 * b2 * b2;
    case 7:
        b2 = base * base;
        return b2 * b2 * b2 * base;
    case 8:
        b2 = base * base;
        b2 = b2 * b2;
        return b2 * b2;
    default:
        break;
    }
    double res = 1.0;
    while (1) {
        if (exp & 1) {
            res *= base;
        }
        exp >>= 1;
        if (!exp) {
            break;
        }
        base *= base;
    }
    return res;
}
 
static inline double complex inverse_imaginary_int_pow(unsigned int exp) {
    static const double complex powers[4] = {1.0, -I, -1.0, I};
    return powers[exp & 3];
}
 
static inline double negative_one_pow(unsigned int exp) {
    if (exp & 1) {
        return -1.;
    }
    return 1.;
}
 
unsigned int mult_abs(unsigned int dim, const unsigned int *alpha);
 
unsigned long long binom(unsigned long long n, unsigned long long k);
 
static inline void kahan_add_c(double complex *restrict sum,
                               double complex *restrict epsilon, double complex x) {
    double complex y = x - *epsilon;
    double complex t = *sum + y;
    *epsilon = (t - *sum) - y;
    *sum = t;
}
 
static inline void kahan_add_r(double *restrict sum, double *restrict epsilon,
                               double x) {
    double y = x - *epsilon;
    double t = *sum + y;
    *epsilon = (t - *sum) - y;
    *sum = t;
}
 
static inline double dot(unsigned int dim, const double *v1, const double *v2) {
    double r = 0;
    for (int i = 0; i < dim; i++) {
        r += v1[i] * v2[i];
    }
    return r;
}
 
static inline void matrix_intVector(unsigned int dim, const double *m, const int *v,
                                    double *res, bool diag) {
    if (diag) {
        for (int i = 0; i < dim; i++) {
            res[i] = m[(i * dim) + i] * v[i];
        }
    } else {
        for (int i = 0; i < dim; i++) {
            res[i] = 0;
            for (int j = 0; j < dim; j++) {
                res[i] += m[(i * dim) + j] * v[j];
            }
        }
    }
}
 
void transpose(unsigned int dim, double *m);
 
bool equals(unsigned int dim, const double *v1, const double *v2);
 
void invert(unsigned int dim, double *m, int *p, double *r);
 
double inf_norm(unsigned int dim, const double *m);
 
double *vectorProj(unsigned int dim, const double *m, const double *m_invt,
                   const double *v);
#endif