1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
/****************************************************************************
* MODULE: R-Tree library
*
* AUTHOR(S): Antonin Guttman - original code
* Daniel Green (green@superliminal.com) - major clean-up
* and implementation of bounding spheres
*
* PURPOSE: Multidimensional index
*
* COPYRIGHT: (C) 2001 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*****************************************************************************/
/*
* SPHERE VOLUME
* by Daniel Green
* dgreen@superliminal.com
*
* Calculates and prints the volumes of the unit hyperspheres for
* dimensions zero through the given value, or 9 by default.
* Prints in the form of a C array of double called sphere_volumes.
*
* From formule in "Regular Polytopes" by H.S.M Coxeter, the volume
* of a hypersphere of dimension d is:
* Pi^(d/2) / gamma(d/2 + 1)
*
* This implementation works by first computing the log of the above
* function and then returning the exp of that value in order to avoid
* instabilities due to the huge values that the real gamma function
* would return.
*
* Multiply the output volumes by R^n to get the volume of an n
* dimensional sphere of radius R.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#ifndef M_PI
# define M_PI 3.1415926535
#endif
static void print_volume(int dimension, double volume)
{
fprintf (stdout, "\t%.6f, /* dimension %3d */\n", volume, dimension);
}
static double sphere_volume(double dimension)
{
/* static const double log_pi = log(M_PI); */
static const double log_pi = 0.49714987268;
double log_gamma, log_volume;
log_gamma = gamma(dimension/2.0 + 1);
log_volume = dimension/2.0 * log_pi - log_gamma;
return exp(log_volume);
}
extern int main(int argc, char *argv[])
{
int dim, max_dims=9;
if(2 == argc)
max_dims = atoi(argv[1]);
fprintf (stdout, "static const double sphere_volumes[] = {\n");
for(dim=0; dim<max_dims+1; dim++)
print_volume(dim, sphere_volume(dim));
fprintf (stdout, "};\n");
return 0;
}
|
