Riemann's Zeta Function - another approach (RPL)
06-17-2017, 01:37 AM (This post was last modified: 06-17-2017 01:41 AM by Gerson W. Barbosa.)
Post: #1
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
Riemann's Zeta Function - another approach (RPL)
Code:
 %%HP: T(3)A(D)F(.); \<< DUP DUP TYPE { C\->R MIN } IFT -1.3 ^ 178. * 1. + 2. / IP DUP + \-> s n   \<< '\GS(k=0.,n-1.,(-1.)^k/(k+1.)^s)' EVAL n .5 + s 1. + 8. n * / + s DUP + 1. + 27. n SQ * / - s ^ DUP + INV + 2. s ^ DUP 2. - / *   \>> \>>

0.5 -> -1.46035450880 (8.37 s)
1.0001 -> 10000.5772771 (3.53 s)
1.27 -> 4.30022020082 (2.64 s)
1.5 -> 2.61237534865 (2.18 s)
2 -> 1.64493406683 (1.54 s)
3 -> 1.20205690315 (1.00 s)
4 -> 1.08232323371 (0.78 s)
5 -> 1.03692775514 (0.64 s)
6 -> 1.01734306198 (0.57 s)
7 -> 1.00834927738 (0.49 s)
19.99 -> 1.0000009606 (0.31 s)
(2,3) -> (0.798021985125,-0.113744308033) (4.89 s)

(Tested on the HP 50g only).

This is based on an alternate series and two correction terms, the latter of which I am not so sure of. If a third correction term is found, both speed and accuracy for arguments close to 1 can be improved. The formula can be extracted from the listing, but I may included it later. (This is an afternoon's work and is still very immature - my original intention is a solution that would work on the HP-42S).

For a more complete and faster solution, with extended range, please refer to Riemann's Zeta Function update (HP-28S, HP-48G/GX/G+, HP-49G/G+/50g) (complex arguments will return accurate results only on a narrow strip, though).

Gerson.
06-17-2017, 06:07 PM (This post was last modified: 06-19-2017 04:47 PM by Gerson W. Barbosa.)
Post: #2
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
$\zeta (s)= \frac{2^{s}}{2^{s}-2}\sum_{k=0}^{\infty} \frac{(-1)^{k}}{(k+1)^{s}}$

$\zeta (s)\approx \frac{2^{s}}{2^{s}-2}\left ( \frac{1}{2\cdot \left (n+\frac{1}{2}+\frac{s+1}{8\cdot n}-\frac{2\cdot s+1}{27\cdot n^{2}} \right )^{s}} + \sum_{k=1}^{n} \frac{(-1)^{k}}{(n-k+1)^{s}} \right )$

If I were to guess, I'd say the third term (n and 1/2 are too trivial to count) is (3*s + 1)/(64*n^3), but I haven't tried it yet.

Here is a Q&D stack-only implementation:

Code:
 %%HP: T(3)A(D)F(.); \<< DUPDUP TYPE { C\->R MIN } IFT -1.3 ^ 178. * 1. + 2. / IP DUP + DUP2 0. SWAP 2. / 1.   FOR k OVER k DUP + DUP 1. - PICK3 ^ INV UNROT SWAP ^ INV SWAP UNROT - + -1.   STEP 4. ROLLD DUP + 1. + OVER SQ 27. * / PICK3 1. + PICK3 8. * / SWAP - SWAP .5 + + OVER ^ DUP + INV ROT + 2. ROT ^ DUP 2. - / * \>>

0.5 -> -1.46035450880 (6.17 s)
1.0001 -> 10000.5772771 (2.54 s)
1.27 -> 4.30022020082 (1.88 s)
1.5 -> 2.61237534865 (1.53 s)
2 -> 1.64493406684 (1.04 s)
3 -> 1.20205690315 (0.64 s)
4 -> 1.08232323371 (0.48 s)
5 -> 1.03692775514 (0.37 s)
6 -> 1.01734306199 (0.32 s)
7 -> 1.00834927738 (0.27 s)
12 -> 1.00024608656 (0.18 s)
19.99 -> 1.0000009606 (0.12 s)
(2,3) -> (0.798021985140,-0.113744308037) (4.22 s)

Update: it looks like the denominator of the second term is 24, not 27. The previous term, (s+1)/(8*n) is correct.

Now using

$\zeta (s)\approx \frac{2^{s}}{2^{s}-2}\left ( \frac{1}{2\cdot \left (n+\frac{1}{2}+\frac{s+1}{8\cdot n}-\frac{2\cdot s+1}{24\cdot n^{2}}+\frac{3\cdot s+1}{30\cdot n^{3}} \right )^{s}} + \sum_{k=1}^{n} \frac{(-1)^{k}}{(n-k+1)^{s}} \right )$

n should always be even.

Not sure about the last two correction terms, though. I'd need to do some calculations using 30+ SD in order to determine them.

Code:
 %%HP: T(3)A(D)F(.); \<< DUPDUP TYPE { C\->R MIN } IFT -1.27 ^ 157. * 1. + 2. / IP DUP + DUP2 0. SWAP 2. / 1.   FOR k OVER k DUP + DUP 1. - PICK3 ^ INV UNROT SWAP ^ INV SWAP UNROT - + -1.   STEP 4. ROLLD 3. * 1. + OVER 3. ^ 30. * / PICK3 DUP + 1. + PICK3 SQ 24. * / - PICK3 1. + PICK3 8. * / + .5 + + OVER ^ DUP + INV ROT + 2. ROT ^ DUP 2. - / * \>>

0.5 -> -1.46035450880 (5.36) [378]
1.0001 -> 10000.5772772 (2.24) [156]
1.27 -> 4.30022020086 (1.70) [116]
1.5 -> 2.61237534869 (1.39) [94]
2 -> 1.64493406686 (0.98) [66]
3 -> 1.20205690316 (0.60) [38]
4 -> 1.08232323371 (0.44) [26]
5 -> 1.03692775514 (0.36) [20]
6 -> 1.01734306198 (0.31) [16]
7 -> 1.00834927738 (0.28) [14]
12 -> 1.00024608655 (0.17) [6]
19.99 -> 1.0000009606 (0.14) [4]
(2,3) -> (0.798021985137,-0.113744308000) (3.90) [66]

(t): time in seconds
[n]: evaluated terms in the alternating series.

06-19-2017, 09:17 PM
Post: #3
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-17-2017 06:07 PM)Gerson W. Barbosa Wrote:  Here is a Q&D stack-only implementation:

