C, a derivative with eps

Question

this is a function to find a gradient from some function f.and everything works with it correctly

void grad(double* y,   
    double* x,    
    int n,        
    double d,    
    double (*f)(double*, int))
{ 
    double f0 = f(x, n);
    for (int i = 0; i < n; ++i) {
        double save = x[i];
        x[i] += d;
        double f1 = f(x, n);
        x[i] = save;
        y[i] = (f1 - f0) / d;
    }
}

but then I wanted to enter the precision (eps) and converted the code a little bit. But now zeros have started to be output. What exactly am I doing wrong to achieve the desired accuracy of calculations?

   void grad(double* y,     
        double* x,     
        int n,        
        double d, 
        double eps;   
        double (*f)(double*, int)) { 
        double f0 = f(x, n);
        double d1, d2;
        for (int i = 0; i < n; ++i) {
            double save = x[i];
            x[i] += d;
            double f1 = f(x, n);
            x[i] = save;
            do {
                y[i] = (f1 - f0) / d;
                  d1 = y[i];
                d = d / 10.0;
                y[i] = (f1 - f0) / d;
                 d2 = y[i];
            } while (fabs(d1 - d2) >= eps);
            
        }
    }

full program Before I changed it

void grad(double* y,    
    double* x,   
    int n,        
    double d,     
    double (*f)(double*, int)) { 
    double f0 = f(x, n);
    for (int i = 0; i < n; ++i) {
        double save = x[i];
        x[i] += d;
        double f1 = f(x, n);
        x[i] = save;
        y[i] = (f1 - f0) / d;
    }
}

double f(double* x, int n) {
    int i; double sum=0;
    for (i = 0; i < n; i++) {
        sum += x[i] * x[i]* x[i];
    }
    return sum;
}

int main() {
    int n; int i;
    scanf_s("%d", &n);
    double* x;
    x = (double*)malloc(n * sizeof(double));
    double* y;
    y = (double*)malloc(n * sizeof(double));
    for (i = 0; i < n; i++)
    {
        printf("x[%d] = ", i);
        scanf_s("%lf", &x[i]);
    }
    grad(y, x, n, 0.001, f);
    printf("grad f(x0,x1,..xn)\n");
    for (i = 0; i < n; i++) {
        printf(" %lf ",y[i]);
    }
   // printf("grad f(1,2) = %lf i + %lf j\n", y[0], y[1]);
}

Answer 1

What exactly am I doing wrong to achieve the desired accuracy of calculations?

Uninitialized object

In fabs(d1 - d2) >= eps, eps is uninitialized.

---

Fixed offset

x[i] += d; (a fixed offset) takes the floating out of floating point math.

x[i] += d is the same as x[i] += 0.0 when x[i] is much larger than d.

x[i] += d is the same as x[i] = d when x[i] is much smaller than d.

Both extremes demo that x[i] += d is not a good way to find a suitable offset for a multi-variable function slope calculation. Easy to incur divide by 0 in later code.

---

Better: when fabs(x[i]) > DBL_MIN, scale instead: x[i] *= e where e is 1.0001, 1.0000000001, 1.000000000001, etc. and adjust slope code calculation accordingly. This too does have limitations and more code to handle wee x[i].

---

Even better: IMO, offset x[i] by N floating point numbers, where N is some modest integer. E.g. N==1 --> x[i] = nextafter(x[i], x[i]*2);.

---

Best: need more info as best is case dependent.

---

FP debug tip

But now zeros have started to be output

Use "%e" or "%g" instead of "%f". It is more informative. "%f" prints zero for about half of all double values.

// printf(" %lf ",y[i]);
// Alternatives
printf(" %e ", y[i]);
printf(" %.*e ", 16, y[i]);
printf(" %g ", y[i]);
printf(" %.*g ", 17, y[i]);
printf(" %a ", y[i]);