# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include "ball_integrals.h"
/******************************************************************************/
double ball01_monomial_integral ( int e[3] )
/******************************************************************************/
/*
Purpose:
BALL01_MONOMIAL_INTEGRAL returns monomial integrals in the unit ball.
Discussion:
The integration region is
X^2 + Y^2 + Z^2 <= 1.
The monomial is F(X,Y,Z) = X^E(1) * Y^E(2) * Z^E(3).
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
10 May 2014
Author:
John Burkardt
Reference:
Philip Davis, Philip Rabinowitz,
Methods of Numerical Integration,
Second Edition,
Academic Press, 1984, page 263.
Parameters:
Input, int E[3], the exponents of X, Y and Z in the
monomial. Each exponent must be nonnegative.
Output, double BALL01_MONOMIAL_INTEGRAL, the integral.
*/
{
int i;
double integral;
const double r = 1.0;
double s;
if ( e[0] < 0 || e[1] < 0 || e[2] < 0 )
{
fprintf ( stderr, "\n" );
fprintf ( stderr, "BALL01_MONOMIAL_INTEGRAL - Fatal error!\n" );
fprintf ( stderr, " All exponents must be nonnegative.\n" );
fprintf ( stderr, " E[0] = %d\n", e[0] );
fprintf ( stderr, " E[1] = %d\n", e[1] );
fprintf ( stderr, " E[2] = %d\n", e[2] );
exit ( 1 );
}
if ( ( e[0] % 2 ) == 1 ||
( e[1] % 2 ) == 1 ||
( e[2] % 2 ) == 1 )
{
integral = 0.0;
}
else
{
integral = 2.0;
for ( i = 0; i < 3; i++ )
{
integral = integral * tgamma ( 0.5 * ( double ) ( e[i] + 1 ) );
}
integral = integral
/ tgamma ( 0.5 * ( double ) ( e[0] + e[1] + e[2] + 3 ) );
}
/*
The surface integral is now adjusted to give the volume integral.
*/
s = e[0] + e[1] + e[2] + 3;
integral = integral * pow ( r, s ) / ( double ) ( s );
return integral;
}
/******************************************************************************/
double *ball01_sample ( int n, int *seed )
/******************************************************************************/
/*
Purpose:
BALL01_SAMPLE uniformly samples points from the unit ball.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
02 January 2014
Author:
John Burkardt
Reference:
Russell Cheng,
Random Variate Generation,
in Handbook of Simulation,
edited by Jerry Banks,
Wiley, 1998, pages 168.
Reuven Rubinstein,
Monte Carlo Optimization, Simulation, and Sensitivity
of Queueing Networks,
Krieger, 1992,
ISBN: 0894647644,
LC: QA298.R79.
Parameters:
Input, int N, the number of points.
Input/output, int *SEED, a seed for the random
number generator.
Output, double X[3*N], the points.
*/
{
const double exponent = 1.0 / 3.0;
int i;
int j;
double norm;
double r;
double *x;
x = r8mat_normal_01_new ( 3, n, seed );
for ( j = 0; j < n; j++ )
{
/*
Compute the length of the vector.
*/
norm = sqrt ( pow ( x[0+j*3], 2 )
+ pow ( x[1+j*3], 2 )
+ pow ( x[2+j*3], 2 ) );
/*
Normalize the vector.
*/
for ( i = 0; i < 3; i++ )
{
x[i+j*3] = x[i+j*3] / norm;
}
/*
Now compute a value to map the point ON the sphere INTO the sphere.
*/
r = r8_uniform_01 ( seed );
for ( i = 0; i < 3; i++ )
{
x[i+j*3] = pow ( r, exponent ) * x[i+j*3];
}
}
return x;
}
/******************************************************************************/
double ball01_volume ( )
/******************************************************************************/
/*
Purpose:
BALL01_VOLUME returns the volume of the unit ball.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
02 January 2014
Author:
John Burkardt
Parameters:
Output, double BALL01_VOLUME, the volume of the unit ball.
*/
{
const double r = 1.0;
const double r8_pi = 3.141592653589793;
double volume;
volume = 4.0 * r8_pi * pow ( r, 3 ) / 3.0;
return volume;
}
/******************************************************************************/
int *i4vec_uniform_ab_new ( int n, int a, int b, int *seed )
/******************************************************************************/
/*
Purpose:
I4VEC_UNIFORM_AB_NEW returns a scaled pseudorandom I4VEC.
Discussion:
The pseudorandom numbers should be uniformly distributed
between A and B.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
06 January 2014
Author:
John Burkardt
Reference:
Paul Bratley, Bennett Fox, Linus Schrage,
A Guide to Simulation,
Second Edition,
Springer, 1987,
ISBN: 0387964673,
LC: QA76.9.C65.B73.
Bennett Fox,
Algorithm 647:
Implementation and Relative Efficiency of Quasirandom
Sequence Generators,
ACM Transactions on Mathematical Software,
Volume 12, Number 4, December 1986, pages 362-376.
Pierre L'Ecuyer,
Random Number Generation,
in Handbook of Simulation,
edited by Jerry Banks,
Wiley, 1998,
ISBN: 0471134031,
LC: T57.62.H37.
Peter Lewis, Allen Goodman, James Miller,
A Pseudo-Random Number Generator for the System/360,
IBM Systems Journal,
Volume 8, Number 2, 1969, pages 136-143.
Parameters:
Input, integer N, the dimension of the vector.
Input, int A, B, the limits of the interval.
Input/output, int *SEED, the "seed" value, which should NOT be 0.
On output, SEED has been updated.
Output, int I4VEC_UNIFORM_AB_NEW[N], a vector of random values
between A and B.
*/
{
int c;
int i;
const int i4_huge = 2147483647;
int k;
float r;
int value;
int *x;
if ( *seed == 0 )
{
fprintf ( stderr, "\n" );
fprintf ( stderr, "I4VEC_UNIFORM_AB_NEW - Fatal error!\n" );
fprintf ( stderr, " Input value of SEED = 0.\n" );
exit ( 1 );
}
/*
Guaranteee A <= B.
*/
if ( b < a )
{
c = a;
a = b;
b = c;
}
x = ( int * ) malloc ( n * sizeof ( int ) );
for ( i = 0; i < n; i++ )
{
k = *seed / 127773;
*seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
if ( *seed < 0 )
{
*seed = *seed + i4_huge;
}
r = ( float ) ( *seed ) * 4.656612875E-10;
/*
Scale R to lie between A-0.5 and B+0.5.
*/
r = ( 1.0 - r ) * ( ( float ) a - 0.5 )
+ r * ( ( float ) b + 0.5 );
/*
Use rounding to convert R to an integer between A and B.
*/
value = round ( r );
/*
Guarantee A <= VALUE <= B.
*/
if ( value < a )
{
value = a;
}
if ( b < value )
{
value = b;
}
x[i] = value;
}
return x;
}
/******************************************************************************/
double *monomial_value ( int m, int n, int e[], double x[] )
/******************************************************************************/
/*
Purpose:
MONOMIAL_VALUE evaluates a monomial.
Discussion:
This routine evaluates a monomial of the form
product ( 1 <= i <= m ) x(i)^e(i)
where the exponents are nonnegative integers. Note that
if the combination 0^0 is encountered, it should be treated
as 1.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
08 May 2014
Author:
John Burkardt
Parameters:
Input, int M, the spatial dimension.
Input, int N, the number of points at which the
monomial is to be evaluated.
Input, int E[M], the exponents.
Input, double X[M*N], the point coordinates.
Output, double MONOMIAL_VALUE[N], the value of the monomial.
*/
{
int i;
int j;
double *v;
v = ( double * ) malloc ( n * sizeof ( double ) );
for ( j = 0; j < n; j++ )
{
v[j] = 1.0;
}
for ( i = 0; i < m; i++ )
{
if ( 0 != e[i] )
{
for ( j =