Gerson, I just want to say: although until now there were no replies to your post this does not mean that it has not been noticed. ;-)

I assume you calculate the value for n with a, say, heuristic method (this –1,3 and 178 thing at the beginnig). Results with errors only in the last place are as good as it gets if you cannot use any additional guard digits. For s very close to 1, maybe you can switch to a simple implemenation, just with 1/s and the Euler-Masceroni constant.

Dieter
06-19-2017, 10:13 PM (This post was last modified: 06-19-2017 10:14 PM by Gerson W. Barbosa.)
Post: #4
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-19-2017 09:17 PM)Dieter Wrote:
(06-17-2017 06:07 PM)Gerson W. Barbosa Wrote:  Here is a Q&D stack-only implementation:

I assume you calculate the value for n with a, say, heuristic method (this –1,3 and 178 thing at the beginnig). Results with errors only in the last place are as good as it gets if you cannot use any additional guard digits. For s very close to 1, maybe you can switch to a simple implemenation, just with 1/s and the Euler-Masceroni constant.

Dieter, thank you for your comment and suggestion!

Yes, I prepared a table with twelve arguments ranging from 1.5 to 25 and the respective number of terms required for the maximum possible accuracy then I simply did a curve fitting of the data on the HP-42S. My goal is something that can work on slow calculators like the HP-11C and HP-15C. Complex arguments on the latter might be a bonus, albeit limited to a narrow strip. Zeta(2) requires the evaluation of only 66 terms of the alternating series (please take a look at the updated formula and program in my previous post) and even less on 10-digit calculators, also on these the last two correction terms, of which I am not sure about, might not be necessary. I think the equivalent program for classic scientific Voyagers would run in acceptable time for all arguments.

Gerson.
06-20-2017, 02:04 AM (This post was last modified: 06-20-2017 02:53 AM by Gerson W. Barbosa.)
Post: #5
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-19-2017 09:17 PM)Dieter Wrote:
(06-17-2017 06:07 PM)Gerson W. Barbosa Wrote:  Here is a Q&D stack-only implementation:
For s very close to 1, maybe you can switch to a simple implemenation, just with 1/s and the Euler-Masceroni constant.

Code:
 C01 LBL C C02 FP C03 x! C04 LASTx C05 x<>y C06 2 C07 x<>y C08 - C09 x<>y C10 / C11 LASTx C12 + C13 RTN

for 1 < x <= 1.00001 ?

This saves the eleven steps required by the constant on the HP-11C, for instance.

Gerson.

Edited to save four steps in the HP-32S II program.
06-20-2017, 05:04 PM (This post was last modified: 06-20-2017 05:55 PM by Dieter.)
Post: #6
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-19-2017 10:13 PM)Gerson W. Barbosa Wrote:  Yes, I prepared a table with twelve arguments ranging from 1.5 to 25 and the respective number of terms required for the maximum possible accuracy then I simply did a curve fitting of the data on the HP-42S.

