FFT Multiplication (HP-48G/GX/G+)

[Gerson W. Barbosa]

This is an RPL adaptation of Valentin Albillo's FFTMULT program for the HP-71B:

  DatafileVA013. (HP-71B Fantastic FOUR).pdf (Size: 109.47 KB / Downloads: 49)

That's a straightforward adaptation, except that I've used IFFT which the HP-71B lacks.

Code:


%%HP: T(3)A(D)F(.);
\<< DUP2 SIZE SWAP
SIZE + ROT ROT DUP2
SIZE SWAP SIZE DUP2
MAX LN 2 LN / DUP
FP NOT NOT SWAP IP
+ 2 SWAP ^ DUP +
ROT OVER - NEG ROT
ROT - NEG ROT ROT 1
SWAP
  START 0 +
  NEXT SWAP ROT 1
ROT
  START 0 +
  NEXT OBJ\-> 1 \->LIST
\->ARRY FFT OBJ\-> 1
GET \->LIST SWAP OBJ\->
1 \->LIST \->ARRY FFT
OBJ\-> 1 GET \->LIST
  \<< *
  \>> DOLIST OBJ\-> 1
\->LIST \->ARRY IFFT
OBJ\-> 1 GET \->LIST 1
ROT SUB RE IP
REVLIST DUP SIZE 1
- 1 SWAP
  FOR i i GETI ROT
ROT GETI SWAP DROP
ROT OVER 10000 MOD
+ OVER 10000 / IP
OVER 10000 / IP +
SWAP 10000 MOD ROT
DROP SWAP ROT ROT i
SWAP PUTI ROT PUTI
DROP
  NEXT REVLIST
\>>

Here is an example from Valentin's article:

{ 5476 1407 } { 8613 2724 } FFTMULT -->

{ 4716 7491 5498 2668 }

This matches Valentin's result. However, I've found a discrepancy in this example:

{ 1234 5678 9012 } { 3456 7890 } FFTMULT -->

{ 426 7639 7023 2002 4680 }

This should have been { 426 7640 7023 2002 4680 }, as we would get using Valentin's algorithm:

10 OPTION BASE 1 @ DIM X(3),Y(2),Z(5)
15 READ X,Y
20 CALL FFTMULT(X,Y,Z) @ MAT DISP Z
100 END
200 DATA 1234,5678,9012,3456,7890
900 SUB FFTMULT(A(),B(),C()) @ I=UBND(A,1) @ J=UBND(B,1)
910 K=I+J @ COMPLEX X(I),Y(J) @ MAT X=A @ MAT Y=B
920 N=LOG2(MAX(I,J)) @ E=10000 @ N=2^(INT(N)+(FP(N)#0))
930 M=2*N @ COMPLEX X(M),Y(M),Z(M) @ MAT X=FOUR(X)
940 MAT Y=FOUR(Y) @ FOR I=1 TO M @ Z(I)=X(I)*Y(I) @ NEXT I
950 COMPLEX Z(M,1) @ MAT Z=TRN(Z) @ MAT Z=FOUR(Z)
960 MAT Z=TRN(Z) @ MAT Z=(M*M)*Z @ COMPLEX Z(M)
965 MAT DISP Z @ DISP
970 FOR I=1 TO K @ C(I)=IROUND(REPT(Z(I))) @ NEXT I
980 FOR I=K-1 TO 1 STEP -1 @ J=I+1 @ C(J)=C(J)+MOD(C(I),E)
990 C(I)=C(I) DIV E+C(J) DIV E @ C(J)=MOD(C(J),E) @ NEXT I

>RUN
(4264704.00002,-0)
(29359428.0003,.000000128)
(75944892,-0)
(71104680,.000000512)
(-.000016,-0)
(-.00025024,-.000000128)
(0,-0)
(-.0002352,-.000000512)

426
7640
7023
2002
4680

As a comparison, the first element of the inverse Fourier transform is computed as (4264703.99998, 0) by the HP-48 calculator, which suggests is should be actually (4264704, 0). Perhaps the addition of a suitable quantity to the real part of these values should fix the problem and allow for chunks of 5 digits, but I haven't tried it yet.
Another issue is speed. My HP-48GX is about 50% slower than the HP-71B when multiplying two 512-digit numbers (173 seconds versus 117 seconds). But of course the RPL program can be improved, as you can see by the code.
Anyway, when these issues are solved (all of you who are interested are invited to try) this can make for a nice tool to compute pi to a couple of thousand digits, for instance.