Riemann's Zeta Function - another approach (RPL)
08-07-2017, 08:52 AM (This post was last modified: 08-08-2017 06:48 PM by Dieter.)
Post: #81
 Dieter Senior Member Posts: 1,640 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(08-04-2017 07:26 PM)Dieter Wrote:  Edit: here is the current version. It addresses the x~0 problem in a ...err... "creative" way in that the smallest |x| is set to 1E–20. #-)

There is a better solution.

Since $$\zeta '(0) = -\ln\sqrt{2 \pi}$$, for x close to zero  $$\zeta(x)\ \approx -(0.5+x\cdot\ln\sqrt{2 \pi})$$. So the program may simply switch to this as soon as |x| is sufficiently small. For the given program a good threshold is at 2 E–11. So if |x| is larger than this, use the regular method, and all |x| up to this limit may use the simple formula above. This is what I now have implemented in my Free42 program.

Here is the modified listing. The changes are in the first lines before the XEQ 77 call, and the STO 00 after LBL 77 has been omitted.
Edit: changed the threshold in line 08 from 1E–11 to 2E–11.

Code:
 00 { 434-Byte Prgm }  01>LBL "ZETA"  02 STO 00  03 0.5  04 X<>Y  05 X>=Y?  06 GTO 77  07 ABS  08 2E-11  09 X<Y?  10 GTO 00  11 PI  12 STO+ ST X  13 SQRT  14 LN  15 RCL× 00  16 0.5  17 +  18 +/-  19 GTO 99  20>LBL 00  21 1  22 RCL- 00  23 STO 00  24 XEQ 77  25 PI  26 STO+ ST X  27 RCL 00  28 Y^X  29 ÷  30 RCL 00  31 GAMMA  32 ×  33 1  34 ASIN  35 RCL× 00  36 COS  37 ×  38 STO+ ST X  39 1  40 RCL- 00  41 STO 00  42 Rv  43 GTO 99  44>LBL 77  45 -1  46 RCL+ 00  47 1/X  48 LASTX  49 X<0?  50 GTO 97  51 2  52 RCL 00  53 X>Y?  54 GTO 96  55 LASTX  56 LASTX  57 LASTX  58 -1.276E-8  59 ×  60 7.05133E-6  61 -  62 ×  63 9.721157E-5  64 +  65 ×  66 3.4243368E-4  67 -  68 ×  69 0.00484515482  70 -  71 ×  72 0.07281584288  73 +  74 ×  75 7.215664988E-3  76 +  77 GTO 98  78>LBL 96  79 100  80 RCL 00  81 SQRT  82 RCL× 00  83 ÷  84 5  85 +  86 IP  87 STO+ ST X  88 STO 01  89 RCL 00  90 +/-  91 STO 00  92 CLX  93>LBL 95  94 RCL ST Y  95 RCL 00  96 Y^X  97 -  98 +/-  99 DSE ST Y 100 GTO 95 101 RCL 00 102 STO+ ST X 103 1 104 - 105 RCL 01 106 X^2 107 24 108 × 109 ÷ 110 1 111 RCL- 00 112 8 113 ÷ 114 RCL÷ 01 115 + 116 0.5 117 + 118 RCL+ 01 119 RCL 00 120 Y^X 121 2 122 ÷ 123 + 124 1 125 RCL 00 126 +/- 127 STO 00 128 - 129 2 130 LN 131 × 132 E^X-1 133 +/- 134 ÷ 135 4 136 RCL- 00 137 44E-13 138 × 139 X<0? 140 CLX 141 - 142 GTO 99 143>LBL 97 144 ENTER 145 ENTER 146 ENTER 147 -5.872675E-6 148 × 149 9.740462E-5 150 + 151 × 152 3.4213951E-4 153 - 154 × 155 0.00484515533 156 - 157 × 158 0.07281584734 159 + 160 × 161 7.21566494432E-3 162 + 163>LBL 98 164 0.57 165 + 166 X<>Y 167 1/X 168 + 169>LBL 99 170 RCL 00 171 X<>Y 172 END