Hehe... I used a similar approach a few years ago for the continued fraction expansion of the Normal distribution CDF, where the number of required terms was computed with a simple formula. This is the preferred method for large x, and the number of terms increases rapidly as x drops below 2...1,5. That's why for smaller arguments a series expansion is the better and faster method. So in a way this is quite similar to what we are discussing here.

(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  What about [snip code sample] ?

This can be done with two steps less: ;-)

Code:
C01 LBL C C02 FP C03 ENTER C04 1/x C05 2 C06 LastX C07 x! C08 - C09 * C10 + C11 RTN

(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  ...for 1 < x <= 1.00001 ?

I would do the switch at the point where both methods, the regular and the simplified one, have the same error. The simplified version should be close to 12 digits with x < 1,00001, or 10 digits for x < 1,0001.

(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  This saves the eleven steps required by the constant on the HP-11C, for instance.

On a ten digit calculator and 1 < x < 1,0001 the result has at most five decimals. So the EM-constant is not required to have more either.

Code:
01 LBL C 02 FP 03 1/x 04 , 05 5 06 7 07 7 08 2 09 2 10 + 11 RTN

Or has this been simplified a bit too much ?-)

Dieter
06-21-2017, 02:46 AM
Post: #7
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-20-2017 05:04 PM)Dieter Wrote:
(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  What about [snip code sample] ?

This can be done with two steps less: ;-)

Code:
C01 LBL C C02 FP C03 ENTER C04 1/x C05 2 C06 LastX C07 x! C08 - C09 * C10 + C11 RTN

I'm glad you missed my 17-step version :-)

(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  ...for 1 < x <= 1.00001 ?

I would do the switch at the point where both methods, the regular and the simplified one, have the same error. The simplified version should be close to 12 digits with x < 1,00001, or 10 digits for x < 1,0001.

(06-20-2017 02:04 AM)Gerson W. Barbosa Wrote:  This saves the eleven steps required by the constant on the HP-11C, for instance.

On a ten digit calculator and 1 < x < 1,0001 the result has at most five decimals. So the EM-constant is not required to have more either.

Code:
01 LBL C 02 FP 03 1/x 04 , 05 5 06 7 07 7 08 2 09 2 10 + 11 RTN

Same number of steps and faster. Sometimes I forget the good old K.I.S.S. rule. But that might be an alternative for 12-digit calculators.

I would like to see how an HP-41C implementation of what we have so far fares against existing solutions, but I'm not willing to do it myself. Maybe later.

Gerson.
06-22-2017, 01:20 AM
Post: #8
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-20-2017 05:04 PM)Dieter Wrote:  On a ten digit calculator and 1 < x < 1,0001 the result has at most five decimals. So the EM-constant is not required to have more either.

Code:
01 LBL C 02 FP 03 1/x 04 , 05 5 06 7 07 7 08 2 09 2 10 + 11 RTN

Or has this been simplified a bit too much ?-)

For 1 < x < 1.001:

Code:
 01 LBL C 02 FRAC 03 1/x 04 LASTx 05 13.73328 13 / 14 .5772157 22 + 23 + 24 RTN

1.009999999 GSB C -> 100.5779539
1.003456789 GSB C -> 289.8632757

(Source)

Gerson.
06-25-2017, 01:16 PM (This post was last modified: 06-27-2017 10:11 PM by Dieter.)
Post: #9
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  For 1 < x < 1.001:

I assume this is 1 < x < 1,01.

(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  1.009999999 GSB C -> 100.5779539
1.003456789 GSB C -> 289.8632757

I think we can top this. ;-)

Edit: the following suggestions are not yet optimized and can be improved further, please see my later post in this thread.

First of all, replace the 13,73328 with 13,7418. Evaluated exactly, this should keep the error within about ±0,7 ULP of a ten-digit result.

But we can do even better:

Zeta(x)  ~  1/u + u/(u + 13,733) + 0,57721568
where u = x–1 and 1 < x ≤ 1,01

Code:
01 LBL C 02 FRAC 03 1/x 04 LastX 05 LastX 06 13,733 07 + 08 / 09 ,57721568 10 + 11 + 12 RTN

1,009999999 => 100,5779533
1,009876543 => 101,8279365
1,003456789 => 289,8632756

According to some quick-and-dirty tests, within the given domain and with exact evaluation (!) the result should have ten significant digits ±0,2 ULP.

If you prefer the result to be either correctly rounded or truncated (error always negative) replace 0,57721568 with ...66:

Zeta (1,001001001) = 999,5772895 48

1,001001001=> 999,5772896 with ...68
1,001001001=> 999,5772895 with ...66

