Need help fixing an algorithm that approximates pi

I'm trying to write the C code for an algorithm that approximates pi. It's supposed to get the volume of a cube and the volume of a sphere inside that cube (the sphere's radius is 1/2 of the cube's side). Then I am supposed to divide the cube's volume by the sphere's and multiply by 6 to get pi.

It's working but it's doing something weird in the part that is supposed to get the volumes. I figure it's something to do the with delta I chose for the approximations. With a cube of side 4 instead of giving me a volume of 64 it's giving me 6400. With the sphere instead of 33 it's giving me 3334. something.

Can someone figure it out? Here is the code (I commented the relevant parts):

#include <stdio.h>      

int in_esfera(double x, double y, double z, double r_esfera){
    double dist = (x-r_esfera)*(x-r_esfera) + (y-r_esfera)*(y-r_esfera) + (z-r_esfera)*(z-r_esfera);

    return  dist <= (r_esfera)*(r_esfera) ? 1 : 0;   
}   

double get_pi(double l_cubo){   
    double r_esfera = l_cubo/2;   
    double total = 0;
    double esfera = 0;    
//this is delta, for the precision. If I set it to 1E anything less than -1 the program continues endlessly. Is this normal?
    double delta = (1E-1);   

    for(double x = 0; x < l_cubo; x+=delta){
        printf("x => %f; delta => %.6f\n",x,delta);
        for(double y = 0; y <l_cubo; y+=delta){
            printf("y => %f; delta => %.6f\n",y,delta);
            for(double z = 0; z < l_cubo; z+=delta){
                printf("z => %f; delta => %.6f\n",z,delta);
                total+=delta;
                if(in_esfera(x,y,z,r_esfera))
                    esfera+=delta;
            }
        }
    }

    //attempt at fixing this
        //esfera/=delta;
        //total/=delta;
    //

//This printf displays the volumes. Notice how the place of the point is off. If delta isn't a power of 10 the values are completely wrong.   
    printf("v_sphere = %.8f; v_cube = %.8f\n",esfera,total);   

    return (esfera)/(total)*6;
}   

void teste_pi(){        
    double l_cubo = 4;    
    double pi = get_pi(l_cubo);

    printf("%.8f\n",pi);
}   

int main(){   
    teste_pi();
}

Solution

The thing is that multiplication over integers like a * b * c is the same as adding 1 + 1 + 1 + 1 + ... + 1 a * b * c times, right?

You're adding delta + delta + ... (x / delta) * (y / delta) * (z / delta) times. Or, in other words, (x * y * z) / (delta ** 3) times.

Now, that sum of deltas is the same as this:

delta * (1 + 1 + 1 + 1 + ...)
         ^^^^^^^^^^^^^^^^^^^^ (x * y * z) / (delta**3) times

So, if delta is a power of 10, (x * y * z) / (delta**3) will be an integer, and it'll be equal to the sum of 1's in parentheses (because it's the same as the product x * y * (z / (delta**3)), where the last term is an integer - see the very first sentence of this answer). Thus, your result will be the following:

delta * ( (x * y * z) / (delta ** 3) ) == (x * y * z) / (delta**2)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the sum of ones

That's how you ended up calculating the product divided by delta squared.

To solve this, multiply all volumes by delta * delta.

However, I don't think it's possible to use this logic for deltas that aren't a power of 10. And indeed, the code will go all kinds of haywire for delta == 0.21 and l_cubo == 2, for example: you'll get 9.261000000000061 instead of 8.