Solving non-linear equation in C using Newton's method?

Question

I have a non-linear equation defined as p^(n-1) * [1-(1-p)^n] * r - c = 0, where n, r, and c are given integers. I want to solve for p in the C programming language and decided to use Newton's method. I am wondering if this is the most suitable approach, and if there are other ways to solve this. (One important note that p stands for probability here, so it ranges from 0 to 1!)

Here is my current implementation, and it seems to work for some test cases but fails for others:

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

#define ld long double
#define TOLERANCE 1e-7
#define ITERATIONS 1000

ld function(ld p, int n, int c, int r) {
    return pow(p, n - 1) * (1 - pow(1 - p, n)) * r - c;
}

ld derivative(ld p, int n, int r) {
    return r * (n - 1) * pow(p, n - 2) * (1 - pow(1 - p, n)) + r * n * pow(p, n - 1) * pow(1 - p, n - 1);
}

ld pFinder(int n, ld c, ld r) {
    ld p = 0.5;

    if (n == 0) {
        return 0;
    }

    for (int i = 0; i < ITERATIONS; i++) {
        ld fValue = function(p, n, c, r);
        ld fPrime = derivative(p, n, r);

        p = p - fValue / fPrime;

        if (fabs(fValue) < TOLERANCE) {
            return p;
        }
    }
    
    printf("Didn't converge after %d iterations", ITERATIONS);
    return -1;
}

int main() {
    int n, c, r;

    scanf("%d %d %d", &n, &c, &r);

    ld p = pFinder(n, c, r);

    printf("%.15Lf\n", p);
    return 0;
}

The issue arises with inputs like 18 30 2, where the result is 1.999993641483941, which is outside the valid range of [0, 1]. I do not understand why this is happening.

The constraints are 1 <= n <= 20, 1 <= c <= 10000, 1 <= r <= 10000, and the absolute error of p is up to 10^-6.

I appreciate any insights into what might be going wrong and how to improve my code!

P.S. This is my first time posting something here, and I am pretty new to coding, so I hope you don't mind that I'm not proficient in Markdown syntax!

Answer 1

For p to be in [0,1], c/r must also be in [0,1]. In that case, if you start with p = 1 (instead of p = 0.5), Newton's method will be well-behaved and will very quickly arrive at the desired precision. Much faster than bisection. The number of iterations will be small. For your tolerance, I can't get them to be more than 16. Usually just a handful.

In your example "10 30 2", c/r is 15, which is not in [0,1]. Check for c <= r before proceeding. (Your c >= 1 and r >= 1 already assures that c/r >= 0.)

Your use of pow() eliminates any precision advantage of using long doubles. Use powl() instead. Also fabsl().

To update, you only need two powl()'s, which are slow. I would:

#include <math.h>
#include <assert.h>

long double root(int n, int c, int r) {
    assert(n >= 1 && c >= 0 && r >= 1 && c <= r);
    int k = n - 1;
    long double x = c / (long double)r, p = 1;
    for (;;) {
        long double q = 1 - p, s = powl(q, k), t = powl(p, k),
                    d = t * (1 - q * s) - x;
        if (fabsl(d) < 1e-7L)
            return p;
        p -= p * d / (t * (k + s * (p * (2 * k + 1) - k)));
    }
}

If I had to do this a lot, I would use an integer power algorithm. It's about three times as fast as powl(), which handles non-integer exponents that do not occur here.

long double pownl(long double x, unsigned n) {
    long double p = 1;
    for (;;) {
        if (n & 1)
            p *= x;
        n >>= 1;
        if (n == 0)
            return p;
        x *= x;
    }
}