#include <algorithm>

extern "C" {
    int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
    int dscal_(int* n, double* da, double* dx, int* incx);
    int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
    void openblas_set_num_threads(int);
    int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda,
               const double* x, int* incx, double* beta, double* y, int* incy);
    double ddot_(int* n, const double* dx, int* incx, const double* dy, int* incy);
}

double shockprob(double p, double rho, double Z, bool give_log);
double dshockprob(double p, double rho, double Z);

//simple class which wraps a dynamically allocated array
//but allows it to be externally allocated
struct array {
    array(int n) : n(n), extern_alloc(false) {
        arr = new double[n];
    }
    array(int n, double* const ptr) :n(n), extern_alloc(true) {
        arr = new(ptr) double[n];
    }
    ~array() {
        if(!extern_alloc) {
            delete[] arr;
        }
    }
    double* begin() {
        return arr;
    }
    double* begin() const {
        return arr;
    }
    double* end() {
        return arr+n;
    }
    double* end() const {
        return arr+n;
    }
private:
    const int n;
    double* arr;
    const bool extern_alloc;
};

struct matrix {
    matrix(int m, int n) : m(m), n(n), mat(m*n) {};
    matrix(int m, int n, double x) : m(m), n(n), mat(m*n) {
        std::fill(mat.begin(), mat.end(), x);
    }
    matrix(int m, int n, double* ptr) : m(m), n(n), mat(m*n, ptr) {};
    double& operator()(int i, int j) {
        return *(mat.begin()+j*m+i);
    }
    double operator()(int i, int j) const {
        return *(mat.begin()+j*m+i);
    }
    double* operator()(int i) {
        return mat.begin()+i*m;
    }
    double* data() {
        return mat.begin();
    }
    double& operator[](int i) {
        return *(mat.begin()+i);
    }
    double operator[](int i) const {
        return *(mat.begin()+i);
    }
private:
    int m, n;
    array mat;
};