Solution of Assignment 3

**[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!

**[Six points inside a circle]**

Suppose we want to compute the circum-radius of the three points`P`,_{1}=(x_{1},y_{1})`P`,_{2}=(x_{2},y_{2})`P`. The perpendicular bisector of the line segment P_{3}=(x_{3},y_{3})_{1}P_{2}contains precisely those points that are equidistant from P_{1}and P_{2}and hence has the equation(x-x

Cancelling the quadratic terms from the two sides of this equation gives a linear equation of the form_{1})^{2}+ (y-y_{1})^{2}= (x-x_{2})^{2}+ (y-y_{2})^{2}.a

Similarly the perpendicular bisector of the segment P_{1}x + b_{1}y = c_{1}._{1}P_{3}is an equation of the forma

Now the three points_{2}x + b_{2}y = c_{2}.`P`are collinear (or repeated) if and only if_{1},P_{2},P_{3}`a`equals zero. If not, we solve these two linear equations to get the circum-center_{1}b_{2}-a_{2}b_{1}`(u,v)`of the three points`P`and the corresponding circum-radius is the distance between_{1},P_{2},P_{3}`(u,v)`and`P`._{1}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]