(Comp.sources.hp48)
Item: 140 by kcooke@milton.u.washington.edu [Ken Cooke]
Subj: Machine Code FFT source code
Date: 14 Jul 1992

I'm glad to hear that someone liked the machine code FFT programs.
In response to email requests, here is the SASM source code.  It can
be assembled with the HP tools.  I didn't include FFT, because it is a
subset of IFFT (I commented the lines to remove).

I apologize for the shortage of higher-level comments (most of this is
still on paper).  Also, I traded some readability to avoid subroutine calls
(hardware stack levels are very precious) and keep variables in registers.
A few of the ROM calls are "unsupported entry points", which may need
fixing in the future.  With rumors of new 48's in the works, ones with
machine code Eqn Writers, re-written plotter user interfaces, and maybe
even faster Saturn processors...well, lets just say that it may be wise to
stick to entries listed in ENTRIES.A.

I hope this is of use to someone interested in ML math software.  While I
like playing games on my 48 as much as the next guy, it would be nice to 
see more math tools using the capabilities of sys-RPL and ML.  (I do, 
occasionally, use my 48 as a calculating device :)  Unfortunately,
there is no documentation on using the low-level floating point routines 
(MULTF, RADDF, etc) so a little hacking with Voyager or MLDL is required.

I would love to hear from other math-hackers out there...

Ken Cooke
N7VFE
kcooke@u.washington.edu

----------

* ANGLE.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
* Copyright 1992 Ken Cooke
* Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
* same size, that contains the angle of each complex element,
* according to the current trig mode.

         TITLE ANGLE
ASSEMBLE
         NIBASC    /HPHP48-A/
RPL
::
CK1NoBlame
CK&DISPATCH1
FOUR
::
* IF REAL ARRAY, EXIT UNCHANGED
         DUP TYPECARRY? NOT?SEMI                         ( ORIG )
* MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
         DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY         ( ORIG ANS )
CODE
****************** UNLISTED ENTRY POINTS *****************
ATTN?Lp  = #0CA88
ARG15    = #2B6D7
SetTrigMode        = #2AEF6
************************ MCODE ANGLE ***********************
* ENTRY:
*       STK2: ORIGINAL ARRY1 (COMPLEX)
*        STK1: ANSWER ARRY (REAL)

* EXIT: STK1: ANSWER ARRAY (REAL)

         NIBHEX  823        ;CLEAR SB AND XM IN HST
         D1=D1+    5        ;DROP STK1
         D=D+1     A        ;RETURN MEM
         GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK

         A=DAT1    A        ;A->ORIG
         D1=D1-    5        ;POINT BACK TO ANS
         C=DAT1    A        ;C->ANS
         R3=C.F    A        ;SAVE ->ANS IN R3
         D0=C
         D0=D0+    15
         D0=D0+    10       ;D0->ANS DATA (OR DIM2 IF ARRAY)

         D1=A               ;D1->ORIG
         D1=D1+    15       ;D1->#DIMS
         A=DAT1    A        ;A=#DIMS
         D1=D1+    5        ;D1->DIM1
         C=DAT1    A        ;C=DIM1
         D1=D1+    5        ;D1->ORIG DATA (OR DIM2 IF ARRAY)
         A=A-1     A        ;CHECK FOR VECTOR
         ?A=0      A
         GOYES     ISVECT
***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
         A=DAT1    A        ;A=DIM2
         D1=D1+    5        ;D1->ORIG DATA
         D0=D0+    5        ;D0->ANS DATA
         SETHEX
         GOSBVL    =MUL#    ;B=A*C
         C=B       A        ;C=DIM1*DIM2

ISVECT  C=C-1      A        ;SUB 1 SO CAN USE BORROW FOR LOOPING
         R4=C.F    A        ;SAVE #ELEMENTS-1 IN R4

ABSLOOP GOSBVL     =ATTN?Lp ;ABORT IF ATTN (ON) WAS PRESSED
         A=DAT1    W        ;A=RE
         D1=D1+    16       ;D1->IM
         CD0EX
         RSTK=C             ;PUSH D0
         C=DAT1    W        ;C=IM
         D1=D1+    16       ;D1->NEXT RE
         SETDEC
         GOSBVL    =SPLTAC  ;CONVERT BOTH TO LONGS
         GOSBVL    =SetTrigMode       ;SET ST ACCORDING TO TRIG MODE FLAGS
         GOSBVL    =ARG15   ;AB=ANGLE (TRASHES A,B,C,D,D0,R0,R1)
         GOSBVL    =PACK    ;A=ABSVAL
         C=RSTK
         D0=C               ;POP D0
         DAT0=A    W        ;STORE ABSVAL IN ANSWER ARRAY
         D0=D0+    16       ;POINT AT NEXT ANSWER
