/***************************************************************************** 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; }