computation challenge



#2

here's one that i thought would be easy.

a lot of people know that 1+1/2+1/3+1/4+1/5 diverges and that 1+1/2^2+1/3^2+1/4^2+... converges to pi^2/6 in fact, and there are simpleish formulas for the limit sums of 1/n^(2k), ie even powers. but no one knows if there are easy formulae for odd powers.

lets look at the most simple case, cubes. write a program to find the limit sum of 1+1/2^3+1/3^3+1/4^3+... if you're really smart, try other odd powers.

thinking this was a five minute job, i wrote a simple loop to do the obvious adding, but after some 600 iterations (some time on a hp41c) it had only managed 5 decimal places. can you do 10?


#3

Hugh,

It's gonna take a lot of terms until you get the value to the tenth decimal place (for the cube). A quickie Fortran program on my PC indicates that you'll need somewhere between 50,000 and 100,000 terms until there is no more change in the tenth place, and you'll need almost 250k terms if you are worried about round-off (i.e. at 250k terms, the eleventh and beyond digits indicate you need to round up by one in the tenth place). However, I'm not sure the calculator has enough precision - you need 11 digits: the answer is 1.202 056 903 15 (to 11 decimal places).

I don't know for sure how the '41 calculates powers, but I suspect you can increase your speed by multiplying 1/n three times rather than taking the third power. On the PC (in Lahey fortran), multiplying was about three times faster than taking the third power of 1/n .

I hate to make your calculator seem like a slowpoke, but on a modest PC (a 200 MHz pentium pro, accelerated to 333 MHz), summing 10 million terms (using multiplication and not exponentiation) takes about 2.5 seconds with an old (10+ years) Lahey Fortran compiler, using double precision.


#4

<just kidding mode>
Yes, Pentium is fast.
If only it were accurate...
</just kidding mode>

#5

I am always trying to improve my own knowledge (and skills, if possible) when we're dealing with subjects like this one.

I don't know what sort of method you use in your program to compute how many terms are needed till the tenth digit does not change. I'm curious because the HP41 uses 12-digit internal precision and 10-digit mantissa results. I know it's not necessary to relate one thing to another, but you mentioned you used a double-precision compiler to compute a 10-million terms summation in 2.5 seconds, right?

I'm afraid you'll think my question is irrelevant, but I'm curious: does your 50,000 to 100,000 estimation take 12-digit to 10-digit rounding after each loop or it is based on FORTRAN's double precision? Maybe it's based in another approach, I don't know.

I'm curious because after reading the HP15C Advanced Functions Handbook and understanding that internal limits may also not change external results, I always try to understand the rounding consequences of series computations, like the one we're talking about here.

Please, I'm not asking you to explain your conclusions, I just want to know how did you get to them. I think they are clever and understanding, and if they can be expanded to systems with different digit numbers (like Pioneers and the RPL models), maybe we can significantly reduce computation time based on the difference achieved between terms instead of wait till they are steady.

Thank you.

Luiz (Brazil)

Edited: 9 Aug 2003, 7:06 p.m.


#6

hi guys,

dave is right. you need too many terms to sum the series directly. this is what i had discovered in my "5 minutes" answer.

so a viable solution involves some kind of series acceleration. one discovery is that a series of the alternating terms (ie 1-1/2^3+1/3^3-1/4^3+...) converges faster.

let z(s) = 1+1/2^s+1/3^s+1/4^s+... we'd like, at first, z(3). now z(s) = (1-1/2^3+1/3^3-1/4^3+...)/(1-2^(1-s))

this is better. i get a good answer in around 2000 terms. but still way too many for the 41c.

there are further accelations possible. here is my new program to compute z(s).

LBL "Z"
STO 00
1
STO 01
STO 02
STO 05
0
STO 06
STO 07
10
STO 09
STO 08
LBL 00
RCL 01
1/X
RCL 00
Y**X
RCL 02
*
ST+ 06
1
ST+ 01
RCL 02
CHS
STO 02
DSE 08
GTO 00
2
RCL 09
Y**X
STO 03
STO 04
LBL 01
RCL 05
ST- 04
RCL 01
1/X
RCL 00
Y**X
RCL 04
*
RCL 02
*
ST+ 07
1
ST+ 01
RCL 02
CHS
STO 02
RCL 09
RCL 08
-
RCL 08
1
+
STO 08
/
ST* 05
RCL 09
RCL 08
X<Y?
GTO 01
RCL 07
RCL 03
/
RCL 06
+
.5
RCL 00
1
-
Y**X
1
-
CHS
/
END

run as 3 ENTER XEQ Z to give 1.202056903 in a few seconds. other values of s work too, for s >= 2.

#7

Luiz,

Here is my fortran program (with proper line breaks, I hope!):

      PROGRAM POWERSUM
C
C CALCULATES SUM OF POWERS OF RECIPROCAL INTEGERS
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C DO A LOT OF NUMBERS TO GET TIMING
C TRY IT TWO WAYS - MULTIPLY & EPONENTIATION
C
SUM = 0D0
DO 100 N = 1,10000000
X = 1D0/DFLOAT(N)
C
C RUN IT ONE WAY OR THE OTHER, BY MOVING THE COMMENT
C
SUM = SUM + X*X*X
C SUM = SUM + X^3
C IF(MOD(N,50000) .EQ. 0) WRITE (*,2000) N,SUM
100 CONTINUE
write (*,2000) n,sum
2000 FORMAT(I9,F25.16)
END

Nothing fancy. The line (commented out here, because I was doing the timing tests) just before the CONTINUE statement gives me the value after every 50,000 terms. I was NOT rounding off; just summing the results to the precision of the compiler/computer.

