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

[Gerson W. Barbosa]

Ever been curious about all digits of 253! ? Out of 500, the HP-48GX gives the first 12 only: 517346099264. Just run the program below with single argument 253 and wait for about one hour and 45 minutes. This is by no means an efficient way to do this, but might be a way to check the accuracy of the algorithm after 253 multiplications, even though one on the operands is small. (Better use Emu48 and get the result under 60 seconds or less).

Code:


%%HP: T(3)A(D)F(,);
\<< 1 -
[ 1 ]
1 ROT ROT 1 ROT
  START OVER 1 +
DUP 1 \->ARRY ROT
FFTMLTS2 TRIM ROT
DROP
  NEXT SWAP DROP
\>>

where FFTMLTS2 is Raymond del Tondo's SysRPL program above and TRIM is the following program, slightly adapted from TRIM, by Wayne H. Scott, part of his POLY package:

Code:


%%HP: T(3)A(D)F(,);                 ;  
\<< OBJ\-> OBJ\-> DROP \->          ; This trims off the leading zeros of the array.
n                                   ; For example,
  \<< n                             ;
    WHILE ROLL DUP                  ; [ 0 0 0 1234 0 5789 ]  -->  [ 1234 0 5789 ]
ABS NOT n 1 - AND                   ;
    REPEAT DROP 'n'
DECR
    END n ROLLD n 1
\->LIST \->ARRY
  \>>
\>>

-->

[ 5173 4609 9264 78 9218 433 899 7295 2746 9542 3561 2720 6639 9607 4846 3616 3134 3029 313 41 2383 1443 7828 1112 1374 4932 5428 7661 7296 3169 484 977 8527 4435 4743 3640 9654 4589 6311 9980 576 3521 219 7345 934 790 1685 4446 6163 7384 4451 7144 4589 2498 2615 9309 2898 1062 2514 4818 9873 9824 3499 6567 2944 9381 9909 5203 1087 3152 8570 9655 6175 4517 6766 2603 4976 5427 6777 1987 6267 959 7099 9373 2257 7683 9082 7849 7279 3284 6880 6763 5727 3110 3332 7966 9572 6049 2114 9638 6749 6804 5622 1513 5307 5201 4396 1440 1249 2800 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

As far as I have checked, all digits are correct, considering the hidden leading zeros in some 4-digit elements.

-----

Update: I've run the User-RPL program on the emulator (50g) and formatted the output to integer. The result was exactly the same I got when computing 253! on the hp 50g in exact mode.

Also, an exact match was obtained for 997! (2559 digits). This was computed using arrays with 5-digit elements. Here they are, in case someone wants to check them out :-)

Code:


40359724461645390234292652765290740211090335246139389143030797373519663190106842​37263758833853587102137006631984661977193944113181265519616824478081989249680516​43330792790545975658652366984953410102994729193397927862707860663312211428139657​33919921028883982924596589308458677218884794980135443761645075224506666559889800​94175577967376951675213432494794136314145342021847264214793926157307811731645263​93982880263279118925406206180689438683308644696334133955184235540077242451165903​81101827719832180031595827989994138156615149091737998105454985248322329275243898​11980802708882543991975745364605704734303715958724031694867571661542949412580453​11241382930836862005052391967478429035362983199050663230586866257612402804942403​44233166394434168335073220412356534986944621623211159899567872446218256850113174​63838577067904001075072667390026316129311241122279096729357421049685332780747960​00335855930432060517447195226436187301231195091058916141500005034486568847599649​00494067769318523221837865944485464570390882493401514455003570460531797737862031​18550953567694888922171302000112504911516415310901200837651592219697553144378802​09281708574493693840125338722070514029362985801732618715060934298236579096167095​85950405331060872571119845720022654435044594115773486340742853243112648568667878​84661486819750191740104532976390040068195207044638407735286912244552652299854897​64356909675383800245969276679872407757924211918488179598530382266647547907226165​47980297654789693965688881325682653906791569527887851625739692098351138902956310​11123253723954647397831433613628798725785501475711681360833913542427351428039887​35616917749898060073075542403509536490539404444972668319521415425667918323473675​96656633239099325959195904942407038086186468220698646372928155733874746654662785​92062875719964916067979790641428194695892008126790265612881240871363598309598670​34513441434850212864818601504529520195828528045600869420646442863720485414968365​31269052383502650854577265971210516113769359526291937135884001947338380202834453​11816794177165630135012424772911390424228141663696011522232935969575275309346520​46662174154235850073391729650007182794396630407081318880947107940245036774649857​42937922077663735689021159654000934909225598804790941759477837570572384191816766​30262770090339396547856717150451221853157302493936160447379021701169807360000000​00000000000000000000000000000000000000000000000000000000000000000000000000000000​00000000000000000000000000000000000000000000000000000000000000000000000000000000​0000000000000000000000000000000000000000000000000000000000000000000000000000000

Update:

It appears 5-digit group arrays yield reliable results up to 128 elements only (1280-digit results, 5*256). 4^2^11 has 1234 digits and is correctly computed using the program below. The next doubling (4^2^12) is not.

Code:



11 

%%HP: T(3)A(D)F(.);
\<< [ 4 ] 1 ROT
  START DUP FFTMULT TRIM
  NEXT A2I
\>>

EVAL

'4^2^11'              

EVAL  -

-->   0 




FFTMULT:

%%HP: T(3)A(D)F(.);
\<< DUP2 SIZE OBJ\-> DROP SWAP SIZE OBJ\-> DROP DUP2 MAX UNROT + 4.
ROLLD LN 2. LN / DUP FP NOT NOT SWAP IP + 2. SWAP ^ DUP + DUP ROT SWAP
1. \->LIST RDM UNROT 1. \->LIST RDM FFT OBJ\-> 1. GET \->LIST SWAP FFT
OBJ\-> 1. GET \->LIST
  \<< *
  \>> DOLIST OBJ\-> 1. \->LIST \->ARRY IFFT 1. ROT SUB RE 0. RND DUP
SIZE OBJ\-> DROP 1. - 1.
  FOR i i GETI 100000. MOD PICK3 ROT DUP UNROT GET ROT + PUT i GETI
100000. / IP PICK3 ROT GET 100000. / IP + i SWAP PUTI OVER OVER GET
100000. MOD PUT -1.
  STEP
\>>

TRIM:

%%HP: T(3)A(D)F(.);
\<< OBJ\-> OBJ\-> DROP \-> n
  \<< n
    WHILE ROLL DUP ABS NOT n 1. - AND
    REPEAT DROP 'n' DECR
    END n ROLLD n 1. \->LIST \->ARRY
  \>>
\>>

A2I:

%%HP: T(3)A(D)F(.);
\<< DUP "" UNROT SIZE OBJ\-> SWAP
  FOR i DUP i GET R\->I \->STR FMT ROT SWAP + SWAP
  NEXT DROP OBJ\->
\>>

FMT:

%%HP: T(3)A(D)F(.);
\<< DUP SIZE 5 - NOT NOT
  IF
  THEN 5. OVER SIZE - 1. SWAP
    START "0" SWAP +
    NEXT
  END
\>>