(49g 50g) Ramanujan Tau Function

[John Keith]

These programs compute the Ramanujan tau function (A000594) for positive integers. All three programs use Gerald Hilliers SUMDIVISORS command from hpcalc.org, and must be run in exact mode. The last two programs also require the ListExt Library. Neither program is very fast, and use of an emulator is recommended for large integers.

Edited 10/13/2018 to replace the first two programs with slightly smaller and faster versions.

The first program is the most general, and can calculate tau(n) for any size integer, given enough time. This is essentially the same algorithm as program #2 here.

Code:


\<<
  IF DUP 1 >
  THEN DUPDUP * \-> n n2
    \<< 0 1 n 1 -
      FOR k k DUPDUP 35 * 52 n * - OVER * 18 n2 * + * * n k DUP 1 SUMDIVISORS UNROT - 1 SUMDIVISORS * * +
      NEXT 24 * n DUP 1 SUMDIVISORS SWAP 4 ^ * SWAP -
    \>>
  END
\>>

The next program is about 25% faster than the previous one. It is limited to integers less than 2000 or so because of the memory required by the precomputed sigma values.

Code:


\<<
  IF DUP 1 >
  THEN DUPDUP * \-> n n2
    \<< n 1 - DUP LSEQ 1 SUMDIVISORS DUP REV * 1 ROT
      FOR k k DUPDUP 35 * 52 n * - OVER * 18 n2 * + * *
      NEXT n 1 - \->LIST * LSUM 24 * n DUP 1 SUMDIVISORS SWAP 4 ^ * SWAP -
    \>>
  END
\>>

The last program returns a list of tau values from 1 through n. It is significantly faster than computing the individual values with the programs above, but still quite slow for large integers.

Code:


\<< DUP 1 - DUP LSEQ 1 SUMDIVISORS \-> n n1 s
  \<< 1 1 n1
    FOR k k DUPN k \->LIST REV s 1 k SUB * LSUM -24 * k /
    NEXT n \->LIST
  \>>
\>>

EDIT: The program above can be used to compute many other sequences if the -24 in line 3 is replaced by another number. Details and a more general program here.

[Gerald H]

If you only want the sum of 1st power of divisors this programme is more economical:

Size: 117.

CkSum: # 87h

Code:

::
  CK1&Dispatch
  # FF
  ::
    FPTR2 ^ZAbs
    FPTR2 ^DupQIsZero?
    caseSIZEERR
    FPTR2 ^DupZIsOne?
    ?SEMI
    FPTR2 ^MZSQFF
    #2/
    ZINT 1
    SWAP
    ZERO_DO
    3PICK
    ROT
    COERCE
    #1+
    FPTR2 ^PPow#
    ZINT 1
    FPTR2 ^QSub
    ROT
    ZINT 1
    FPTR2 ^QSub
    FPTR2 ^ZQUOText
    FPTR2 ^QMul
    LOOP
  ;
;

[John Keith]

Thanks, Gerald, I just now noticed your reply. I'm not very knowledgeable about SysRPL but I will study your program carefully. I would like to port the Ramanujan tau program to the Prime but it has only the Divisors function (like the 50g), not a fast way of calculating sums of divisor powers.

[Gerald H]

Here's my shortest & fastest effort at a stand alone programme, only good for input up to 1,048,575 but that's no great handicap as the programme is slow, taking 9 sec to process input 100 & 18 sec for input 200, the time characteristic being linear on the input.

Size: 350.0000

CkSum: # 59690d

Code:

::
  CK1&Dispatch
  # FF
  ::
    DUP
    ZINT 2
    Z<
    ?SEMI
    DUPDUP
    FPTR2 ^QMul
    '
    ::
      FPTR2 ^DupZIsOne?
      ?SEMI
      FPTR2 ^MZSQFF
      #2/
      ZINT 1
      SWAP
      ZERO_DO
      3PICK
      ROT
      COERCE
      #1+
      FPTR2 ^PPow#
      ZINT 1
      FPTR2 ^QSub
      ROT
      ZINT 1
      FPTR2 ^QSub
      FPTR2 ^ZQUOText
      FPTR2 ^QMul
      LOOP
    ;
    '
    NULLLAM
    BINT3
    NDUPN
    DOBIND
    ZINT 0
    3GETLAM
    FPTR2 ^Z2BIN
    ONE_DO
    INDEX@
    FPTR2 ^#>Z
    DUPDUP
    ZINT 35
    FPTR2 ^QMul
    ZINT 52
    3GETLAM
    FPTR2 ^QMul
    FPTR2 ^QSub
    OVER
    FPTR2 ^QMul
    ZINT 18
    2GETLAM
    FPTR2 ^QMul
    FPTR2 ^QAdd
    FPTR2 ^QMul
    FPTR2 ^QMul
    3GETLAM
    INDEX@
    FPTR2 ^#>Z
    DUP
    1GETLAM
    EVAL
    3UNROLL
    FPTR2 ^QSub
    1GETLAM
    EVAL
    FPTR2 ^QMul
    FPTR2 ^QMul
    FPTR2 ^QAdd
    LOOP
    ZINT 24
    FPTR2 ^QMul
    3GETLAM
    DUP
    1GETLAM
    EVAL
    SWAP
    BINT4
    FPTR2 ^RP#
    FPTR2 ^QMul
    SWAP
    FPTR2 ^QSub
    ABND
  ;
;

[John Keith]

Thanks, really neat!