I just looked at those values to make my statement that it was between 50k and 100k terms that the sum settled down in the tenth decimal place. (I was too lazy to try to find exactly how many terms were needed!) It was looking at the printout that also indicated that somewhere around 200k-250k terms, the tenth decimal place should be rounded up (IF you were interested in only the "correct" properly rounded-off value to ten decimal places).

Double precision fortran should be good to 15 or 16 significant figures - plenty enough for this test. If amybody cares, there is quadruple precision, too, I think.

This doesn't offer much insight into HP calculator precision or round-off. Conversely, if you can define the algorithm that the calculator uses, it would be easy to check it this way.

In any case, Hugh's solution with alternating terms is far more elegant!


#8

Hello, Dave;

thank you very much twice: allowing me to remember FORTRAN and showing me (us) the listing. About fanciness: I like to see programs that run, like yours; it's fancy enough for me.

I did not dedicate enough time to understand FORTRAN in deep, because I learnt it enough to accomplish all tasks related to the language at the university (1979 & 1980). But I want to know if it is possible to add a small modification at this line:

      X = 1D0/DFLOAT(N)
such that 1D0/DFLOAT(N) is rounded to ten significant digits. Is it easy to do, if possible? I think that, given that specific rounding propagation errors do not occur so often, rounding a 16-digit mantissa or a 12-digit mantissa to ten places will give closer but equal approximated results. I thought about writing a program in PASCAL, but I see that yours is ready and we can, at least, see if rounding to ten decimal places will reduce the number of cycles.

If you think it's worth doing and you have the time to (on Sunday...), I'd just like to know what happens.

Thanks again, Dave. It's been enlightening, though.

Luiz (Brazil)

Edited: 10 Aug 2003, 8:04 p.m.


#9

Luiz,

I don't recall a truncation or rounding function in fortran, but it is fairly easy to accomplish.

After the statment X = 1D0/DFLOAT(N) , you would multiply X by 1e10 (where 10 is the number of decimal places you want), then set the floating point value equal to an integer, and then refloat the integer and divide by 1e10.

Because IX = X truncates the value in X, for proper rounding you should add 0.5 to X before integerizing it. (You have to be a bit careful about the sign of X, because setting the floating point value equal to an integer TRUNCATES the floating point number, causing a rounding in the wrong direction for a negative value. You can use the SIGN and ABS functions to take care of the general case.)

Thus a few steps after X = 1D0/DFLOAT(N) that go like this:

(I omit FLOAT and INT operations, which are often inserted automatically by the compiler, although it is good practice to do it yourself, so you KNOW what is happening.)

     X = X * 1E10
IX = X + 0.5
X = IX/1E10

will round X off to 10 decimal places. You need to worry about the overall precision of fortran. You'd probably want to use double precision integers! A regular integer (4 byte) has about 9 or 10 digits of precision (2^31 or thereabouts, you need a bit for the sign - although "twos-complement" style numbers may give this bit back). I never worried to much about that. Astronomers seldom do anything that needs so much precision! The joke is that we worry about significant figures IN THE EXPONENT!

Doing this could falsely indicate how fast a sum converges - you'd be throwing away higher order digits that otherwise add up. However, it should show what a calculator would do.

#10

The value of the limit to 36 digits after the decimal is:

1.202 056 903 159 594 285 399 738 161 511 449 990
^
The 11th digit | after the decimal point should be rounded
up to 6, if what we want is a result correct to 11 digits
after the decimal point.

Summing 329894 terms gives:

1.202 056 903 154 999 980 070

Summing 329895 terms gives:

1.202 056 903 155 000 007 923

So, we would have to sum 329895 terms to get a correctly
rounded result to 11 digits after the decimal point.

Isn't this fun?!! Thanks to Mathematica for all these
results.

#11

Hello!

You can approximate the minimal numbers of terms: Nmin=ceil((1e10)^(1/3))=2155. But this serie is not alternate, the numbers of terms you must to add is N=70711. (with MAPLE7)

If you want more precision, begin the summation with 1/70711^3, then decrease the number in the denominator.

Csaba


#12

70711... Its suspicious... I try it with some numbers...

Yes! The numbers of terms you must to add is N=1/SQRT(2*eps), where eps is the desirable precision.

Cs.


#13

thats remarkable. you do indeed get 70711 with that calculation.

the program i wrote was based on further series acceleration. various techniques here,
http://numbers.computation.free.fr/Constants/Miscellaneous/seriesacceleration.html

so that the number of terms for eps is given by -log((1-2^(1-s))*eps)/log(8). which for eps=1e-10, rounds up to 12.
in fact i use 10 in the program which happens to be enough.


#14

Try it in MAPLE, and play with 'n':

restart;with(plots):
Digits:=20;
a:=sum(1/x^3,x=1..infinity);
b:=sum(1/x^3,x=1..n);
d:=a-b;
n:=424;evalf(d);

And if you plot 'd' in loglog plot, you'll see an linear function:

loglogplot(evalf(d),n=1..100000);

The equation of linear function is:

log(N)=-1/2*log(2)-1/2*log(eps), where base of 'log'-s are anything.

('d' are equal 'eps')

(I know, its not a very elegant solution, I just use it for a fast 'research'...)

Csaba


Possibly Related Threads…
Thread Author Replies Views Last Post
  RE: 35s sorting routine challenge - Gene's Challenge Miguel Toro 4 1,637 08-01-2007, 08:36 AM
Last Post: Miguel Toro
  Internal Computation of Cumulative Distributions Les Wright 15 3,640 06-23-2006, 02:26 PM
Last Post: Les Wright
  Prelude to a 49g+ computation Eric Lundgren 6 1,655 10-06-2003, 12:50 PM
Last Post: Tom (UK)

Forum Jump: