## CS13002 PDS Lab, Spring 2003, Section 1/ASolution of Assignment 3

1. [Linear sequence]

This is straightforward as the following snippet shows:
```#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 us
```
The 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!

2. [Six points inside a circle]

Suppose we want to compute the circum-radius of the three points P1=(x1,y1), P2=(x2,y2), P3=(x3,y3). The perpendicular bisector of the line segment P1P2 contains precisely those points that are equidistant from P1 and P2 and hence has the equation
```(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]