* CHECK FOR DONE
         C=R4.F    A        ;GET COUNTER
         SETHEX
         C=C-1     A
         R4=C.F    A        ;DEC COUNTER
         GONC    ABSLOOP    ;LOOP UNTIL BORROW

***NOW ANSWER CONTAINS ABS VALUES OF ORIG
         A=R3.F    A
         GOSBVL    =GPOverWrALp       ;GETPTR, REPLACE TOS WITH ANS, LOOP

ENDCODE
;
;
END_SRC


BEGIN_SRC ifft.s
******************************************************************
* MACHINE CODE IFFT by Ken Cooke
* Copyright 1992
******************************************************************

* RADIX-2 DIT IFFT
* Input: real or complex vector on STK1, with size=2^n
* Ouput: complex vector containing IFFT
* Note: uses the definition that FFT uses exp(-xxx) and IFFT uses exp(+xxx)

         TITLE IFFT
ASSEMBLE
         NIBASC    /HPHP48-A/
RPL
::
CK1NoBlame
CK&DISPATCH1
FOUR
::
         DUP MDIMS                              ( ARRY #N flag )
* CHECK FOR MATRIX                              (NcaseDIMERR = 37DF6 )
         NOT NcaseDIMERR            ( ARRY #N )
* CHECK FOR SIZE=2^N
         DUP                                    ( ARRY #N #N )

CODE
****************** CHECK THAT SIZE=POWER OF TWO **************
*CHKSIZE2 REPLACES BINT (ON TOS) WITH TRUE IF BINT IS A
* POWER OF 2, FALSE OTHERWISE.
*NOTE: WORKS UP TO 2^19, AND RETURNS TRUE FOR BINT=0
WrTLoop  EQU       #03B1A
WrFLoop  EQU       #03B06
         C=DAT1    A        ;C->BINT
         CD1EX              ;SAVE D1 IN C.A
         D1=D1+    5        ;D1->DATA
         A=DAT1    A        ;A.A=BINT
         D1=C               ;RESTORE D1

         SETHEX
         C=A       A
         A=-A-1    A        ;A=ONES COMP
         C=-C      A        ;C=TWOS COMP
         A=A!C     A        ;NO ZEROS IF BINT=2^N
         A=A+1     A        ;WILL SET CARRY IF BINT=2^N
         GOC       DoTrue
         GOVLNG    WrFLoop  ;OVERWRITE BINT WITH FALSE AND CONT
DoTrue   GOVLNG    WrTLoop  ;OVERWRITE BINT WITH TRUE AND CONT
ENDCODE

*                  ( ARRY #N FLAG )
NcaseDIMERR        ( IF FALSE THEN "INVALID DIMENSION" ERROR )

* IF SIZE=1, RETURN WITH ARRY ON STACK
         DUP#1= ITE DROP    ( ARRY #N )

::
* CREATE TEMP STORAGE (MUST BE 172 NIBS BIGGER THAN UNPACKED A FOR VARS)

         FORTYTWO OVER #* 172 #+ MAKEHEX$ SWAP  ( ARRY HEX$ #N )
         ONE{}N C%0 MAKEARRY                             ( ARRY HEX$ ANSARRY )

* DO THE FFT

CODE
****************** UNLISTED ENTRY POINTS *****************
PUTAB1   = #2C066
GETCD1   = #2C017
aSIN15   = #2B6E0
ATTN?Lp  = #0CA88
longpi   = #2A45D
varnibs  = 172
************************ MCODE FFT ***********************
* ENTRY:
*        STK3: (orig) real/complex data array
*       STK2: (temp) hex$ (to hold vars and unpacked reals during FFT)
*        STK1: (ans) complex answer array

* EXIT: STK1: complex FFT array (ans)

         NIBHEX  823        ;CLEAR SB AND XM IN HST
         D1=D1+    10       ;DROP temp AND ans
         D=D+1     A        ;RETURN MEM
         D=D+1     A
         GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK

* SAVE NEEDED POINTERS IN SCRATCH REGS
         D1=D1-    5        ;POINT BACK TO temp
         A=DAT1    A        ;A->temp
         A=A+CON   A,10     ;SKIP HEADER TO vars
         R4=A.F    A        ;R4->vars
         P=        0        ;*****REMEMBER:ALWAYS HAVE P=0 BEFORE LC(5)*****
         LC(5)     varnibs  ;
         SETHEX
         A=A+C     A        ;SKIP OVER VARS
         R1=A.F    A        ;R1->data

         D1=D1+    5        ;POINT TO orig
         A=DAT1    A        ;A->orig
         D0=A
         D0=D0+    10       ;D0->OBJ PROLOG

* SET FLAG IF REAL, FOR USE WHEN UNPACKING
* USE D.15 INSTEAD OF ST TO AVOID TROUBLE
         D=0       S        ;USE D.15 AS FLAG (0=COMPLEX,1=REAL)
         A=DAT0    A        ;A=OB PROLOG
         LC(5)     #02977   ;C=COMPLEX PROLOG
         ?A=C      A
         GOYES     CMPLEX
         D=D+1     S        ;IF REAL, D.15=1

CMPLEX   D0=D0+    10       ;D0->SIZE
         A=DAT0    A        ;A.A=SIZE
         R2=A.F    A        ;STORE SIZE IN R2
         D0=D0+    5        ;D0->DATA
         CD0EX
         R0=C.F    A        ;R0->origdata

*FIND LOG2 OF SIZE
         C=0       S        ;C.S WILL HOLD SIZE2=LOGbase2(SIZE)
LOG2LP   ASRB.F    A
         C=C+1     S
         ?ABIT=0   0
         GOYES     LOG2LP
         R2=C.F    S        ;STORE SIZE2 IN R2.15

*************************** BITREV AND UNPACK ****************************
*NOW, CONVERT REALS IN orig INTO EXT.REALS IN data
* AT THE SAME TIME, REARRANGE INTO BIT-REVERSED ORDER

         A=R0.F    A
         D0=A               ;D0->origdata
         D=0       A        ;D HOLDS INDEX INTO orig (START AT 0)

BITREV   CPEX      15       ;PUT BIT-COUNTER IN P
         P=P-1           ;SUB ONE SO CAN USE BORROW FOR BITLOOP
         C=D       A        ;C GETS INDEX INTO orig
         B=0       A        ;B WILL BE BIT-REV INDEX

*BIT REVERSAL
BITLOOP B=B+B      A        ;SHIFT B LEFT ONE BIT
         ?CBIT=0   0        ;TEST LSB
         GOYES     NOCARRY  ;SKIP IF ZERO
         B=B+1     A        ;ELSE FEED INTO B
NOCARRY  CSRB.F    A       ;SHIFT C RIGHT ONE BIT
         P=P-1              ;DEC BIT-COUNTER (BORROW WHEN DONE)
         GONC      BITLOOP  ;LOOP IF MORE

*MULT B (BITREV INDEX) BY 21 AND ADD TO ->data
         A=B       A        ;COULD USE =MUL# BUT THIS IS FASTER/SHORTER
         B=B+B     A
         B=B+B     A        ;B=4B
         A=A+B     A        ;A=5B
         B=B+B     A
         B=B+B     A        ;B=16B
         B=A+B     A        ;B=21B
         B=B+B     A        ;B=42B

         A=R1.F    A        ;A->data
         A=A+B     A        ;
         D1=A               ;D1 POINTS INTO data AT BITREV POSITION
         A=DAT0    W        ;GET RE(Ai)
         D0=D0+    16       ;POINT AT NEXT REAL
         GOSBVL    =SPLITA  ;NO NEED TO SETDEC (DONE IN SPLITA)
         GOSBVL  =PUTAB1    ;PUT LONG AT DAT1, D1+=21

* IF INPUT ARRAY WAS REAL, USE IMAG=0
         ?D=0      S        ;IF COMPLEX
         GOYES     GETIM    ; GET IMAG PART
         A=0       W        ;ELSE USE IMAG=0
         B=0       W
         GONC      NOGETIM  ;GO ALWAYS

GETIM    A=DAT0    W        ;GET IM(Ai)
         D0=D0+    16
         GOSBVL    =SPLITA

NOGETIM  GOSBVL  =PUTAB1    ;STORE IMAG PART
         SETHEX
         D=D+1     A        ;INC INDEX
         C=R2               ;LOAD SIZE AND SIZE2 INTO C
         ?D=C      A        ;COMPARE INDEX TO SIZE
         GOYES     BITDONE
         GOTO      BITREV   ;REPEAT IF NOT DONE
BITDONE

********** data NOW HOLDS LONG BITREV ARRAY **********


******************** BUTTERFLIES *******************
* ENTRY: temp CONTAINS BIT-REV LONG COMPLEX DATA
*        A IS NO LONGER NEEDED
* Define template for var storage: (var names from Num.Recipes)

istep    = 0
n        = (istep)+5
mmax     = (n)+5
data0    = (mmax)+5
wr       = (data0)+5
wi       = (wr)+21
tempr    = (wi)+21
wtemp    = (tempr)+21
theta    = (wtemp)+21
wpi      = (theta)+21
wpr      = (wpi)+21
m        = (wpr)+21

*FIRST, STORE SCRATCH REGS IN VARS

         SETHEX
         A=R4.F    A        ;A->istep
         D1=A
         D1=D1+    5        ;D1->n
         C=R2.F    A        ;C=size
         C=C+C     A        ;C=2*size
         DAT1=C    A        ;n=2*size

         D1=D1+    5        ;D1->mmax
         P=        0
         LC(5)     #2       ;C=2
         DAT1=C    A        ;mmax=2

         D1=D1+    5        ;D1->data0          ***STORE POINTER TO data[0] HERE
         C=R1.F    A        ;***THIS CORRECTS FOR "FORTRAN INDEXING"
         C=C-CON   A,16     ;*** SINCE ALGORITHM ASSUMES INDEXES START AT 1
         C=C-CON   A,5      ;*** I.E. data[1] = data0+21 = FIRST SPOT
         DAT1=C    A

**REDUCED TO ONE SINE/MAINLOOP BY USING wtemp TO HOLD LAST SIN(theta)
** AND CALCULATING SIN(theta/2)
         LC(5)     wtemp   ;INIT wtemp TO AVOID SIN(PI) SINCE NOT RETURN ZERO
         C=C+A     A        ;A STILL =R4
         D0=C               ;D0->wtemp
         A=0       W
         B=0       W        ;AB=0.0
         GOSBVL    =PUTAB0  ;STORE wtemp=0.0, D0->theta

**INIT theta TO -pi FOR FFT, +pi FOR IFFT
         D1=(5)  =longpi
         GOSBVL    =GETAB1  ;get pi from ROM
*        SETDEC                                 *********enable for FFT
*        A=-A-1    S        ;A=-PI              *********enable for FFT
         GOSBVL    =PUTAB0  ;STORE IN theta

************************ BEGIN MAIN LOOP ***********************
MAINLOOP
         C=R4.F    A        ;C->vars
         D1=C
         D1=D1+    5        ;D1->n
         C=DAT1    A        ;C=n
         D1=D1+    5        ;POINT AT mmax
         A=DAT1    A        ;A=mmax
* MAINLOOP TEST
         ?A<C      A        ;IF (mmax>=n)
         GOYES     MAINCONT ; THEN DONE
         GOTO      MAINDONE
** A=mmax, D1->mmax
MAINCONT
         SETHEX
         A=A+A     A        ;A=2*mmax
         D1=D1-    10       ;istep
         DAT1=A    A        ;STORE istep=2*mmax

         D1=D1+    10
         D1=D1+    10       ;D1->wr
         A=0       W
         B=0       W
         P=        14       ;CREATE AB=1.0
         B=B+1     P
         GOSBVL    =PUTAB1  ;STORE wr=1.0, D1->wi
         B=0       P
         GOSBVL    =PUTAB1  ;STORE wi=0.0, D1->tempr

         D1=D1+    16
         D1=D1+    5        ;D1->wtemp
         CD1EX
         RSTK=C             ;PUSH ->wtemp
         D0=C               ;D0->wtemp
         D1=C               ;D1->wtemp
         D1=D1+    16
         D1=D1+    16
         D1=D1+    10       ;D1->wpi

         SETDEC             :***NOTE: SETDEC BEFORE ALL FP MATH (WILL STAY IN DEC)
         GOSBVL    =GETAB0  ;AB=wtemp, D0->theta
         GOSBVL    =PUTAB1  ;STORE wpi=wtemp, D1->wpr
         GOSBVL    =GETAB0  ;AB=theta, D0->wpi
         GOSBVL    =DIV2    ;AB=theta/2
         D0=D0-    16
         D0=D0-    5        ;D0->theta
         GOSBVL    =PUTAB0  ;STORE theta=theta/2
         ST=1      9        ;SET UP FOR SIN (IN RADIANS)
         ST=0      4
         GOSBVL    =aSIN15  ;**NOTE:TRASHES A,B,C,D,P,D0,R0,R1!!!
         C=RSTK             ;POP ->wtemp
         D0=C               ;D0->wtemp
         GOSBVL    =PUTAB0  ;STORE wtemp=SIN(theta)

         C=B       W
         D=C       W        ;CD=AB
         C=A       W
         GOSBVL    =MULTF   ;SQUARE wtemp
         C=B       W
         D=C       W        ;CD=AB
         C=A       W
         GOSBVL    =RADDF   ;DOUBLE AB (D0 TRASHED)
         A=-A-1    S        ;NEGATE AB
         GOSBVL    =PUTAB1  ;STORE wpr=-2*SQR(wtemp), D1->m

         A=0       A
         A=A+1     B        ;A=1
         DAT1=A    A        ;STORE m=1

************************ BEGIN OUTER LOOP ***********************
OUTLOOP  SETHEX
         A=R4.F    A
         D0=A
         D0=D0+    10       ;D0->mmax
         P=        0
         LC(5)     m
         C=C+A     A
         D1=C               ;D1->m
         A=DAT0    A        ;A=mmax
         C=DAT1    A        ;C=m
*OUTLOOP TEST
         ?C<A      A        ;IF (m>=mmax)
         GOYES     NOTDONE  ; THEN DONE
         GOTO      OUTDONE
NOTDONE  RSTK=C             ;PUSH i=m (i KEPT ON TOS)
         A=R4.F    A
         D1=A               ;D1->n-5

************************ BEGIN INNER LOOP ***********************
INLOOP
         GOSBVL    =ATTN?Lp ;IF ATTN(ON) HAS BEEN PRESSED, GETPTR AND LOOP
         D1=D1+    5        ;D1->n
         A=DAT1    A        ;A=n
         C=RSTK             ;C=i
*INLOOP TEST
         ?C<=A     A        ;IF (i>n)
         GOYES     INCONT   ; THEN DONE
         GOTO      INDONE
INCONT   RSTK=C             ;PUSH i BACK
         D1=D1+    5        ;D1->mmax
         A=DAT1    A        ;A=mmax
         SETHEX
         A=A+C     A        ;A=j
         B=A       A
         B=B+B     A
         B=B+B     A
         A=A+B     A
         B=B+B     A
         B=B+B     A
         A=A+B     A        ;A=21*j

         B=C       A
         C=C+C     A
         C=C+C     A
         B=B+C     A
         C=C+C     A
         C=C+C     A
         B=C+B     A        ;B=21*i

         D1=D1+    5        ;D1->data0
         C=DAT1    A        ;C->data
         A=A+C     A        ;A->data[j]
         C=B+C     A        ;C->data[i]
         RSTK=C             ;PUSH ->data[i]
         C=A       A
         RSTK=C             ;PUSH ->data[j]

         D1=D1+    5        ;D1->wr
         D0=C               ;D0->data[j]
         GOSBVL    =GETAB1  ;AB=wr, D1->wi
         GOSBVL    =STAB2   ;SAVE wr IN R2/R3
         GOSBVL    =GETCD0  ;CD=data[j], D0->data[j+1]
         SETDEC
         GOSBVL    =MULTF   ;AB=wr*data[j]
         GOSBVL    =STAB0   ;SAVE IN R0/R1

         GOSBVL    =GETAB1  ;AB=wi, D1->tempr
         GOSBVL    =GETCD0  ;CD=data[j+1]
         GOSBVL    =MULTF   ;AB=wi*data[j+1]
         GOSBVL    =RCCD0   ;CD=wr*data[j]
         A=-A-1    S        ;NEGATE AB
         GOSBVL    =RADDF   ;AB=RESULT (D0 TRASHED)

         GOSBVL    =PUTAB1  ;tempr=RESULT, D1->tempi
         C=RSTK
         RSTK=C             ;COPY ->data[j] FROM TOS
         D0=C               ;D0->data[j]
         GOSBVL    =GETCD0  ;CD=data[j], D0->data[j+1]
         D1=D1-    16
         D1=D1-    16
         D1=D1-    10       ;D1->wi
         GOSBVL    =GETAB1  ;AB=wi

         GOSBVL    =MULTF
         GOSBVL    =STAB0   ;STORE wi*data[j] IN R0/R1
         GOSBVL    =GETCD0  ;CD=data[j+1]
         GOSBVL    =RCAB2   ;AB=wr
         GOSBVL    =MULTF   ;AB=wr*data[j+1]
         GOSBVL    =RCCD0   ;CD=wi*data[j]
         GOSBVL    =RADDF   ;AB=tempi (D0 TRASHED)

         GOSBVL    =STAB0   ;SAVE tempi IN R0/R1
         C=RSTK
***STORE D1=>data[j]
         D1=C               ;D1->data[j]
         C=RSTK
         RSTK=C             ;COPY ->data[i] FROM TOS
         C=C+CON   A,16
         C=C+CON   A,5      ;C->data[i+1]
         RSTK=C             ;PUSH ->data[i+1]
         D0=C
         GOSBVL    =GETCD0  ;CD=data[i+1]
         GOSBVL    =STCD2   ;SAVE data[i+1] IN R2/R3

         A=-A-1    S        ;NEGATE tempi
         GOSBVL    =RADDF   ;AB=data[i+1] - tempi
         CD1EX
         D1=C               ;C->data[j]
         D0=C
         D0=D0+    16
         D0=D0+    5        ;D0->data[j+1]

         GOSBVL    =PUTAB0  ;STORE data[j+1]
         GOSBVL    =RCAB0   ;AB=tempi
         GOSBVL    =RCCD2   ;CD=data[i+1]
         GOSBVL    =RADDF   ;AB=tempi+(data[i+1])
         C=RSTK             ;POP ->data[i+1]
         D0=C               ;D0->data[i+1]
         GOSBVL    =PUTAB0  ;STORE data[i+1]

         A=R4.F    A
         P=        0
         LC(5)     tempr
         SETHEX
         A=A+C     A
         D0=A               ;D0->tempr
         GOSBVL    =GETAB0  ;AB=tempr
         GOSBVL    =STAB0   ;STORE tempr IN R0/R1
         C=RSTK
         RSTK=C             ;COPY ->data[i] FROM TOS
         D0=C               ;D0->data[i]
         GOSBVL    =GETCD0  ;CD=data[i]
         GOSBVL    =STCD2   ;SAVE data[i] IN R2/R3
         SETDEC
         A=-A-1    S        ;NEGATE tempr
         GOSBVL    =RADDF   ;AB=data[i]-tempr
**USE D1=>data[j]
         GOSBVL    =PUTAB1  ;STORE data[j]

         GOSBVL    =RCAB0   ;AB=tempr
         GOSBVL    =RCCD2   ;CD=data[i]
         GOSBVL    =RADDF   ;AB=data[i]+tempr
         C=RSTK             ;POP ->data[i]
         D0=C               ;D0->data[i]
         GOSBVL    =PUTAB0  ;STORE data[i]

         A=R4.F    A
         D1=A               ;D1->istep
         A=DAT1    A        ;A=istep
         C=RSTK             ;POP i
         SETHEX
         C=C+A     A        ;i=i+istep
         RSTK=C             ;PUSH NEW i

         GOTO      INLOOP   ;NEXT ITERATION
************************ END INNER LOOP ***********************
INDONE
         A=R4.F    A
         P=        0
         LC(5)     wr
         SETHEX
         C=C+A     A
         RSTK=C             ;PUSH *wr
         D0=C               ;D0->wr
         LC(5)     wpr
         C=C+A     A
         D1=C               ;D1->wpr

         SETDEC
         GOSBVL    =GETAB0  ;AB=wr
         GOSBVL    =STAB0   ;SAVE OLDwr IN R0/R1
         GOSBVL    =GETCD1  ;CD=wpr, D1->wpr+21
         GOSBVL    =MULTF   ;AB=wr*wpr
         GOSBVL    =STAB2   ;SAVE IN R2/R3
         GOSBVL    =GETCD0  ;CD=wi
         D1=D1-    16
         D1=D1-    16
         D1=D1-    10       ;D1->wpi
         GOSBVL    =GETAB1  ;AB=wpi, D1->wpr
         GOSBVL    =MULTF   ;AB=wi*wpi
         A=-A-1    S        ;NEGATE AB
         GOSBVL    =RCCD2   ;CD=wr*wpr
         GOSBVL    =RADDF   ;AB=wr*wpr - wi*wpi (D0 TRASHED)
         GOSBVL    =RCCD0   ;CD=OLDwr
         GOSBVL    =RADDF   ;AB=RESULT (D0 TRASHED)

         C=RSTK             ;POP *wr
         D0=C               ;D0->wr
         GOSBVL    =PUTAB0  ;STORE NEW wr, D0->wi
         D1=D1-    16
         D1=D1-    5        ;D1->wpi
         GOSBVL    =GETAB1  ;AB=wpi, D1->wpr
         GOSBVL    =RCCD0   ;CD=OLDwr
         GOSBVL    =MULTF   ;AB=OLDwr*wpi
         GOSBVL    =STAB0   ;SAVE IN R0/R1

         GOSBVL    =GETCD0  ;CD=wi
         GOSBVL    =STCD2   ;SAVE wi IN R2/R3
         GOSBVL    =GETAB1  ;AB=wpr, D1->m
         GOSBVL    =MULTF   ;AB=wpr*wi
         GOSBVL    =RCCD0   ;CD=OLDwr*wpi
         GOSBVL    =RADDF   ;AB=wpr*wi + OLDwr*wpi (D0 TRASHED)
         GOSBVL    =RCCD2   ;CD=wi
         GOSBVL    =RADDF   ;AB=RESULT

         C=R4.F    A
         D0=C
         D0=D0+    16
         D0=D0+    16
         D0=D0+    9        ;D0->wi
         GOSBVL    =PUTAB0  ;STORE NEW wi
         C=DAT1    A        ;C=m
         SETHEX
         C=C+1     A
         C=C+1     A        ;m=m+2
         DAT1=C    A        ;STORE NEW m

         GOTO      OUTLOOP  ;NEXT ITERATION
************************ END OUTER LOOP ***********************

OUTDONE  SETHEX             ;D0 STILL POINTS TO mmax, A=mmax
         A=A+A     A        ;A=2*mmax
         DAT0=A    A        ;STORE NEW mmax

         GOTO      MAINLOOP
************************ END MAIN LOOP ***********************

MAINDONE
* temp NOW HOLDS vars (FIRST 172 NIBS) AND LONG COMPLEX ARRAY OF ANSWERS

***************** PACK AND COPY INTO ans *****************

         A=R4.F    A
         D1=A               ;D1->vars
         D1=D1+    5        ;D1->n
         A=DAT1    A        ;A=n
         R2=A.F    A        ;SAVE n IN R2 (COUNTER IN PACKLOOP)

***REMOVE FOR FFT
         ASRB.F    A        ;A=size
************* CONVERT A.A INTO EXT.REAL ************
* USING HORNER'S RULE, 10100 = 2^4 + 2^2 = ((((1)2+0)2+1)2+0)2+0
*  SINCE THE *2 AND ADDITION IS DONE IN DEC, RESULT IN DEC
* ENTER: BINT IN A.A (ANY)
* EXIT: EXT.REAL IN A/B (DEC)
* TRASHES: A,B,C,P

         P=        0
         LCHEX  19          ;REPEAT FOR EACH BIT IN A.A
         B=0    W  ;B WILL HOLD MANTISSA

BITLP    B=B+B  W  ;B*=2
         SETHEX
         A=A+A  A  ;SHIFT MSB OF BINT INTO CARRY
         SETDEC
         GONC   SKIP
         B=B+1  W
SKIP     C=C-1  B  ;DEC BIT COUNTER
         GONC   BITLP

         P=     5  ;A.A IS 0 FROM SHIFTING IN ZEROS
         A=A-1  X
NORMLP   A=A+1  X  ;A ACCUM'S EXPONENT BY COUNTING DIGITS
         BSRC      ;ROTATE UNTIL MANT IS LEFT-JUSTIFIED
         ?B#0   WP
         GOYES  NORMLP

         BSR    W  ;SHIFT MANT INTO PLACE (B.S=0)
         A=0    S  ;SIGN IS POSITIVE
* AB NOW HOLDS EXT.REAL
**************** END OF CONVERSION ****************
         GOSBVL    =STAB0   ;SAVE SIZE IN R0/R1
*******REMOVE FOR FFT

         SETHEX
         A=R4.F    A        ;**THIS IS THE LAST TIME WE USE R4(temp)
        P=         0
         LC(5)     varnibs  ;SKIP OVER VARS TO DATA
         C=A+C     A        ;C->data
         D0=C               ;D0->data (SOURCE OF COPY)

         D1=(5)    =DSKTOP  ;RECALL SAVED D1
         C=DAT1    A
         D1=C               ;D1->STK3 (orig)
         D1=D1-    10       ;D1->STK1 (ans)
         C=DAT1    A
         R3=C.F    A        ;SAVE ->ans IN R3
         D1=C               ;D1->ans
         D1=D1+    15      ;SKIP HEADER
         D1=D1+    10       ;D1->ansdata (DEST OF COPY)

* NOW PACK AND COPY
* NOTE: OVER/UNDERFLOW ERRORS ARE TRAPPED IN =PACK. orig IS SAVED AS TOS
PAKLOOP
         GOSBVL    =GETAB0  ;GET NEXT LONG
         GOSBVL    =RCCD0   ;RECALL SIZE                 *REMOVE FOR FFT
         SETDEC
         GOSBVL    =DIVF    ;DIVIDE BY SIZE              *REMOVE FOR FFT
         GOSBVL    =PACK    ;EXT.REAL->REAL

         DAT1=A    W        ;PUT REAL IN PLACE
         D1=D1+    16       ;POINT AT NEXT REAL

         C=R2.F    A        ;GET COUNTER
         SETHEX
         C=C-1     A        ;DEC COUNTER
         R2=C.F    A
         ?C#0      A
         GOYES     PAKLOOP


* NOW RESULT IS PACKED INTO ans ARRAY
* AT LAST, REPLACE STK1 (ORIGINAL ARRAY) WITH RESULT ARRAY
         A=R3.F    A                  ;A->ans
         GOVLNG    =GPOverWrALp       ;RESTORE REGS, PUT @A ON TOS, CONT RPL
ENDCODE
;
;
;
END_SRC



BEGIN_SRC absv.s
* ABSV.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
* Copyright 1992 Ken Cooke
* Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
* same size, that contains the element-by-element absolute values.

         TITLE ABSV
ASSEMBLE
         NIBASC    /HPHP48-A/
RPL
::
CK1NoBlame
CK&DISPATCH1
FOUR
::
* IF REAL ARRAY, EXIT UNCHANGED
         DUP TYPECARRY? NOT?SEMI                         ( ORIG )
* MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
         DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY         ( ORIG ANS )
CODE
****************** UNLISTED ENTRY POINTS *****************
ATTN?Lp  = #0CA88
SQR15    = #2B9F3

************************ MCODE ABSV ***********************
* ENTRY:
*       STK2: ORIGINAL ARRY1 (COMPLEX)
*        STK1: ANSWER ARRY (REAL)

* EXIT: STK1: ANSWER ARRAY (REAL)

         NIBHEX  823        ;CLEAR SB AND XM IN HST
         D1=D1+    5        ;DROP STK1
         D=D+1     A        ;RETURN MEM
         GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK

* NOW IF USER HITS "ON" KEY, WILL LEAVE ORIGINAL STACK
         A=DAT1    A        ;A->ORIG
         D1=D1-    5        ;POINT BACK TO ANS
         C=DAT1    A        ;C->ANS
         R0=C.F    A        ;SAVE ->ANS IN R0
         D0=C
         D0=D0+    15
         D0=D0+    10       ;D0->ANS DATA (OR DIM2 IF ARRAY)

         D1=A               ;D1->ORIG
         D1=D1+    15       ;D1->#DIMS
         A=DAT1    A        ;A=#DIMS
         D1=D1+    5        ;D1->DIM1
         C=DAT1    A        ;C=DIM1
         D1=D1+    5        ;D1->ORIG DATA (OR DIM2 IF ARRAY)
         A=A-1     A        ;CHECK FOR VECTOR
         ?A=0      A
         GOYES     ISVECT
***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
         A=DAT1    A        ;A=DIM2
         D1=D1+    5        ;D1->ORIG DATA
         D0=D0+    5        ;D0->ANS DATA
         SETHEX
         GOSBVL    =MUL#    ;B=A*C
         C=B       A        ;C=DIM1*DIM2

ISVECT  C=C-1      A        ;SUB 1 SO CAN USE BORROW FOR LOOPING
         R4=C.F    A        ;SAVE #ELEMENTS-1 IN R4

ABSLOOP GOSBVL     =ATTN?Lp ;ABORT IF ATTN (ON) WAS PRESSED
         SETDEC
         GOSUB     SQDAT1   ;SQUARE RE PART
         GOSBVL    =STAB2   ;SAVE IN R2/R3
         GOSUB     SQDAT1   ;SQUARE IM PART

         CD0EX
         RSTK=C             ;PUSH D0
         GOSBVL    =RCCD2   ;CD=SQR(RE)
         GOSBVL    =RADDF   ;AB=SQR(RE)+SQR(IM) (D0 TRASHED)
         GOSBVL    =SQR15   ;AB=LONG ABSVAL
         GOSBVL    =PACK    ;A=ABSVAL
         C=RSTK
         D0=C               ;POP D0
         DAT0=A    W        ;STORE ABSVAL IN ANSWER ARRAY
         D0=D0+    16       ;POINT AT NEXT ANSWER
* CHECK FOR DONE
         C=R4.F    A        ;GET COUNTER
         SETHEX
         C=C-1     A
         R4=C.F    A        ;DEC COUNTER
         GONC    ABSLOOP    ;LOOP UNTIL BORROW

***NOW ANSWER CONTAINS ABS VALUES OF ORIG
         GOVLNG    =GPOverWrR0Lp      ;GETPTR, REPLACE TOS WITH ANS, LOOP

************************ SQDAT1 ***********************************
* GET REAL FROM DAT1, UNPACK AND SQUARE IT, INC D1

SQDAT1   A=DAT1    W        ;A=REAL
         D1=D1+    16       ;D1->NEXT
         GOSBVL    =SPLITA  ;AB=REAL
         C=B       W
         D=C       W
         C=A       W        ;COPY INTO CD
         GOSBVL    =MULTF   ;AB=SQR(REAL)
         RTN

ENDCODE
;
;
