# Newton Raphson and polynomials in C

I have the following code, that defines:

1. A polynomial struct with some useful functions.
2. The newton Raphson algorithm for polynomials.

and calculates sqrt(2). What can I improve?

There are some things I am unsure about myself:

1. Is size_t a good datatype for my purposes? I needed to be very careful when writing the loop condition in the eval function. Is this a good place to switch to signed arithmetics?

2. Is my definition of my_nan a good way to ensure portability when switching to another floating point type?

Things I am aware of: I know, that the Newton Raphson for polynomials does not require the explicit construction of a derived polynomial and I am a bit wasteful on the memory.

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>

typedef double real;
const real my_nan = ((real) 0.0) / ((real) 0.0);

typedef struct Polynomial {
real* coeffs;
size_t capacity;
} Polynomial;

Polynomial* create_poly(size_t capacity)
{
real* coeffs;

coeffs = malloc(capacity * sizeof(real));
if (coeffs == NULL) return NULL;

Polynomial* poly = malloc(sizeof(Polynomial));

if (poly == NULL) {
free(coeffs);
return NULL;
}
poly->coeffs = coeffs;
poly->capacity = capacity;

return poly;
}

void delete_poly(Polynomial* p)
{
free(p->coeffs);
free(p);
p = NULL;
}

size_t deg(const Polynomial*const p)
{
for (size_t i = p->capacity - 1; i > 0; i--) {
/* Here we actually want to compare reals exactly instead of |a - b| < eps */
if (p->coeffs[i] != 0.0) return i;
}
return 0;
}

void print(const Polynomial*const p)
{
size_t i;
for (i = 0; i < deg(p); ++i) {
printf("%f * x^%zu + ", p->coeffs[i], i);
}
printf("%f * x^%zu\n", p->coeffs[i], i);
}

real eval(const Polynomial*const p, real x)
{
/* Use Horner Scheme for evaluation */
size_t i = deg(p);
real res = p->coeffs[i];

for (; i-- > 0;) {
res = res * x + p->coeffs[i];
}
return res;
}

Polynomial* derive(const Polynomial*const p)
{
Polynomial* Dp = create_poly(p->capacity);
if (Dp == NULL) return NULL;

for (size_t i = 1; i < p->capacity; ++i) {
Dp->coeffs[i - 1] = ((real) i) * p->coeffs[i];
}
return Dp;
}

real newton_raphson_poly(const Polynomial*const p, real x0, real eps)
{
Polynomial* Dp = derive(p);
real x, prev = x0;
const int max_iter = 100;

for (int i = 0; i < max_iter; ++i) {
x = prev - eval(p, prev) / eval(Dp, prev);
if (fabs(x - prev) < eps) {
return x;
} else {
prev = x;
}
}

return my_nan;
}

int main()
{
Polynomial* p;
const real EPS = pow(10, -7);
p = create_poly(3);

p->coeffs[0] = -2;
p->coeffs[1] = 0;
p->coeffs[2] = 1;

printf("The result of sqrt(2) is given by the root of\n");
print(p);
printf("Its value is: %f \n", newton_raphson_poly(p, 1.0, EPS));

delete_poly(p);
}
$$$$


## Overall Impression

The code displays some good programming practices or habits. It is generally well structured, with well thought out functions. Good definitions of typedefs. The type size_t is used to index arrays instead of int.

There are several potential bugs in the code. They are potential because they could happen, not that they would happen as the code is written currently.

## Error Checking

While the return values of malloc() are checked, there should be some additional error checking in 2 places in the code. The first is that in Polynomial* create_poly(size_t capacity) before any memory allocation takes place, the value of the parameter capacity should be checked that is larger than zero, if not the memory allocation should not take place.

The second place for additional error checking is in main(). If p after the call to create_poly() is NULL, assignments to p will fail, possibly causing a catastrophic error. This is one potential bug.

## Memory Allocation

It might be better to use calloc() rather than malloc() in the following statement:

    coeffs = malloc(capacity * sizeof(real));


