1. Assignment One: Evaluation Efficiency

The code to evaluate the geometric sum in two ways is given by two Maple procedures

>	S1 := proc(x,n)  description `sums consecutive powers of x, ending with x^n`:  local s,k,y:  s := 1+x:  y := x:  for k from 2 to n do    y := x*y:    s := s+y:  end do:  return s: end proc:

The loop in the procedure S1 is executed n-1 times.  Every step in the loop has one addition and one multiplication.  Notice the one addition before the loop.  Therefore, the number of additions is n and the number of multiplications is n-1 in the execution of S1.

>	S2 := proc(x,n)  description `explicit formula for the geometric sum`:  local s,k,y,xm1:  y := x:  for k from 1 to n do    y := x*y:  end do:  xm1 := x-1:  s := y/xm1 - 1/xm1:  return s: end proc:

We see n multiplications in the loop to compute x^(n+1).  There are two subtractions and two divisions.  Thus the cost to execute S2 is n multiplications, two divisions, and two subtractions.

Let us first see whether the two procedures compute the same thing at a random sample:

>	S1(2312,23);

>	S2(2312,23);

For exact data (except for x = 1), as expected, we see that S1 and S2 are equivalent.

Comparing the execution efficiency of S1 and S2, let us have Maple summarize our earlier observations:

>	codegen[cost](S1);

>	codegen[cost](S2);

Looking only at the arithmetical operations (ignoring assignments and storage), and taking into account that n is typically large, S1 has n-2 additions more than S2.  The cost of these n-2 additions will quickly outgrow the constant cost of executing the one multiplication and two divisions which S2 has to execute more.

Thus, S2 is more efficient than S1.

2. Assignment Two: Evaluation at Random Points

We will evaluate S1 and S2 for n = 20, at random numbers, uniformly taken between 0 and 1, for increasing values of Digits, each time computing the relative error between the values returned by S1 and S2.

>

for k from 2 to 16 do  Digits := k;  r := stats[random,uniform[0,1]](1): rr := convert(r,rational):  y := S1(rr,20): y1 := S1(r,20): y2 := S2(r,20):  err_pls1 := abs(y1-y)/abs(y); err_pls2 := abs(y2-y)/abs(y):  y := S1(-rr,20): y1 := S1(-r,20): y2 := S2(-r,20):  err_min1 := abs(y1-y)/abs(y): err_min2 := abs(y2-y)/abs(y):

>	printf("Digits = %2d: er1+ = %.2e, er2+ = %.2e, er1- = %.2e, er2- =%.2e\n",k,err_pls1,err_pls2,err_min1,err_min2); end do:

Digits =  2: er1+ = 0.00e+00, er2+ = 5.70e-02, er1- = 0.00e+00, er2- =1.40e-02
Digits =  3: er1+ = 6.80e-03, er2+ = 0.00e+00, er1- = 0.00e+00, er2- =0.00e+00
Digits =  4: er1+ = 6.56e-04, er2+ = 6.56e-04, er1- = 1.34e-04, er2- =2.69e-04
Digits =  5: er1+ = 1.05e-04, er2+ = 5.26e-05, er1- = 1.47e-05, er2- =4.42e-05
Digits =  6: er1+ = 4.42e-06, er2+ = 0.00e+00, er1- = 1.56e-06, er2- =0.00e+00
Digits =  7: er1+ = 2.54e-06, er2+ = 2.28e-06, er1- = 3.49e-07, er2- =3.49e-07
Digits =  8: er1+ = 0.00e+00, er2+ = 0.00e+00, er1- = 1.03e-08, er2- =2.06e-08
Digits =  9: er1+ = 8.32e-09, er2+ = 5.55e-09, er1- = 1.72e-09, er2- =0.00e+00
Digits = 10: er1+ = 0.00e+00, er2+ = 3.96e-10, er1- = 1.60e-10, er2- =0.00e+00
Digits = 11: er1+ = 2.55e-11, er2+ = 0.00e+00, er1- = 1.74e-11, er2- =1.74e-11
Digits = 12: er1+ = 7.40e-12, er2+ = 7.40e-12, er1- = 2.52e-12, er2- =6.30e-12
Digits = 13: er1+ = 6.90e-13, er2+ = 6.90e-13, er1- = 0.00e+00, er2- =1.31e-13
Digits = 14: er1+ = 3.68e-13, er2+ = 3.27e-13, er1- = 0.00e+00, er2- =1.78e-14
Digits = 15: er1+ = 0.00e+00, er2+ = 0.00e+00, er1- = 1.04e-15, er2- =0.00e+00
Digits = 16: er1+ = 9.12e-16, er2+ = 0.00e+00, er1- = 2.18e-16, er2- =0.00e+00

We see that the relative error is bounded by the working precision 10^(-Digits), for both positive and negative values.

While we may expect to see more roundoff with negative values, this trend is not obvious.

For random numbers in [-1,+1], both formulas give the same accuracy.

Let us now see what happens if we let the degree grow larger:

>	Digits := 4:

>	r1 := 0.9511; r2 := 0.01342;

If we give S1 a rational number on input, then Maple will use exact arithmetic to execute the sum.

>	exact_r1 := convert(r1,rational); exact_r2 := convert(r2,rational);

>	abs(S1(r1,2000) - S1(exact_r1,2000))/abs(S1(r1,2000));

>	abs(S1(r2,2000) - S1(exact_r2,2000))/abs(S1(r2,2000));

The relative error tends to grow as x is closer to 1.  Let us compare with S2.

>	abs(S2(r1,2000) - S2(exact_r1,2000))/abs(S1(r1,2000));

>	abs(S2(r2,2000) - S2(exact_r2,2000))/abs(S2(r2,2000));

For larger values of n, we see that S2 is more accurate than S1.

3. Assignment Three: Evaluation around One

>	Digits := 4:

>	one := 0.99999:

>	S1(one,11) = S2(one,11);

>	Digits := 5:

>	S1(one,11) = S2(one,11);

S2 fails dramatically when x comes too close to 1, "close"  is relative to the working precision:

>	printf("distance of one to 1 = %e",one - 1);

distance of one to 1 = -1.000000e-05

If Digits = 5, our working precision is 10^-5, and we just have enough precision to capture the difference between "one" and 1.  If we use fewer than five Digits, the x - 1 in S2 evaluates to zero and the division by zero leads to an undefined float.  In general, S2 will always fail when x-1 becomes smaller than the working precision.

Conclusion

The symbolic formula S2 requires less operations and is as accurate as the geometric sum S1 for random numbers, uniformly distributed in [-1,+1-eps[.  As the degree gets larger, we observed a higher accuracy of the result returned by S2, in particular as x gets larger, and takes values around 0.9.  S2 cannot be used when the difference of x to 1 is less than the machine precision.  In case x gets too close to one, either S2 must be executed with sufficiently high precision, or S1 must be invoked.