Double precision

[helios]

After running this:

#include <iostream>
#include <iomanip>
#include <cmath>

const double log_of_2=log(2.0);

inline double log2(double x){
	return log(x)/log_of_2;
}

int main(){
	double n=0xFFFFFFFF,
		accum=32.0;
	for (;n>0;--n)
		accum+=log2(n);
	std::cout <<std::setprecision(16)<<accum;
	return 0;
}

I got 131242625471.8555 as output, which is within 2^-32 of error from 131242625454.7232, which is what I got as the definite integral of log2(x) between 1 and 2^32, which I had estimated to be equal or very close to the exact value of the sum.
Considering one computation took billions of floating point operations more than the other, should I trust the value of the integral to be the exact value, or the operations I performed in the program conserve enough precision for the result to be trusted?
The program was compiled with g++ test.cpp -o test -O3 -s and ran on a Core 2 Duo.

[hamsterman]

That sum is not equal to the integral. Consider both methods for range (0; 1].
In fact I was just going to complain how your method is wrong as the integral is smaller than the sum..
^{wait, is the correct expression "to be just about to" ?..}

[helios]

I know, but it's definitely possible that the integral is closer to the actual value of the sum than what the program returned. The error being such a special number is certainly suspicious.
The question is more about the precision of floating point operations, particularly in the FPU->memory transition, than about math.

Consider both methods for range (0; 1].

Is it possible to integrate over open intervals? _{EDIT: Ah, I see. Improper integrals are limits.} In any case, the point is moot because there are no natural numbers in (0;1), so the sum doesn't make sense, either.

wait, is the correct expression "to be just about to" ?..

Both are fine.

[kev82]

I've been playing with this and haven't really got anywhere because the difference between the sum and the integral is

\sum_{i=1}^{2^32} 1 + i\ln(\frac{i+1}{i})

The i'th term is strictly increasing with i so the difference gets worse the further you go.

The best I can do is bound the sum as follows. This does show the integral isn't far off though.

Let s_n = \sum_{i=1}^{n} \ln(i)
Let i_n = \int_{x=1}^{n} \ln(x) dx

Simple numerical integration gives s_{2^32-1} < i_{2^32} < s_{2^32}

s_{2^32} - s_{2^32 - 1} = \ln{2^32}

Substituting gives

s_{2^32} < i_{2^32} + \ln{2^32}
s_{2^32} > i_{2^32}

To convert ln to log2 just divide by ln 2 everywhere. That gives a lower bound of the integral and an upper bound of the integral+32.

[helios]

Blast these darn large numbers!

[kev82]

Not that it matters as the thread is about precision (which I know nothing about), but I have improved my bound slightly.

Consider triangles that go through the points (x,log(x)), (x+1,log(x)), (x+1,log(x+1)). They have gradient 1/(x+1) and go through the point (x,log(x)). Because log(x) is convex, they can't cross log(x), so if we add up their area we still stay below the integral. Their total area is 1/2 (H_n - 3/2) where H_n is the nth harmonic number.

Substituting this in the inequality gives

s_n - log_2 (n) + 1/2 (H_n - 3/2) < i_n

Because H_n appears with positive coefficient on the less than side of the inequality, we can replace it by a lower bound. It just so happens that ln(n) is a lower bound for the harmonic function, so substituting and tidying up gives.

s_n < i_n + log2(n)(1-ln(2)/2) + 3/4

So the difference between your sum and the integral can be at most 21.66.

You can do exactly the same trick above log(x) and derive the difference must be at least 10.59.

[helios]

If that's true, then it proves that whatever the actual value of the sum, the directly computed one will always be closer than the integral's. Which, after thinking about it for a while, is to be expected. There's a whole chunk of surface it misses at every step, and it looks a bit like this:

*************************************
*******************
********
****
**
*

Well, it was fun. Did you have fun, kev82?

[kev82]

I really wanted a series expansion where I could bound the error term, and hence calculate the result to a desired accuracy, but the integral for the Chebyshev coefficients looked horrible, and Maple confirmed it.

I do think it's quite neat that with relatively simple maths, and a little picture, you can derive

s-i<21.66
s-i>10.59

A little bit of fun, yeah. But anything is fun when you've still got 50 exams to mark on the independent samples t test. I think I just miss real maths. I don't do anything like this anymore.

[helios]

Would you happen to know why when I give sum(i=1,2^32,log2(i)) to a CAS, it gives back log2(gamma(2^32+1)), when I arrived at the summation from log2((2^32)!)? It just doesn't make sense.

What do you do now?

[kev82]

The gamma function is a generalization of factorial, but for natural n is equivalent to (n-1)! So your CAS has basically undone the manipulation that you did.

I am supposed to be finishing my PhD, but I tend to rock climb and run most of my time. I also work part time for a consultancy that looks at the scheduling of home care nurses/staff.

[helios]

for natural n is equivalent to (n-1)!

Ah. I didn't remember that. I probably should have looked it up.

I mean what do you do that would provoke the "I miss real math" part.

[kev82]

My masters degree is in maths, with a hefty emphasis on approximation theory and numerical analysis.

My consultancy work/PhD is metaheuristic optimization. There is very little `proper' maths in it.

[helios]

That's applied math, isn't it?

Actually, some of the parts that most interest me in CS overlap with mathematics (e.g. language theory and anything to do with symbolic manipulations; computer graphics, which applies many tools from linear algebra and calculus; data compression), so I guess I'll end up moving to a similar region as you from a different direction. Hopefully I won't then complain about a lack of proper programming.

[kev82]

Approx theory and numerical analysis are pretty much just applied analysis. Metaheuristic optimization, at least what I do, is simulation and graph theory.