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
|
|||
|
|||
Riemann's Zeta Function - another approach (RPL)
Code:
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
|
|||
|
|||
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:
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:
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
|
|||
|
|||
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
|
|||
|
|||
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: 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
|
|||
|
|||
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. What about Code:
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
|
|||
|
|||
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 (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 Or has this been simplified a bit too much ?-) Dieter |
|||
06-21-2017, 02:46 AM
Post: #7
|
|||
|
|||
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] ? 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
|
|||
|
|||
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. For 1 < x < 1.001: Code:
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
|
|||
|
|||
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 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 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
|
|||
|
|||
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: 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:
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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
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: Maybe even Gerson is afraid you will almost instantly find a way to shave a byte or speed it up... ![]() --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
|
|||
|
|||
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 ?-) 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
|
|||
|
|||
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. 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
|
|||
|
|||
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: Here it is: Code:
(06-27-2017 05:47 PM)Dieter Wrote: Finally, here are some optimized simple approximations for 1 < x ≤ 1,01. 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
|
|||
|
|||
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 :-) Indeed. An example in my listing above: Code:
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
|
|||
|
|||
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 ?-) 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
|
|||
|
|||
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: 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
|
|||
|
|||
RE: Riemann's Zeta Function - another approach (RPL)
HP-15C program update:
Code:
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)