#include <stdio.h>
#include <time.h>
/* The recursive procedure */
int recursive ( int n )
{
if (n == 1) return (1);
if (n == 2) return (1);
return (6*recursive(n-2)-recursive(n-1));
}
/* The iterative procedure */
int iterative ( int n )
{
int i, a, prev, temp;
if (n == 1) return (1);
if (n == 2) return (1);
a = prev = 1;
for (i=3; i<=n; i++) {
temp = a;
a = 6 * prev - a;
prev = temp;
}
return(a);
}
/* The formula */
int formula ( int n )
{
int i, twoPower, threePower;
if (n == 1) return (1);
twoPower = 4;
threePower = 1;
for (i=2; i<=n; i++) {
twoPower *= 2;
threePower *= -3;
}
return((twoPower + threePower) / 5);
}
int main ()
{
clock_t t1, t2;
int a20;
int i, repcntrec = 10000, repcntiter = 10000000, repcntform = 10000000;
printf("Recursive procedure : \n");
t1 = clock();
for (i=0; i<repcntrec; i++) a20 = recursive(20);
t2 = clock();
printf("a_20 = %d.\n",a20);
printf("Time taken = %f us\n",
1000000.*(float)(t2-t1)/(float)CLOCKS_PER_SEC/(float)repcntrec);
printf("\n");
printf("Iterative procedure : \n");
t1 = clock();
for (i=0; i<repcntiter; i++) a20 = iterative(20);
t2 = clock();
printf("a_20 = %d.\n",a20);
printf("Time taken = %f us\n",
1000000.*(float)(t2-t1)/(float)CLOCKS_PER_SEC/(float)repcntiter);
printf("\n");
printf("Formula : \n");
t1 = clock();
for (i=0; i<repcntform; i++) a20 = formula(20);
t2 = clock();
printf("a_20 = %d.\n",a20);
printf("Time taken = %f us\n",
1000000.*(float)(t2-t1)/(float)CLOCKS_PER_SEC/(float)repcntform);
}
The above code reports times in micro-seconds (us). Typical timing
figures on my machine (1.7GHz Pentium processor) are given below:
Recursive procedure : a_20 = -232032863. Time taken = 183.000000 us Iterative procedure : a_20 = -232032863. Time taken = 0.145000 us Formula : a_20 = -232032863. Time taken = 0.215000 usThe whole point in measuring the timings is to highlight the dangers associated with indiscriminate uses of recursion. Recursion has additional overheads. Moreover a blind use of recursion as above makes many repeated calls (call of recursive on the same argument value). The motto is: Avoid recursion as far as possible. There are specialized situations, however, where using recursion simplifies coding massively. If you encounter such a situation and are sure that you are not making multiple calls on the same argument values, go ahead with recursion!
(x-x1)2 + (y-y1)2 = (x-x2)2 + (y-y2)2.Cancelling the quadratic terms from the two sides of this equation gives a linear equation of the form
a1x + b1y = c1.Similarly the perpendicular bisector of the segment P1P3 is an equation of the form
a2x + b2y = c2.Now the three points P1,P2,P3 are collinear (or repeated) if and only if a1b2-a2b1 equals zero. If not, we solve these two linear equations to get the circum-center (u,v) of the three points P1,P2,P3 and the corresponding circum-radius is the distance between (u,v) and P1.
The routine circumRadius below is based on the above facts.
#include <stdio.h>
#include <math.h>
#define NPTS 6
double circumRadius ( double x1 , double y1 , double x2 , double y2 , double x3 , double y3 )
{
double a1, b1, c1, a2, b2, c2;
double u, v, w;
a1 = 2 * (x1 - x2);
b1 = 2 * (y1 - y2);
c1 = x1 * x1 + y1 * y1 - x2 * x2 - y2 * y2;
a2 = 2 * (x1 - x3);
b2 = 2 * (y1 - y3);
c2 = x1 * x1 + y1 * y1 - x3 * x3 - y3 * y3;
w = a1 * b2 - a2 * b1;
if (w == 0) return (-1); /* The three points are collinear */
u = ( c1 * b2 - c2 * b1 ) / w;
v = ( a1 * c2 - a2 * c1 ) / w;
return (sqrt((u-x1)*(u-x1)+(v-y1)*(v-y1)));
}
int main ()
{
int i, j, k;
double x[NPTS], y[NPTS], crad;
for (i=0; i<NPTS; i++) {
printf("Input the coordinates of Point %d : ",i);
scanf("%lf%lf",&x[i],&y[i]);
printf("Point %d : (%lf,%lf)\n",i,x[i],y[i]);
}
for (i=0; i<NPTS; i++) for (j=i+1; j<NPTS; j++) for (k=j+1; k<NPTS; k++) {
crad = circumRadius(x[i],y[i],x[j],y[j],x[k],y[k]);
printf("Circum-radius of points %d, %d and %d is %lf ", i, j, k, crad);
if (crad < 0) printf("(COLLINEAR POINTS)");
else if (crad > 1) printf("(BIG CIRCUM-RADIUS)");
printf("\n");
}
}
A typical run of the program is shown below.
Input the coordinates of Point 0 : 0 0 Point 0 : (0.000000,0.000000) Input the coordinates of Point 1 : -0.3 -0.4 Point 1 : (-0.300000,-0.400000) Input the coordinates of Point 2 : 0.3 -0.4 Point 2 : (0.300000,-0.400000) Input the coordinates of Point 3 : -0.3 0.4 Point 3 : (-0.300000,0.400000) Input the coordinates of Point 4 : -0.6 0.8 Point 4 : (-0.600000,0.800000) Input the coordinates of Point 5 : -0.6 -0.8 Point 5 : (-0.600000,-0.800000) Circum-radius of points 0, 1 and 2 is 0.312500 Circum-radius of points 0, 1 and 3 is 0.416667 Circum-radius of points 0, 1 and 4 is 0.644235 Circum-radius of points 0, 1 and 5 is -1.000000 (COLLINEAR POINTS) Circum-radius of points 0, 2 and 3 is -1.000000 (COLLINEAR POINTS) Circum-radius of points 0, 2 and 4 is -1.000000 (COLLINEAR POINTS) Circum-radius of points 0, 2 and 5 is 0.512961 Circum-radius of points 0, 3 and 4 is -1.000000 (COLLINEAR POINTS) Circum-radius of points 0, 3 and 5 is 0.644235 Circum-radius of points 0, 4 and 5 is 0.833333 Circum-radius of points 1, 2 and 3 is 0.500000 Circum-radius of points 1, 2 and 4 is 0.773082 Circum-radius of points 1, 2 and 5 is 0.615554 Circum-radius of points 1, 3 and 4 is 1.030776 (BIG CIRCUM-RADIUS) Circum-radius of points 1, 3 and 5 is 1.030776 (BIG CIRCUM-RADIUS) Circum-radius of points 1, 4 and 5 is 1.030776 (BIG CIRCUM-RADIUS) Circum-radius of points 2, 3 and 4 is 3377699720527872.500000 (BIG CIRCUM-RADIUS) Circum-radius of points 2, 3 and 5 is 0.634498 Circum-radius of points 2, 4 and 5 is 0.820738 Circum-radius of points 3, 4 and 5 is 1.030776 (BIG CIRCUM-RADIUS)The input/output behavior looks normal and one can check for the correctness of the values computed by the code. The only problem is with that of the triple (2,3,4) of input points. These three points are actually collinear, but the program does not report that, it issues an enormously big circum-radius instead. Where is the problem?
The problem stems from floating point approximations. Look at Point 2. Its x coordinate is stored in the compiler's double format. The value stored is very close to 0.3 but not exactly equal to it. Similar is the case for (most of) the other coordinates. These representation errors accumulate in a non-zero (albeit small) value of w in the function circumRadius. So instead of returning -1 it goes on computing u and v and subsequently the circum-radius.
What is the remedy then? Frankly speaking, there is no remedy at all. One possibility is to treat very small values of w as zero. But then three points that are actually not collinear may be certified so by this strategy.
The errors arising out of approximate (fixed-precision) representations of real numbers and the propagations of these errors pose severe difficulties in numerical computations. Using clever tricks one can sometimes minimize their effects, but one cannot eliminate them altogether.
When you are working with integers, use int, long or even long long variables. Resorting to double and/or deploying routines (like pow) involving double computations can be perilous.
Take care!
[Course home] [Home]