Why does using double instead of float gives a wrong result in this double integration code?

Question

I found the following code in this page to compute a double integral. Whenever I run it with all variables being declared as float, it gives the right result for the example integral, which is 3.91905. However, if I just change all float variables to double, the program gives a completely wrong result (2.461486) for this integral.

Could you help me undestanding why this happens? I expected to have a better result using double precision, but that's evidently not the case here.

Below is the code pasted from the aforementioned website.

// C++ program to calculate
// double integral value

#include <bits/stdc++.h>
using namespace std;

// Change the function according to your need
float givenFunction(float x, float y)
{
    return pow(pow(x, 4) + pow(y, 5), 0.5);
}

// Function to find the double integral value
float doubleIntegral(float h, float k,
                    float lx, float ux,
                    float ly, float uy)
{
    int nx, ny;

    // z stores the table
    // ax[] stores the integral wrt y
    // for all x points considered
    float z[50][50], ax[50], answer;

    // Calculating the number of points
    // in x and y integral
    nx = (ux - lx) / h + 1;
    ny = (uy - ly) / k + 1;

    // Calculating the values of the table
    for (int i = 0; i < nx; ++i) {
        for (int j = 0; j < ny; ++j) {
            z[i][j] = givenFunction(lx + i * h,
                                    ly + j * k);
        }
    }

    // Calculating the integral value
    // wrt y at each point for x
    for (int i = 0; i < nx; ++i) {
        ax[i] = 0;
        for (int j = 0; j < ny; ++j) {
            if (j == 0 || j == ny - 1)
                ax[i] += z[i][j];
            else if (j % 2 == 0)
                ax[i] += 2 * z[i][j];
            else
                ax[i] += 4 * z[i][j];
        }
        ax[i] *= (k / 3);
    }

    answer = 0;

    // Calculating the final integral value
    // using the integral obtained in the above step
    for (int i = 0; i < nx; ++i) {
        if (i == 0 || i == nx - 1)
            answer += ax[i];
        else if (i % 2 == 0)
            answer += 2 * ax[i];
        else
            answer += 4 * ax[i];
    }
    answer *= (h / 3);

    return answer;
}

// Driver Code
int main()
{
    // lx and ux are upper and lower limit of x integral
    // ly and uy are upper and lower limit of y integral
    // h is the step size for integration wrt x
    // k is the step size for integration wrt y
    float h, k, lx, ux, ly, uy;

    lx = 2.3, ux = 2.5, ly = 3.7,
    uy = 4.3, h = 0.1, k = 0.15;

    printf("%f", doubleIntegral(h, k, lx, ux, ly, uy));
    return 0;
}

Thanks in advance for your help!

Answer 1

Due to numeric imprecisions, this line:

ny = (uy - ly) / k + 1;  // 'ny' is an int.

Evaluates to 5 when the types of uy, ly and k are float. When the type is double, it yields 4.

You may use std::round((uy - ly) / k) or a different formula (I haven't checked the mathematical correctness of the whole program).