Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/math/math_tools.cc
8ae23a91b711lammps_2020

Last Change
Raw File

math_tools.cc
View Options

/**
* @file math_tools.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
*
* @date Thu Jun 27 14:49:39 2013
*
* @brief This provides access to a set of math functions for various usage
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#include "lm_common.hh"
#include <cmath>
#include <cstdlib>
#include <vector>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
namespace MathTools {
Real ipow(Real x, int y) {
int i;
Real res = 1.0;
for (i = y; i > 0; --i) {
res *= x;
}
return res;
}
/* -------------------------------------------------------------------------- */
Real zeta7[20], zeta13[20], zeta12[20], zeta6[20];
/* -------------------------------------------------------------------------- */
Real zeta(UInt puis, UInt ordre) {
Real res = 0;
UInt i, pos;
pos = ordre;
if (ordre > 19)
pos = 19;
if (puis == 7) {
if (zeta7[pos] == 0.0) {
for (i = 1; i <= ordre; ++i) {
zeta7[pos] += 1 / ipow(i, puis);
}
}
return zeta7[pos];
}
if (puis == 13) {
if (zeta13[pos] == 0.0) {
for (i = 1; i <= ordre; ++i) {
zeta13[pos] += 1 / ipow(i, puis);
}
}
return zeta13[pos];
}
if (puis == 12) {
if (zeta12[pos] == 0.0) {
for (i = 1; i <= ordre; ++i) {
zeta12[pos] += 1 / ipow(i, puis);
}
}
return zeta12[pos];
}
if (puis == 6) {
if (zeta6[pos] == 0.0) {
for (i = 1; i <= ordre; ++i) {
zeta6[pos] += 1 / ipow(i, puis);
}
}
return zeta6[pos];
}
for (i = 1; i <= ordre; ++i) {
res += 1 / ipow(i, puis);
}
return res;
}
/* ------------------------------------------------------------------------ */
static long seed;
/* -------------------------------------------------------------------------- */
void setSeed(long s) {
seed = s;
srand48(seed);
}
/* ------------------------------------------------------------------------- */
/* A uniform distribution (between a and b) */
Real uniformRandom(Real a, Real b) { return drand48() * (b - a) + a; }
/* -------------------------------------------------------------------------- */
/* A uniform distribution (between 0 and 1) */
Real uniform() { return drand48(); }
/* -------------------------------------------------------------------------- */
/* A normal (Gaussian) distribution (with mean mu and stdev sig) */
Real gaussian() {
// static int iset=0;
Real fac, rsq, v1, v2;
/* If selec=0 we use zero mean and variance 1 */
/* else we use the specified mean and variance */
do {
v1 = (2.0 * uniform() - 1.0);
v2 = (2.0 * uniform() - 1.0);
rsq = v1 * v1 + v2 * v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log(rsq) / rsq);
return (v2 * fac);
}
/* -------------------------------------------------------------------------- */
//! function to compute standard deviation of a set of values stored in an array
Real stdev(Real *x, Real mu, UInt n) {
UInt i;
Real std = 0.0;
Real tmp = 0.0;
for (i = 0; i < n; i++) {
tmp = (*(x + i) - mu);
std += (tmp * tmp);
}
return sqrt(std / n);
}
Real stdev(std::vector<Real> &x, Real mu) { return stdev(&x[0], mu, x.size()); }
/* ------------------------------------------------------------------------ */
//! function to compute the kurtosis of a set of values stored in an array
Real kurtosis(Real *x, Real mu, Real std, int n) {
int i;
Real kur = 0.0;
Real tmp = 0.0;
for (i = 0; i < n; i++) {
tmp = (*(x + i) - mu);
kur += (tmp * tmp * tmp * tmp);
}
kur /= (n * std * std * std * std);
return kur - 3;
}
/* -------------------------------------------------------------------------- */
//! function to compute the skewness of a set of values stored in an array
Real skewness(Real *x, Real mu, Real std, int n) {
int i;
Real skw = 0.0;
Real tmp = 0.0;
for (i = 0; i < n; i++) {
tmp = (*(x + i) - mu);
skw += (tmp * tmp * tmp);
}
skw /= (n * std * std * std);
return skw;
}
/* -------------------------------------------------------------------------- */
//! function to compute the average of a set of values stored in an array
Real average(Real *x, UInt n) {
UInt i;
Real avg = 0.0;
for (i = 0; i < n; i++) {
avg += *(x + i);
}
return avg / n;
}
Real average(std::vector<Real> &x) { return average(&x[0], x.size()); }
/* -------------------------------------------------------------------------- */
}
__END_LIBMULTISCALE__
c4science · Help