! File: fft31t.i90
! Public domain 2004 James Van Buskirk
subroutine fft31t(h)
! 658 adds, 234 muls
   complex(wp) h(0:30)
   real(wp) xr0, xi0, xr1, xi1, xr2, xi2, xr3, xi3, xr4, xi4, xr5, xi5
   real(wp) xr6, xi6, xr7, xi7, xr8, xi8, xr9, xi9, xra, xia, xrb, xib
   real(wp) xrc, xic, xrd, xid, xre, xie, xrf, xif
   real(wp) yr1, yi1, yr2, yi2, yr3, yi3, yr4, yi4, yr5, yi5, yr6, yi6
   real(wp) yr7, yi7, yr8, yi8, yr9, yi9, yra, yia, yrb, yib, yrc, yic
   real(wp) yrd, yid, yre, yie, yrf, yif
   real(wp) xr16b, xi16b, xr49e, xi49e, xr27c, xi27c
   real(wp) xr5af, xi5af, xr38d, xi38d, xr7_2, xi7_2
   real(wp) xr8_3, xi8_3, xra_5, xia_5, xrb_6, xib_6
   real(wp) xrd_8, xid_8, xre_9, xie_9, xr1_b, xi1_b
   real(wp) xr2_c, xi2_c, xr4_e, xi4_e, xr5_f, xi5_f
   real(wp) xr16b49e, xi16b49e, xr27c5af, xi27c5af
   real(wp) xrsigma, xisigma, ursigma, uisigma
   real(wp) sr0, si0, ur1234, ui1234, ur12, ui12, ur34, ui34
   real(wp) ur1, ui1, ur2, ui2
   real(wp) ur3, ui3, ur4, ui4, ur5, ui5
   real(wp) xra4_5e, xia4_5e, xrae_54, xiae_54
   real(wp) xrd1_8b, xid1_8b, xrdb_81, xidb_81
   real(wp) xra4d1_5e8b, xia4d1_5e8b, pr6, pi6, pr8a, pi8a
   real(wp) pr8, pi8, pra, pia, prc, pic, pre, pie
   real(wp) xrb5_6f, xib5_6f, xrbf_65, xibf_65
   real(wp) xre2_9c, xie2_9c, xrec_92, xiec_92
   real(wp) xrb5e2_6f9c, xib5e2_6f9c, pr7, pi7, pr9b, pi9b
   real(wp) pr9, pi9, prb, pib, prd, pid, prf, pif
   real(wp) ur6, ui6, ur7, ui7
   real(wp) pr89, pi89, prcd, picd
   real(wp) ur89, ui89, ur8, ui8, ur9, ui9
   real(wp) urcd, uicd, urc, uic, urd, uid
   real(wp) prab, piab, pref, pief
   real(wp) urab, uiab, ura, uia, urb, uib
   real(wp) uref, uief, ure, uie, urf, uif
   real(wp) ur8a, ui8a, wr6, wi6, wr8eac, wi8eac
   real(wp) wr8e, wi8e, wrac, wiac
   real(wp) wr8_e, wi8_e, wra_c, wia_c
   real(wp) wr8, wi8, wre, wie, wra, wia, wrc, wic
   real(wp) ur9b, ui9b, wr7, wi7, wr9fbd, wi9fbd
   real(wp) wr9f, wi9f, wrbd, wibd
   real(wp) wr9_f, wi9_f, wrb_d, wib_d
   real(wp) wr9, wi9, wrf, wif, wrb, wib, wrd, wid
   real(wp) sr1, si1, sr2, si2, sr3, si3, sr4, si4, sr5, si5
   real(wp) sr6, si6, sr7, si7, sr8, si8, sr9, si9, sra, sia
   real(wp) srb, sib, src, sic, srd, sid, sre, sie, srf, sif
   real(wp) yr16b, yi16b, yr49e, yi49e, yr27c, yi27c
   real(wp) yr5af, yi5af, yr38d, yi38d, yr7_2, yi7_2
   real(wp) yr8_3, yi8_3, yra_5, yia_5, yrb_6, yib_6
   real(wp) yrd_8, yid_8, yre_9, yie_9, yr1_b, yi1_b
   real(wp) yr2_c, yi2_c, yr4_e, yi4_e, yr5_f, yi5_f
   real(wp) yr16b49e, yi16b49e, yr27c5af, yi27c5af
   real(wp) yrsigma, yisigma, vrsigma, visigma
   real(wp) vr1234, vi1234, vr12, vi12, vr34, vi34
   real(wp) vr1, vi1, vr2, vi2
   real(wp) vr3, vi3, vr4, vi4, vr5, vi5
   real(wp) yra4_5e, yia4_5e, yrae_54, yiae_54
   real(wp) yrd1_8b, yid1_8b, yrdb_81, yidb_81
   real(wp) yra4d1_5e8b, yia4d1_5e8b, qr6, qi6, qr8a, qi8a
   real(wp) qr8, qi8, qra, qia, qrc, qic, qre, qie
   real(wp) yrb5_6f, yib5_6f, yrbf_65, yibf_65
   real(wp) yre2_9c, yie2_9c, yrec_92, yiec_92
   real(wp) yrb5e2_6f9c, yib5e2_6f9c, qr7, qi7, qr9b, qi9b
   real(wp) qr9, qi9, qrb, qib, qrd, qid, qrf, qif
   real(wp) vr6, vi6, vr7, vi7
   real(wp) qr89, qi89, qrcd, qicd
   real(wp) vr89, vi89, vr8, vi8, vr9, vi9
   real(wp) vrcd, vicd, vrc, vic, vrd, vid
   real(wp) qrab, qiab, qref, qief
   real(wp) vrab, viab, vra, via, vrb, vib
   real(wp) vref, vief, vre, vie, vrf, vif
   real(wp) vr8a, vi8a, zr6, zi6, zr8eac, zi8eac
   real(wp) zr8e, zi8e, zrac, ziac
   real(wp) zr8_e, zi8_e, zra_c, zia_c
   real(wp) zr8, zi8, zre, zie, zra, zia, zrc, zic
   real(wp) vr9b, vi9b, zr7, zi7, zr9fbd, zi9fbd
   real(wp) zr9f, zi9f, zrbd, zibd
   real(wp) zr9_f, zi9_f, zrb_d, zib_d
   real(wp) zr9, zi9, zrf, zif, zrb, zib, zrd, zid
   real(wp) tr1, ti1, tr2, ti2, tr3, ti3, tr4, ti4, tr5, ti5
   real(wp) tr6, ti6, tr7, ti7, tr8, ti8, tr9, ti9, tra, tia
   real(wp) trb, tib, trc, tic, trd, tid, tre, tie, trf, tif
   real(wp), parameter :: dc1 = -.41848031916582736768981687098748646_wp
   real(wp), parameter :: dc2 = 1.03708005397391977997160126838915860_wp
   real(wp), parameter :: dc3 = .06899615305609372142805240007781954_wp
   real(wp), parameter :: dc4 = -.71199788901815727304919706300135852_wp
   real(wp), parameter :: dc5 = .49457655717669416346035981441598679_wp
   real(wp), parameter :: dc6 = -.96808390091782605854354886831133920_wp
   real(wp), parameter :: dc7 = .32508216495576250692240420538780022_wp
   real(wp), parameter :: dc8 = .16949439222093165653795560902818660_wp
   real(wp), parameter :: dc9 = -.29637372110299413755460095857226913_wp
   real(wp), parameter :: dc10 = .04534684817389996223192362526889384_wp
   real(wp), parameter :: dc11 = -.34172056927689409978652458384116289_wp
   real(wp), parameter :: dc12 = .07310633583025898405268087482211124_wp
   real(wp), parameter :: dc13 = -.20167486904767495859952810665520680_wp
   real(wp), parameter :: dc14 = -.84201778467603236410459766557180865_wp
   real(wp), parameter :: dc15 = -.13041239332209095847634591191071144_wp
   real(wp), parameter :: dc16 = -.45935884940861915763712570904324634_wp
   real(wp), parameter :: dc17 = .48562853299381401878246608937483484_wp
   real(wp), parameter :: dc18 = .05290024156536825355103499921807512_wp
   real(wp), parameter :: dc19 = -.14593330617823442758015434230035912_wp
   real(wp), parameter :: dc20 = -.60928979281795509222872426261479928_wp
   real(wp), parameter :: dc21 = .50714554918236722335933886804513386_wp
   real(wp), parameter :: dc22 = -.09614247372823604317748361539732969_wp
   real(wp), parameter :: dc23 = -.57652834602834558094251739609128330_wp
   real(wp), parameter :: dc24 = .06326947458204131524249727162375559_wp
   real(wp), parameter :: dc25 = -.76932220806904896036979558584064568_wp
   real(wp), parameter :: dc26 = -.41013896184380627585461186107586186_wp
   real(wp), parameter :: dc27 = .36697396683700721180660373150156373_wp
   real(wp), parameter :: dc28 = -.06956934754225036504904405641510024_wp
   real(wp), parameter :: dc29 = -.41717983028166295195804014713757581_wp
   real(wp), parameter :: csigma = -.16666666666666666666666666666666700_wp
   real(wp), parameter :: ds1 = .75538558300353835543581057551249341_wp
   real(wp), parameter :: ds2 = -.71798870119142816504365819783899618_wp
   real(wp), parameter :: ds3 = -1.14129841620797351615515640500369610_wp
   real(wp), parameter :: ds4 = .17619908884183581953599075099065958_wp
   real(wp), parameter :: ds5 = -.44789045813036076345265894471834662_wp
   real(wp), parameter :: ds6 = -.42330971501654535111149820716469997_wp
   real(wp), parameter :: ds7 = -.54178961234959234550766744684833662_wp
   real(wp), parameter :: ds8 = .09389915421923158205500850212998998_wp
   real(wp), parameter :: ds9 = .10209749786491606368824206751661151_wp
   real(wp), parameter :: ds10 = .36010442196019251577804178188101286_wp
   real(wp), parameter :: ds11 = -.25800692409527645208979971436440136_wp
   real(wp), parameter :: ds12 = -.68362751420088060523602271315212364_wp
   real(wp), parameter :: ds13 = 1.42448046407360845685461718962915170_wp
   real(wp), parameter :: ds14 = .71705834223306922870900492103008074_wp
   real(wp), parameter :: ds15 = -.43581585016667742322916851390840979_wp
   real(wp), parameter :: ds16 = .44515550282226738565114480689187057_wp
   real(wp), parameter :: ds17 = .52796329829794124367930049003364155_wp
   real(wp), parameter :: ds18 = -.49467751640467748805281786601391635_wp
   real(wp), parameter :: ds19 = 1.03076374706570777841622291802646680_wp
   real(wp), parameter :: ds20 = .51886829082317972874576505841504815_wp
   real(wp), parameter :: ds21 = .82187006169806796277221650036735597_wp
   real(wp), parameter :: ds22 = -.79391727157204814808670251421336608_wp
   real(wp), parameter :: ds23 = -.75874209018537673384980278323642215_wp
   real(wp), parameter :: ds24 = -.19833215631171570317346941647854694_wp
   real(wp), parameter :: ds25 = -.43149551708674121082911820080227901_wp
   real(wp), parameter :: ds26 = .27843331019209505198390355841028348_wp
   real(wp), parameter :: ds27 = .59471076351191660168095301828858968_wp
   real(wp), parameter :: ds28 = -.57448393456065017248053353113326010_wp
   real(wp), parameter :: ds29 = -.54903093419716620563980211536607990_wp
   real(wp), parameter :: ssigma = .92796072713833698701991188315309154_wp
   real(wp), parameter :: c50 = -.89442719099991587856366946749251043_wp
   real(wp), parameter :: c51 = -.27639320225002103035908263312687237_wp
   real(wp), parameter :: s52 = .61803398874989484820458683436563806_wp

   xr0 = real(h(0))
   xi0 = aimag(h(0))
   xr1 = real(h(1))+real(h(30))
   xi1 = aimag(h(1))+aimag(h(30))
   yr1 = real(h(1))-real(h(30))
   yi1 = aimag(h(1))-aimag(h(30))
   xr2 = real(h(7))+real(h(24))
   xi2 = aimag(h(7))+aimag(h(24))
   yr2 = real(h(7))-real(h(24))
   yi2 = aimag(h(7))-aimag(h(24))
   xr3 = real(h(18))+real(h(13))
   xi3 = aimag(h(18))+aimag(h(13))
   yr3 = real(h(18))-real(h(13))
   yi3 = aimag(h(18))-aimag(h(13))
   xr4 = real(h(2))+real(h(29))
   xi4 = aimag(h(2))+aimag(h(29))
   yr4 = real(h(2))-real(h(29))
   yi4 = aimag(h(2))-aimag(h(29))
   xr5 = real(h(14))+real(h(17))
   xi5 = aimag(h(14))+aimag(h(17))
   yr5 = real(h(14))-real(h(17))
   yi5 = aimag(h(14))-aimag(h(17))
   xr6 = real(h(5))+real(h(26))
   xi6 = aimag(h(5))+aimag(h(26))
   yr6 = real(h(5))-real(h(26))
   yi6 = aimag(h(5))-aimag(h(26))
   xr7 = real(h(4))+real(h(27))
   xi7 = aimag(h(4))+aimag(h(27))
   yr7 = real(h(4))-real(h(27))
   yi7 = aimag(h(4))-aimag(h(27))
   xr8 = real(h(28))+real(h(3))
   xi8 = aimag(h(28))+aimag(h(3))
   yr8 = real(h(28))-real(h(3))
   yi8 = aimag(h(28))-aimag(h(3))
   xr9 = real(h(10))+real(h(21))
   xi9 = aimag(h(10))+aimag(h(21))
   yr9 = real(h(10))-real(h(21))
   yi9 = aimag(h(10))-aimag(h(21))
   xra = real(h(8))+real(h(23))
   xia = aimag(h(8))+aimag(h(23))
   yra = real(h(8))-real(h(23))
   yia = aimag(h(8))-aimag(h(23))
   xrb = real(h(25))+real(h(6))
   xib = aimag(h(25))+aimag(h(6))
   yrb = real(h(25))-real(h(6))
   yib = aimag(h(25))-aimag(h(6))
   xrc = real(h(20))+real(h(11))
   xic = aimag(h(20))+aimag(h(11))
   yrc = real(h(20))-real(h(11))
   yic = aimag(h(20))-aimag(h(11))
   xrd = real(h(16))+real(h(15))
   xid = aimag(h(16))+aimag(h(15))
   yrd = real(h(16))-real(h(15))
   yid = aimag(h(16))-aimag(h(15))
   xre = real(h(19))+real(h(12))
   xie = aimag(h(19))+aimag(h(12))
   yre = real(h(19))-real(h(12))
   yie = aimag(h(19))-aimag(h(12))
   xrf = real(h(9))+real(h(22))
   xif = aimag(h(9))+aimag(h(22))
   yrf = real(h(9))-real(h(22))
   yif = aimag(h(9))-aimag(h(22))
   xr16b = xr1+xr6+xrb
   xi16b = xi1+xi6+xib
   xr49e = xr4+xr9+xre
   xi49e = xi4+xi9+xie
   xr27c = xr2+xr7+xrc
   xi27c = xi2+xi7+xic
   xr5af = xr5+xra+xrf
   xi5af = xi5+xia+xif
   xr38d = xr3+xr8+xrd
   xi38d = xi3+xi8+xid
   xr7_2 = xr7-xr2
   xi7_2 = xi7-xi2
   xr8_3 = xr8-xr3
   xi8_3 = xi8-xi3
   xra_5 = xra-xr5
   xia_5 = xia-xi5
   xrb_6 = xrb-xr6
   xib_6 = xib-xi6
   xrd_8 = xrd-xr8
   xid_8 = xid-xi8
   xre_9 = xre-xr9
   xie_9 = xie-xi9
   xr1_b = xr1-xrb
   xi1_b = xi1-xib
   xr2_c = xr2-xrc
   xi2_c = xi2-xic
   xr4_e = xr4-xre
   xi4_e = xi4-xie
   xr5_f = xr5-xrf
   xi5_f = xi5-xif
   xr16b = xr16b-xr38d
   xi16b = xi16b-xi38d
   xr49e = xr49e-xr38d
   xi49e = xi49e-xi38d
   xr27c = xr27c-xr38d
   xi27c = xi27c-xi38d
   xr5af = xr5af-xr38d
   xi5af = xi5af-xi38d
   xr16b49e = xr16b+xr49e
   xi16b49e = xi16b+xi49e
   xr27c5af = xr27c+xr5af
   xi27c5af = xi27c+xi5af
   xrsigma = xr16b49e+xr27c5af
   xisigma = xi16b49e+xi27c5af
   ursigma = xrsigma+5*xr38d
   uisigma = xisigma+5*xi38d
   sr0 = xr0+ursigma
   si0 = xi0+uisigma
   ur1234 = xr0+csigma*xr38d+dc1*xrsigma
   ui1234 = xi0+csigma*xi38d+dc1*xisigma
   ur12 = ur1234+dc2*xr27c5af
   ui12 = ui1234+dc2*xi27c5af
   ur34 = ur1234+dc7*xr16b49e
   ui34 = ui1234+dc7*xi16b49e
   ur1 = ur12+dc3*xr49e+dc4*xr5af
   ui1 = ui12+dc3*xi49e+dc4*xi5af
   ur2 = ur12+dc5*xr16b+dc6*xr27c
   ui2 = ui12+dc5*xi16b+dc6*xi27c
   ur3 = ur34+dc8*xr49e+dc3*xr5af
   ui3 = ui34+dc8*xi49e+dc3*xi5af
   ur4 = ur34-dc4*xr16b+dc5*xr27c
   ui4 = ui34-dc4*xi16b+dc5*xi27c
   ur5 = csigma*ursigma-(ur1+ur2+ur3+ur4)+5*xr0
   ui5 = csigma*uisigma-(ui1+ui2+ui3+ui4)+5*xi0
   xra4_5e = xra_5+xr4_e
   xia4_5e = xia_5+xi4_e
   xrae_54 = xra_5-xr4_e
   xiae_54 = xia_5-xi4_e
   xrd1_8b = xrd_8+xr1_b
   xid1_8b = xid_8+xi1_b
   xrdb_81 = xrd_8-xr1_b
   xidb_81 = xid_8-xi1_b
   xra4d1_5e8b = xra4_5e+xrd1_8b
   xia4d1_5e8b = xia4_5e+xid1_8b
   pr6 = xr7_2+xra4d1_5e8b
   pi6 = xi7_2+xia4d1_5e8b
   pr8a = c50*xr7_2+c51*xra4d1_5e8b
   pi8a = c50*xi7_2+c51*xia4d1_5e8b
   pr8 = pr8a+xrd1_8b
   pi8 = pi8a+xid1_8b
   pra = pr8a+xra4_5e
   pia = pi8a+xia4_5e
   prc = xrae_54+s52*xrdb_81
   pic = xiae_54+s52*xidb_81
   pre = s52*xrae_54-xrdb_81
   pie = s52*xiae_54-xidb_81
   xrb5_6f = xrb_6+xr5_f
   xib5_6f = xib_6+xi5_f
   xrbf_65 = xrb_6-xr5_f
   xibf_65 = xib_6-xi5_f
   xre2_9c = xre_9+xr2_c
   xie2_9c = xie_9+xi2_c
   xrec_92 = xre_9-xr2_c
   xiec_92 = xie_9-xi2_c
   xrb5e2_6f9c = xrb5_6f+xre2_9c
   xib5e2_6f9c = xib5_6f+xie2_9c
   pr7 = xr8_3+xrb5e2_6f9c
   pi7 = xi8_3+xib5e2_6f9c
   pr9b = c50*xr8_3+c51*xrb5e2_6f9c
   pi9b = c50*xi8_3+c51*xib5e2_6f9c
   pr9 = pr9b+xre2_9c
   pi9 = pi9b+xie2_9c
   prb = pr9b+xrb5_6f
   pib = pi9b+xib5_6f
   prd = xrbf_65+s52*xrec_92
   pid = xibf_65+s52*xiec_92
   prf = s52*xrbf_65-xrec_92
   pif = s52*xibf_65-xiec_92
   ur6 = dc9*pr6+dc10*pr7
   ui6 = dc9*pi6+dc10*pi7
   ur7 = dc11*pr6+dc9*pr7
   ui7 = dc11*pi6+dc9*pi7
   pr89 = pr8+pr9
   pi89 = pi8+pi9
   prcd = prc+prd
   picd = pic+pid
   ur89 = dc12*pr89+dc15*prcd
   ui89 = dc12*pi89+dc15*picd
   ur8 = ur89+dc13*pr9+dc16*prd
   ui8 = ui89+dc13*pi9+dc16*pid
   ur9 = ur89+dc14*pr8+dc17*prc
   ui9 = ui89+dc14*pi8+dc17*pic
   urcd = dc18*prcd-dc15*pr89
   uicd = dc18*picd-dc15*pi89
   urc = urcd-dc16*pr9+dc19*prd
   uic = uicd-dc16*pi9+dc19*pid
   urd = urcd-dc17*pr8+dc20*prc
   uid = uicd-dc17*pi8+dc20*pic
   prab = pra+prb
   piab = pia+pib
   pref = pre+prf
   pief = pie+pif
   urab = dc21*prab+dc24*pref
   uiab = dc21*piab+dc24*pief
   ura = urab+dc22*prb+dc25*prf
   uia = uiab+dc22*pib+dc25*pif
   urb = urab+dc23*pra+dc26*pre
   uib = uiab+dc23*pia+dc26*pie
   uref = dc27*pref-dc24*prab
   uief = dc27*pief-dc24*piab
   ure = uref-dc25*prb+dc28*prf
   uie = uief-dc25*pib+dc28*pif
   urf = uref-dc26*pra+dc29*pre
   uif = uief-dc26*pia+dc29*pie
   ur8a = ur8+ura
   ui8a = ui8+uia
   wr6 = ur6+c50*ur8a
   wi6 = ui6+c50*ui8a
   wr8eac = ur6+c51*ur8a
   wi8eac = ui6+c51*ui8a
   wr8e = wr8eac+ura
   wi8e = wi8eac+uia
   wrac = wr8eac+ur8
   wiac = wi8eac+ui8
   wr8_e = urc+s52*ure
   wi8_e = uic+s52*uie
   wra_c = s52*urc-ure
   wia_c = s52*uic-uie
   wr8 = wr8e+wr8_e
   wi8 = wi8e+wi8_e
   wre = wr8e-wr8_e
   wie = wi8e-wi8_e
   wra = wrac+wra_c
   wia = wiac+wia_c
   wrc = wrac-wra_c
   wic = wiac-wia_c
   ur9b = ur9+urb
   ui9b = ui9+uib
   wr7 = ur7+c50*ur9b
   wi7 = ui7+c50*ui9b
   wr9fbd = ur7+c51*ur9b
   wi9fbd = ui7+c51*ui9b
   wr9f = wr9fbd+urb
   wi9f = wi9fbd+uib
   wrbd = wr9fbd+ur9
   wibd = wi9fbd+ui9
   wr9_f = urd+s52*urf
   wi9_f = uid+s52*uif
   wrb_d = s52*urd-urf
   wib_d = s52*uid-uif
   wr9 = wr9f+wr9_f
   wi9 = wi9f+wi9_f
   wrf = wr9f-wr9_f
   wif = wi9f-wi9_f
   wrb = wrbd+wrb_d
   wib = wibd+wib_d
   wrd = wrbd-wrb_d
   wid = wibd-wib_d
   sr1 = ur1+wrc
   si1 = ui1+wic
   sr2 = ur3-wr6+wrd
   si2 = ui3-wi6+wid
   sr3 = ur5-wr7
   si3 = ui5-wi7
   sr4 = ur2+wre
   si4 = ui2+wie
   sr5 = ur4-wr8+wrf
   si5 = ui4-wi8+wif
   sr6 = ur1-wr9
   si6 = ui1-wi9
   sr7 = ur3+wr6
   si7 = ui3+wi6
   sr8 = ur5+wr7-wra
   si8 = ui5+wi7-wia
   sr9 = ur2-wrb
   si9 = ui2-wib
   sra = ur4+wr8
   sia = ui4+wi8
   srb = ur1+wr9-wrc
   sib = ui1+wi9-wic
   src = ur3-wrd
   sic = ui3-wid
   srd = ur5+wra
   sid = ui5+wia
   sre = ur2+wrb-wre
   sie = ui2+wib-wie
   srf = ur4-wrf
   sif = ui4-wif

   yr16b = yr1+yr6+yrb
   yi16b = yi1+yi6+yib
   yr49e = yr4+yr9+yre
   yi49e = yi4+yi9+yie
   yr27c = yr2+yr7+yrc
   yi27c = yi2+yi7+yic
   yr5af = yr5+yra+yrf
   yi5af = yi5+yia+yif
   yr38d = yr3+yr8+yrd
   yi38d = yi3+yi8+yid
   yr7_2 = yr7-yr2
   yi7_2 = yi7-yi2
   yr8_3 = yr8-yr3
   yi8_3 = yi8-yi3
   yra_5 = yra-yr5
   yia_5 = yia-yi5
   yrb_6 = yrb-yr6
   yib_6 = yib-yi6
   yrd_8 = yrd-yr8
   yid_8 = yid-yi8
   yre_9 = yre-yr9
   yie_9 = yie-yi9
   yr1_b = yr1-yrb
   yi1_b = yi1-yib
   yr2_c = yr2-yrc
   yi2_c = yi2-yic
   yr4_e = yr4-yre
   yi4_e = yi4-yie
   yr5_f = yr5-yrf
   yi5_f = yi5-yif
   yr16b = yr16b-yr38d
   yi16b = yi16b-yi38d
   yr49e = yr49e-yr38d
   yi49e = yi49e-yi38d
   yr27c = yr27c-yr38d
   yi27c = yi27c-yi38d
   yr5af = yr5af-yr38d
   yi5af = yi5af-yi38d
   yr16b49e = yr16b+yr49e
   yi16b49e = yi16b+yi49e
   yr27c5af = yr27c+yr5af
   yi27c5af = yi27c+yi5af
   yrsigma = yr16b49e+yr27c5af
   yisigma = yi16b49e+yi27c5af
   vrsigma = yrsigma+5*yr38d
   visigma = yisigma+5*yi38d
   vr1234 = ssigma*yr38d+ds1*yrsigma
   vi1234 = ssigma*yi38d+ds1*yisigma
   vr12 = vr1234+ds2*yr27c5af
   vi12 = vi1234+ds2*yi27c5af
   vr34 = vr1234+ds7*yr16b49e
   vi34 = vi1234+ds7*yi16b49e
   vr1 = vr12+ds3*yr49e+ds4*yr5af
   vi1 = vi12+ds3*yi49e+ds4*yi5af
   vr2 = vr12+ds5*yr16b+ds6*yr27c
   vi2 = vi12+ds5*yi16b+ds6*yi27c
   vr3 = vr34+ds8*yr49e+ds3*yr5af
   vi3 = vi34+ds8*yi49e+ds3*yi5af
   vr4 = vr34-ds4*yr16b+ds5*yr27c
   vi4 = vi34-ds4*yi16b+ds5*yi27c
   vr5 = ssigma*vrsigma-(vr1+vr2+vr3+vr4)
   vi5 = ssigma*visigma-(vi1+vi2+vi3+vi4)
   yra4_5e = yra_5+yr4_e
   yia4_5e = yia_5+yi4_e
   yrae_54 = yra_5-yr4_e
   yiae_54 = yia_5-yi4_e
   yrd1_8b = yrd_8+yr1_b
   yid1_8b = yid_8+yi1_b
   yrdb_81 = yrd_8-yr1_b
   yidb_81 = yid_8-yi1_b
   yra4d1_5e8b = yra4_5e+yrd1_8b
   yia4d1_5e8b = yia4_5e+yid1_8b
   qr6 = yr7_2+yra4d1_5e8b
   qi6 = yi7_2+yia4d1_5e8b
   qr8a = c50*yr7_2+c51*yra4d1_5e8b
   qi8a = c50*yi7_2+c51*yia4d1_5e8b
   qr8 = qr8a+yrd1_8b
   qi8 = qi8a+yid1_8b
   qra = qr8a+yra4_5e
   qia = qi8a+yia4_5e
   qrc = yrae_54+s52*yrdb_81
   qic = yiae_54+s52*yidb_81
   qre = s52*yrae_54-yrdb_81
   qie = s52*yiae_54-yidb_81
   yrb5_6f = yrb_6+yr5_f
   yib5_6f = yib_6+yi5_f
   yrbf_65 = yrb_6-yr5_f
   yibf_65 = yib_6-yi5_f
   yre2_9c = yre_9+yr2_c
   yie2_9c = yie_9+yi2_c
   yrec_92 = yre_9-yr2_c
   yiec_92 = yie_9-yi2_c
   yrb5e2_6f9c = yrb5_6f+yre2_9c
   yib5e2_6f9c = yib5_6f+yie2_9c
   qr7 = yr8_3+yrb5e2_6f9c
   qi7 = yi8_3+yib5e2_6f9c
   qr9b = c50*yr8_3+c51*yrb5e2_6f9c
   qi9b = c50*yi8_3+c51*yib5e2_6f9c
   qr9 = qr9b+yre2_9c
   qi9 = qi9b+yie2_9c
   qrb = qr9b+yrb5_6f
   qib = qi9b+yib5_6f
   qrd = yrbf_65+s52*yrec_92
   qid = yibf_65+s52*yiec_92
   qrf = s52*yrbf_65-yrec_92
   qif = s52*yibf_65-yiec_92
   vr6 = ds9*qr6+ds10*qr7
   vi6 = ds9*qi6+ds10*qi7
   vr7 = ds11*qr6+ds9*qr7
   vi7 = ds11*qi6+ds9*qi7
   qr89 = qr8+qr9
   qi89 = qi8+qi9
   qrcd = qrc+qrd
   qicd = qic+qid
   vr89 = ds12*qr89+ds15*qrcd
   vi89 = ds12*qi89+ds15*qicd
   vr8 = vr89+ds13*qr9+ds16*qrd
   vi8 = vi89+ds13*qi9+ds16*qid
   vr9 = vr89+ds14*qr8+ds17*qrc
   vi9 = vi89+ds14*qi8+ds17*qic
   vrcd = ds18*qrcd-ds15*qr89
   vicd = ds18*qicd-ds15*qi89
   vrc = vrcd-ds16*qr9+ds19*qrd
   vic = vicd-ds16*qi9+ds19*qid
   vrd = vrcd-ds17*qr8+ds20*qrc
   vid = vicd-ds17*qi8+ds20*qic
   qrab = qra+qrb
   qiab = qia+qib
   qref = qre+qrf
   qief = qie+qif
   vrab = ds21*qrab+ds24*qref
   viab = ds21*qiab+ds24*qief
   vra = vrab+ds22*qrb+ds25*qrf
   via = viab+ds22*qib+ds25*qif
   vrb = vrab+ds23*qra+ds26*qre
   vib = viab+ds23*qia+ds26*qie
   vref = ds27*qref-ds24*qrab
   vief = ds27*qief-ds24*qiab
   vre = vref-ds25*qrb+ds28*qrf
   vie = vief-ds25*qib+ds28*qif
   vrf = vref-ds26*qra+ds29*qre
   vif = vief-ds26*qia+ds29*qie
   vr8a = vr8+vra
   vi8a = vi8+via
   zr6 = vr6+c50*vr8a
   zi6 = vi6+c50*vi8a
   zr8eac = vr6+c51*vr8a
   zi8eac = vi6+c51*vi8a
   zr8e = zr8eac+vra
   zi8e = zi8eac+via
   zrac = zr8eac+vr8
   ziac = zi8eac+vi8
   zr8_e = vrc+s52*vre
   zi8_e = vic+s52*vie
   zra_c = s52*vrc-vre
   zia_c = s52*vic-vie
   zr8 = zr8e+zr8_e
   zi8 = zi8e+zi8_e
   zre = zr8e-zr8_e
   zie = zi8e-zi8_e
   zra = zrac+zra_c
   zia = ziac+zia_c
   zrc = zrac-zra_c
   zic = ziac-zia_c
   vr9b = vr9+vrb
   vi9b = vi9+vib
   zr7 = vr7+c50*vr9b
   zi7 = vi7+c50*vi9b
   zr9fbd = vr7+c51*vr9b
   zi9fbd = vi7+c51*vi9b
   zr9f = zr9fbd+vrb
   zi9f = zi9fbd+vib
   zrbd = zr9fbd+vr9
   zibd = zi9fbd+vi9
   zr9_f = vrd+s52*vrf
   zi9_f = vid+s52*vif
   zrb_d = s52*vrd-vrf
   zib_d = s52*vid-vif
   zr9 = zr9f+zr9_f
   zi9 = zi9f+zi9_f
   zrf = zr9f-zr9_f
   zif = zi9f-zi9_f
   zrb = zrbd+zrb_d
   zib = zibd+zib_d
   zrd = zrbd-zrb_d
   zid = zibd-zib_d
   tr1 = vr1+zrc
   ti1 = vi1+zic
   tr2 = vr3-zr6+zrd
   ti2 = vi3-zi6+zid
   tr3 = vr5-zr7
   ti3 = vi5-zi7
   tr4 = vr2+zre
   ti4 = vi2+zie
   tr5 = vr4-zr8+zrf
   ti5 = vi4-zi8+zif
   tr6 = vr1-zr9
   ti6 = vi1-zi9
   tr7 = vr3+zr6
   ti7 = vi3+zi6
   tr8 = vr5+zr7-zra
   ti8 = vi5+zi7-zia
   tr9 = vr2-zrb
   ti9 = vi2-zib
   tra = vr4+zr8
   tia = vi4+zi8
   trb = vr1+zr9-zrc
   tib = vi1+zi9-zic
   trc = vr3-zrd
   tic = vi3-zid
   trd = vr5+zra
   tid = vi5+zia
   tre = vr2+zrb-zre
   tie = vi2+zib-zie
   trf = vr4-zrf
   tif = vi4-zif

   h(0) = cmplx(sr0,si0,wp)
   h(1) = cmplx(srf+tif,sif-trf,wp)
   h(2) = cmplx(src+tic,sic-trc,wp)
   h(3) = cmplx(sr8-ti8,si8+tr8,wp)
   h(4) = cmplx(sr9+ti9,si9-tr9,wp)
   h(5) = cmplx(sra+tia,sia-tra,wp)
   h(6) = cmplx(sr5-ti5,si5+tr5,wp)
   h(7) = cmplx(sre+tie,sie-tre,wp)
   h(8) = cmplx(sr6+ti6,si6-tr6,wp)
   h(9) = cmplx(sr1+ti1,si1-tr1,wp)
   h(10) = cmplx(sr7+ti7,si7-tr7,wp)
   h(11) = cmplx(sr4-ti4,si4+tr4,wp)
   h(12) = cmplx(sr2-ti2,si2+tr2,wp)
   h(13) = cmplx(srd-tid,sid+trd,wp)
   h(14) = cmplx(srb+tib,sib-trb,wp)
   h(15) = cmplx(sr3-ti3,si3+tr3,wp)
   h(16) = cmplx(sr3+ti3,si3-tr3,wp)
   h(17) = cmplx(srb-tib,sib+trb,wp)
   h(18) = cmplx(srd+tid,sid-trd,wp)
   h(19) = cmplx(sr2+ti2,si2-tr2,wp)
   h(20) = cmplx(sr4+ti4,si4-tr4,wp)
   h(21) = cmplx(sr7-ti7,si7+tr7,wp)
   h(22) = cmplx(sr1-ti1,si1+tr1,wp)
   h(23) = cmplx(sr6-ti6,si6+tr6,wp)
   h(24) = cmplx(sre-tie,sie+tre,wp)
   h(25) = cmplx(sr5+ti5,si5-tr5,wp)
   h(26) = cmplx(sra-tia,sia+tra,wp)
   h(27) = cmplx(sr9-ti9,si9+tr9,wp)
   h(28) = cmplx(sr8+ti8,si8-tr8,wp)
   h(29) = cmplx(src-tic,sic+trc,wp)
   h(30) = cmplx(srf-tif,sif+trf,wp)
end subroutine fft31t