If that's still too much, try
Zeta(x) ~ 1/u + u/(0,9 · u + 13,7335) + 0,577215668

(resp. ...666, see above)

Here the largest error should be about ±0,02 ULP. If evaluated exactly, that is. But this improvement does not always show up on a ten-digit calculator. Actually the limiting factor seems to be the built-in 1/x function with its inherent error of 0,5 ULP.

Dieter
06-25-2017, 11:03 PM
Post: #10
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-25-2017 01:16 PM)Dieter Wrote:
(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  For 1 < x < 1.001:

I assume this is 1 < x < 1,01.

(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  1.009999999 GSB C -> 100.5779539
1.003456789 GSB C -> 289.8632757

I think we can top this. ;-)

First of all, replace the 13,73328 with 13,7418. Evaluated exactly, this should keep the error within about ±0,7 ULP of a ten-digit result.

But we can do even better:

Zeta(x)  ~  1/u + u/(u + 13,733) + 0,57721568
where u = x–1 and 1 < x ≤ 1,01

Code:
01 LBL C 02 FRAC 03 1/x 04 LastX 05 LastX 06 13,733 07 + 08 / 09 ,57721568 10 + 11 + 12 RTN

1,009999999 => 100,5779533
1,009876543 => 101,8279365
1,003456789 => 289,8632756

According to some quick-and-dirty tests, within the given domain and with exact evaluation (!) the result should have ten significant digits ±0,2 ULP.

If you prefer the result to be either correctly rounded or truncated (error always negative) replace 0,57721568 with ...66:

Zeta (1,001001001) = 999,5772895 48

1,001001001=> 999,5772896 with ...68
1,001001001=> 999,5772895 with ...66

If that's still too much, try
Zeta(x) ~ 1/u + u/(0,9 · u + 13,7335) + 0,577215668

(resp. ...666, see above)

Here the largest error should be about ±0,02 ULP. If evaluated exactly, that is. But this improvement does not always show up on a ten-digit calculator. Actually the limiting factor seems to be the built-in 1/x function with its inherent error of 0,5 ULP.

Thank you for your analysis and improvement. That's what I'll use if I ever try this on the HP-41C.

On the HP 50g I have used 5 Stieltjes constants for 1 < x <= 1.3. Since the size of the numbers doesn't matter much I've kept all constants in the equivalent Horner expression with twelve significant digits.

Code:
 %%HP: T(3)A(D)F(.); DIR   ZetaX   \<< RCLF SWAP RAD DUP NOT { .000000000001 + } IFT DUP .5 < { DUP \pi * 2. / SIN DUP + 1. ROT - SWAP OVER GAMMA * \pi DUP + PICK3 ^ / SWAP Zeta * } { Zeta } IFTE SWAP STOF   \>>   Zeta   \<< DUP DUP 1. - ABS .3 > { -1.27 ^ 157. * 1. + 2. / IP DUP + DUP2 0. SWAP 2. / 1.     FOR k OVER k DUP + DUP 1. - PICK3 ^ INV UNROT SWAP ^ INV SWAP UNROT - + -1.     STEP 4. ROLLD 3. * 1. + OVER 3. ^ 30. * / PICK3 DUP + 1. + PICK3 SQ 24. * / - PICK3 1. + PICK3 8. * / + .5 + + OVER ^ DUP + INV ROT + 2. ROT ^ DUP 2. - / * } { DUP DUP2 DUP2 -6.66328326918E-6 * 1.3020683574E-4 + * 7.96500246987E-4 - * 3.17028903723E-3 - * 8.10584133725E-2 + * .500000497261 + SWAP 1. - INV + NIP } IFTE   \>> END

Zeta: x > 1;

ZetaX: x <> 1.

Real arguments only.

We can use this to check whether 1 + 2 + 3 + 4 + 5 + ... = -1/12    :-)     (YouTube video)

1 +/- ZetaX -> -8.33333333338E-2

1/x -> -11.9999999999

Gerson.
06-27-2017, 01:26 AM
Post: #11
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
On the classic HP-15C:

1 CHS GSB B -> -0.083333333(28) (103 s)
0.02 GSB B -> -0.518788211(0) (11 s)
0 GSB -> -0.500000000 (1 s)
1.02 GSB A -> 50.5786700(5) (4 s)
1.04 GSB A -> 25.580120(64) (4 s)
2 GSB A -> 1.64493406(6) (99 s)
3 GSB A -> 1.202056902 (60 s)
10 GSB A -> 1.000994575 (19 s)
06-27-2017, 05:47 PM (This post was last modified: 06-27-2017 06:17 PM by Dieter.)
Post: #12
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 01:26 AM)Gerson W. Barbosa Wrote:  On the classic HP-15C:

