Riemann's Zeta Function - another approach (RPL)

[Gerson W. Barbosa]

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.

[Gerson W. Barbosa]

\[\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)

Edited to correct a sign in the second formula.

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.

[Dieter]

(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

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

(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.

What about

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.

[Dieter]

(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

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

(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.

[Dieter]

(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

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

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)

[Dieter]

(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

[rprosperi]

(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.

[Gerson W. Barbosa]

(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

[Dieter]

(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

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

(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.

[Dieter]

(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

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

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}}\]