Title: Student Code Critique Competition 12

Let us firstly consider the functionality of the code before we address any imperfections in semantics or design. The program begins by getting an initial approximation for the root from the user (I will later mention some other improvements which might be employed). The program then enters the main calculation loop. Here, the current best guess for the root location is used to calculate a new value for x, which is, at least in most cases, a closer approximation to the root sought. When the best guess is precise enough, in this case within 0.000001 of the root, checked by comparing the deviate of the new and previous best guess, the estimation is printed and the program terminates.

The bulk of the work, however, is performed within three auxiliary functions. Firstly, NewtonRough() calculates a new best guess for the root, based on a previous estimate, x0, passed as an argument. It does this by calling function func() to calculate the value of f(x) for x = x0, and derivative() to calculate the value of the derivative for which x = x0. The quotient of the function and derivative (f(x)/f'(x)) is subtracted from the previous estimate x0 to evaluate the new estimate, which is returned to main(). If we go further and examine func(), when we put the various components together we find the following function is evaluated:

f(x) = 25/7 * ln(3750 * (|x| ^ 0.5)) - 1/(|x| ^ 0.5) - 6.

The same process for derivative() shows that

f'(x) = 1/2x * (25/7 + 1/(|x| ^ 0.5)).

To improve upon the code as set, I will present two programs; the first will correct some of the immediate problems with the code, and the second will show a more robust and flexible version to give a taste of what might be achieved with a little further work. I will also state up front, that the code is (I hope) compliant with C89, but doesn't attempt to use any new C99 features. I only have a passing knowledge of C99 having not worked with C in anger for a couple of years.

I have not, by and large, commented on spacing issues as it seems clear that this is mostly to fit the code onto the page, and in any case is largely a matter of preference. I will mention, however, a couple of details where I am not certain if it is just publication issue. So onto the code and some specifics, the first program is given in newton1.h:

#ifndef NEWTON1_H
#define NEWTON1_H
#define PRECISION 1E-6
#define MAX_TRIES 30
#define INPUT_BUFFER_LEN 50
#define CONSTANT1 25.0 / 7.0
#define CONSTANT2 3750.0
#define CONSTANT3 6.0
#define CONSTANT4 2.0
double Func(const double current_gues );
double Derivative(const double current_guess);
double Reciprocal(const double val_to_recip);
void NewtonMethod(double current_guess );
double GetNumber(void );
#endif

and newton1.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include "newton1.h"
int main(void){
  double current_guess = 0;
  printf("Enter an estimate of the root\n> ");
  current_guess = GetDouble();
  NewtonMethod( current_guess );
  return (EXIT_SUCCESS);
}
void NewtonMethod(double current_guess){
  int tries = 1;
  double next_guess = 0;
  for(tries = 1; tries <= MAX_TRIES; tries++){
    errno = 0;
    next_guess = current_guess - (Func(current_guess)/Derivative(current_guess));
    if(errno == EDOM) {
      perror("Domain error encountered calculating next root estimate\n");
      break;
    } else if(errno == ERANGE) {
      perror("Range error encountered calculating next root estimate\n");
      break;
    } else if (errno != 0 ) {
      perror("Unexpected error encountered calculating next root estimate\n");
      break;
    }
    if(fabs(next_guess - current_guess) 
                  < PRECISION) {
      printf("Estimation of root %f\n", next_guess);
      break;
    } else if( tries == MAX_TRIES ) {
      printf("Failed to reach convergence after %d attempts\n", MAX_TRIES);
      break;
    }
    current_guess = next_guess;
  }
}
double Func(const double current_guess){
  double positive_sq_root = 
              sqrt(fabs(current_guess));
  double ln_eval = log(CONSTANT2 
                  * positive_sq_root);
  return ((CONSTANT1 * ln_eval ) -
  (Reciprocal(positive_sq_root)) - CONSTANT3);
}
double Derivative(const double current_guess){
  double factor = Reciprocal(CONSTANT4 
                    * current_guess);
  double reciprocal_root = 
    Reciprocal(sqrt(fabs(current_guess)));
  return((factor * CONSTANT1) + reciprocal_root);
}
double Reciprocal(const double val_to_recip){
  if(val_to_recip == 0) {
    errno = EDOM;
    return ( HUGE_VAL );
  }
  return (1.0/val_to_recip);
}
double GetDouble(void){
  char number_buf[INPUT_BUFFER_LEN];
  char * end_ptr = number_buf;
  double user_value = 0;
  while (1) {
    fgets(number_buf, INPUT_BUFFER_LEN, stdin);
    errno = 0;
    user_value = strtod(number_buf, &end_ptr);
    if(user_value == 0 && end_ptr == number_buf){
      printf( "Invalid number format, 
non-numeric input or underflow\n> " );
    } else if ( errno == ERANGE ) {
      printf("Overflow in number parsing\n> ");
    } else break;
  }
  return (user_value);
}

Let us consider the overall structure of the code first of all. The main changes for the first program are just to farm out some of the work to extra functions. Most significantly I have moved the main calculation loop from main() to NewtonMethod(). Personally I tend not to have anything but the most basic of startup code in main(), and in the original code, NewtonRough() was hardly earning its keep. The other new functions are Reciprocal() and GetNumber(), which are discussed further below.

There are also a couple of minor details that are related to my personal preference which are worth noting. I have capitalised the first letter of each function (except, of course, main()) for consistency. As it happens, I tend to capitalise function names, but whichever way you prefer to name functions, I feel it ought to be consistent. I have also made all the function returns consistent. Again, this is a matter of preference; I always put the return in parentheses, but put a space between the return keyword and the first parenthesis to indicate it is not a function. This is one of those places where I am not sure if the inconsistency in the original code was a space issue, but I felt it worth flagging.

So onto the more significant details, roughly in the order they appear as you follow the program execution. The first issue is in the definition of main() itself. To my knowledge, unless it is different in C99 (with which I am as yet not totally familiar), in ISO C, main should be defined as one of:

int main(void)

or

int main(int argc, char *argv[])

of which the first form is preferred here as we are not using command line arguments.

I have added declarations for the user-defined functions as part of a header file. These are not required, and the code works fine without them, with the compiler treating the definitions as prototypes. That said, once larger software is considered, it is wise to have a declaration for every function (even if you do not intend to call the function outside of its translation unit in my opinion), so this is a good habit to develop. I like to give the compiler every chance to spot my mistakes.

As well as function declarations, I have declared the parameters as const where they are not modified. Once again, in a small piece of code, this makes little difference, but on a larger scale, it enables you to see the intent of the designer, and gives the compiler a chance to stop you doing something in error.

Throughout the code, I have renamed many of the variables to make them somewhat more meaningful. This adds clarity to the designer's intentions and makes the life of future maintainers (the original designer included) much simpler. Another potential benefit of using better names is revealed by a couple of typographical errors in the code. Those who only glanced over it, may have missed the upper case identifiers used by mistake in a few places (I wonder if these were deliberate or accidental on Francis' part?). To my mind, it is less easy to make this kind of mistake with more expressive names, as you tend to think about them more as you type, and once you make the mistake, it is visually more obvious.

The user input facility has been removed into a separate function, and made more robust. Trying to make decent input code with scanf() is notoriously difficult, and I guess many would say, impossible. The first issue with the scanf() statement as published, is that %f is the wrong conversion specifier for type double. Unintuitively, while %f is correct for double in printf(), %lf is used in scanf(). This is just one of the problems with scanf(). The second problem which might be encountered here is when some kind of sanity check is placed on the input as I have done. Unfortunately scanf() doesn't remove unexpected non-numeric data from the input stream, so any attempt to call scanf() again if an error is detected results in scanf() picking up the same bad character. Anyway, if anyone wants to know more about the problems, I am sure a browse through the accu-general archives will probably help, or those of comp.lang.c if it doesn't.

Furthermore, by rewriting the user input code to use strtod() from an input string, we can handle any numeric input (at least within the overflow and underflow boundaries for a double type), and easily implement a basic sanity check. The only drawback being, at least as far as I understand, that it is not possible to distinguish between an underflow and an invalid numeric format in a strictly ISO conforming library. Others of course, may, and probably will, know better?

Throughout the code, I have replaced 'magic' numbers with preprocessor macros. This is probably the most idiomatic way to deal with magic numbers in C, so is the method I have chosen. Not only is this more elegant and less error prone during later maintenance, but in some cases faster as it avoids a function call. In the case of the constants in the equation, the solution is still unsatisfying, however, as I couldn't think of great names for them. When we look at a more advanced version of the program, we will see a better solution. I have also factored out the taking of reciprocals into a function (it could equally well be a macro) so the multiple occurrences of the 'magic' number 1.0 are avoided.

Another problem can be identified from a consideration of Newton's method. It is not always foolproof in various circumstances, and no method is bulletproof in the face of a rootless function. With my maths being a little rusty I am not certain, but I believe the function implemented in the sample doesn't have any roots, only a singularity at x = 0. As a safety net, therefore, the calculation loop in main() is limited to a maximum number of tries, defined in this case to be 30. This is sufficient for most purposes as Newton's method converges rapidly in general.

Whilst the nature of the function under consideration removes the risk of encountering many of the possible errors in the standard C maths library, some can still occur. In this code, the log() function can produce a range error if the value is close to zero. By checking the external errno variable each loop cycle, we can handle this error, and also implement error checking in the Reciprocal() function.

A final change in main() is the return value to EXIT_SUCCESS. Whilst a return of 0 is perfectly valid ISO C, I prefer the standard macros, especially if the code is for the benefit of people relatively inexperienced. After all, how many people, inexperienced or otherwise, have occasionally parsed code such as the idiomatic (although perhaps unwise)

if(!(strcmp(s1, s2))) { ... }

as meaning 'if string s1 is not equal to string s2'? Although the return value of 0 for success in many situations makes sense to the experienced C coder, it is perhaps a little unintuitive for the unwary.

So this code deals with the basic aspects of the problem. We have tidied the code up so that it can be compiled, and is reasonably idiomatic standard C (well standard C89 anyway). The overall design is still very limited and lightweight, however.

For those who have managed to read this far, I have also written a somewhat larger and more flexible version of the software. This has slightly fewer rough edges, and can be adapted to wide range of functions with ease. It still maintains a pretty lightweight approach, but I hope is developed enough to be usable in real root finding problems. If a student really wanted to pursue the issue, they could add a fully-fledged parser, e.g. a recursive descent system, to the code I present below, but that is really beyond the scope of this article. The modified code is presented in newton2.c and newton2.h:

(For reasons of space the code has not been printed but is available on our web site, www.accu.org. - ed.)

All I have done there is break the concept of a mathematical function down into components separated by + and -, as one tends to while differentiating by hand. I have illustrated the mechanism with a simple quadratic equation, f(x) = 6x^2 + 9x - 6, which has two roots at -2 and -0.5, rather than the more complex function used before. I have created three calculation functions for the quadratic, linear and constant components of this equation, but it is a simple task to add further functions for other standard differentials. Each calculation function needs to be able to handle the calculation for both f(x) and f'(x) based on an enumerated mode parameter passed as an argument. With a library of standard calculation functions in place, to work with different equations all the user has to do is alter the static array of structures in CalculateBetterApprox() to use the right constants and calculation function. The op parameter of the structure determines whether the component is to be added to or subtracted from the evaluation, or in the case of op_end, indicates that the end of the equation has been reached. To create functions to handle more complex standard differentials we might need to add further constants to the structure, but again, this is a simple task.

This still isn't the most flexible system, but represents a significant improvement on the single function program originally outlined. This code can at least be used reasonably successfully in a real world situation, and could be extended to be friendlier if required. For example, it would be relatively trivial to add functionality to read the equation details from a file into a dynamically allocated list rather than having to recompile each time.

Finally we can briefly mention a few other additions that could be used to make the system more flexible. Adding a search facility to locate an initial estimate for a root rather than relying on the user might be useful. Other improvements might revolve around combining Newton's method with other root finding algorithms to improve global convergence and handle cases where one particular method runs into difficulties.

func(): Returns the value of a function for argument x. The function is

(25/7)(ln(3750*sqrt(abs(x)) - 1/sqrt(abs(x)) - 6

where abs is the absolute value function. The function can be simplified slightly using the properties of logarithms to separate the constant factor. The function is symmetric about the "y-axis" and has a singularity at x = 0

derivative(): Returns the value of the derivative (with respect to x) of func() at x.

The (derivative) function is

(1/(2x))*( (25/7) + 1/(sqrt(abs(x))) )

Note that the factor of x outside the sum correctly allows for the change in sign of the derivative with the sign of x. Note also that the value of the derivative is undefined for x = 0 (that is the limit from both sides does not exist). If the derivative function is called with a value of zero, then it will fail with a divide by zero error.

NewtonRough(): Returns one iteration of the Newton-Raphson algorithm. The choice of name is poor. Something like doNewtonRaphsonIteration() would have been more descriptive.

The Newton-Raphson algorithm is not guaranteed to converge to a steady value. The convergence depends critically on the nature of the function and the initial starting point. In some pathological cases NR may appear to converge but will not converge to a genuine zero of the function. The student should not be surprised by this behaviour, but all too frequently such routines are used without the users being aware of the restrictions. When it works the Newton-Raphson method is quadratically convergent, but when it doesn't work it can rapidly diverge or enter a loop of values. so, the answer to the direct question is that the student's understanding of numerical analysis and algorithms is lacking. The student should be made aware of other root-finding methods to enrich his/her armoury of techniques.

No iterative algorithm should have a potentially infinite loop. If it hasn't converged after a reasonable number of iterations then a new starting point or algorithm needs to be considered. If the number of iterations is set too high then consideration will have to be given to numerical stability and accumulated rounding errors.

There has been much discussion in C Vu in the past about how convergence tests should be done. In this instance the absolute difference in the value of x0 returned by two consecutive iterations is a reasonable measure since graphing the function shows that the zeros lie at around (+/-) 0.005.The maximum number of iterations and the convergence criterion need to be made more visible and stop being magic numbers. The lack of function prototypes forces an order of function definition that is restrictive for future devel

Notes: 

More fields may be available via dynamicdata ..