Gerson, what about a program listing ?-)

Finally, here are some optimized simple approximations for 1 < x ≤ 1,01.

With u = x–1:

Zeta(x) ~ 1/u + u/13,7433 + 0,57721576
(error less than ±0,5 units in the 10th significant digit)

Zeta(x) ~ 1/u + u/(u + 13,73234) + 0,577215656
(error less than ±0,5 units in the 11th significant digit)

Zeta(x) ~ 1/u + u/(u + 13,73234) + 0,577215651
(error between 0 and –1 unit in the 11th significant digit)

Zeta(x) ~ 1/u + u/(0,9 · u + 13,733437) + 0,5772156664
(error less than ±1 unit in the 12th significant digit)

The mentioned error bounds assume exact evaluation, i.e. with more digits than the target accuracy. Otherwise the resulting errors may be slightly larger.

I tried the last approximation on the 35s and indeed in the results I got only the last digit was off by one here and there. Which does not mean that larger errors may occur due to roundoff errors in intermediate results.

Dieter
06-27-2017, 07:12 PM
Post: #13
 rprosperi Senior Member Posts: 2,915 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 05:47 PM)Dieter Wrote:
(06-27-2017 01:26 AM)Gerson W. Barbosa Wrote:  On the classic HP-15C:

Gerson, what about a program listing ?-)

Maybe even Gerson is afraid you will almost instantly find a way to shave a byte or speed it up... But please keep it up Dieter, I almost always learn something from your elegant improvements.

--Bob Prosperi
06-27-2017, 07:54 PM (This post was last modified: 06-27-2017 08:47 PM by Gerson W. Barbosa.)
Post: #14
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 07:12 PM)rprosperi Wrote:
(06-27-2017 05:47 PM)Dieter Wrote:  Gerson, what about a program listing ?-)

Maybe even Gerson is afraid you will almost instantly find a way to shave a byte or speed it up... But please keep it up Dieter, I almost always learn something from your elegant improvements.

Surely anyone will be able to shave not one, but lots of bytes :-)
Speed and size optimizations are not my concern for the time being. I was primarily interested to check whether my first HP calculator, the HP-15C (2343B75099) I bought back in the day, would be able to do it in finite time. 100 seconds to compute Zeta(2) is not fast enough, but I still can use my 15C, at least while speed optimizations are not available :-)

The 160-step listing isn't anything I should be proud of, but I'll publish it just the same. Rather than do it manually, I intend to obtain it from a simulator or emulator. I know there's at least one around with automatic listing capability. I only have to find it.

Dieter's optimizations – and anyone else's – are always welcome.

Gerson.

-------

PS: That's it: A Simulator for Windows, Linux and Mac OS X, by Torsten Manz
06-27-2017, 10:23 PM
Post: #15
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 07:54 PM)Gerson W. Barbosa Wrote:  Surely anyone will be able to shave not one, but lots of bytes :-)

That's the fun of it. ;-)

(06-27-2017 07:54 PM)Gerson W. Barbosa Wrote:  The 160-step listing isn't anything I should be proud of, but I'll publish it just the same. Rather than do it manually, I intend to obtain it from a simulator or emulator. I know there's at least one around with automatic listing capability. I only have to find it.
...
PS: That's it: A Simulator for Windows, Linux and Mac OS X, by Torsten Manz

I just realized that I also got Torsten's simulator. Here the help screen says "Version 8.5.9" while the website mentions a current version 3.4 – ?!?

Anyway, this is a really nice simulator. But it doesn't seem to be a true emulator of the original 15C microcode, so the results may not exactly match those of your hardware 15C.

(06-27-2017 07:54 PM)Gerson W. Barbosa Wrote:  Dieter's optimizations – and anyone else's – are always welcome.

Thank you very much.

Dieter
06-27-2017, 10:28 PM (This post was last modified: 06-28-2017 04:44 PM by Gerson W. Barbosa.)
Post: #16
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 05:47 PM)Dieter Wrote:
(06-27-2017 01:26 AM)Gerson W. Barbosa Wrote:  On the classic HP-15C:

Gerson, what about a program listing ?-)

