Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88433677
math_tools.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 18, 18:50
Size
6 KB
Mime Type
text/x-c
Expires
Sun, Oct 20, 18:50 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21766996
Attached To
rLIBMULTISCALE LibMultiScale
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__
Event Timeline
Log In to Comment