|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 7680 (0x1e00)
Types: TextFile
Names: »tfftipow«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦09b4e9619⟧ »thcømat«
└─⟦this⟧
message fftipow
(fftipow=slang fpnames type.yes insertproc entry.no
fftipow)
; 8-10-74 ojh/rm
; powerspectrum of integer data
b. g1,e3 w.
k=10000
s. g6,j40,a30,b10
h.
g0=0
e2:
g1: g3,g2
j13: g0+13,0 ; last used
j30: g0+30,0 ; saved sr,w3
j8: g0+8,0 ; end address expr
j4: g0+4,0 ; take expr
j18: g0+18,0 ; zone alarm (index)
g2=k-2-g1, g3=g2 ; end of abs words,points
w.
a13: e0: g0
a14: 0
a11: 5 10 74
a12: 16 00 00
e1: rl. w2 (j13.)
ds. w3 (j30.)
dl w1 x2+8 ; take m
so w0 16
jl. w3 (j4.)
ds. w3 (j30.)
dl w1 x1
rl w3 x2+6
sz w3 1
cf w1 0 ; w1=m
sh w1 13 ; if m>13 or
sh w1 0 ; m<=0 then
jl. w3 (j18.) ; alarm index <m>
ac w3 x1-24
al w0 -1
ls w0 x3
rs. w0 a3. ; a3: m ones, 24-m zeroes
al w3 1
as w3 x1
al w1 x3+1 ; w1:=2**m+1
rl w3 x2+12 ; take array
rl w0 x3
ba w3 x2+10
wa w0 x3
wa. w0 a25.
rs. w0 a0. ; a0:=address R(0)
ba. w0 a26.
rs. w0 a1. ; a1:=address R(1)
rl w0 x3-2 ; index check
ws w0 x3
as w0 -1
sh w0 x1 ;if not 2**m+1 elements
jl. w3 (j18.) ; alarm index <2**m+1>
al w1 x1-2 ; w1:=2**m-1
rs. w1 a4. ; a4:=2**m-1
wa w1 3 ; w1:=2*(2**m-1)
wa. w1 a0.
rs. w1 a2. ; a2:=addr R(2**m-1)
al w1 1
b7: rs. w1 a5. ; a5:=l:=1
rl. w3 a3.
b0: ld w2 -1
al w0 x2
ld w0 1
sh w3 -1
jl. b0.
rl. w2 a5. ; w3:=bitreverse(l); w2:=l
sl w2 x3 ; if l>=bitreverse(l) then
jl. b1. ; no exchange
ld w3 1
wa. w2 a0.
wa. w3 a0.
rl w0 x3
rx w0 x2
rs w0 x3
b1: rl. w1 a5.
al w1 x1+1
sh. w1 (a4.)
jl. b7.
al. w2 a22.
al. w1 a23.
ds. w2 a20. ; init s2v addr, c2v addr
al w1 1
al w2 2
a26=k-1
ds. w2 a5. ; init p,q (1,1/2)
al w2 4
al w3 4 ; init pp (2)
b6: ds. w3 a3. ; store pp
dl. w3 a0. ; w2:=addr R(p),w3:=addr R(0)
b5: rl w1 x3 ; w1:=R(i)
al w0 x1 ; w0:=R(i)
wa w0 x2 ; w0:=R(i)+R(i+p)
ws w1 x2 ; w1:=R(i)-R(i+p)
rs w0 x3 ; R(i):=w0
rs w1 x2 ; R(i+p):=w1
aa. w3 a3. ; "i":="i"+pp, "i+p":="i+p"+pp
sh. w2 (a2.) ; if addr R(i+p)<addr R(2**m-1)
jl. b5. ; then next step
dl. w1 (a21.)
ds. w1 a18. ; a17 a18:=c2v
dl. w1 (a20.)
ds. w1 a16. ; a15 a16:=s2v
al w2 2 ; n:=1
b3: sl. w2 (a19.) ; if n=q then
jl. b4. ; then goto next l-step
rs. w2 a4. ; store n
ac w1 x2 ; w1:=-n
aa. w2 a0. ; a7:=addr R(p+n)
ds. w2 a7. ; a8:=addr R(p-n)
wa. w1 a5. ; a9:=addr R(p+n)
wa. w2 a5. ; a10:=addr R(p+p-n)
b2: ds. w2 a9. ; store addr R(i+p) R(j+p)
rl w1 x1
ci w1 ; float R(i+p)
ds. w1 a14.
rl w3 x2
ci w3 ; float R(j+p)
ds. w3 a12.
fm. w3 a18. ; w2w3:=R(i+p)*c
fm. w1 a16. ; w0w1:=R(j+p)*s
fs w3 3 ; w2w3:=re:=R(i+p)*c-R(j+p)*s
cf w3
dl. w1 a14.
rs. w3 a14. ; a14:=re
fm. w1 a18. ; w0w1:=R(j+p)*c
dl. w3 a12.
fm. w3 a16. ; w2w3:=R(i+p)*s
fa w3 3 ; w2w3:=im:=R(i+p)*s+R(j+p)*c
cf w3
al w0 x3
ws. w0 (a8.)
wa. w3 (a8.)
rs. w0 (a9.) ; R(i+p):=-R(j)+im
rs. w3 (a10.) ; R(j+p):=R(j)+im
rl. w1 a14.
ac w0 x1
wa. w0 (a7.)
wa. w1 (a7.)
rs. w0 (a8.) ; R(j):=R(i)-re
rs. w1 (a7.) ; R(i):=R(i)+re
dl. w2 a7.
aa. w2 a3.
ds. w2 a7. ; "i":="i"+pp; "j":="j"+pp
dl. w2 a9.
aa. w2 a3. ; "i+p":="i+p"+pp; "j+p":="j+p"+pp
sh. w1 (a2.) ; if addr R(j+p)<=addr R(2**m-1) then
jl. b2. ; goto next step
dl. w1 a16.
fm. w1 (a20.) ; w0w1:=s*s2v
dl. w3 a18.
fm. w3 (a21.) ; w2w3:=c*c2v
fs w3 3
rx. w3 a18.
rx. w2 a17. ; c:=c*c2v-s*s2v
fm. w3 (a20.) ; w2w3:=c*s2v
dl. w1 a16.
fm. w1 (a21.) ; w0w1:=s*c2v
fa w3 3
ds. w3 a16. ; s:=c*s2v+s*c2v
rl. w2 a4.
al w2 x2+2
jl. b3.
b4: dl. w3 a20.
aa. w3 a25.
ds. w3 a20. ; new addr s2v c2v
rl. w1 a5.
dl. w3 a3.
ds. w2 a5. ; new p,q
ad w3 1 ; new pp
wa. w1 a1.
rs. w1 a1. ; new addr R(p)
sh. w1 (a2.) ; if addr R(p)<=addr R(2**m-1)
jl. b6. ; then next step
al w0 0
rl. w2 a0.
wa. w2 a19.
al w2 x2-2 ; w2:=addr R(2**(m-1)-1)
al w3 x2+6 ; w3:=addr R(2**(m-1)+2)
b9: rx w0 x2
rx w0 x3
aa. w3 a25.
sl. w2 (a0.) ; if w2>addr R(0) then
jl. b9. ; next interchange
al w0 0
rs w0 x2 ; R(-1):=0
rl. w3 a0.
b8: dl w2 x3
ci w1
fm w1 3
ds. w1 a12.
ci w2
fm w2 5
fa. w2 a12.
ds w2 x3
al w3 x3+4
sh. w3 (a1.) ; if w3<=addr R(2**m) then
jl. b8. ; next convert
jl. (j8.) ; exit
a1: 0 ; addr R(p)
a0: 0 ; addr R(0)
a2: 0 ; addr R(2**m-1)
a6: 0 ; pp
a3: 0 ; pp
a4: 0 ; n
a19: 0 ; q
a5: 0 ; p
a8: 0 ; addr R(j)
a7: 0 ; addr R(i)
a10: 0 ; addr R(j+p)
a9: 0 ; addr R(i+p)
a15: 0
a16: 0 ; s
a17: 0
a18: 0 ; c
a21: 0 ; addr c2v
a20: 0 ; addr s2v
a24: -4
a25: 4
f.
0.999999705792
0.999998823456
0.999995293760
0.999981175328
0.999924701856
0.999698818720
0.998795456192
0.995184726720
a22: 0.980785280416
0.923879532576
0.707106781152
0.382683432368
a23: 0.195090322008
0.980171403296'-1
0.490676743248'-1
0.245412285224'-1
0.122715382856'-1
0.613588464896'-2
0.306795676288'-2
0.153398018624'-2
0.766990318720'-3
w.
0
0
<:fftipow <0>:>
e.
g0:
g1: 1
0,0,0,0
1<23+e1-e2
1<18+26<12+13<6+0 ; procedure fftipow(m,R);
0 ; value m; integer m; array R;
4<12+e0-e2
1<12+0
n.
▶EOF◀