FFT Multiplication (HP-48G/GX/G+)
|
02-10-2014, 09:25 PM
Post: #1
|
|||
|
|||
FFT Multiplication (HP-48G/GX/G+)
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: 51) That's a straightforward adaptation, except that I've used IFFT which the HP-71B lacks. Code:
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. |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)