And again: please note that this version is intended for high-precision environments like Free42 and may not perform well on a regular 42s.

Dieter
08-08-2017, 04:05 AM
Post: #82
 Gerson W. Barbosa Senior Member Posts: 874 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(08-07-2017 08:52 AM)Dieter Wrote:
(08-04-2017 07:26 PM)Dieter Wrote:  Edit: here is the current version. It addresses the x~0 problem in a ...err... "creative" way in that the smallest |x| is set to 1E–20. #-)

There is a better solution.

Since $$\zeta '(0) = -\ln\sqrt{2 \pi}$$, for x close to zero  $$\zeta(x)\ \approx -(0.5+x\cdot\ln\sqrt{2 \pi})$$. So the program may simply switch to this as soon as |x| is sufficiently small. For the given program a good threshold is at about 1E–11. So if |x| is larger than this, use the regular method, and all |x| up to this limit may use the simple formula above. This is what I now have implemented in my Free42 program.

Here is the modified listing. The changes are in the first lines before the XEQ 77 call, and the STO 00 after LBL 77 has been omitted.

Code:
 00 { 434-Byte Prgm }  01>LBL "ZETA"  02 STO 00  03 0.5  04 X<>Y  05 X>=Y?  06 GTO 77  07 ABS  08 1E-11  09 X<Y?  10 GTO 00  11 PI  12 STO+ ST X  13 SQRT  14 LN  15 RCL× 00  16 0.5  17 +  18 +/-  19 GTO 99  20>LBL 00  21 1  22 RCL- 00  23 STO 00  24 XEQ 77  25 PI  26 STO+ ST X  27 RCL 00  28 Y^X  29 ÷  30 RCL 00  31 GAMMA  32 ×  33 1  34 ASIN  35 RCL× 00  36 COS  37 ×  38 STO+ ST X  39 1  40 RCL- 00  41 STO 00  42 Rv  43 GTO 99  44>LBL 77  45 -1  46 RCL+ 00  47 1/X  48 LASTX  49 X<0?  50 GTO 97  51 2  52 RCL 00  53 X>Y?  54 GTO 96  55 LASTX  56 LASTX  57 LASTX  58 -1.276E-8  59 ×  60 7.05133E-6  61 -  62 ×  63 9.721157E-5  64 +  65 ×  66 3.4243368E-4  67 -  68 ×  69 0.00484515482  70 -  71 ×  72 0.07281584288  73 +  74 ×  75 7.215664988E-3  76 +  77 GTO 98  78>LBL 96  79 100  80 RCL 00  81 SQRT  82 RCL× 00  83 ÷  84 5  85 +  86 IP  87 STO+ ST X  88 STO 01  89 RCL 00  90 +/-  91 STO 00  92 CLX  93>LBL 95  94 RCL ST Y  95 RCL 00  96 Y^X  97 -  98 +/-  99 DSE ST Y 100 GTO 95 101 RCL 00 102 STO+ ST X 103 1 104 - 105 RCL 01 106 X^2 107 24 108 × 109 ÷ 110 1 111 RCL- 00 112 8 113 ÷ 114 RCL÷ 01 115 + 116 0.5 117 + 118 RCL+ 01 119 RCL 00 120 Y^X 121 2 122 ÷ 123 + 124 1 125 RCL 00 126 +/- 127 STO 00 128 - 129 2 130 LN 131 × 132 E^X-1 133 +/- 134 ÷ 135 4 136 RCL- 00 137 44E-13 138 × 139 X<0? 140 CLX 141 - 142 GTO 99 143>LBL 97 144 ENTER 145 ENTER 146 ENTER 147 -5.872675E-6 148 × 149 9.740462E-5 150 + 151 × 152 3.4213951E-4 153 - 154 × 155 0.00484515533 156 - 157 × 158 0.07281584734 159 + 160 × 161 7.21566494432E-3 162 + 163>LBL 98 164 0.57 165 + 166 X<>Y 167 1/X 168 + 169>LBL 99 170 RCL 00 171 X<>Y 172 END

Very nice!

I would only make a minor modification. Not that one less square root evaluation really might make a difference in execution time in Free42:

Code:
 … 13 LN 14 RCL× 00 15 RCL× ST T 16 RCL+ ST T 17 +/- …

Gerson.
08-08-2017, 06:46 PM (This post was last modified: 08-08-2017 06:58 PM by Dieter.)
Post: #83
 Dieter Senior Member Posts: 1,640 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(08-08-2017 04:05 AM)Gerson W. Barbosa Wrote:  Very nice!

Wait, it even gets a tiny bit nicer. ;-)

(08-08-2017 04:05 AM)Gerson W. Barbosa Wrote:  I would only make a minor modification. Not that one less square root evaluation really might make a difference in execution time in Free42:

Code:
… 13 LN 14 RCL× 00 15 RCL× ST T 16 RCL+ ST T 17 +/- …

Ah, great. It's not beause of the sqrt, but I love this elegant stack content reuse of the 0.5 that still sits in T. ;-)

But there even is one more tweak. The second term of the Zeta Taylor series is quite exactly –x², i.e.  $$\zeta(x)\ \approx -(0.5+x^2+x\cdot\ln\sqrt{2 \pi})$$.
And this can be implemented with just one more step:

Code:
… 13 LN 14 RCL× ST T 15 RCL+ 00 16 RCL× 00 17 RCL+ ST T 18 +/- …

This way the threshold in line 08 should be somewhere near 1 E–8. The perfect value differs slightly for negative resp. positive x, but 8 E–9 works very nicely for both.

BTW, if you want to stick to the original implementation the perfect switchpoint seems to be close to 2 E–11, so this value should be used there. I edited my original post accordingly. But we're talking about differences beyond the 20th digit here. ;-)

Dieter
08-12-2017, 07:03 PM
Post: #84
 Gerson W. Barbosa Senior Member Posts: 874 Joined: Dec 2013
RE: Riemann's Zeta Function - another approach (RPL)
(08-04-2017 07:26 PM)Dieter Wrote:  [quote='Gerson W. Barbosa' pid='76831' dateline='1501861065']
Back to the HP-41:
Code:
 01▸LBL "ZETA"  02 X>0?  03 GTO 00  04 X=0?  05 GTO 00  ...

Just one remark: Testing condition1 OR condition2 can be done with a simple trick from the olden days. Simply have the inverse of condition1 followed by condition2.

Code:
 01▸LBL "ZETA"  02 X≠0?  03 X>0?  04 GTO 00  ...

And the calculation of the number of terms should be replaced by this:

Code:
 64▸LBL 96  65 24  66 RCL 00  67 /  68 4  69 +  70 INT  71 ST+ X  72 STO 01

Just for the record, here is the new HP-41 listing with your suggestions:

