The Museum of HP Calculators


Ramanujan Tau Numbers for the HP-41

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

-Ramanujan numbers Tau(n) are defined by

   x.[ (1-x).(1-x2).(1-x3)....... ]24  = Tau(1).x + Tau(2).x2 + Tau(3).x3 + ..... + Tau(n).xn + ......

-"TAU" calculates Tau(n) by the recurrence relation:

    (n-1).Tau(n) = SUM1<=m<=b(n) (-1)m+1 (2m+1).(n-1-9m(m+1)/2).Tau(n-m(m+1)/2)     where   b(n) = INT ((8n+1)1/2-1)/2

-However, this algorithm is very sensitive to roundoff errors and results are exact up to n = 43 only.
-Therefore, another approach is necessary to compute Tau(n) for larger arguments.
-"TAU2" employs the formula:

    Tau(n) = n4.s(n) - 24 SUM0<k<n k2.(35.k2-52.k.n+18.n2).s(k).s(n-k)     ( n > 1 )  where  s(k) is the sum of the divisors of k.

-Multiple precision is used and Tau(n) may be exactly computed up to n = 13868.
 

1°) Program#1   ( small arguments )
 

-This program calculates and stores Tau(1) , Tau(2) , ...... , Tau(n) into registers R01 , R02 , ....... , Rnn

Data Registers:  R00 = 0 , R01 = Tau(1) , ........... , Rnn = Tau(n)
Flags: /
Subroutines: /

-Synthetic registers M , N , O may be replaced by any unused data registers,
 for instance R97 , R98 , R99  if  n < 97
 

01  LBL "TAU"
02  STO O
03  .1
04  %
05  2
06  +
07  STO M
08  CLX
09  STO 00
10  SIGN
11  STO 01
12  LBL 01
13  RCL M
14  INT
15  8
16  *
17  RCL 01
18  +
19  SQRT
20  RCL 01
21  -
22  2
23  /
24  INT
25  STO N
26  CLX
27  FIX 0
28  LBL 02
29  RCL M
30  RCL 01
31  RCL N
32  ST+ Y
33  *
34  2
35  /
36  -
37  X<>Y
38  RCL IND Y
39  LASTX
40  9
41  *
42  RCL 01
43  +
44  RCL M
45  INT
46  -
47  *
48  RCL N
49  ST+ X
50  RCL 01
51  +
52  *
53  RCL 01
54  CHS
55  RCL N
56  Y^X
57  *
58  +
59  DSE N
60  GTO 02
61  RCL M
62  INT
63  RCL 01
64  -
65  /
66  RND
67  STO IND M
68  ISG M
69  GTO 01
70  RCL O
71  SIGN
72  RCL IND L
73  CLA
74  FIX 4
75  END

( 106 bytes / SIZE nnn+1, but at least 003 )
 
 
      STACK        INPUTS      OUTPUTS
           X             n          tau(n)
           L             /             n

Example:    Evaluate Tau(10)

   10  XEQ "TAU"  >>>>  Tau(10) = -115920  in X-register and in R10  ( execution time = 32 seconds )

-We also have   1    , -24 ,  252 , -1472 , 4830 , -6048 , -16744 , 84480 , -113643  ( i-e Tau(1) , .......... , Tau(9) )
    in registers   R01 , R02 , R03 ,   R04  ,  R05 ,    R06 ,     R07  ,    R08  ,     R09       respectively

Notes:   -Without lines 27 and 66  ( FIX 0 and RND ) , Tau(33) would already be wrong!
            -This program yields Tau(84) = 6,188,510,373 whereas the correct value is Tau(84) = 6,211,086,336
      and 300 XEQ "TAU" produces a completely meaningless result!
-Obviously, roundoff errors overwhelm the wanted function.
 

2°) Program#2   ( 1 < n < 13869 )
 

-The following program computes Tau(n) by groups of 4 digits in registers R15 R16 R17 R18 R19 R20
 

Data Registers:        R01 = n   ;   R00 and R02 thru R20 are used for temporary data storage and when the program stops:
                                   R03 = Tau(n) = X-register ( approximate value )
                          and   R15 to R20 contain the digits of Tau(n) ( exact value )
Flags: /
Subroutines: /
 
 

  01  LBL "TAU2"
  02  STO 00
  03  STO 01
  04   E4
  05  STO 02
  06  15.02
  07  STO 14
  08  CLRGX             if you don't have an HP-41CX , replace line 08 by  0  LBL 00  STO IND Y  ISG Y  GTO 00   or add  CLRG  after line 01
  09  DSE 00
  10  LBL 01
  11  RCL 00
  12  X^2
  13  ENTER^
  14  STO Z
  15  RCL 02
  16  MOD
  17  STO 04
  18  -
  19  RCL 02
  20  /
  21  STO 03
  22  CLX
  23  35
  24  *
  25  RCL 00
  26  RCL 01
  27  *
  28  52
  29  *
  30  -
  31  RCL 01
  32  X^2
  33  18
  34  *
  35  +
  36  STO Y
  37  RCL 02
  38  MOD
  39  STO 07
  40  -
  41  RCL 02
  42  /
  43  STO Y
  44  RCL 02
  45  MOD
  46  STO 06
  47  -
  48  RCL 02
  49  /
  50  STO 05
  51  RCL 00
  52  XEQ 05
  53  RCL 01
  54  RCL 00
  55  -
  56  XEQ 05
  57  *
  58  STO Y
  59  RCL 02
  60  MOD
  61  STO 09
  62  -
  63  RCL 02
  64  /
  65  STO 08
  66  XEQ 08
  67  DSE 00
  68  GTO 01
  69  RCL 14
  70  24
  71  CHS
  72  LBL 02
  73  ST* IND Y
  74  ISG Y
  75  GTO 02
  76  RCL 01
  77  X^2
  78  STO Y
  79  RCL 02
  80  MOD
  81  STO 04
  82  STO 07
  83  -
  84  RCL 02
  85  /
  86  STO 03
  87  STO 06
  88  CLX
  89  STO 05
  90  RCL 01
  91  XEQ 05
  92  STO Y
  93  RCL 02
  94  MOD
  95  STO 09
  96  -
  97  RCL 02
  98  /
  99  STO 08
