/*****************************************************************************
Treibergs Mar. 10, 2006
Problem 6a.
Probability that a random triangle is acute.
Generate triangles with random vertices (x1,y1), (x2,y2), (x3,y3)
in the unit square, so in principle,
x1 = rand() / randmaxpo;
where randmaxpo = RAND_MAX + 1.0 which is the largest of possible nonnegative
random integers generated by the computer. Then
(S1) (x2-x1)*(y2-y1) + (x3-x1)*(y3-y1) > 0.0
iff the vertex at (x1,y1) has an acute angle. Similar expressions (S2) and
(S3) test if the other vertices are acute.
Then we loop from 1 to n. we generate three random coordinates from 0 to 1.
Then if (S1) && (S2) && (S3) the triangle is acute, so we increment a
counter j, so that after the loop is done, ther fraction of acute triangles
is j/n which is the approximate Monte Carlo probability that the triangle
is acute.
We can speed up the computation first by using a little probability theory.
First, triangles that have right angles occur with probability zero, so they
may be ignored in the considerations. Then, since there can be at most one
obtuse angle in a triangle (for otherwise the sum of the angles at the
vertices would exceed 180 degrees) a random triangle has either (A) no obtuse
angles or (B) the first vertex is obtuse and the others are not, or (C) the
second is the obtuse angle or (D) the third is the obtuse angle. This means that
the set of all triangles consist of four mutually exclusive possibilities
(A) or (B) or (C) or (D). This means that the probabilities of these events
satisfy p(A) + p(B) + p(C) + p(D) = 1. Furthermore, by symmetry one angle is
as good as another so that the probability of the first vertex being obtuse is
the same as the probability of the second being obtuse is the same as the
probability of the third being obtuse or p(B)=p(C)=p(D). it follows that if
we knew p(B) then the probability that a random triangle is acute is
p(A) = 1 - 3p(B). Now estimating the probability that vertex one is acute
can be done faster in the loop because the test has only one part
(x2-x1)*(y2-y1) + (x3-x1)*(y3-y1) <= 0.0
Two things speed up the test. First is that we can multiply this inequality by
the square of randmaxpo so that this is true whether or not we divide the \
coordinates by randmaxpo. Thus the test is done by forming products of
large integers. Since this may overflow, we do the computation in double
precision. This can be further sped up by not bothering to store x2,y2,x3,y3
because those quantities occur only once in the test. Thus the main
loop that computes p(B) is
for ( i=1; i <= n; i=i+1 )
{
x1 = rand ();
y1 = rand ();
if ( (double)( rand () - x1 ) * ( rand () - x1 ) +
(double)( rand () - y1 ) * ( rand () - y1 ) <= 0.0 )
j++;
}
p(B) is approximately j/n.
Some other programming nortes: We call the time once to initializre the random
number seed. (Otherwise, each time we ran the program would yield the same
random numbers.) We call time at the start and end of the loop to time the
computation, and compute the difference in start and end times. The header
file time.h has the time functions as well as the time type called
time_t. It takes about 17 sec for a hundred million triangles!
If INT_MAX = 32767 on your machine, you will have to decrease n or use
long integers.
M2160s6a.c
*****************************************************************************/
# include
# include
# include
int
main ( void )
{
register int x1,y1;
double p ,s=0.0, diff;
int i, j, k, n = 10000000;
time_t start_time, end_time;
srand ( (unsigned int)time ( NULL ) ); /*seed for rand is time of day*/
printf ( "\n Monte Carlo Probability that Random Triangle in Unit Square is Acute\n\n");
start_time=time(NULL);
for ( k=1; k<= 10; k++ )
{
j=0;
for ( i=1; i <= n; i=i+1 )
{
x1 = rand ();
y1 = rand () ;
if ( (double)( rand () - x1 ) * ( rand () - x1 ) +
(double)( rand () - y1 ) * ( rand () - y1 ) <= 0.0 )
j++;
}
p = 1.0-3.0* j / n;
s = s + p;
printf ( "%4d Approximate probability triangle is acute = %21.15f\n", k, p );
}
end_time = time(NULL);
diff = difftime( end_time, start_time );
p = s / 10.0;
printf ( " Average approximate probability is %21.15f\n\n", p );
printf ( " Number of triangles simulated = %d\n", n * 10);
printf ( " Total time taken is about %d seconds.\n", (int)diff);
return EXIT_SUCCESS;
}
/***********************************************************************************
Treibergs April 5, 2006
Homework 6b. Define structures point, segment, circle. Enter segment and circle.
Determine whether the segment is inside the circle, outside the circle or neither.
Let S denote the segment and C the circle. If we let D to denote the open
disk strictly inside the circle and E the open region outside the circle not
including the circle itself, we shall say that S is OUTSIDE if S is a subset of
E, and S is INSIDE if S is a subset of D. In all other cases we shall say
NEITHER.
Here is the logic in outline form. Let p1 and p2 denote the endpoints of the
segment. Let r be the radius of the circle and c its center. Let r1 be the
distance p1 to the center and r2 the denote the distance from p2 to the
center. If p1 and p2 are distinct, then T = p2 - p1 is a vector tangent to
the line determined by the segment. Let n be a unit vector perpendicular to T
and d = (p1 - C) . n the signed distance from the center to the line. (This is
the length of the projection of the p1 - c vector on the line in the n
direction.) Then the point on the line closest to the center is w = c + d n.
If all we wanted was to see if the segment is INSIDE, OUTSIDE or NEITHER, here is
the decision outline. This program is given first.
A.) If both r1 < r AND r2 < r. Then the segment is INSIDE.
B.) If either (r1 <= r AND r2 >= r) OR (r1 >= r AND r2 <= r). Then the segment
is NEITHER.
C.) If both r1 > r AND r2 > r (DEFAULT)
1.) If p1 = p2 then OUTSIDE.
2.) p1 != p2 (DEFAULT)
a.) If |d| > r then OUTSIDE.
b.) If |d| <= r (DEFAULT) then
i.) If (p1 - w) . (p2 - w) <= 0 then NEITHER
ii.) If (p1 - w) . (p2 - w) > 0 (DEFAULT) then OUTSIDE.
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
Here is a second version if you want to also count the number of intersection
points. It is a bit more involved. This program is given second.
A.) If both r1 < r AND r2 < r. Then the segment is INSIDE and 0 int. pt.
B.) If either (r1 < r AND r2 >= r) OR (r1 >= r AND r2 < r). Then the segment
is NEITHER inside or ourside and there is one intersection point.
C.) If both r1 = r and r2 = r then the segment is NETHER in or out
1.) If p1 = p2 then there is 1 intersection point
2.) if (p1 != p2 (DEFAULT) there are two intersection points.
D.) If ( r1 > r AND r2 = r ) OR ( r1 = r AND r2 > r ) then segment is NEITHER
inside nor outside.
1.) If |r1^2 - r2^2| >= |p1 - p2|^2 then there is one intersection point.
2.) If |r1^2 - r2^2| < |p1 - p2|^2 (DEFAULT) then there are two intersection pts.
E.) If ( r1 > r AND r2 > r) (DEFAULT)
1.) If p1 = p2 then OUTSIDE and there are no intersection points.
2.) If p1 != p2 (DEFAULT) then
a.) If |d| > r OUTSIDE and there are no intersection points.
b.) If |d| = r then
i.) If (p1 - w) . (p2 - w) <= 0 NEITHER and there is one intersection.
ii.) If (p1 - w) . (p2 - w) > 0 (DEFAULT) OUTSIDE with no intersections.
c.) If |d| < r (DEFAULT) then
i.) If (p1 - w) . (p2 - w) <= 0 NEITHER and there are two intersections
ii.) If (p1 - w) . (p2 - w) > 0 (DEFAULT) OUTSIDE with no intersections.
Programming notes. The default cases correspond to else.
Because roundoff errors may cause two theoretically equal doubles to be unequal, we
always test equality up to a fudge as in fabs( a - b ) < TOL, where TOL is a small
number.
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
Here is third version of the program using another interpretation of INSIDE and
OUTSIDE. This time we analyze the relative positions from the point of view of the
line by considering the square of the distance function to the center restricted to
the line. This function is a quadratic function.
The first structure is a pair of doubles which represents a point.
The second structure is a pair of points which is re[presenmts a segment.
The third structure is a point and double that represents a circle.
For this program, INSIDE, OUTSIDE and NEITHER have a different meaning:
We shall say that the segment is outside if no point of the segment has a
distance to the center less than the radius. Thus an outside segment may touch
the circle at one point or at no points. Similarly a segment is inside the
circle if none of its points are at a distance to the center greater than the
radius. Thus an inside segment may touch the circle at no points, one endpoint
or two endpointspoints.
We first deal with the case the endpoints are equal. Then
|p1-center|^2 is either less than, equal to or greater than r^2,
corresponding to inside, on or outside the circle.
If there are two points then a paramentric version of the line can be given by
y(t) = p2 (1 + t)/2 + p1 (1 - t)/2
so that y(-1)=p1 and y(1) = p2 are the endpoints. Then the function
|y(t) - center|^2 - radius^2 = at^2 + bt + c has a>0 and is
is a quadratic function along the line. The function is zero at
the points where the circle intersects the extended line along the segment.
The case that the discriminant d = b^2-4ac is negative means no zeros: whole
extended line is outside the circle. If d=0 then there is one point of the
extended line on the circle. Then the program determines if the zero, tz is in
the segment |tz|<=1 or not. Finally, if d>0 then the extended line intersects
the circle at two points t1 and t2. If either t1 or t2 are in [-1,1] then
that's where the circle intersects the segment. If both points t1 and t2 are
on the same side outside the interval [-1 1], then the segment is outside the
circle. Otherwise if t1 and t2 are outside but on opposite sides, then the
segment is inside the circle. In the remaining case, the segment is neither
inside or outside the circle. The intersection points correspond to zeros and
can be reconstructed using y(t). We define a point valued function of two
points and two doubles that gives the linear combination of points.
Finally, there is a possibility that roundoff errors may mean that the lengths
are not exact, so that equality can fail. To avoid this, we consider a point
on the circle if its distance to the center is "almost" the radius. Here,
almost means |distance - radius| < TOL, where TOL is a preset tolerance
value, which I take to be 5e-7.
*************************************************************************************/
# include
# include
# include
# define TOL 0.5e-7
typedef struct
{
double x;
double y;
}
vector;
typedef struct
{
vector p1;
vector p2;
}
segment;
typedef struct
{
vector c;
double r;
}
circle;
vector
lincomb( vector, vector, double );
double
dot( vector, vector );
vector
unitperp ( vector );
vector
enterv( void );
void
printm( int, int );
int
main ( void )
{
segment s;
vector t, n, w;
circle q;
double rr, d, d1, d2;
printf ( "First endpoint of the segment,\n");
s.p1 = enterv();
printf ( "\nSecond endpoint of the segment\n");
s.p2 = enterv();
printf ( "\nCenter of the circle\n" );
q.c = enterv();
printf ( "\nEnter the radius of the circle\n");
scanf ( "%lf",&q.r);
rr = q.r * q.r;
d1 = dot( lincomb(s.p1, q.c, -1.0), lincomb(s.p1, q.c, -1.0) );
d2 = dot( lincomb(s.p2, q.c, -1.0), lincomb(s.p2, q.c, -1.0) );
if ( (d1 <= rr - TOL) && (d2 <= rr - TOL) )
printm ( 1 , 0);
else if( ( (d1 < rr+TOL) && (d2 > rr-TOL)) || ( (d1 > rr-TOL) && (d2 < rr+TOL)) )
printm ( 2 );
else
{
if ( fabs( s.p1.x - s.p2.x ) + fabs ( s.p1.y - s.p2.y ) < TOL )
printm( 3 );
else
{
t = lincomb( s.p2, s.p1, -1.0);
n = unitperp( t );
d = dot ( lincomb( s.p1, q.c, -1.0), n);
if ( d*d >= rr + TOL )
printm( 3 );
else
{
w = lincomb ( q.c, n, d);
if( dot( lincomb( s.p1, w, -1.0), lincomb( s.p2, w, -1.0) ) < TOL)
printm( 2 );
else
printm( 3 );
}
}
}
printf ( "Oops! An error has occurred.\n");
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
lincomb( vector u, vector v, double t )
{
vector w;
w.x = u.x + t * v.x;
w.y = u.y + t * v.y;
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
double
dot( vector u, vector v)
{
return u.x * v.x + u.y * v.y;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
unitperp ( vector u)
{
vector w;
double s;
s = sqrt( u.x * u.x + u.y * u.y );
w.x = - u.y / s;
w.y = u.x / s;
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
enterv( void )
{
vector w;
printf("Enter the x and y components :");
scanf( "%lf%lf", &w.x, &w.y );
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm( int k, int n )
{
char mes[4][30] = { "oops","INSIDE","NEITHER INSIDE NOR OUTSIDE","OUTSIDE" };
printf("\nThe segment is %s the circle.\n",
"The circle and segtment intersect in %d points.\n", mes[k],n);
exit (0);
}
/*************************************************************************************
Second version - - - with point counting
**************************************************************************************/
# include
# include
# include
# define TOL 0.5e-7
typedef struct
{
double x;
double y;
}
vector;
typedef struct
{
vector p1;
vector p2;
}
segment;
typedef struct
{
vector c;
double r;
}
circle;
vector
lincomb( vector, vector, double );
double
dot( vector, vector );
vector
unitperp ( vector );
vector
enterv( void );
void
printm( int, int );
int
main ( void )
{
segment s;
vector t, n, w;
circle q;
double rr, d, d1, d2, len2;
int k;
printf ( "First endpoint of the segment,\n");
s.p1 = enterv();
printf ( "\nSecond endpoint of the segment\n");
s.p2 = enterv();
printf ( "\nCenter of the circle\n" );
q.c = enterv();
printf ( "\nEnter the radius of the circle : ");
scanf ( "%lf",&q.r);
rr = q.r * q.r;
d1 = dot( lincomb(s.p1, q.c, -1.0), lincomb(s.p1, q.c, -1.0) );
d2 = dot( lincomb(s.p2, q.c, -1.0), lincomb(s.p2, q.c, -1.0) );
if ( (d1 <= rr - TOL) && (d2 <= rr - TOL) )
printm ( 1 , 0);
else if( ( (d1 <= rr-TOL) && (d2 > rr-TOL)) || ( (d1 > rr-TOL) && (d2 <= rr-TOL)) )
printm ( 2, 1 );
else if( ( fabs(d1 - rr ) < TOL ) && ( fabs(d2 - rr) < TOL ) )
{
if ( fabs( s.p1.x - s.p2.x ) + fabs ( s.p1.y - s.p2.y ) < TOL )
printm( 2, 1);
else
printm( 2, 2);
}
else if( ((d1 >= rr + TOL ) && ( fabs(d2 - rr) < TOL )) || (( fabs(d1 - rr) < TOL ) && (d2 >= rr + TOL )) )
{
len2 = dot( lincomb(s.p1, s.p2, -1.0), lincomb(s.p1, s.p2, -1.0) );
if ( fabs ( d1 - d2 ) > len2 - TOL )
printm( 2, 1);
else
printm( 2, 2);
}
else
{
if ( fabs( s.p1.x - s.p2.x ) + fabs ( s.p1.y - s.p2.y ) < TOL )
printm( 3, 0 );
else
{
t = lincomb( s.p2, s.p1, -1.0);
n = unitperp( t );
d = dot ( lincomb( s.p1, q.c, -1.0), n);
if ( d*d >= rr + TOL )
printm( 3, 0 );
else
{
if (fabs (d*d - rr) < TOL )
k = 1;
else
k = 2;
w = lincomb ( q.c, n, d);
if( dot( lincomb( s.p1, w, -1.0), lincomb( s.p2, w, -1.0) ) < TOL)
printm( 2 , k);
else
printm( 3 ,0);
}
}
}
printf ( "Oops! An error has occurred.\n");
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
lincomb( vector u, vector v, double t )
{
vector w;
w.x = u.x + t * v.x;
w.y = u.y + t * v.y;
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
double
dot( vector u, vector v)
{
return u.x * v.x + u.y * v.y;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
unitperp ( vector u)
{
vector w;
double s;
s = sqrt( u.x * u.x + u.y * u.y );
w.x = - u.y / s;
w.y = u.x / s;
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
vector
enterv( void )
{
vector w;
printf(" Enter the x and y components : ");
scanf( "%lf%lf", &w.x, &w.y );
return w;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm( int k, int n )
{
char mes[4][30] = { "oops","INSIDE","NEITHER INSIDE NOR OUTSIDE","OUTSIDE" };
printf("\nThe segment is %s the circle.\n"
"The circle and segment intersect in %d points.\n", mes[k],n);
exit (0);
}
/*************************************************************************************
Third version - - - analysis from the point of view of the line
**************************************************************************************/
# include
# include
# include
# include
# define TOL (0.5e-7)
typedef struct /* Structure consists of two doubles */
{ /* Every type "point" has two vars called .x and .y */
double x;
double y;
} point;
typedef struct /* Represents a circle with center and radius */
{
point center;
double radius;
} circle;
typedef struct /* Nested structures: */
{ /* consists of two "point" structures */
point p1;
point p2;
} segment;
double
lengthsq( point v);
void
printv( point t);
void
printmn( int i, int m );
point
linco( point v, double cv, point w, double cw );
int
main( void )
{
double r,r1,rho1,rho2,rho3,a,b,c,d,w1,w2,t1,t2,tz,t;
int n;
point middle;
segment seg;
circle cir;
printf ( "Circle and Segment \n\n"
"For this program, a segment is INSIDE if none of its points have a\n"
"distance to the center greater than the radius. A segment is OUTSIDE\n"
"if none of its points have a distance to the center less than the radius.\n"
"Thus inside or outside segments may touch the circle at a point.\n\n" );
printf ( "Enter the center of the circle x y : " );
scanf ( "%lf%lf", &cir.center.x, &cir.center.y );
printf ( "Enter the radius of the circle r :");
scanf("%lf",&cir.radius);
r = cir.radius*cir.radius;
printf ( "Enter the segment's starting point x y : " );
scanf ( "%lf%lf", &seg.p1.x, &seg.p1.y );
printf ( "Enter the segment's ending point x y : " );
scanf ( "%lf%lf", &seg.p2.x, &seg.p2.y );
if ( fabs(seg.p1.x - seg.p2.x)+fabs(seg.p1.y-seg.p2.y) < TOL ) /* degenerate case p1=p2 */
{
r1 = lengthsq(linco(seg.p1,1.0,cir.center,-1.0))-r;
if ( fabs(r1)< TOL)
{
printv( seg.p1);
printmn(1,1);
}
else if (r1>0)
{
printmn(2,0);
}
else
{
printmn(0,0);
}
}
else
{
middle = linco ( seg.p1, 0.5, seg.p2, 0.5);
rho1 = lengthsq ( linco (seg.p1, 1.0, cir.center,-1.0) );
rho2 = lengthsq ( linco (seg.p2, 1.0, cir.center,-1.0) );
rho3 = lengthsq ( linco (middle, 1.0, cir.center,-1.0) );
a = (rho1+rho2)/2.0 - rho3;
b = (rho2-rho1)/2.0;
c = rho3-r;
d = b*b-4.0*a*c;
tz = -0.5*b/a;
if( fabs(d) < TOL) /* case extended line along segment is tangent to the circle */
{
if ( fabs(tz)<1.0+ TOL )
{
printv ( linco ( seg.p2,(1.0+tz)/2.0,seg.p1,(1.0- tz)/2.0));
printmn(2,1);
}
else
{
printmn(2,0);
}
}
else if( d< 0.0) /* case extended line segment is outside the circle */
{
printmn(2,0);
}
else /* case extended segment passes through the circle */
{
t = sqrt(d)/(2.0*a);
t1 = tz-t; /* quadratic formula for zeros of quadratic */
t2 = tz+t;
n=0; /*n counts number of zeros (intersections with circle) in segment */
if( fabs(t1)<1.0+ TOL ) /* if first root is in segment */
{
n=1;
printv(linco ( seg.p2,(1.0+t1)/2.0,seg.p1,(1.0- t1)/2.0));
}
if( fabs(t2) < 1.0+TOL ) /* if second root is in segment */
{
printv(linco ( seg.p2,(1.0+t2)/2.0,seg.p1,(1.0- t2)/2.0));
n = n + 1;
}
if( (t1>1.0-TOL) && ( t2 >1.0-TOL) ) /* if both roots are on positive side of segment */
printmn(2,n);
else if( (t1< TOL - 1.0) && ( t2 < TOL - 1.0 ) ) /* else if both on neg side */
printmn(2,n);
else if( ( fabs(t1)> 1.0-TOL ) && (fabs(t2) > 1.0-TOL) ) /* else if one on each side */
printmn(0,n);
else /* remaining casi is that one of the roots is an interior point of segment */
printmn(3,n);
}
}
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* Example of a double valued function that takes a "struct coords' argument. */
double
lengthsq( point v)
{
return v.x * v.x + v.y * v.y ;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
point
linco ( point v, double cv, point w, double cw )
{
point temp;
temp.x = cv * v.x + cw * w.x;
temp.y = cv * v.y + cw * w.y;
return temp;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printv( point t)
{
printf(" A point of the intersection is ( %f, %f )\n",t.x,t.y);
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printmn( int i, int m )
{
static char *mess[4]={"inside", "on","outside","neither inside nor outside"};
printf(" The segment is %s the circle.\n"
" The number of intersection points is %d \n\n", mess[i], m);
return;
}