diff --git a/include/numerics/detail/dot_serial.h b/include/numerics/detail/dot_serial.h new file mode 100644 index 0000000..6e3a806 --- /dev/null +++ b/include/numerics/detail/dot_serial.h @@ -0,0 +1,42 @@ +#pragma once + +#include //uint64_t +#include // std::runtime_error + +#include "../utils/vector.h" +#include "../utils/matrix.h" + +namespace numerics::detail{ + + // ---------------- Vector * Vector ---------------- + template + inline T dot_serial(const utils::Vector& v, const utils::Vector& p) { + const uint64_t N = v.size(); + if (N != p.size()) { + throw std::runtime_error("dot_serial: dimension mismatch"); + } + T dot_product = T{0}; + for (uint64_t i = 0; i < N; ++i){ + dot_product += v[i]*p[i]; + } + return dot_product; + } + + // ---------------- Matrix * Vector ---------------- + template + inline utils::Vector dot_serial(const utils::Matrix& A, const utils::Vector& v) { + const uint64_t rows = A.rows(); + const uint64_t cols = A.cols(); + if (cols != v.size()) { + throw std::runtime_error("dot_serial: dimension mismatch"); + } + utils::Vector dot_product(rows, T{0}); + for (uint64_t j = 0; j < cols; ++j){ + for (uint64_t i = 0; i < rows; ++i){ + dot_product[i] += A(i,j)*v[j]; + } + } + return dot_product; + } +} // namespace numerics + diff --git a/include/numerics/detail/equal_serial.h b/include/numerics/detail/equal_serial.h new file mode 100644 index 0000000..197e46b --- /dev/null +++ b/include/numerics/detail/equal_serial.h @@ -0,0 +1,47 @@ +#pragma once + +#include //uint64_t +//#include // std::runtime_error + +#include "../utils/vector.h" +#include "../utils/matrix.h" + +namespace numerics::detail{ + + // ---------------- Matrix ---------------- + template + inline bool equal_serial(const utils::Matrix& A, const utils::Matrix & B) { + const uint64_t rows = A.rows(); + const uint64_t cols = A.cols(); + + if ((rows != B.rows()) || (cols != B.cols())){ + return false; + } + + for (uint64_t i = 0; i < rows; ++i){ + for (uint64_t j = 0; j < cols; ++j){ + if (A(i,j) != B(i,j)){ + return false; + } + } + } + return true; + } + + // ---------------- Vector ---------------- + template + inline bool equal_serial(const utils::Vector& v, const utils::Vector& p) { + const uint64_t N = v.size(); + if (N != p.size()){ + return false; + } + for (uint64_t i = 0; i < N; ++i){ + if ((v[i] != p[i])){ + return false; + } + } + return true; + } + +} // namespace numerics + diff --git a/include/numerics/detail/isclose_serial.h b/include/numerics/detail/isclose_serial.h new file mode 100644 index 0000000..6dac405 --- /dev/null +++ b/include/numerics/detail/isclose_serial.h @@ -0,0 +1,49 @@ +#pragma once + +#include //uint64_t +//#include // std::runtime_error + +#include "../utils/vector.h" +#include "../utils/matrix.h" +#include // std::abs + + +namespace numerics::detail{ + + // ---------------- Matrix ---------------- + template + inline bool isclose_serial(const utils::Matrix& A, const utils::Matrix & B, const T tol = T{1e-5}) { + const uint64_t rows = A.rows(); + const uint64_t cols = A.cols(); + + if ((rows != B.rows()) || (cols != B.cols())){ + return false; + } + + for (uint64_t i = 0; i < rows; ++i){ + for (uint64_t j = 0; j < cols; ++j){ + if (static_cast(std::abs(A(i,j)-B(i,j))) > tol){ + return false; + } + } + } + return true; + } + + // ---------------- Vector ---------------- + template + inline bool isclose_serial(const utils::Vector& v, const utils::Vector& p, const T tol = T{1e-5}) { + const uint64_t N = v.size(); + if (N != p.size()){ + return false; + } + for (uint64_t i = 0; i < N; ++i){ + if (static_cast(std::abs((v[i]-p[i]))) > tol){ + return false; + } + } + return true; + } + +} // namespace numerics + diff --git a/include/numerics/detail/matmul_serial.h b/include/numerics/detail/matmul_serial.h new file mode 100644 index 0000000..86b9033 --- /dev/null +++ b/include/numerics/detail/matmul_serial.h @@ -0,0 +1,35 @@ +#pragma once + +#include //uint64_t +#include // std::runtime_error + +#include "../utils/vector.h" +#include "../utils/matrix.h" + +namespace numerics::detail{ + +// ---------------- Matrix * Matrix ---------------- + template + inline utils::Matrix matmul_serial(const utils::Matrix& A, const utils::Matrix& B){ + const uint64_t m = A.rows(); + const uint64_t n = A.cols(); // also B.rows() + const uint64_t p = B.cols(); + if(n != B.rows()){ + throw std::runtime_error("matmul: dimension mismatch"); + } + T tmp; + utils::Matrix C(m, p, T{0}); + for (uint64_t i = 0; i < m; ++i){ + for (uint64_t j = 0; j < n; ++j){ + tmp = A(i,j); + for (uint64_t k = 0; k < p; ++k){ + C(i,k) += tmp * B(j,k); + } + } + } + return C; + } + + +} // namespace numerics + diff --git a/include/numerics/detail/random_serial.h b/include/numerics/detail/random_serial.h new file mode 100644 index 0000000..61be467 --- /dev/null +++ b/include/numerics/detail/random_serial.h @@ -0,0 +1,63 @@ +#pragma once + +#include //uint64_t +//#include // std::runtime_error +#include +#include + + +#include "../utils/vector.h" +#include "../utils/matrix.h" + +namespace numerics::detail{ + + // Shared engine + inline std::mt19937& rng() { + static std::random_device rd; + static std::mt19937 gen(rd()); + return gen; + } + + // Integral overload + template < + typename T, + typename std::enable_if::value, int>::type = 0 + > + inline T random(const T low, const T high) { + std::uniform_int_distribution dist(low, high); + return dist(rng()); + } + + // Floating-point overload + template < + typename T, + typename std::enable_if::value, int>::type = 0 + > + inline T random(const T low, const T high) { + std::uniform_real_distribution dist(low, high); + return dist(rng()); + } + + // ---------------- Matrix ---------------- + template + inline void inplace_random_serial(utils::Matrix& A, const T low, const T high) { + const uint64_t rows = A.rows(); + const uint64_t cols = A.cols(); + for (uint64_t i = 0; i < rows; ++i){ + for (uint64_t j = 0; j < cols; ++j){ + A(i,j) = random(low, high); + } + } + } + + // ---------------- Vector ---------------- + template + inline void inplace_random_serial(utils::Vector& v, const T low, const T high) { + const uint64_t N = v.size(); + for (uint64_t i = 0; i < N; ++i){ + v[i] = random(low, high); + } + } + +} // namespace numerics + diff --git a/include/numerics/dot.h b/include/numerics/dot.h new file mode 100644 index 0000000..3cee66a --- /dev/null +++ b/include/numerics/dot.h @@ -0,0 +1,21 @@ +#pragma once + +#include "./core/omp_config.h" +#include "detail/dot_serial.h" + + +namespace numerics{ + + // ---------------- Vector * vector ---------------- + template + inline T dot(const utils::Vector& v, const utils::Vector& p) { + return detail::dot_serial(v,p); + } + + // ---------------- Matrix * vector ---------------- + template + inline utils::Vector dot(const utils::Matrix& A, const utils::Vector& v) { + return detail::dot_serial(A,v); + } + +} \ No newline at end of file diff --git a/include/numerics/equal.h b/include/numerics/equal.h new file mode 100644 index 0000000..c132844 --- /dev/null +++ b/include/numerics/equal.h @@ -0,0 +1,22 @@ +#pragma once + +#include "./core/omp_config.h" +#include "detail/equal_serial.h" + + + +namespace numerics{ + + // ---------------- equal ---------------- + template + inline bool equal(const utils::Vector& v, const utils::Vector& p) { + return detail::equal_serial(v, p); + } + + template + inline bool equal(const utils::Matrix& A, const utils::Matrix& B) { + return detail::equal_serial(A, B); + } + + +} \ No newline at end of file diff --git a/include/numerics/isclose.h b/include/numerics/isclose.h new file mode 100644 index 0000000..e69de29 diff --git a/include/numerics/matdot.h b/include/numerics/matdot.h deleted file mode 100644 index 5de22b6..0000000 --- a/include/numerics/matdot.h +++ /dev/null @@ -1,63 +0,0 @@ -#pragma once - -#include "./core/omp_config.h" - -#include "./utils/matrix.h" -#include "./utils/vector.h" - - -namespace numerics{ - - - template - utils::Vector matdot_row(const utils::Matrix& A, const utils::Matrix& B){ - - if (A.rows() != B.rows() || A.cols() != B.cols()){ - throw std::runtime_error("matmul: dimension mismatch"); - } - - const uint64_t rows = A.rows(); - const uint64_t cols = A.cols(); - - - utils::Vector c(rows, T{0}); - - for (uint64_t i = 0; i < rows; ++i){ - T sum = T{0}; - for (uint64_t j = 0; j < cols; ++j){ - sum += A(i,j) * A(i,j); - } - c[i] = sum; - } - return c; - } - - template - utils::Vector matdot_col(const utils::Matrix& A, const utils::Matrix& B){ - - if (A.rows() != B.rows() || A.cols() != B.cols()){ - throw std::runtime_error("matmul: dimension mismatch"); - } - - const uint64_t rows = A.rows(); - const uint64_t cols = A.cols(); - - - utils::Vector c(cols, T{0}); - - for (uint64_t j = 0; j < cols; ++j){ - T sum = T{0}; - for (uint64_t i = 0; i < rows; ++i){ - sum += A(i,j) * A(i,j); - } - c[j] = sum; - } - return c; - } - - - - - -} // namespace numerics - diff --git a/include/numerics/matequal.h b/include/numerics/matequal.h deleted file mode 100644 index 8e1aaf7..0000000 --- a/include/numerics/matequal.h +++ /dev/null @@ -1,85 +0,0 @@ -#pragma once - -#include "./core/omp_config.h" - -#include "./utils/matrix.h" -#include "./numerics/abs.h" - - -namespace numerics{ - - // -------------- Serial ---------------- - template - bool matequal(const utils::Matrix& A, const utils::Matrix& B, double tol = 1e-9) { - - if (A.rows() != B.rows() || A.cols() != B.cols()) { - return false; - } - - bool decimal = std::is_floating_point::value; - const uint64_t rows=A.rows(), cols=A.cols(); - - for (uint64_t i = 0; i < rows; ++i) - for (uint64_t j = 0; j < cols; ++j) - if (decimal) { - if (numerics::abs(A(i,j) - B(i,j)) > static_cast(tol)){ - return false; - } - } else { - if (A(i,j) != B(i,j)){ - return false; - } - } - return true; - } - - // -------------- Parallel ---------------- - template - bool matequal_omp(const utils::Matrix& A, const utils::Matrix& B, double tol = 1e-9) { - - if (A.rows() != B.rows() || A.cols() != B.cols()) { - return false; - } - - bool decimal = std::is_floating_point::value; - bool eq = true; - const uint64_t rows=A.rows(), cols=A.cols(); - - #pragma omp parallel for collapse(2) schedule(static) reduction(&&:eq) - for (uint64_t i = 0; i < rows; ++i) - for (uint64_t j = 0; j < cols; ++j) - if (decimal) { - eq = eq && (numerics::abs(A(i,j) - B(i,j)) <= static_cast(tol)); - } else { - eq = eq && (A(i,j) == B(i,j)); - } - return eq; - } - - // -------------- Auto OpenMP ---------------- - template - bool matequal_auto(const utils::Matrix& A, const utils::Matrix& B, double tol = 1e-9) { - - if (A.rows() != B.rows() || A.cols() != B.cols()) { - return false; - } - - uint64_t work = A.rows() * A.cols(); - - #ifdef _OPENMP - bool can_parallel = omp_config::omp_parallel_allowed(); - uint64_t threads = static_cast(omp_get_max_threads()); - #else - bool can_parallel = false; - uint64_t threads = 1; - #endif - - if (can_parallel || work > threads * 4ull) { - return matequal_omp(A,B,tol); - } - else{ - // Safe fallback - return matequal(A,B,tol); - } - } -} // namespace numerics diff --git a/include/numerics/matmul.h b/include/numerics/matmul.h index c2a25a7..52ee2c4 100644 --- a/include/numerics/matmul.h +++ b/include/numerics/matmul.h @@ -1,124 +1,15 @@ -#ifndef _matmul_n_ -#define _matmul_n_ +#pragma once - -#include "./utils/matrix.h" #include "./core/omp_config.h" +#include "detail/matmul_serial.h" namespace numerics{ -// ---------------- Serial baseline ---------------- +// ---------------- matmul ---------------- template - utils::Matrix matmul(const utils::Matrix& A, const utils::Matrix& B){ - - if(A.cols() != B.rows()){ - throw std::runtime_error("matmul: dimension mismatch"); - } - - const uint64_t m = A.rows(); - const uint64_t n = A.cols(); // also B.rows() - const uint64_t p = B.cols(); - T tmp; - - utils::Matrix C(m, p, T{0}); - - for (uint64_t i = 0; i < m; ++i){ - for (uint64_t j = 0; j < n; ++j){ - tmp = A(i,j); - for (uint64_t k = 0; k < p; ++k){ - C(i,k) += tmp * B(j,k); - } - } - } - return C; - } - -// ---------------- Rows-only OpenMP ---------------- -template -utils::Matrix matmul_rows_omp(const utils::Matrix& A, - const utils::Matrix& B) { - if (A.cols() != B.rows()) throw std::runtime_error("matmul_rows_omp: dim mismatch"); - const uint64_t m=A.rows(), n=A.cols(), p=B.cols(); - - utils::Matrix C(m, p, T{0}); - - #pragma omp parallel for schedule(static) - for (uint64_t i=0;i -utils::Matrix matmul_collapse_omp(const utils::Matrix& A, - const utils::Matrix& B) { - if (A.cols() != B.rows()) throw std::runtime_error("matmul_collapse_omp: dim mismatch"); - const uint64_t m=A.rows(), n=A.cols(), p=B.cols(); - utils::Matrix C(m, p, T{0}); - - #pragma omp parallel for collapse(2) schedule(static) - for (uint64_t i=0;i -utils::Matrix matmul_auto(const utils::Matrix& A, - const utils::Matrix& B) { - const uint64_t m=A.rows(), p=B.cols(); - const uint64_t work = m * p; - - - - #ifdef _OPENMP - bool can_parallel = omp_config::omp_parallel_allowed(); - uint64_t threads = static_cast(omp_get_max_threads()); - #else - bool can_parallel = false; - uint64_t threads = 1; - #endif - - - // Tiny problems: serial is cheapest. - if (!can_parallel || work < threads*4ull) { - - return matmul(A,B); - } - // Plenty of (i,j) work → collapse(2) is a great default. - else if (work >= 8ull * threads) { - return matmul_collapse_omp(A,B); - } - // Many rows and very few columns → rows-only cheaper overhead. - else if (m >= static_cast(threads) && p <= 4) { - return matmul_rows_omp(A,B); - } - else{ - // Safe fallback - return matmul(A,B); - } -} - - - - - + inline utils::Matrix matmul(const utils::Matrix& A, const utils::Matrix& B){ + return detail::matmul_serial(A, B); + } + } // namespace numerics - -#endif // _matmul_n_ \ No newline at end of file diff --git a/include/numerics/matrandom.h b/include/numerics/matrandom.h deleted file mode 100644 index b64093d..0000000 --- a/include/numerics/matrandom.h +++ /dev/null @@ -1,55 +0,0 @@ -#ifndef _matrandom_n_ -#define _matrandom_n_ - -#include "./utils/vector.h" -#include "./utils/matrix.h" -#include "./core/omp_config.h" - -namespace numerics{ - - template - void inplace_matrandom_add(utils::Matrix& A, const T lower, const T higher) { - - const uint64_t rows = A.rows(); - const uint64_t cols = A.cols(); - - utils::Matrix B; - B.random(rows,cols, lower, higher); - - numerics::inplace_matadd_mat(A, B); - } - - template - void inplace_matrandom_mul(utils::Matrix& A, const T lower, const T higher) { - - const uint64_t rows = A.rows(); - const uint64_t cols = A.cols(); - - utils::Matrix B; - B.random(rows,cols, lower, higher); - - for (uint64_t i = 0; i < rows; ++i){ - for (uint64_t j = 0; j < cols; ++j){ - A(i,j) *= B(i,j); - } - } - } - - - template - utils::Matrix matrandom_add(const utils::Matrix& A, const T lower, const T higher) { - - const uint64_t rows = A.rows(); - const uint64_t cols = A.cols(); - - utils::Matrix B = A; - - numerics::inplace_matadd_mat(B, lower, higher); - - return B; - } - - -} // namespace numerics - -#endif // _matrandom_n_ \ No newline at end of file diff --git a/include/numerics/matvec.h b/include/numerics/matvec.h deleted file mode 100644 index b6d944b..0000000 --- a/include/numerics/matvec.h +++ /dev/null @@ -1,156 +0,0 @@ -#ifndef _matvec_n_ -#define _matvec_n_ - - -#include "./utils/matrix.h" -#include "./core/omp_config.h" - -namespace numerics{ - -// ================================================= -// y = A * x (Matrix–Vector product) -// ================================================= - template - utils::Vector matvec(const utils::Matrix& A, const utils::Vector& x) { - if (A.cols() != x.size()) { - throw std::runtime_error("matvec: dimension mismatch"); - } - - const uint64_t m = A.rows(); - const uint64_t n = A.cols(); - - utils::Vector y(m, T{0}); - - for (uint64_t i = 0; i < m; ++i) { - for (uint64_t j = 0; j < n; ++j) { - y[i] += A(i, j) * x[j]; - } - } - return y; - } - // -------------- Collapse(2) OpenMP ---------------- - template - utils::Vector matvec_omp(const utils::Matrix& A, const utils::Vector& x) { - if (A.cols() != x.size()) { - throw std::runtime_error("matvec: dimension mismatch"); - } - - const uint64_t m = A.rows(); - const uint64_t n = A.cols(); - - utils::Vector y(m, T{0}); // <-- y has length m (rows) - - - const T* xptr = x.data(); - const T* Aptr = A.data(); // row-major: A(i,j) == Aptr[i*n + j] - - // Each row i is an independent dot product: y[i] = dot(A[i,*], x) - #pragma omp parallel for schedule(static) - for (uint64_t i = 0; i < m; ++i) { - const T* row = Aptr + i * n; // contiguous row i - T acc = T{0}; - #pragma omp simd reduction(+:acc) - for (uint64_t j = 0; j < n; ++j) { - acc += row[j] * xptr[j]; - } - y[i] = acc; - } - - return y; - } - - // -------------- Auto OpenMP ---------------- - template - utils::Vector matvec_auto(const utils::Matrix& A, - const utils::Vector& x) { - - - uint64_t work = A.rows() * A.cols(); - - bool can_parallel = omp_config::omp_parallel_allowed(); - #ifdef _OPENMP - int threads = omp_get_max_threads(); - #else - int threads = 1; - #endif - - if (can_parallel || work > static_cast(threads) * 4ull) { - return matvec_omp(A,x); - } - else{ - // Safe fallback - return matvec(A,x); - } - - } - -// ================================================= -// y = x * A (Vector–Matrix product) -// ================================================= - template - utils::Vector vecmat(const utils::Vector& x, const utils::Matrix& A) { - if (x.size() != A.rows()) { - throw std::runtime_error("vecmat: dimension mismatch"); - } - const uint64_t m = A.rows(); - const uint64_t n = A.cols(); - - utils::Vector y(n, T{0}); - - for (uint64_t j = 0; j < n; ++j) { - for (uint64_t i = 0; i < m; ++i) { - y[j] += x[i] * A(i, j); - } - } - return y; - } - - // -------------- Collapse(2) OpenMP ---------------- - template - utils::Vector vecmat_omp(const utils::Vector& x, const utils::Matrix& A) { - if (x.size() != A.rows()) { - throw std::runtime_error("vecmat: dimension mismatch"); - } - const uint64_t m = A.rows(); - const uint64_t n = A.cols(); - - utils::Vector y(n, T{0}); - #pragma omp parallel for schedule(static) - for (uint64_t j = 0; j < n; ++j) { - T acc = T{0}; - for (uint64_t i = 0; i < m; ++i) { - acc += x[i] * A(i, j); - } - y[j] = acc; - } - return y; - } - - // -------------- Auto OpenMP ---------------- - template - utils::Vector vecmat_auto(const utils::Vector& x, - const utils::Matrix& A) { - - uint64_t work = A.rows() * A.cols(); - - bool can_parallel = omp_config::omp_parallel_allowed(); - #ifdef _OPENMP - int threads = omp_get_max_threads(); - #else - int threads = 1; - #endif - - if (can_parallel || work > static_cast(threads) * 4ull) { - return vecmat_omp(x,A); - } - else{ - // Safe fallback - return vecmat(x,A); - } - - } - - -} // namespace numerics - -#endif // _matvec_n_ \ No newline at end of file diff --git a/include/numerics/random.h b/include/numerics/random.h new file mode 100644 index 0000000..c010eae --- /dev/null +++ b/include/numerics/random.h @@ -0,0 +1,97 @@ +#pragma once + +#include "./core/omp_config.h" +#include "detail/random_serial.h" +#include "add.h" +#include "mul.h" + + + +namespace numerics{ + + // ---------------- inplace_random ---------------- + template + inline void inplace_random(utils::Vector& v, const T low, const T high) { + detail::inplace_random_serial(v, low, high); + } + + template + inline void inplace_random(utils::Matrix& A, const T low, const T high) { + detail::inplace_random_serial(A, low, high); + } + + + // ---------------- random ---------------- + template + inline utils::Vector random_vector(const uint64_t n, const T low, const T high) { + utils::Vector v(n); + inplace_random(v, low, high); + return v; + } + + template + inline utils::Matrix random_matrix(const uint64_t n, const uint64_t m, const T low, const T high) { + utils::Matrix A(n, m); + inplace_random(A, low, high); + return A; + } + + // ---------------- inplace random add---------------- + template + inline void inplace_random_add(utils::Vector& v, const T low, const T high){ + utils::Vector noise = random_vector(v.size(), low, high); + inplace_add(v, noise); + } + + template + inline void inplace_random_add(utils::Matrix& A, const T low, const T high){ + utils::Matrix noise = random_matrix(A.rows(), A.cols(), low, high); + inplace_add(A, noise); + } + + // ---------------- random add---------------- + template + inline utils::Vector random_add(const utils::Vector& v, const T low, const T high){ + utils::Vector out = v; + inplace_random_add(out, low, high); + return out; + } + + template + inline utils::Matrix random_add(const utils::Matrix& A, const T low, const T high){ + utils::Matrix out = A; + inplace_random_add(out, low, high); + return out; + } + + + // ---------------- inplace random mul---------------- + template + inline void inplace_random_mul(utils::Vector& v, const T low, const T high){ + utils::Vector noise = random_vector(v.size(), low, high); + inplace_mul(v, noise); + } + + template + inline void inplace_random_mul(utils::Matrix& A, const T low, const T high){ + utils::Matrix noise = random_matrix(A.rows(), A.cols(), low, high); + inplace_mul(A, noise); + } + + // ---------------- random mul---------------- + template + inline utils::Vector random_mul(const utils::Vector& v, const T low, const T high){ + utils::Vector out = v; + inplace_random_mul(out, low, high); + return out; + } + + template + inline utils::Matrix random_mul(const utils::Matrix& A, const T low, const T high){ + utils::Matrix out = A; + inplace_random_mul(out, low, high); + return out; + } + + +} \ No newline at end of file diff --git a/include/numerics/vecrandom.h b/include/numerics/vecrandom.h deleted file mode 100644 index 75e60b5..0000000 --- a/include/numerics/vecrandom.h +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef _vecrandom_n_ -#define _vecrandom_n_ - -#include "./utils/vector.h" -#include "./utils/matrix.h" -#include "./core/omp_config.h" - -namespace numerics{ - - template - void inplace_vecrandom_add(utils::Vector& a, const T lower, const T higher) { - - const uint64_t N = a.size(); - - utils::Vector b; - b.random(N, lower, higher); - - a.inplace_add(b); - } - - template - void inplace_vecrandom_mul(utils::Vector& a, const T lower, const T higher) { - - const uint64_t N = a.size(); - - utils::Vector b; - b.random(N, lower, higher); - - a.inplace_multiply(b); - - } - - - template - utils::Vector vecrandom_add(const utils::Vector& a, const T lower, const T higher) { - - const uint64_t N = a.size(); - - utils::Vector b; - inplace_vecrandom_add(b, lower, higher); - - return b; - } - - -} // namespace numerics - -#endif // _vecrandom_n_ \ No newline at end of file