ECE 8473 - Assignment #3 - Due: 6 October 2022

Programs must compile with no warnings or errors using:

  gcc -std=c11 -pedantic -Wall

1. [60] a3/p1.c - testing math.h functions

   Each student has been randomly assigned 6 of the functions from math.h in the form:

     double func( double)
   

   The assigned functions are listed in a comment at the top of your initial a3p1.c file in the Assignment Results.

   For each function, write your own version of the function, named differently from the original, e.g.

     double my_func( double)
   

   so that you can invoke both your version of the function and the original for comparison.

   Implement your version of the function by using other functions from math.h, do not use any loops or series expansions. For example, sqrt(x) can be implemented using pow(x,0.5):

     double my_sqrt( double x)
     {
       return pow( x, 0.5);
     }
     ...
     double x, y1, y2;
   
     x = 2.0;
     y1 = sqrt(x);
     y2 = my_sqrt(x);
     ...
   

   If you need the value of π use 4*atan(1), and for e use exp(1).

   Compare your function result values with those from the standard library at two non-trivial points, i.e. points where the input and/or output are not just 0.0 or 1.0. Try to use one point where you think the function results should agree closely and one where you think they may differ more significantly.

   Display the function names, inputs, results, and differences in results. A format like "%.18g" may be suitable to see all of the significant digits.

   ---

2. [20] a3/p2.c - Newton square root

   Write a function in the form:

     double nsqrt( double a)
   

   which performs 5 iterations of the Newton algorithm for square root, displaying the result after each iteration (%.18f format), and returning the result of the final iteration.

   The details of the Newton algorithm are shown below.

   In your main() program, compare results of nsqrt(a) with sqrt(a) for a = 2 and a = 10, and display the difference (%.18f format).

   ---

3. [20] a3/p3.c - Goldschmidt square root

   Write a function in the form:

     double gsqrt( double a)
   

   which performs 5 iterations of the Goldschmidt algorithm for square root, displaying the result after each iteration (%.18f format), and returning the result of the final iteration.

   The details of the Goldschmidt algorithm are shown below.

   In your main() program, compare results of gsqrt(a) with sqrt(a) for a = 2 and a = 10, and display the difference (%.18f format).

   Example output for a = 2

Newton said:

  f(x) ~= f(x0) + (x-x0)*f'(x0) = 0

so:

  x = x0 - f(x0)/f'(x0)

For sqrt(a):

  f(x) = x2 - a, and f'(x) = 2*x

so:

  x = x0 - (x02 - a)/(2*x0)

    = (x0 + a/x0) / 2

Starting with an initial estimate, e.g. x = a/2, each iteration of the Newton algorithm produces a better estimate using:

  x = (x + a/x) / 2

Goldschmidt said:

  sqrt(a) ~= x0 = a*k0,  if  k0 ~= 1/sqrt(a)

Define the error e as:

  e = 1 - a*k0*k0 

then:

  a*k0*k0 = 1 - e 

To reduce the error, multiply both sides by (1 + e)

  a*k0*k0*(1 + e) = 1 - e2 

and let k1*k1 = (1 + e)

  a*k0*k0*k1*k1 = a*(k0*k1)*(k0*k1) = 1 - e2 

so:

  x = a*(k0*k1)

is better than x0 as an estimate of sqrt(a).

For k1 we use the approximation:

  k1 = sqrt(1 + e) ~= (1 + e/2)

The Goldschmidt algorithm is preferred over Newton's algorithm for hardware implementations because it does not use division, it only uses multiplication and addition (the apparent division e/2 can be implemented as a bit-shift right in hardware). In fact, Goldschmidt's basic algorithm is for division, and it accomplishes division using only multiplication and addition. The Goldschmidt sqrt algorithm is just an extension of the Goldschmidt algorithm for division.

Starting with an initial estimate, e.g. k = 1/2, x = a*k, each iteration of the Goldschmidt algorithm produces a better estimate using:

  e = 1 - a*k*k
  k = k*(1 + e/2)
  x = a*k