# include # include # include # include # include "polygon_monte_carlo.h" /******************************************************************************/ int between ( double xa, double ya, double xb, double yb, double xc, double yc ) /******************************************************************************/ /* Purpose: BETWEEN is TRUE if vertex C is between vertices A and B. Discussion: The points must be (numerically) collinear. Given that condition, we take the greater of XA - XB and YA - YB as a "scale" and check where C's value lies. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, double XA, YA, XB, YB, XC, YC, the coordinates of the vertices. Output, int BETWEEN, is TRUE if C is between A and B. */ { int value; double xmax; double xmin; double ymax; double ymin; if ( ! collinear ( xa, ya, xb, yb, xc, yc ) ) { value = 0; } else if ( fabs ( ya - yb ) < fabs ( xa - xb ) ) { xmax = fmax ( xa, xb ); xmin = fmin ( xa, xb ); value = ( xmin <= xc && xc <= xmax ); } else { ymax = fmax ( ya, yb ); ymin = fmin ( ya, yb ); value = ( ymin <= yc && yc <= ymax ); } return value; } /******************************************************************************/ int collinear ( double xa, double ya, double xb, double yb, double xc, double yc ) /******************************************************************************/ /* Purpose: COLLINEAR returns a measure of collinearity for three points. Discussion: In order to deal with collinear points whose coordinates are not numerically exact, we compare the area of the largest square that can be created by the line segment between two of the points to (twice) the area of the triangle formed by the points. If the points are collinear, their triangle has zero area. If the points are close to collinear, then the area of this triangle will be small relative to the square of the longest segment. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, double XA, YA, XB, YB, XC, YC, the coordinates of the vertices. Output, int COLLINEAR, is TRUE if the points are judged to be collinear. */ { double area; const double r8_eps = 2.220446049250313E-016; double side_ab_sq; double side_bc_sq; double side_ca_sq; double side_max_sq; int value; area = 0.5 * ( ( xb - xa ) * ( yc - ya ) - ( xc - xa ) * ( yb - ya ) ); side_ab_sq = pow ( xa - xb, 2 ) + pow ( ya - yb, 2 ); side_bc_sq = pow ( xb - xc, 2 ) + pow ( yb - yc, 2 ); side_ca_sq = pow ( xc - xa, 2 ) + pow ( yc - ya, 2 ); side_max_sq = fmax ( side_ab_sq, fmax ( side_bc_sq, side_ca_sq ) ); if ( side_max_sq <= r8_eps ) { value = 1; } else if ( 2.0 * fabs ( area ) <= r8_eps * side_max_sq ) { value = 1; } else { value = 0; } return value; } /******************************************************************************/ int diagonal ( int im1, int ip1, int n, int prev[], int next[], double x[], double y[] ) /******************************************************************************/ /* Purpose: DIAGONAL: VERTEX(IM1) to VERTEX(IP1) is a proper internal diagonal. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, int IM1, IP1, the indices of two vertices. Input, int N, the number of vertices. Input, int PREV[N], the previous neighbor of each vertex. Input, int NEXT[N], the next neighbor of each vertex. Input, double X[N], Y[N], the coordinates of each vertex. Output, int DIAGONAL, the value of the test. */ { int value; int value1; int value2; int value3; value1 = in_cone ( im1, ip1, n, prev, next, x, y ); value2 = in_cone ( ip1, im1, n, prev, next, x, y ); value3 = diagonalie ( im1, ip1, n, next, x, y ); value = ( value1 && value2 && value3 ); return value; } /******************************************************************************/ int diagonalie ( int im1, int ip1, int n, int next[], double x[], double y[] ) /******************************************************************************/ /* Purpose: DIAGONALIE is true if VERTEX(IM1):VERTEX(IP1) is a proper diagonal. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, int IM1, IP1, the indices of two vertices. Input, int N, the number of vertices. Input, int NEXT[N], the next neighbor of each vertex. Input, double X[N], Y[N], the coordinates of each vertex. Output, int DIAGONALIE, the value of the test. */ { int first; int j; int jp1; int value; int value2; first = im1; j = first; jp1 = next[first]; value = 1; /* For each edge VERTEX(J):VERTEX(JP1) of the polygon: */ while ( 1 ) { /* Skip any edge that includes vertex IM1 or IP1. */ if ( j == im1 || j == ip1 || jp1 == im1 || jp1 == ip1 ) { } else { value2 = intersect ( x[im1], y[im1], x[ip1], y[ip1], x[j], y[j], x[jp1], y[jp1] ); if ( value2 ) { value = 0; break; } } j = jp1; jp1 = next[j]; if ( j == first ) { break; } } return value; } /******************************************************************************/ int in_cone ( int im1, int ip1, int n, int prev[], int next[], double x[], double y[] ) /******************************************************************************/ /* Purpose: IN_CONE is TRUE if the diagonal VERTEX(IM1):VERTEX(IP1) is strictly internal. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, int IM1, IP1, the indices of two vertices. Input, int N, the number of vertices. Input, int PREV[N], the previous neighbor of each vertex. Input, int NEXT[N], the next neighbor of each vertex. Input, double X[N], Y[N], the coordinates of each vertex. Output, int IN_CONE, the value of the test. */ { int i; int im2; double t1; double t2; double t3; double t4; double t5; int value; im2 = prev[im1]; i = next[im1]; t1 = triangle_area ( x[im1], y[im1], x[i], y[i], x[im2], y[im2] ); if ( 0.0 <= t1 ) { t2 = triangle_area ( x[im1], y[im1], x[ip1], y[ip1], x[im2], y[im2] ); t3 = triangle_area ( x[ip1], y[ip1], x[im1], y[im1], x[i], y[i] ); value = ( ( 0.0 < t2 ) && ( 0.0 < t3 ) ); } else { t4 = triangle_area ( x[im1], y[im1], x[ip1], y[ip1], x[i], y[i] ); t5 = triangle_area ( x[ip1], y[ip1], x[im1], y[im1], x[im2], y[im2] ); value = ! ( ( 0.0 <= t4 ) && ( 0.0 <= t5 ) ); } return value; } /******************************************************************************/ int intersect ( double xa, double ya, double xb, double yb, double xc, double yc, double xd, double yd ) /******************************************************************************/ /* Purpose: INTERSECT is true if lines VA:VB and VC:VD intersect. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, double XA, YA, XB, YB, XC, YC, XD, YD, the X and Y coordinates of the four vertices. Output, int INTERSECT, the value of the test. */ { int value; if ( intersect_prop ( xa, ya, xb, yb, xc, yc, xc, yd ) ) { value = 1; } else if ( between ( xa, ya, xb, yb, xc, yc ) ) { value = 1; } else if ( between ( xa, ya, xb, yb, xd, yd ) ) { value = 1; } else if ( between ( xc, yc, xd, yd, xa, ya ) ) { value = 1; } else if ( between ( xc, yc, xd, yd, xb, yb ) ) { value = 1; } else { value = 0; } return value; } /******************************************************************************/ int intersect_prop ( double xa, double ya, double xb, double yb, double xc, double yc, double xd, double yd ) /******************************************************************************/ /* Purpose: INTERSECT_PROP is TRUE if lines VA:VB and VC:VD have a proper intersection. Licensing: This code is distributed under the MIT license. Modified: 04 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, double XA, YA, XB, YB, XC, YC, XD, YD, the X and Y coordinates of the four vertices. Output, int INTERSECT_PROP, the result of the test. */ { double t1; double t2; double t3; double t4; int value; int value1; int value2; int value3; int value4; if ( collinear ( xa, ya, xb, yb, xc, yc ) ) { value = 0; } else if ( collinear ( xa, ya, xb, yb, xd, yd ) ) { value = 0; } else if ( collinear ( xc, yc, xd, yd, xa, ya ) ) { value = 0; } else if ( collinear ( xc, yc, xd, yd, xb, yb ) ) { value = 0; } else { t1 = triangle_area ( xa, ya, xb, yb, xc, yc ); t2 = triangle_area ( xa, ya, xb, yb, xd, yd ); t3 = triangle_area ( xc, yc, xd, yd, xa, ya ); t4 = triangle_area ( xc, yc, xd, yd, xb, yb ); value1 = ( 0.0 < t1 ); value2 = ( 0.0 < t2 ); value3 = ( 0.0 < t3 ); value4 = ( 0.0 < t4 ); value = ( l4_xor ( value1, value2 ) ) && ( l4_xor ( value3, value4 ) ); } return value; } /******************************************************************************/ int l4_xor ( int l1, int l2 ) /******************************************************************************/ /* Purpose: L4_XOR returns the exclusive OR of two L4's. Discussion: An L4 is a logical value. Licensing: This code is distributed under the MIT license. Modified: 04 May 2014 Author: John Burkardt Parameters: Input, int L1, L2, two values whose exclusive OR is needed. Output, int L4_XOR, the exclusive OR of L1 and L2. */ { int value; int value1; int value2; value1 = ( l1 && ( ! l2 ) ); value2 = ( ( ! l1 ) && l2 ); value = ( value1 || value2 ); return value; } /******************************************************************************/ 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 MIT 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 = 0; j < n; j++ ) { v[j] = v[j] * pow ( x[i+j*m], e[i] ); } } } return v; } /******************************************************************************/ double polygon_area ( int nv, double v[] ) /******************************************************************************/ /* Purpose: POLYGON_AREA determines the area of a polygon. Licensing: This code is distributed under the MIT license. Modified: 09 May 2014 Author: John Burkardt Parameters: Input, int NV, the number of vertices of the polygon. Input, double V[2*N], the vertex coordinates. Output, double POLYGON_AREA, the area of the polygon. */ { double area; int e[2]; e[0] = 0; e[1] = 0; area = polygon_monomial_integral ( nv, v, e ); return area; } /******************************************************************************/ double polygon_monomial_integral ( int nv, double v[], int e[] ) /******************************************************************************/ /* Purpose: POLYGON_INTEGRAL integrates a monomial over a polygon. Discussion: Nu(P,Q) = Integral ( x, y in polygon ) x^e(1) y^e(2) dx dy Licensing: This code is distributed under the MIT license. Modified: 09 May 2014 Author: John Burkardt Reference: Carsten Steger, On the calculation of arbitrary moments of polygons, Technical Report FGBV-96-05, Forschungsgruppe Bildverstehen, Informatik IX, Technische Universitaet Muenchen, October 1996. Parameters: Input, int NV, the number of vertices of the polygon. Input, double V[2*N], the vertex coordinates. Input, int E[2], the exponents of the monomial. Output, double POLYGON_MONOMIAL_INTEGRAL, the unnormalized moment Nu(P,Q). */ { int i; int k; int l; double nu_pq; int p; int q; double s_pq; double xi; double xj; double yi; double yj; p = e[0]; q = e[1]; nu_pq = 0.0; xj = v[0+(nv-1)*2]; yj = v[1+(nv-1)*2]; for ( i = 0; i < nv; i++ ) { xi = v[0+i*2]; yi = v[1+i*2]; s_pq = 0.0; for ( k = 0; k <= p; k++ ) { for ( l = 0; l <= q; l++ ) { s_pq = s_pq + r8_choose ( k + l, l ) * r8_choose ( p + q - k - l, q - l ) * pow ( xi, k ) * pow ( xj, p - k ) * pow ( yi, l ) * pow ( yj, q - l ); } } nu_pq = nu_pq + ( xj * yi - xi * yj ) * s_pq; xj = xi; yj = yi; } nu_pq = nu_pq / ( double ) ( p + q + 2 ) / ( double ) ( p + q + 1 ) / r8_choose ( p + q, p ); return nu_pq; } /******************************************************************************/ double *polygon_sample ( int nv, double v[], int n, int *seed ) /******************************************************************************/ /* Purpose: POLYGON_SAMPLE uniformly samples a polygon. Licensing: This code is distributed under the MIT license. Modified: 07 May 2014 Author: John Burkardt Parameters: Input, int NV, the number of vertices. Input, double V[2*NV], the vertices of the polygon, listed in counterclockwise order. Input, int N, the number of points to create. Input/output, int *SEED, a seed for the random number generator. Output, double POLYGON_SAMPLE[2*N], the points. */ { double *area_cumulative; double area_polygon; double *area_relative; double *area_triangle; double area_percent; int i; int j; int k; double *r; double *s; int *triangles; double *x; double *y; /* Triangulate the polygon. */ x = ( double * ) malloc ( nv * sizeof ( double ) ); y = ( double * ) malloc ( nv * sizeof ( double ) ); for ( i = 0; i < nv; i++ ) { x[i] = v[0+i*2]; y[i] = v[1+i*2]; } triangles = polygon_triangulate ( nv, x, y ); /* Determine the areas of each triangle. */ area_triangle = ( double * ) malloc ( ( nv - 2 ) * sizeof ( double ) ); for ( i = 0; i < nv - 2; i++ ) { area_triangle[i] = triangle_area ( v[0+triangles[0+i*3]*2], v[1+triangles[0+i*3]*2], v[0+triangles[1+i*3]*2], v[1+triangles[1+i*3]*2], v[0+triangles[2+i*3]*2], v[1+triangles[2+i*3]*2] ); } /* Normalize the areas. */ area_polygon = r8vec_sum ( nv - 2, area_triangle ); area_relative = ( double * ) malloc ( ( nv - 2 ) * sizeof ( double ) ); for ( i = 0; i < nv - 2; i++ ) { area_relative[i] = area_triangle[i] / area_polygon; } /* Replace each area by the sum of itself and all previous ones. */ area_cumulative = ( double * ) malloc ( ( nv - 2 ) * sizeof ( double ) ); area_cumulative[0] = area_relative[0]; for ( i = 1; i < nv - 2; i++ ) { area_cumulative[i] = area_relative[i] + area_cumulative[i-1]; } s = ( double * ) malloc ( 2 * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { /* Choose triangle I at random, based on areas. */ area_percent = r8_uniform_01 ( seed ); for ( k = 0; k < nv - 2; k++ ) { i = k; if ( area_percent <= area_cumulative[k] ) { break; } } /* Now choose a point at random in triangle I. */ r = r8vec_uniform_01_new ( 2, seed ); if ( 1.0 < r[0] + r[1] ) { r[0] = 1.0 - r[0]; r[1] = 1.0 - r[1]; } s[0+j*2] = ( 1.0 - r[0] - r[1] ) * v[0+triangles[0+i*3]*2] + r[0] * v[0+triangles[1+i*3]*2] + r[1] * v[0+triangles[2+i*3]*2]; s[1+j*2] = ( 1.0 - r[0] - r[1] ) * v[1+triangles[0+i*3]*2] + r[0] * v[1+triangles[1+i*3]*2] + r[1] * v[1+triangles[2+i*3]*2]; free ( r ); } free ( area_cumulative ); free ( area_relative ); free ( area_triangle ); free ( triangles ); free ( x ); free ( y ); return s; } /******************************************************************************/ int *polygon_triangulate ( int n, double x[], double y[] ) /******************************************************************************/ /* Purpose: POLYGON_TRIANGULATE determines a triangulation of a polygon. Discussion: There are N-3 triangles in the triangulation. For the first N-2 triangles, the first edge listed is always an internal diagonal. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: Original C version by Joseph ORourke. This C version by John Burkardt. Reference: Joseph ORourke, Computational Geometry in C, Cambridge, 1998, ISBN: 0521649765, LC: QA448.D38. Parameters: Input, int N, the number of vertices. Input, double X[N], Y[N], the coordinates of each vertex. Output, int TRIANGLES[3*(N-2)], the triangles of the triangulation. */ { double area; int *ear; int i; int i0; int i1; int i2; int i3; int i4; int *next; int node; int node_m1; int *prev; int triangle_num; int *triangles; /* We must have at least 3 vertices. */ if ( n < 3 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "POLYGON_TRIANGULATE - Fatal error!\n" ); fprintf ( stderr, " N < 3.\n" ); exit ( 1 ); } /* Consecutive vertices cannot be equal. */ node_m1 = n - 1; for ( node = 0; node < n; node++ ) { if ( x[node_m1] == x[node] && y[node_m1] == y[node] ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "POLYGON_TRIANGULATE - Fatal error!\n" ); fprintf ( stderr, " Two consecutive nodes are identical.\n" ); exit ( 1 ); } node_m1 = node; } /* Area must be positive. */ area = 0.0; for ( node = 0; node < n - 2; node++ ) { area = area + 0.5 * ( ( x[node+1] - x[node] ) * ( y[node+2] - y[node] ) - ( x[node+2] - x[node] ) * ( y[node+1] - y[node] ) ); } if ( area <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "POLYGON_TRIANGULATE - Fatal error!\n" ); fprintf ( stderr, " Polygon has zero or negative area.\n" ); exit ( 1 ); } triangles = ( int * ) malloc ( 3 * ( n - 2 ) * sizeof ( int ) ); /* PREV and NEXT point to the previous and next nodes. */ prev = ( int * ) malloc ( n * sizeof ( int ) ); next = ( int * ) malloc ( n * sizeof ( int ) ); i = 0; prev[i] = n - 1; next[i] = i + 1; for ( i = 1; i < n - 1; i++ ) { prev[i] = i - 1; next[i] = i + 1; } i = n - 1; prev[i] = i - 1; next[i] = 0; /* EAR indicates whether the node and its immediate neighbors form an ear that can be sliced off immediately. */ ear = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ear[i] = diagonal ( prev[i], next[i], n, prev, next, x, y ); } triangle_num = 0; i2 = 0; while ( triangle_num < n - 3 ) { /* If I2 is an ear, gather information necessary to carry out the slicing operation and subsequent "healing". */ if ( ear[i2] ) { i3 = next[i2]; i4 = next[i3]; i1 = prev[i2]; i0 = prev[i1]; /* Make vertex I2 disappear. */ next[i1] = i3; prev[i3] = i1; /* Update the earity of I1 and I3, because I2 disappeared. */ ear[i1] = diagonal ( i0, i3, n, prev, next, x, y ); ear[i3] = diagonal ( i1, i4, n, prev, next, x, y ); /* Add the diagonal [I3, I1, I2] to the list. */ triangles[0+triangle_num*3] = i3; triangles[1+triangle_num*3] = i1; triangles[2+triangle_num*3] = i2; triangle_num = triangle_num + 1; } /* Try the next vertex. */ i2 = next[i2]; } /* The last triangle is formed from the three remaining vertices. */ i3 = next[i2]; i1 = prev[i2]; triangles[0+triangle_num*3] = i3; triangles[1+triangle_num*3] = i1; triangles[2+triangle_num*3] = i2; triangle_num = triangle_num + 1; free ( ear ); free ( next ); free ( prev ); return triangles; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: R8_CHOOSE computes the binomial coefficient C(N,K) as an R8. Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in R8 arithmetic. The formula used is: C(N,K) = N! / ( K! * (N-K)! ) Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Parameters: Input, int N, K, the values of N and K. Output, double R8_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; double value; if ( k < n - k ) { mn = k; mx = n - k; } else { mn = n - k; mx = k; } if ( mn < 0 ) { value = 0.0; } else if ( mn == 0 ) { value = 1.0; } else { value = ( double ) ( mx + 1 ); for ( i = 2; i <= mn; i++ ) { value = ( value * ( double ) ( mx + i ) ) / ( double ) i; } } return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { const int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 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, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ double triangle_area ( double xa, double ya, double xb, double yb, double xc, double yc ) /******************************************************************************/ /* Purpose: TRIANGLE_AREA computes the signed area of a triangle. Licensing: This code is distributed under the MIT license. Modified: 05 May 2014 Author: John Burkardt Parameters: Input, double XA, YA, XB, YB, XC, YC, the coordinates of the vertices of the triangle, given in counterclockwise order. Output, double TRIANGLE_AREA, the signed area of the triangle. */ { double value; value = 0.5 * ( ( xb - xa ) * ( yc - ya ) - ( xc - xa ) * ( yb - ya ) ); return value; }