Multiplying small numbers and getting Inf

[jk22]

At approx. step 46000 the following code puts "Inf" out, but I can't understand why since it multiplies by itself a value smaller than 1. The sin is ok, but the loop for the power fails :

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

double power(double x, int n)
{
	double result=1.0;
	for(int i=0;i<n;i++)
		result*=x;
	return(result);
}

int main(void)
{
	int nsub=800;
	int iter=50000;
	double sum=0.0;

	int *distr=new int[nsub];

	for(int i=1;i<iter;i++)
	{
		double val=pow(sin((double)i),i*i);
		int subn=(int)(val*(double)nsub/2.0);

		distr[subn+nsub/2]++;
		if(fabs(val)>1.0) printf("%d %lf %lf\n",i,sin((double)i),val);
		sum+=val;
	}
	printf("%lf\n",sum);

	FILE *seriesout=fopen("series.out","w");

	for(int i=0;i<nsub;i++)
		fprintf(seriesout,"%lf %lf\n",2*(double)i/(double)nsub-1.0,(double)distr[i]/(double)iter);

}

Is it because the i squared is too big for the int type ?

 

[Mark44]

When you post code here, please use [code[/color]] and [/code[/color]] tags. These tags preserve your indentation, which makes your code easier to understand. I have done this below.

As to your problem, I suspect that when i = 46000, i*i is too large for a signed int. When this happens, the result becomes negative. Raising a number between -1 and 1 to a negative power results in a very large number.

An int is typically stored in four bytes. For your system, I'll bet that "int" means "signed int." The largest signed int would be 2³¹ - 1 = 2,147,483,647. When i is 46,000, i² is 2,116,000,000, which is pretty close to the maximum signed integer.

Changing from int to long might fix your problem.

At approx. step 46000 the following code puts "Inf" out, but I can't understand why since it multiplies by itself a value smaller than 1. The sin is ok, but the loop for the power fails :

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

double power(double x, int n)
{
	double result=1.0;
	for(int i=0;i<n;i++)
		result*=x;
	return(result);
}

int main(void)
{
	int nsub=800;
	int iter=50000;
	double sum=0.0;

	int *distr=new int[nsub];

	for(int i=1;i<iter;i++)
	{
		double val=pow(sin((double)i),i*i);
		int subn=(int)(val*(double)nsub/2.0);

		distr[subn+nsub/2]++;
		if(fabs(val)>1.0) printf("%d %lf %lf\n",i,sin((double)i),val);
		sum+=val;
	}
	printf("%lf\n",sum);

	FILE *seriesout=fopen("series.out","w");

	for(int i=0;i<nsub;i++)
		fprintf(seriesout,"%lf %lf\n",2*(double)i/(double)nsub-1.0,(double)distr[i]/(double)iter);

}

Is it because the i squared is too big for the int type ?

 

[jk22]

As to your problem, I suspect that when i = 46000, i*i is too large for a signed int. When this happens, the result becomes negative. Raising a number between -1 and 1 to a negative power results in a very large number.

Thanks, I'll try with unsigned int.

But if let say i*i becomes a negative number, the power subroutine does not even enter the loop so shouldn't the result be 1.0 ? This subroutine power does not compute negative powers.

 

[Mark44]

The loop is controlled by the value of i, not i*i.

... the power subroutine does not even enter the loop

?

The call to pow() is inside the loop. The loop runs 49,999 times. When i is large enough (a little above 44,000) i * i becomes negative due to overflow. Your compiler should have some documentation about the pow() function.

BTW, C and C++ programmers pretty much never use the term "subroutine."

 

[D H]

You are calling the system function pow() rather than your power(). That's a good thing. You could fix your program so it calls your function instead of pow(), fix your calculation so it computes i*i correctly, and you fix your function so it takes a long instead of an int. Don't do that. Your program would take an hour or so to run if you were calling your function.

There is a fast way to compute xⁿ that doesn't use pow(). x⁸ can be done with only three multiplies rather than eight. x¹⁶, with only four multiplies. x³², with only five multiplies. Done properly, the number of multiplications will grow logarithmically rather than linearly.

 

[jk22]

Thanks I misregarding that I used the system function.

As I can see in this programm the sin(n)^{n^2} seems to tend towards 0, but mathematically we saw in another thread that the limit should not exist : https://www.physicsforums.com/showthread.php?t=715871&page=2

Hence, for sin power n squared, there should exist a subsequence not tending towards 0.

I suppose however the numerical result to be correct, hence there should be an error in the mathematical derivation. I don't know in which section I should write this since it's a comparison between numerical results and analytical ones.

 

[Mark44]

Thanks I misregarding that I used the system function.

As I can see in this programm the sin(n)^{n^2} seems to tend towards 0,

What does this mean? Is n being raised to the power n², or is sin(n) being raised to that power?

I could probably figure out your intent by looking at your code again, but I shouldn't need to do that.

but mathematically we saw in another thread that the limit should not exist : https://www.physicsforums.com/showthread.php?t=715871&page=2

Hence, for sin power n squared, there should exist a subsequence not tending towards 0.

I suppose however the numerical result to be correct, hence there should be an error in the mathematical derivation. I don't know in which section I should write this since it's a comparison between numerical results and analytical ones.

 

[D H]

Thanks I misregarding that I used the system function.

As I can see in this programm the sin(n)^{n^2} seems to tend towards 0, but mathematically we saw in another thread that the limit should not exist : https://www.physicsforums.com/showthread.php?t=715871&page=2

Hence, for sin power n squared, there should exist a subsequence not tending towards 0.

I suppose however the numerical result to be correct, hence there should be an error in the mathematical derivation. I don't know in which section I should write this since it's a comparison between numerical results and analytical ones.

The error is yours. You are assuming two key things here, both of them false:

• You are assuming you are going to find it by looking at the sequence of integers from 1 to 50,000. You aren't.
• You are assuming that you are going to find this sequence using a computer's built-in floating point arithmetic. You aren't.

This is a problem of behavior as n grows very large. 50,000 is not a large number. You missed, for example, n=51819, 52174, 156167, and 260515. There are lots of integers where sin(n) is very close to ±1. A lot of those candidates are going to go to near zero upon raising sin(n) to the n² power, but some won't. (I picked four that won't. Try them.)

Here's what happens with the partial sums of $\sum_{n=1}^{\infty} (\sin n)^{(n^2)}$. It's going to bounce around some but then appear to stabilize to about 1.926. Then, blam, you come across n=51819 and n=52174. Then a slight jump when you hit n=156167. Then a huge jump when you hit n=260515. You'll keep getting these seemingly random points where the partial sums jump. That's because you hit a number n that is very close to pi/2 or 3*pi/2, modulo 2*pi.

If n is very close to pi/2 (modulo 2*pi), then sin(n) is close to 1, so sin(n)⁽ⁿ²⁾ is going to be positive. If n is is very close to 3*pi/2 (modulo 2*pi), then sin(n) is close to -1, so sin(n)⁽ⁿ²⁾ is going to be positive if n is even, negative if n is odd. On average, you'll see three significant positive values of sin(n)⁽ⁿ²⁾ to every significant negative value amidst the sea of values that are extremely close to zero. The spacing between these significant values will increase, but they're always going to be there. The series grows without bound because of that 3:1 positive to negative ratio.