100  XEQ 08
101  RCL 02
102  RCL 02
103  RCL 14
104  STO 00
105  CLX
106  LBL 03
107  RCL IND 00
108  +
109  *
110  ISG 00
111  GTO 03
112  LASTX
113  STO 03
114  SIGN
115  STO 04
116  ST* 02
117  20.015
118  STO 00
119  CLX
120  LBL 04
121  RCL IND 00
122  +
123  STO Y
124  RCL 02
125  MOD
126  STO IND 00
127  -
128  RCL 02
129  /
130  RCL 04
131  *
132  DSE 00
133  GTO 04
134  ST+ 15
135  RCL 14
136  RCL 03
137  RTN
138  LBL 05                 lines 138 to 161 calculates the sum of the divisors of the number in X-register.
139  STO 08
140  SQRT
141  INT
142  STO 09
143  CLX
144  LBL 06
145  RCL 08
146  RCL 09
147  /
148  FRC
149  X#0?
150  GTO 07
151  X<> L
152  ST+ Y
153  RCL 09
154  X#Y?
155  ST+ Z
156  RDN
157  LBL 07
158  RDN
159  DSE 09
160  GTO 06
161  RTN
162  LBL 08                 lines 162 to 207 calculates the product P = k2.(35.k2-52.k.n+18.n2) and stores the result in registers R10 thru R13
163  RCL 04
164  RCL 07
165  *
166  STO Y
167  RCL 02
168  MOD
169  STO 13
170  -
171  RCL 02
172  /
173  RCL 03
174  RCL 07
175  *
176  +
177  RCL 04
178  RCL 06
179  *
180  +
181  STO Y
182  RCL 02
183  MOD
184  STO 12
185  -
186  RCL 02
187  /
188  RCL 03
189  RCL 06
190  *
191  +
192  RCL 04
193  RCL 05
194  *
195  +
196  STO Y
197  RCL 02
198  MOD
199  STO 11
200  -
201  RCL 02
202  /
203  RCL 03
204  RCL 05
205  *
206  +
207  STO 10
208  RCL 09                 lines 208 to 275 computes the product of P by s(k).s(n-k) and add the result to the partial sum in registers R15 thru R20.
209  RCL 13
210  *
211  STO Y
212  RCL 02
213  MOD
214  ST+ 20
215  -
216  RCL 02
217  /
218  RCL 08
219  RCL 13
220  *
221  +
222  RCL 09
223  RCL 12
224  *
225  +
226  STO Y
227  RCL 02
228  MOD
229  ST+ 19
230  -
231  RCL 02
232  /
233  RCL 08
234  RCL 12
235  *
236  +
237  RCL 09
238  RCL 11
239  *
240  +
241  STO Y
242  RCL 02
243  MOD
244  ST+ 18
245  -
246  RCL 02
247  /
248  RCL 08
249  RCL 11
250  *
251  +
252  RCL 09
253  RCL 10
254  *
255  +
256  STO Y
257  RCL 02
258  MOD
259  ST+ 17
260  -
261  RCL 02
262  /
263  RCL 08
264  RCL 10
265  *
266  +
267  STO Y
268  RCL 02
269  MOD
270  ST+ 16
271  -
272  RCL 02
273  /
274  ST+ 15
275  RTN
276  END

( 352 bytes / SIZE 021 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /        15.020
           X             n   tau(n) ( approx.)

  15.020  = control number of the solution ( in registers R15 thru R20 )

Example1:     Calculate  Tau(84)

   84  XEQ "TAU2"  >>>>  Tau(84) =  6211086336 ( 13m07s )    In this case, the approximate value is quite exact!

which is confirmed by   R15 = 0 , R16 = 0 , R17 = 0 , R18 = 62 , R19 = 1108 , R20 = 6336

Example2:    Compute Tau(314)

  314  XEQ "TAU2"  >>>>  Tau(314) = -3.156280211 1013  = approximate value ( in 1h07mn20s )

and we have:  R15 = 0 , R16 = 0 , R17 = -31 , R18 = -5628 , R19 = -210 ( which is to be read -0210 ) , R20 = -5744

 whence   Tau(314) = -31562802105744

-Remark that registers R15 to R20 are inferior or equal to zero when Tau(n) < 0

Notes:   -With a good 41-emulator in turbo mode  13868 XEQ "TAU2"  yields:

   R15 = -336 ; R16 = -1948 ; R17 = -0538 ; R18 = -4247 ; R19 = -3535 ; R20 = -1552   ( in about 1h23mn on my PC )

  whence  Tau(13868) = -33619480538424735351552      with a "real" HP-41, execution time is probably of the order of  8 or 9 days!

-For n > 13868 , the calculations involve numbers greater than 1010 and the results are only approximate,
 but if we use more registers to perform multiplications and additions, this limit can be overcome.

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall