Double integral help

I think I'm doing everything correctly but for some reason the program refuses to compute, it just freezes.

I'm given an equation for a curve, and for certain values of x and y along that I want to take the double integral and set that over the double integral of the whole curve and compute that fraction.

To do this I loop through values of x and y and then multiply them to get the double integral of the whole curve.

For the isolated section of the curve I substitute a new variable z for y and then within the x loop I cycle through z values for each value of x and if x and z are within certain values I compute the smaller integral.

I then set the two integrals over each other and compute that fraction.

What am I doing wrong?

Code:

#include <stdio.h>
#include <math.h>
#define STEP 5.556e-5

int main() {
float a=0.3489, b=0.255, R=0.68085, xlimit=0.2799, zlimit=-0.0008; 
float utotalarea=0; 
float vtotalarea=0; 
float x, xarea, xtotalarea=0; 
float y, yarea, ytotalarea=0; 
float z, zarea; 
float result, result2, fraction;
//First loop through x values.
for (x=STEP; x<=10; x=x+STEP){
	xarea=0.5*(exp(-log(2)*pow(x/a,2)-1)+exp(-log(2)*pow((x-STEP)/a,2)-1))*STEP;
	xtotalarea=xtotalarea+xarea;
		//Next, loop through z values for each x value. Note: z is really y, its a different variable here in order to calculate the integral for certain x and y values.  
		for (z=STEP; z<=10; z=z+STEP){
			zarea=0.5*(exp(-log(2)*pow(z/b,2)-1) + exp(-log(2)*pow((z-STEP)/b,2)-1))*STEP;
			//If x and z are within certain values take integral.			
			if((pow(x-xlimit,2)+pow(z-zlimit,2))<=pow(R,2)){
				utotalarea=utotalarea+xarea;				
				vtotalarea=vtotalarea+zarea;
			}
		}
	}
//Now, loop through all y values.
for (y=STEP; y<=10; y=y+STEP){
	yarea=0.5*(exp(-log(2)*pow(y/b,2)-1)+exp(-log(2)*pow((y-STEP)/b,2)-1))*STEP;
	ytotalarea=ytotalarea+yarea;
	}
//Compute double integral for the whole curve.
result=xtotalarea*ytotalarea;
//Compute double integral for the certain x and z values.
result2=utotalarea*vtotalarea;
//Set the smaller integral over the larger integral.
fraction=result2/result;
//print the fraction. 
printf("%.5f\n", fraction);
}

Since both your outer and inner loop are repeating 10 / .00005556 = 180000 times, your inner loop is repeating about 32.4 billion times. It might finish in a two or three hours (days?).

for kicks just increase STEP to 1 and see if it runs. the answer won't be what you want but you will see if it runs to completion.

I would suggest to never use float. Use double instead. float is a relic from decades ago where one would sacrifice numerical accuracy for memory. Nowadays all math.h functions use and return double. It's actually faster to stick with double than implicit casting to float all the time.

Originally Posted by nonoob

I would suggest to never use float. Use double instead. float is a relic from decades ago where one would sacrifice numerical accuracy for memory. Nowadays all math.h functions use and return double. It's actually faster to stick with double than implicit casting to float all the time.

Only when memory bandwidth is not a concern. Afterall, 16-bit halffloats are popular in 3D graphics because they take half the memory of floats, which is great if you have a lot of them, say for a z-buffer.
Last time I checked the difference in speed on floats vs doubles in my software 3D engine, floats were still faster.

Having said that, I agree that doubles are best in this case.

I'd suggest it's the algorithm (the complete logic of how you organise the calculations) that needs refining. Changing from float to double might paper over the gaps, but can equally fail again if some parameters get tweaked.

However, if you're not willing to apply some thought into changing algorithm, you can simplify the code you have......

It is generally considered a no-no to use a floating point variable to control a loop. Floating point variables are NOT real numbers, they are a finite approximation, so all sorts of errors creep in. Particularly if you are trying to iterate (i.e. add a fixed value each time through the loop). Thanks to finite precision of floating point, rounding error can build up, the value of the control variable will be affected by error, and may well be too inaccurate for needs. Do something like.

Code:

   for (int i = 1; i <= some_number; ++i)
   {
          x = STEP * i;
           .....

rather than using x as the control variable and adding a small value of STEP to it.

You might want to look into simplifying various expressions in the loops, to avoid replication of calculations.

x*x is better - computationally - for squaring x, than pow(x,2).

The fact you have exp(-log(2) .... in several expressions suggests the formula you're using actually seeks to raise 2 to various powers. In the calculation of xarea, for example, exp(-log(2)*pow(x/a,2)-1) can be simplified as 1.0/(pow(2.0, x*x/a*a)*exp(1.0)). That can be simplified even further, since a and exp(1.0) are both non-varying values in your code. So you might hoist a fair few repetitive calculations out of the loops (only compute them once, rather than for every nested iteration, and store their values in a temporary variable).

As I said, however, you would be better off looking at whether there is a better algorithm. There are lots of libraries out there for doing numerical integrations of various special functions.

I agree that for this problem, using float doesn't make a lot of sense. But the memory saving aspects of float are certainly not irrelevant even with huge main memories these days. You might have 10 GB of RAM in your box, but the poor little core processing your data still only has 32k of L1 cache (64k on AMD)