Some sort of issue with simple int/double operations-CodePudding

I am new to C and I tried to put together a (only slightly) nontrivial piece of code, in which I use an array to store the values of atan(x) in bins of width dx from 0 to 1, then use the rectangular method to compute the integral of atan(x)dx from 0 to 1. The code should loop while successively making dx smaller to see the integral converge to the analytical result. I have been totally unable to figure out why I get expected outputs for most entries in the loop, but (some number)e 268 for the 7th/8th/15th outputs; I imagine it has something to do with int/double casting but the irregularity of it is very strange to me.

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

int main() {
    int array_len = 20;
    int array_len_new;
    int num_conv = 18;
    double linspace[200];
    double conv_array[200];
    double PI = 3.142857;
    double result;
    int i;
    int j;
            
    for (i = 0; i < num_conv; i  ) {
        array_len_new = array_len   10*i;
        double dx = 1./array_len_new;
        for (j = 0; j < array_len_new; j  ) {
            linspace[j] = j* 1.0/array_len_new;
            conv_array[i]  = atan(linspace[j])*dx;
        }
        printf("Result for %d bins is: %e\n", array_len_new, conv_array[i]);
    }

    printf("Converged result: %e\n", (PI-log(4.))/4.0);
    return 0;
}

Output:

Result for 20 bins is: 4.190854e-001
Result for 30 bins is: 4.256883e-001
Result for 40 bins is: 4.289811e-001
Result for 50 bins is: 4.309539e-001
Result for 60 bins is: 4.322680e-001
Result for 70 bins is: 4.332061e-001
Result for 80 bins is: 2.308177e 268
Result for 90 bins is: 2.308177e 268
Result for 100 bins is: 4.348934e-001
Result for 110 bins is: 4.352511e-001
Result for 120 bins is: 4.355492e-001
Result for 130 bins is: 4.358013e-001
Result for 140 bins is: 4.360175e-001
Result for 150 bins is: 4.362047e-001
Result for 160 bins is: 2.316093e 268
Result for 170 bins is: 4.365131e-001
Result for 180 bins is: 4.366416e-001
Result for 190 bins is: 4.367566e-001
Converged result: 4.391407e-001

EDIT: I have found that the issue resolves itself if I change the length of conv_array (which only needs 18 values, I made it big because I thought it wouldn't matter) from 200 to 18, or even 100. Why could this possibly be the case?

CodePudding user response：

You do not initialize the elements of conv_array. As a non-static local variable, its initial value is indeterminate. In particular, it is not safe to assume that it will be zero-initialized, but your code relies on exactly that assumption. Declaring it with an initializer should help:

    // ...
    double linspace[200];
    double conv_array[200] = {0};
    double PI = 3.14159265358979;
    // ...

The = {0} says that the first element is to be initialized to 0, and if any element is initialized then those without explicit initializers are implicitly default-initialized (to 0).

Personally, however, I would not use an array at all where you use conv_array, nor where you use linspace. Try this on for size:

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

int main() {
    int base_num_bins = 20;
    int num_cycles = 18;
    double PI = 3.14159265358979;
            
    for (int i = 0; i < num_cycles; i  ) {
        int num_bins = base_num_bins   10 * i;
        double dx = 1. / num_bins;
        double integral = 0;

        for (int j = 0; j < num_bins; j  ) {
            integral  = atan((double) j / num_bins) * dx;
        }
        printf("Result for %d bins is: %e\n", num_bins, integral);
    }

    printf("Converged result: %e\n", (PI - log(4.)) / 4.0);
    return 0;
}

CodePudding user response：

conv_array has not been initialised, but you do

conv_array[i]  = atan(linspace[j])*dx;

I suggest

double conv_array[200] = { 0 };

Side issue:

double PI = 3.142857;

might be ok with float but it isn't a very accurate value for π when working with double. I suggest using M_PI instead of double PI; and with MSVC you'll also need

#define _USE_MATH_DEFINES