3

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();
}
Roberto Caboni
  • 7,252
  • 10
  • 25
  • 39
Hypercube
  • 89
  • 6
  • Your `get_pi` function is like `O((l_cubo / delta)**3)` (_cubic_ complexity! No wonder it "continues endlessly"). Is it supposed to calculate the volumes of the cube & the sphere? – ForceBru Apr 10 '20 at 17:48
  • @ForceBru Partly. First that's what it does yes (check the last printf on that function) and then returns the ratio * 6 (or pi). – Hypercube Apr 10 '20 at 17:50
  • Turns out, you need to _multiply_ both volumes by `(delta * delta)`, like `esfera *= delta * delta;` to fix this. It's simple math, but I'm not quite sure how to put it into words lol – ForceBru Apr 10 '20 at 18:02
  • @ForceBru That works even when delta isn't a power of 10. Thanks! I'm curious as to why it works... – Hypercube Apr 10 '20 at 18:06
  • @ForceBru they could have picked square and circle, if they meant to reduce time complexity. Also, seems to slightly lower error. – JunisvaultCo Apr 10 '20 at 19:24
  • Additionally, the precision could be seriously improved by not using double or floating point, but int (or maybe even long long) wherever possible. Instead of having double delta and double l_cubo, use double *only* when returning pi(and the division where pi is calculated). Otherwise you are adding up precision errors. – JunisvaultCo Apr 10 '20 at 19:34

2 Answers2

2
total+=delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta;

total and esfera are three-dimensional volumes whereas delta is a one-dimensional length. If you were tracking units you'd have m3 on the left and m on the right. The units are incompatible.

To fix it, cube delta so that you're conceptually accumulating tiny cubes instead of tiny lines.

total+=delta*delta*delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta*delta*delta;

Doing that fixes the output, and also works for any value of delta:

v_sphere = 33.37400000; v_cube = 64.00000000
3.12881250

Note that this algorithm "works" for arbitrary delta values, but it has severe accuracy issues. It's incredibly prone to rounding problems. It works best when delta is a power of two: 1/64.0 is better than 1/100.0, for example:

v_sphere = 33.50365448; v_cube = 64.00000000
3.14096761

Also, if you want your program to run faster get rid of all those printouts! Or at least the ones in the inner loops...

John Kugelman
  • 349,597
  • 67
  • 533
  • 578
0

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.

ForceBru
  • 43,482
  • 10
  • 63
  • 98