Here it is:
Code:
 [code]# -------------------------------------------- # HEWLETT·PACKARD 15C Simulator program # Created with version 3.4.01 # -------------------------------------------- # --------------------------------------------    000 {             }     001 {    42 21 11 } f LBL A    002 {       44  3 } STO 3    003 {           1 } 1    004 {          30 } -    005 {       43 11 } g x²    006 {          11 } √x̅    007 {          48 } .    008 {           0 } 0    009 {           4 } 4    010 {          34 } x↔y    011 {       43 10 } g x≤y    012 {       22  1 } GTO 1    013 {       45  3 } RCL 3    014 {          36 } ENTER    015 {          36 } ENTER    016 {           1 } 1    017 {          48 } .    018 {           3 } 3    019 {          16 } CHS    020 {          14 } y^x    021 {           7 } 7    022 {           5 } 5    023 {          20 } ×    024 {           1 } 1    025 {          40 } +    026 {       43 44 } g INT    027 {          36 } ENTER    028 {          40 } +    029 {       44 25 } STO I    030 {       44  2 } STO 2    031 {          34 } x↔y    032 {          16 } CHS    033 {       44  1 } STO 1    034 {           0 } 0    035 {       44  0 } STO 0    036 {    42 21  0 } f LBL 0    037 {       45 25 } RCL I    038 {       45  1 } RCL 1    039 {          14 } y^x    040 {    44 30  0 } STO - 0    041 {           1 } 1    042 {    44 30 25 } STO - I    043 {       45 25 } RCL I    044 {       45  1 } RCL 1    045 {          14 } y^x    046 {    44 40  0 } STO + 0    047 {    42  5 25 } f DSE I    048 {       22  0 } GTO 0    049 {       45  1 } RCL 1    050 {          36 } ENTER    051 {          40 } +    052 {           1 } 1    053 {          30 } -    054 {       45  2 } RCL 2    055 {       43 11 } g x²    056 {           2 } 2    057 {           4 } 4    058 {          20 } ×    059 {          10 } ÷    060 {           1 } 1    061 {    45 30  1 } RCL - 1    062 {       45  2 } RCL 2    063 {           8 } 8    064 {          20 } ×    065 {          10 } ÷    066 {          48 } .    067 {           5 } 5    068 {          40 } +    069 {    45 40  2 } RCL + 2    070 {       45  1 } RCL 1    071 {          14 } y^x    072 {           2 } 2    073 {          10 } ÷    074 {    45 40  0 } RCL + 0    075 {           2 } 2    076 {       45  1 } RCL 1    077 {          16 } CHS    078 {          14 } y^x    079 {          36 } ENTER    080 {          36 } ENTER    081 {           2 } 2    082 {          30 } -    083 {          10 } ÷    084 {          20 } ×    085 {       43 32 } g RTN    086 {    42 21  1 } f LBL 1    087 {       45  3 } RCL 3    088 {           1 } 1    089 {          30 } -    090 {          15 } 1/x    091 {       43 36 } g LSTx    092 {          48 } .    093 {           9 } 9    094 {       43 36 } g LSTx    095 {          20 } ×    096 {           1 } 1    097 {           3 } 3    098 {          48 } .    099 {           7 } 7    100 {           3 } 3    101 {           3 } 3    102 {           4 } 4    103 {           4 } 4    104 {          40 } +    105 {          10 } ÷    106 {          48 } .    107 {           5 } 5    108 {           7 } 7    109 {           7 } 7    110 {           2 } 2    111 {           1 } 1    112 {           5 } 5    113 {           6 } 6    114 {           7 } 7    115 {          40 } +    116 {          40 } +    117 {       43 32 } g RTN    118 {    42 21 12 } f LBL B    119 {          48 } .    120 {           5 } 5    121 {          34 } x↔y    122 {    43 30  0 } g TEST x≠0    123 {       22  2 } GTO 2    124 {          34 } x↔y    125 {          16 } CHS    126 {       43 32 } g RTN    127 {    42 21  2 } f LBL 2    128 {       43 10 } g x≤y    129 {       22  3 } GTO 3    130 {       32 11 } GSB A    131 {       43 32 } g RTN    132 {    42 21  3 } f LBL 3    133 {           1 } 1    134 {          34 } x↔y    135 {          30 } -    136 {       44  4 } STO 4    137 {       32 11 } GSB A    138 {       43 26 } g π    139 {          36 } ENTER    140 {          40 } +    141 {       45  4 } RCL 4    142 {          14 } y^x    143 {          10 } ÷    144 {           1 } 1    145 {    45 30  4 } RCL - 4    146 {       43 26 } g π    147 {          20 } ×    148 {           2 } 2    149 {          10 } ÷    150 {       43  8 } g RAD    151 {          23 } SIN    152 {          20 } ×    153 {           1 } 1    154 {          16 } CHS    155 {    45 40  4 } RCL + 4    156 {       42  0 } f x!    157 {          20 } ×    158 {          36 } ENTER    159 {          40 } +    160 {       43 32 } g RTN # --------------------------------------------

(06-27-2017 05:47 PM)Dieter Wrote:  Finally, here are some optimized simple approximations for 1 < x ≤ 1,01.