Code:
  01▸LBL "ZETA"  02 X<=0?  03 X=0?  04 GTO 00  05 CHS  06 1  07 +  08 STO 02  09 XEQ "GAMMA"  10 STO 03  11 RCL 02  12 XEQ 00  13 RCL 03  14 *  15 PI  16 ST+ X  17 RCL 02  18 Y^X  19 /  20 1  21 ASIN  22 RCL 02  23 *  24 COS  25 *  26 ST+ X  27 GTO 99  28▸LBL 00  29 STO 00  30 1  31 -  32 1/X  33 LASTX  34 X<0?  35 GTO 97  36 2  37 RCL 00  38 X>Y?  49 GTO 96  40 LASTX  41 LASTX  42 LASTX  43 -1.276E -8  44 *  45 7.05133E -6  46 -  47 *  48 9.721157E -5  49 +  50 *  51 3.4243368E -4  52 -  53 *  54 4.84515482E -3  55 -  56 *  57 7.281584288E -2  58 +  59 *  60 7.215664988E -3      61 +  62 GTO 98  63▸LBL 96  64 24  65 RCL 00  66 /  67 4  68 +  69 INT  70 ST+ X  71 STO 01  72 RCL 00  73 CHS  74 STO 00  75 CLX  76▸LBL 01  77 RCL Y  78 RCL 00  79 Y^X  80 -  81 CHS  82 DSE Y  83 GTO 01  84 RCL 00  85 ST+ X  86 1  87 -  88 RCL 01  89 X^2  90 24  91 *  92 /  93 1  94 RCL 00  95 -  96 8  97 /  98 RCL 01  99 / 100 + 101 .5 102 + 103 RCL 01 104 + 105 RCL 00 106 Y^X 107 2 108 / 109 + 110 1 111 RCL 00 112 + 113 2 114 LN 115 * 116 E^X-1 117 CHS 118 / 119 GTO 99 120▸LBL 97 121 ENTER^ 122 ENTER^ 123 ENTER^ 124 -8.4715E -7 125 * 126 7.51334E -6 127 - 128 * 129 9.609657E -5 130 + 131 * 132 3.42683396E -4 133 - 134 * 135 4.84527616E -3 136 - 137 * 138 7.281583446E -2 139 + 140 * 141 7.215664464E -3 142 + 143▸LBL 98 144 RDN 145 1/X 146 INT 147 ST* Z 148 SIGN 149 STO/ X 150 STO- Z 151 X<> L 152 RDN 153 / 154 - 155 R^ 156 .57 157 + 158 + 159▸LBL 99 160 END

The primitive Gamma implementation I've been using is slightly faster:

Code:
  01▸LBL "GAMMA"  02 -1  03 X<>Y  04 +  05 1  06 STO 00  07 X>Y?  08 GTO 01  09 ST- L  10 LASTX  11 STO 00  12 INT  13 1  14▸LBL 00  15 RCL 00  16 *  17 DSE 00  18 ABS  19 DSE Y  20 GTO 00  21 STO 00  22▸LBL 01  23 LASTX  24 ENTER^  25 ENTER^  26 ENTER^  27 1.604589926 E-2  28 *  29 2.641400819 E-1  30 -  31 *  32 1.96580515  33 +  34 *  35 8.729327662  36 -  37 *  38 25.69590147  39 +  40 *  41 52.63472435  42 -  43 *  44 76.53826433  45 +  46 *  47 78.92639573  48 -  49 *  50 56.50761084  51 +  52 *  53 26.37243835  54 -  55 *  56 7.203398486  57 +  58 RCL 00  59 *  60 END

And here is the new table:

41 XEQ ZETA  -->   1.000000000           ( 8.24 s)
25     R/S   -->   1.000000030           ( 8.28 s)
3     R/S   -->   1.202056903           (17.89 s)
2.001 R/S   -->   1.643997513       (2) (22.33 s)
2     R/S   -->   1.644934067           ( 4.04 s)
1.5   R/S   -->   2.612375349           ( 4.08 s)
0.5   R/S   -->  -1.460354509           ( 4.74 s)
0     R/S   -->  -0.500000000           ( 4.03 s)
-0.5   R/S   -->  -0.2078862450    (250) ( 9.72 s)
-1     R/S   -->  -0.08333333384    (33) (10.07 s)
-1.001 R/S   -->  -0.08316803696   (723) (28.32 s)
-1.5   R/S   -->  -0.02548520436   (190) (25.38 s)
-2     R/S   -->   0.00000000000         (24.84 s)
-3     R/S   -->   0.008333333384   (33) (22.09 s)
-5     R/S   -->  -0.003968253990   (68) (20.24 s)
-7     R/S   -->   0.004166666686   (67) (19.39 s)
-15.16  R/S   -->   0.4964873534     (85) (18.99 s)
-33.34  R/S   -->  -1.924684098E10  (152) (22.15 s)
-41.42  R/S   -->  -3.506595602E16  (584) (24.03 s)
-48.49  R/S   -->  -3.653091058E22   (22) (25.98 s)
-58.59  R/S   -->   1.136304829E32  (792) (28.54 s)
-67.97  R/S   -->   1.832460467E40 (1182) (31.02 s)

Times obtained from the HP-41CX built-in stopwatch
 « Next Oldest | Next Newest »

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