The memory allocation function calloc(size_t count, size_t size) was written with arrays in mind. In addition to being slightly more readable, calloc() sets all of the values to zero when the array is allocated, which means the values in the array are initialized.

    real* coeffs = calloc(capacity, sizeof(*coeffs));


Note that in the example of calloc() usage the sizeof() argument is what coeffs points to. This allows who ever is maintaining the code to change the type of coeffs without have to modify more than the type itself. If the type real was used in the statement, there would be 2 places to change the code and not one.

## Unnecessary Statement

In the function void delete_poly(Polynomial* p) the statement p = NULL; is unnecessary. Since p was passed in instead of a pointer to p it only affects the local value of p, it does not affect the value of p in main().

## Initialize Variables When They are Declared

In at lest 2 places in the code variables are declared on one line and then initialized on another line like the initialization is an after thought. A better habit to get into is to initialize the variables as they are declared. This can lead to less bugs and less debugging of code.

In main():

    Polynomial* p;
const real EPS = pow(10, -7);
p = create_poly(3);


Versus

    const real EPS = pow(10, -7);
Polynomial* p = create_poly(3);


In create_poly()

    real* coeffs;

coeffs = malloc(capacity * sizeof(real));
if (coeffs == NULL) return NULL;


Versus

    real* coeffs = malloc(capacity * sizeof(real));
if (coeffs == NULL) return NULL;


## Magic Numbers

There are Magic Numbers in the main() function (10 and -7), it might be better to create symbolic constants for them to make the code more readable and easier to maintain. These numbers may be used in many places and being able to change them by editing only one line makes maintenance easier.

Numeric constants in code are sometimes referred to as Magic Numbers, because there is no obvious meaning for them. There is a discussion of this on stackoverflow.

The code already has a symbolic constant defined in real newton_raphson_poly(const Polynomial*const p, real x0, real eps).

    const int max_iter = 100;


Write the code consistently.

## Possible Optimizations

The code in main() could be more flexible or extensible if it modified to create an array of coefficients and get the size of the array as the capacity. The code to copy the coefficients into the Polynomial struct could be an additional function, or it could also be added to create_poly(size_t capacity, real coeffs[]).

int main()
{
real coeffs[] = {-2, 0, 1};
size_t capacity = sizeof(coeffs) / sizeof(*coeffs);

const real EPS = pow(10, -7);
Polynomial* p = create_poly(capacity);

real* poly_coeffs_ptr = &p->coeffs[0];
real* coeffs_ptr = &coeffs[0];
for (size_t i = 0; i < capacity; i++)
{
*poly_coeffs_ptr = *coeffs_ptr;
}

printf("The result of sqrt(2) is given by the root of\n");
print(p);
printf("Its value is: %f \n", newton_raphson_poly(p, 1.0, EPS));

delete_poly(p);
}


Reordering of malloc()s might lead to less code.

Polynomial* create_poly(size_t capacity)
{
if (capacity > 0)
{
Polynomial* poly = malloc(sizeof(Polynomial));
if (poly != NULL)
{
poly->coeffs = calloc(capacity, sizeof(real));
if (poly->coeffs == NULL)
{
free(poly);
poly = NULL;
}
}
return poly;
}

return NULL;
}

• Thank you very much for this thorough review. Regarding the inconsistent initialization in the declaration, I just realized an unconscious pattern in my code. I have to program a lot in Fortran and if you initialize something in the declaration it becomes a static variable (for whatever reason). So I had to hammer into my head to do this only for constants and especially not for local variables in a function. Apparently I took this over to C. Commented Apr 30, 2020 at 7:26

Weak EPS test

The following is only useful for differences on a narrow power-of-2 range. This approach is taking the float out of floating point.

if (fabs(x - prev) < eps)  // weak


When x, prev are small like 10e-10 the result is always true. Not useful to find the of sqrt(10e-20).

When x, prev are large like 10e+12 the result may never be true as the large difference wobbles around 0.0.

Instead with floating point, this nearness test needs to consider the magnitudes of x, prev.

Something like fabs(x - prev)/fabs(x + prev) < eps` can make sense - with additional code to protect against division by zero and overflow.

This is a deep subject and better tests exist. Usually what is best is depends on the situation.