With u = x–1:

...

Zeta(x) ~ 1/u + u/(0,9 · u + 13,733437) + 0,5772156664
(error less than ±1 unit in the 12th significant digit)

The mentioned error bounds assume exact evaluation, i.e. with more digits than the target accuracy. Otherwise the resulting errors may be slightly larger.

That's what I'd been using, with your previous constants. As you've pointed out, we'd need more digits to take full advantage of that. Two extra digits as on the HP-12C Platinum would be nice. HP not truncating intermediate results to the number of digits in the display, at least when running a program, would also help. The simulator gives more accurate results, however.

Thanks again for your valuable suggestions and improvements!

Gerson.

Edited to remove attached file.
06-27-2017, 11:57 PM (This post was last modified: 06-28-2017 12:00 AM by Gerson W. Barbosa.)
Post: #17
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 10:23 PM)Dieter Wrote:
(06-27-2017 07:54 PM)Gerson W. Barbosa Wrote:  Surely anyone will be able to shave not one, but lots of bytes :-)

That's the fun of it. ;-)

Indeed. An example in my listing above:

Code:
 005 {       43 11 } g x² 006 {          11 } √x̅

That's simply ABS, which the HP-15C of course doesn't lack. There are probably more obvious silly mistakes like this one, but I won't spoil the fun :-)

Gerson.
06-28-2017, 12:21 PM
Post: #18
 Dieter Senior Member Posts: 2,177 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-27-2017 10:28 PM)Gerson W. Barbosa Wrote:
(06-27-2017 05:47 PM)Dieter Wrote:  Gerson, what about a program listing ?-)

Here it is:

Lunch break – time for a closer look at your 15C program.

Line 007 ff: the simple approximation is done for 0,96 < s < 1,04. Why did you choose these limits?

Line 016 ff: the number of loops n seems to be calculated with the same or a similar formula like the one you used for the 50g, i.e. for 12 digit precision. For the 10-digit 15C less iterations should be sufficient. What do you think?

Line 065...066: it looks like there is a "+" missing between these lines.

The program stores the argument s in R3 and –s in R1. Afterwards only R1 is used, so R3 can be omitted.

And finally, what does the routine at LBL B do?

Maybe I'll try a 35s version. ;-)

Dieter
06-28-2017, 02:00 PM (This post was last modified: 06-28-2017 02:05 PM by Gerson W. Barbosa.)
Post: #19
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(06-28-2017 12:21 PM)Dieter Wrote:
(06-27-2017 10:28 PM)Gerson W. Barbosa Wrote:  Here it is:

Lunch break – time for a closer look at your 15C program.

...

Line 065...066: it looks like there is a "+" missing between these lines.

...

And finally, what does the routine at LBL B do?

Maybe I'll try a 35s version. ;-)

Indeed there is a "+" missing between lines 065 and 065. I think this might be the answer to all your other questions. When the parameters of the curve fitting expression are properly changed, I estimate a 3x speed-up factor in the computation of Zeta(2). Thank you very much!

A: Zeta(x), x > 1

B: Zeta(x), x<>1

It'll be fast enough on the 35s, despite the need of the evaluation of more terms.

Gerson.
06-28-2017, 06:07 PM (This post was last modified: 06-28-2017 07:19 PM by Gerson W. Barbosa.)
Post: #20
 Gerson W. Barbosa Senior Member Posts: 1,056 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
HP-15C program update:

Code:
 # -------------------------------------------- # HEWLETT·PACKARD 15C Simulator program # Created with version 3.4.01 # -------------------------------------------- # --------------------------------------------    000 {             }     001 {    42 21 11 } f LBL A    002 {       44  3 } STO 3    003 {           1 } 1    004 {          30 } -    005 {       43 16 } g ABS    006 {          48 } .    007 {           0 } 0    008 {           4 } 4    009 {          34 } x↔y    010 {       43 10 } g x≤y    011 {       22  1 } GTO 1    012 {       45  3 } RCL 3    013 {          36 } ENTER    014 {          36 } ENTER    015 {          48 } .    016 {           8 } 8    017 {           2 } 2    018 {          16 } CHS    019 {          14 } y^x    020 {           2 } 2    021 {           5 } 5    022 {          20 } ×    023 {       43 44 } g INT    024 {          36 } ENTER    025 {          40 } +    026 {       44 25 } STO I    027 {       44  2 } STO 2    028 {          34 } x↔y    029 {          16 } CHS    030 {       44  1 } STO 1    031 {           0 } 0    032 {       44  0 } STO 0    033 {    42 21  0 } f LBL 0    034 {       45 25 } RCL I    035 {       45  1 } RCL 1    036 {          14 } y^x    037 {    44 30  0 } STO - 0    038 {           1 } 1    039 {    44 30 25 } STO - I    040 {       45 25 } RCL I    041 {       45  1 } RCL 1    042 {          14 } y^x    043 {    44 40  0 } STO + 0    044 {    42  5 25 } f DSE I    045 {       22  0 } GTO 0    046 {       45  1 } RCL 1    047 {          36 } ENTER    048 {          40 } +    049 {           1 } 1    050 {          30 } -    051 {       45  2 } RCL 2    052 {       43 11 } g x²    053 {           2 } 2    054 {           4 } 4    055 {          20 } ×    056 {          10 } ÷    057 {           1 } 1    058 {    45 30  1 } RCL - 1    059 {           8 } 8    060 {    45 20  2 } RCL × 2    061 {          10 } ÷    062 {          40 } +    063 {          48 } .    064 {           5 } 5    065 {          40 } +    066 {    45 40  2 } RCL + 2    067 {       45  1 } RCL 1    068 {          14 } y^x    069 {           2 } 2    070 {          10 } ÷    071 {    45 40  0 } RCL + 0    072 {           2 } 2    073 {       45  1 } RCL 1    074 {          16 } CHS    075 {          14 } y^x    076 {          36 } ENTER    077 {          36 } ENTER    078 {           2 } 2    079 {          30 } -    080 {          10 } ÷    081 {          20 } ×    082 {       43 32 } g RTN    083 {    42 21  1 } f LBL 1    084 {       45  3 } RCL 3    085 {           1 } 1    086 {          30 } -    087 {          15 } 1/x    088 {       43 36 } g LSTx    089 {          48 } .    090 {           9 } 9    091 {       43 36 } g LSTx    092 {          20 } ×    093 {           1 } 1    094 {           3 } 3    095 {          48 } .    096 {           7 } 7    097 {           3 } 3    098 {           3 } 3    099 {           4 } 4    100 {           4 } 4    101 {          40 } +    102 {          10 } ÷    103 {          48 } .    104 {           5 } 5    105 {           7 } 7    106 {           7 } 7    107 {           2 } 2    108 {           1 } 1    109 {           5 } 5    110 {           6 } 6    111 {           7 } 7    112 {          40 } +    113 {          40 } +    114 {       43 32 } g RTN    115 {    42 21 12 } f LBL B    116 {          48 } .    117 {           5 } 5    118 {          34 } x↔y    119 {    43 30  0 } g TEST x≠0    120 {       22  2 } GTO 2    121 {          34 } x↔y    122 {          16 } CHS    123 {       43 32 } g RTN    124 {    42 21  2 } f LBL 2    125 {       43 10 } g x≤y    126 {       22  3 } GTO 3    127 {       32 11 } GSB A    128 {       43 32 } g RTN    129 {    42 21  3 } f LBL 3    130 {           1 } 1    131 {          34 } x↔y    132 {          30 } -    133 {       44  4 } STO 4    134 {       32 11 } GSB A    135 {       43 26 } g π    136 {          36 } ENTER    137 {          40 } +    138 {       45  4 } RCL 4    139 {          14 } y^x    140 {          10 } ÷    141 {           1 } 1    142 {    45 30  4 } RCL - 4    143 {           1 } 1    144 {          16 } CHS    145 {       43 24 } g COS-¹    146 {          20 } ×    147 {           2 } 2    148 {          10 } ÷    149 {          23 } SIN    150 {          20 } ×    151 {           1 } 1    152 {          16 } CHS    153 {    45 40  4 } RCL + 4    154 {       42  0 } f x!    155 {          20 } ×    156 {          36 } ENTER    157 {          40 } +    158 {       43 32 } g RTN # --------------------------------------------

1 CHS GSB B -> -0.08333333332 (52 s)
0.5 GSB B -> -1.460354506 (143 s)
0.02 GSB B -> -0.5187882110 (11 s)
0.959 GSB B -> -23.81602182 (83 s)
0.961 GSB B -> -25.06665702 (5 s)
0 GSB B -> -0.500000000 (1 s)
1.02 GSB A -> 50.57867005 (4 s)
1.04 GSB A -> 25.580120(64) (4 s)
2 GSB A -> 1.644934067 (48 s)
3 GSB A -> 1.202056903 (36 s)
4 GSB A -> 1.082323234 (31 s)
10 GSB A -> 1.000994575 (16 s)

PS: Edited to include formula used in program B:

$\zeta(x)=\frac{2\cdot (-x)!\cdot \sin \left ( \frac{x\cdot \cos^{-1}\left ( -1 \right )}{2} \right )\zeta (1-x)}{(2 \pi )^{1-x}}$
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)