! File: fft31a.i90
! Public domain 2004 James Van Buskirk
subroutine fft31a(h)
! 632 adds, 344 muls
   complex(wp), intent(inout) :: h(0:30)
   real(wp) xr1, xi1, yr1, yi1
   real(wp) xr2, xi2, yr2, yi2
   real(wp) xr3, xi3, yr3, yi3
   real(wp) xr4, xi4, yr4, yi4
   real(wp) xr5, xi5, yr5, yi5
   real(wp) xr6, xi6, yr6, yi6
   real(wp) xr7, xi7, yr7, yi7
   real(wp) xr8, xi8, yr8, yi8
   real(wp) xr9, xi9, yr9, yi9
   real(wp) xra, xia, yra, yia
   real(wp) xrb, xib, yrb, yib
   real(wp) xrc, xic, yrc, yic
   real(wp) xrd, xid, yrd, yid
   real(wp) xre, xie, yre, yie
   real(wp) xrf, xif, yrf, yif
   real(wp) xr4d, xi4d, xr4_d, xi4_d
   real(wp) xr7a, xi7a, xr7_a, xi7_a
   real(wp) xr47ad, xi47ad, ur1, ui1
   real(wp) ur4a, ui4a, ur4, ui4, ura, uia
   real(wp) ur7, ui7, urd, uid
   real(wp) xr5e, xi5e, xr5_e, xi5_e
   real(wp) xr8b, xi8b, xr8_b, xi8_b
   real(wp) xr58be, xi58be, ur2, ui2
   real(wp) ur5b, ui5b, ur5, ui5, urb, uib
   real(wp) ur8, ui8, ure, uie
   real(wp) xr6f, xi6f, xr6_f, xi6_f
   real(wp) xr9c, xi9c, xr9_c, xi9_c
   real(wp) xr69cf, xi69cf, ur3, ui3
   real(wp) ur6c, ui6c, ur6, ui6, urc, uic
   real(wp) ur9, ui9, urf, uif
   real(wp) ur123, ui123, sr0, si0, pr123, pi123
   real(wp) pr1, pi1, pr2, pi2, pr3, pi3
   real(wp) ur7_4, ui7_4, ur8_5, ui8_5, ur9_6, ui9_6
   real(wp) pr47, pi47, pr4, pi4, pr7, pi7
   real(wp) pr58, pi58, pr5, pi5, pr8, pi8
   real(wp) pr69, pi69, pr6, pi6, pr9, pi9
   real(wp) urd_a, uid_a, ure_b, uie_b, urf_c, uif_c
   real(wp) prad, piad, pra, pia, prd, pid
   real(wp) prbe, pibe, prb, pib, pre, pie
   real(wp) prcf, picf, prc, pic, prf, pif
   real(wp) pr4a, pi4a, sr1, si1, sr47ad, si47ad
   real(wp) sr4d, si4d, sr4_d, si4_d
   real(wp) sr4, si4, srd, sid
   real(wp) sr7a, si7a, sr7_a, si7_a
   real(wp) sr7, si7, sra, sia
   real(wp) pr5b, pi5b, sr2, si2, sr58be, si58be
   real(wp) sr5e, si5e, sr5_e, si5_e
   real(wp) sr5, si5, sre, sie
   real(wp) sr8b, si8b, sr8_b, si8_b
   real(wp) sr8, si8, srb, sib
   real(wp) pr6c, pi6c, sr3, si3, sr69cf, si69cf
   real(wp) sr6f, si6f, sr6_f, si6_f
   real(wp) sr6, si6, srf, sif
   real(wp) sr9c, si9c, sr9_c, si9_c
   real(wp) sr9, si9, src, sic
   real(wp) yr4d, yi4d, yr4_d, yi4_d
   real(wp) yr7a, yi7a, yr7_a, yi7_a
   real(wp) yr47ad, yi47ad, vr1, vi1
   real(wp) vr4a, vi4a, vr4, vi4, vra, via
   real(wp) vr7, vi7, vrd, vid
   real(wp) yr5e, yi5e, yr5_e, yi5_e
   real(wp) yr8b, yi8b, yr8_b, yi8_b
   real(wp) yr58be, yi58be, vr2, vi2
   real(wp) vr5b, vi5b, vr5, vi5, vrb, vib
   real(wp) vr8, vi8, vre, vie
   real(wp) yr6f, yi6f, yr6_f, yi6_f
   real(wp) yr9c, yi9c, yr9_c, yi9_c
   real(wp) yr69cf, yi69cf, vr3, vi3
   real(wp) vr6c, vi6c, vr6, vi6, vrc, vic
   real(wp) vr9, vi9, vrf, vif
   real(wp) qr1, qi1, qr2, qi2, qr3, qi3
   real(wp) vr7_4, vi7_4, vr8_5, vi8_5, vr9_6, vi9_6
   real(wp) qr47, qi47, qr4, qi4, qr7, qi7
   real(wp) qr58, qi58, qr5, qi5, qr8, qi8
   real(wp) qr69, qi69, qr6, qi6, qr9, qi9
   real(wp) vrd_a, vid_a, vre_b, vie_b, vrf_c, vif_c
   real(wp) qrad, qiad, qra, qia, qrd, qid
   real(wp) qrbe, qibe, qrb, qib, qre, qie
   real(wp) qrcf, qicf, qrc, qic, qrf, qif
   real(wp) qr4a, qi4a, tr1, ti1, tr47ad, ti47ad
   real(wp) tr4d, ti4d, tr4_d, ti4_d
   real(wp) tr4, ti4, trd, tid
   real(wp) tr7a, ti7a, tr7_a, ti7_a
   real(wp) tr7, ti7, tra, tia
   real(wp) qr5b, qi5b, tr2, ti2, tr58be, ti58be
   real(wp) tr5e, ti5e, tr5_e, ti5_e
   real(wp) tr5, ti5, tre, tie
   real(wp) tr8b, ti8b, tr8_b, ti8_b
   real(wp) tr8, ti8, trb, tib
   real(wp) qr6c, qi6c, tr3, ti3, tr69cf, ti69cf
   real(wp) tr6f, ti6f, tr6_f, ti6_f
   real(wp) tr6, ti6, trf, tif
   real(wp) tr9c, ti9c, tr9_c, ti9_c
   real(wp) tr9, ti9, trc, tic

   real(wp),parameter :: c1 = -.8944271909999158785636694674925104_wp
   real(wp),parameter :: c2 = -.2763932022500210303590826331268724_wp
   real(wp),parameter :: s1 = .6180339887498948482045868343656382_wp
   real(wp),parameter :: dc1 = -.3297070544363274708879342919056026_wp
   real(wp),parameter :: dc2 = .6380942903798882373411255424134321_wp
   real(wp),parameter :: dc3 = .2510268729290941753226773333033754_wp
   real(wp),parameter :: dc4 = -.3056489870319706225490819189370370_wp
   real(wp),parameter :: dc5 = -.0836204891340032060287324012733157_wp
   real(wp),parameter :: dc6 = .4663103174051567617121751092152637_wp
   real(wp),parameter :: dc7 = .1185198271993044376989431474045266_wp
   real(wp),parameter :: dc8 = 1.2587167895346618298423896206185935_wp
   real(wp),parameter :: dc9 = .7065344455204369076802243236418372_wp
   real(wp),parameter :: dc10 = -.5229597775845896829713923186462800_wp
   real(wp),parameter :: dc11 = -.7713286973491624553605992621309607_wp
   real(wp),parameter :: dc12 = .3432383388931257185335662802706080_wp
   real(wp),parameter :: dc13 = -.5375570547960443575471925003995269_wp
   real(wp),parameter :: dc14 = -.8979444240956453223528385703234650_wp
   real(wp),parameter :: dc15 = .2767783045476210221271537127440506_wp
   real(wp),parameter :: dc16 = -.4545858312766339718632518226423785_wp
   real(wp),parameter :: dc17 = -.8315695740161540356642021423935391_wp
   real(wp),parameter :: dc18 = .1817978323890879742278485954562746_wp
   real(wp),parameter :: dc19 = .6689938540280712956570326234710885_wp
   real(wp),parameter :: dc20 = .0672653473106970601112424697542567_wp
   real(wp),parameter :: dc21 = -1.1044086918887071759621095370958563_wp
   real(wp),parameter :: ds1 = -.2876896432925834610922244441472298_wp
   real(wp),parameter :: ds2 = -.4435990695229438494937820909950198_wp
   real(wp),parameter :: ds3 = .1745122765325251183740594052503946_wp
   real(wp),parameter :: ds4 = .2998145478576453478122750014261742_wp
   real(wp),parameter :: ds5 = .0372577528951103256163666300172535_wp
   real(wp),parameter :: ds6 = -.4898024294879831247820037838217848_wp
   real(wp),parameter :: ds7 = -.4867290865768665735673427785714679_wp
   real(wp),parameter :: ds8 = -.4639301227436081939265253137860307_wp
   real(wp),parameter :: ds9 = .5032265717882682551371923740675969_wp
   real(wp),parameter :: ds10 = -.9968303766175870335965646301934142_wp
   real(wp),parameter :: ds11 = -.7889681380851011225379679528474900_wp
   real(wp),parameter :: ds12 = 1.1472409054152214516271510263778878_wp
   real(wp),parameter :: ds13 = .9977591866678393034097541872801121_wp
   real(wp),parameter :: ds14 = .4976723612826765805499394917955183_wp
   real(wp),parameter :: ds15 = -1.3596254129817497946555050772146859_wp
   real(wp),parameter :: ds16 = .3140391925370850578721469975569951_wp
   real(wp),parameter :: ds17 = .3527942467158582649715245177385254_wp
   real(wp),parameter :: ds18 = -.2859957718861560364921745106621358_wp
   real(wp),parameter :: ds19 = -.1293773399525311412341577663709124_wp
   real(wp),parameter :: ds20 = .2588356792521043090571708180839408_wp
   real(wp),parameter :: ds21 = .4102909196240484871164009213276759_wp
   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(25))+real(h(6))
   xi2 = aimag(h(25))+aimag(h(6))
   yr2 = real(h(25))-real(h(6))
   yi2 = aimag(h(25))-aimag(h(6))
   xr3 = real(h(5))+real(h(26))
   xi3 = aimag(h(5))+aimag(h(26))
   yr3 = real(h(5))-real(h(26))
   yi3 = aimag(h(5))-aimag(h(26))
   xr4 = real(h(16))+real(h(15))
   xi4 = aimag(h(16))+aimag(h(15))
   yr4 = real(h(16))-real(h(15))
   yi4 = aimag(h(16))-aimag(h(15))
   xr5 = real(h(28))+real(h(3))
   xi5 = aimag(h(28))+aimag(h(3))
   yr5 = real(h(28))-real(h(3))
   yi5 = aimag(h(28))-aimag(h(3))
   xr6 = real(h(18))+real(h(13))
   xi6 = aimag(h(18))+aimag(h(13))
   yr6 = real(h(18))-real(h(13))
   yi6 = aimag(h(18))-aimag(h(13))
   xr7 = real(h(8))+real(h(23))
   xi7 = aimag(h(8))+aimag(h(23))
   yr7 = real(h(8))-real(h(23))
   yi7 = aimag(h(8))-aimag(h(23))
   xr8 = real(h(14))+real(h(17))
   xi8 = aimag(h(14))+aimag(h(17))
   yr8 = real(h(14))-real(h(17))
   yi8 = aimag(h(14))-aimag(h(17))
   xr9 = real(h(9))+real(h(22))
   xi9 = aimag(h(9))+aimag(h(22))
   yr9 = real(h(9))-real(h(22))
   yi9 = aimag(h(9))-aimag(h(22))
   xra = real(h(4))+real(h(27))
   xia = aimag(h(4))+aimag(h(27))
   yra = real(h(4))-real(h(27))
   yia = aimag(h(4))-aimag(h(27))
   xrb = real(h(7))+real(h(24))
   xib = aimag(h(7))+aimag(h(24))
   yrb = real(h(7))-real(h(24))
   yib = aimag(h(7))-aimag(h(24))
   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(2))+real(h(29))
   xid = aimag(h(2))+aimag(h(29))
   yrd = real(h(2))-real(h(29))
   yid = aimag(h(2))-aimag(h(29))
   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(10))+real(h(21))
   xif = aimag(h(10))+aimag(h(21))
   yrf = real(h(10))-real(h(21))
   yif = aimag(h(10))-aimag(h(21))
   xr4d = xr4+xrd
   xi4d = xi4+xid
   xr4_d = xr4-xrd
   xi4_d = xi4-xid
   xr7a = xr7+xra
   xi7a = xi7+xia
   xr7_a = xr7-xra
   xi7_a = xi7-xia
   xr47ad = xr4d+xr7a
   xi47ad = xi4d+xi7a
   ur1 = xr1+xr47ad
   ui1 = xi1+xi47ad
   ur4a = c1*xr1+c2*xr47ad
   ui4a = c1*xi1+c2*xi47ad
   ur4 = ur4a+xr7a
   ui4 = ui4a+xi7a
   ura = xr4d+ur4a
   uia = xi4d+ui4a
   ur7 = xr4_d+s1*xr7_a
   ui7 = xi4_d+s1*xi7_a
   urd = s1*xr4_d-xr7_a
   uid = s1*xi4_d-xi7_a
   xr5e = xr5+xre
   xi5e = xi5+xie
   xr5_e = xr5-xre
   xi5_e = xi5-xie
   xr8b = xr8+xrb
   xi8b = xi8+xib
   xr8_b = xr8-xrb
   xi8_b = xi8-xib
   xr58be = xr5e+xr8b
   xi58be = xi5e+xi8b
   ur2 = xr2+xr58be
   ui2 = xi2+xi58be
   ur5b = c1*xr2+c2*xr58be
   ui5b = c1*xi2+c2*xi58be
   ur5 = ur5b+xr8b
   ui5 = ui5b+xi8b
   urb = xr5e+ur5b
   uib = xi5e+ui5b
   ur8 = xr5_e+s1*xr8_b
   ui8 = xi5_e+s1*xi8_b
   ure = s1*xr5_e-xr8_b
   uie = s1*xi5_e-xi8_b
   xr6f = xr6+xrf
   xi6f = xi6+xif
   xr6_f = xr6-xrf
   xi6_f = xi6-xif
   xr9c = xr9+xrc
   xi9c = xi9+xic
   xr9_c = xr9-xrc
   xi9_c = xi9-xic
   xr69cf = xr6f+xr9c
   xi69cf = xi6f+xi9c
   ur3 = xr3+xr69cf
   ui3 = xi3+xi69cf
   ur6c = c1*xr3+c2*xr69cf
   ui6c = c1*xi3+c2*xi69cf
   ur6 = ur6c+xr9c
   ui6 = ui6c+xi9c
   urc = xr6f+ur6c
   uic = xi6f+ui6c
   ur9 = xr6_f+s1*xr9_c
   ui9 = xi6_f+s1*xi9_c
   urf = s1*xr6_f-xr9_c
   uif = s1*xi6_f-xi9_c
   ur123 = ur1+ur2+ur3
   ui123 = ui1+ui2+ui3
   sr0 = real(h(0))+ur123
   si0 = aimag(h(0))+ui123
   pr123 = real(h(0))+dc1*ur123
   pi123 = aimag(h(0))+dc1*ui123
   pr1 = pr123+dc2*ur2+dc3*ur3
   pi1 = pi123+dc2*ui2+dc3*ui3
   pr2 = dc3*ur1+pr123+dc2*ur3
   pi2 = dc3*ui1+pi123+dc2*ui3
   pr3 = dc2*ur1+dc3*ur2+pr123
   pi3 = dc2*ui1+dc3*ui2+pi123
   ur7_4 = ur7-ur4
   ui7_4 = ui7-ui4
   ur8_5 = ur8-ur5
   ui8_5 = ui8-ui5
   ur9_6 = ur9-ur6
   ui9_6 = ui9-ui6
   pr47 = dc4*ur7_4+dc7*ur8_5+dc10*ur9_6
   pi47 = dc4*ui7_4+dc7*ui8_5+dc10*ui9_6
   pr4 = pr47+dc5*ur4+dc8*ur5+dc11*ur6
   pi4 = pi47+dc5*ui4+dc8*ui5+dc11*ui6
   pr7 = pr47+dc6*ur7+dc9*ur8+dc12*ur9
   pi7 = pi47+dc6*ui7+dc9*ui8+dc12*ui9
   pr58 = dc10*ur7_4+dc4*ur8_5+dc7*ur9_6
   pi58 = dc10*ui7_4+dc4*ui8_5+dc7*ui9_6
   pr5 = pr58+dc11*ur4+dc5*ur5+dc8*ur6
   pi5 = pi58+dc11*ui4+dc5*ui5+dc8*ui6
   pr8 = pr58+dc12*ur7+dc6*ur8+dc9*ur9
   pi8 = pi58+dc12*ui7+dc6*ui8+dc9*ui9
   pr69 = dc7*ur7_4+dc10*ur8_5+dc4*ur9_6
   pi69 = dc7*ui7_4+dc10*ui8_5+dc4*ui9_6
   pr6 = pr69+dc8*ur4+dc11*ur5+dc5*ur6
   pi6 = pi69+dc8*ui4+dc11*ui5+dc5*ui6
   pr9 = pr69+dc9*ur7+dc12*ur8+dc6*ur9
   pi9 = pi69+dc9*ui7+dc12*ui8+dc6*ui9
   urd_a = urd-ura
   uid_a = uid-uia
   ure_b = ure-urb
   uie_b = uie-uib
   urf_c = urf-urc
   uif_c = uif-uic
   prad = dc13*urd_a+dc16*ure_b+dc19*urf_c
   piad = dc13*uid_a+dc16*uie_b+dc19*uif_c
   pra = prad+dc14*ura+dc17*urb+dc20*urc
   pia = piad+dc14*uia+dc17*uib+dc20*uic
   prd = prad+dc15*urd+dc18*ure+dc21*urf
   pid = piad+dc15*uid+dc18*uie+dc21*uif
   prbe = dc19*urd_a+dc13*ure_b+dc16*urf_c
   pibe = dc19*uid_a+dc13*uie_b+dc16*uif_c
   prb = prbe+dc20*ura+dc14*urb+dc17*urc
   pib = pibe+dc20*uia+dc14*uib+dc17*uic
   pre = prbe+dc21*urd+dc15*ure+dc18*urf
   pie = pibe+dc21*uid+dc15*uie+dc18*uif
   prcf = dc16*urd_a+dc19*ure_b+dc13*urf_c
   picf = dc16*uid_a+dc19*uie_b+dc13*uif_c
   prc = prcf+dc17*ura+dc20*urb+dc14*urc
   pic = picf+dc17*uia+dc20*uib+dc14*uic
   prf = prcf+dc18*urd+dc21*ure+dc15*urf
   pif = picf+dc18*uid+dc21*uie+dc15*uif
   pr4a = pr4+pra
   pi4a = pi4+pia
   sr1 = pr1+c1*pr4a
   si1 = pi1+c1*pi4a
   sr47ad = pr1+c2*pr4a
   si47ad = pi1+c2*pi4a
   sr4d = sr47ad+pra
   si4d = si47ad+pia
   sr4_d = pr7+s1*prd
   si4_d = pi7+s1*pid
   sr4 = sr4d+sr4_d
   si4 = si4d+si4_d
   srd = sr4d-sr4_d
   sid = si4d-si4_d
   sr7a = sr47ad+pr4
   si7a = si47ad+pi4
   sr7_a = s1*pr7-prd
   si7_a = s1*pi7-pid
   sr7 = sr7a+sr7_a
   si7 = si7a+si7_a
   sra = sr7a-sr7_a
   sia = si7a-si7_a
   pr5b = pr5+prb
   pi5b = pi5+pib
   sr2 = pr2+c1*pr5b
   si2 = pi2+c1*pi5b
   sr58be = pr2+c2*pr5b
   si58be = pi2+c2*pi5b
   sr5e = sr58be+prb
   si5e = si58be+pib
   sr5_e = pr8+s1*pre
   si5_e = pi8+s1*pie
   sr5 = sr5e+sr5_e
   si5 = si5e+si5_e
   sre = sr5e-sr5_e
   sie = si5e-si5_e
   sr8b = sr58be+pr5
   si8b = si58be+pi5
   sr8_b = s1*pr8-pre
   si8_b = s1*pi8-pie
   sr8 = sr8b+sr8_b
   si8 = si8b+si8_b
   srb = sr8b-sr8_b
   sib = si8b-si8_b
   pr6c = pr6+prc
   pi6c = pi6+pic
   sr3 = pr3+c1*pr6c
   si3 = pi3+c1*pi6c
   sr69cf = pr3+c2*pr6c
   si69cf = pi3+c2*pi6c
   sr6f = sr69cf+prc
   si6f = si69cf+pic
   sr6_f = pr9+s1*prf
   si6_f = pi9+s1*pif
   sr6 = sr6f+sr6_f
   si6 = si6f+si6_f
   srf = sr6f-sr6_f
   sif = si6f-si6_f
   sr9c = sr69cf+pr6
   si9c = si69cf+pi6
   sr9_c = s1*pr9-prf
   si9_c = s1*pi9-pif
   sr9 = sr9c+sr9_c
   si9 = si9c+si9_c
   src = sr9c-sr9_c
   sic = si9c-si9_c

   yr4d = yr4+yrd
   yi4d = yi4+yid
   yr4_d = yr4-yrd
   yi4_d = yi4-yid
   yr7a = yr7+yra
   yi7a = yi7+yia
   yr7_a = yr7-yra
   yi7_a = yi7-yia
   yr47ad = yr4d+yr7a
   yi47ad = yi4d+yi7a
   vr1 = yr1+yr47ad
   vi1 = yi1+yi47ad
   vr4a = c1*yr1+c2*yr47ad
   vi4a = c1*yi1+c2*yi47ad
   vr4 = vr4a+yr7a
   vi4 = vi4a+yi7a
   vra = yr4d+vr4a
   via = yi4d+vi4a
   vr7 = yr4_d+s1*yr7_a
   vi7 = yi4_d+s1*yi7_a
   vrd = s1*yr4_d-yr7_a
   vid = s1*yi4_d-yi7_a
   yr5e = yr5+yre
   yi5e = yi5+yie
   yr5_e = yr5-yre
   yi5_e = yi5-yie
   yr8b = yr8+yrb
   yi8b = yi8+yib
   yr8_b = yr8-yrb
   yi8_b = yi8-yib
   yr58be = yr5e+yr8b
   yi58be = yi5e+yi8b
   vr2 = yr2+yr58be
   vi2 = yi2+yi58be
   vr5b = c1*yr2+c2*yr58be
   vi5b = c1*yi2+c2*yi58be
   vr5 = vr5b+yr8b
   vi5 = vi5b+yi8b
   vrb = yr5e+vr5b
   vib = yi5e+vi5b
   vr8 = yr5_e+s1*yr8_b
   vi8 = yi5_e+s1*yi8_b
   vre = s1*yr5_e-yr8_b
   vie = s1*yi5_e-yi8_b
   yr6f = yr6+yrf
   yi6f = yi6+yif
   yr6_f = yr6-yrf
   yi6_f = yi6-yif
   yr9c = yr9+yrc
   yi9c = yi9+yic
   yr9_c = yr9-yrc
   yi9_c = yi9-yic
   yr69cf = yr6f+yr9c
   yi69cf = yi6f+yi9c
   vr3 = yr3+yr69cf
   vi3 = yi3+yi69cf
   vr6c = c1*yr3+c2*yr69cf
   vi6c = c1*yi3+c2*yi69cf
   vr6 = vr6c+yr9c
   vi6 = vi6c+yi9c
   vrc = yr6f+vr6c
   vic = yi6f+vi6c
   vr9 = yr6_f+s1*yr9_c
   vi9 = yi6_f+s1*yi9_c
   vrf = s1*yr6_f-yr9_c
   vif = s1*yi6_f-yi9_c
   qr1 = ds1*vr1+ds2*vr2+ds3*vr3
   qi1 = ds1*vi1+ds2*vi2+ds3*vi3
   qr2 = ds3*vr1+ds1*vr2+ds2*vr3
   qi2 = ds3*vi1+ds1*vi2+ds2*vi3
   qr3 = ds2*vr1+ds3*vr2+ds1*vr3
   qi3 = ds2*vi1+ds3*vi2+ds1*vi3
   vr7_4 = vr7-vr4
   vi7_4 = vi7-vi4
   vr8_5 = vr8-vr5
   vi8_5 = vi8-vi5
   vr9_6 = vr9-vr6
   vi9_6 = vi9-vi6
   qr47 = ds4*vr7_4+ds7*vr8_5+ds10*vr9_6
   qi47 = ds4*vi7_4+ds7*vi8_5+ds10*vi9_6
   qr4 = qr47+ds5*vr4+ds8*vr5+ds11*vr6
   qi4 = qi47+ds5*vi4+ds8*vi5+ds11*vi6
   qr7 = qr47+ds6*vr7+ds9*vr8+ds12*vr9
   qi7 = qi47+ds6*vi7+ds9*vi8+ds12*vi9
   qr58 = ds10*vr7_4+ds4*vr8_5+ds7*vr9_6
   qi58 = ds10*vi7_4+ds4*vi8_5+ds7*vi9_6
   qr5 = qr58+ds11*vr4+ds5*vr5+ds8*vr6
   qi5 = qi58+ds11*vi4+ds5*vi5+ds8*vi6
   qr8 = qr58+ds12*vr7+ds6*vr8+ds9*vr9
   qi8 = qi58+ds12*vi7+ds6*vi8+ds9*vi9
   qr69 = ds7*vr7_4+ds10*vr8_5+ds4*vr9_6
   qi69 = ds7*vi7_4+ds10*vi8_5+ds4*vi9_6
   qr6 = qr69+ds8*vr4+ds11*vr5+ds5*vr6
   qi6 = qi69+ds8*vi4+ds11*vi5+ds5*vi6
   qr9 = qr69+ds9*vr7+ds12*vr8+ds6*vr9
   qi9 = qi69+ds9*vi7+ds12*vi8+ds6*vi9
   vrd_a = vrd-vra
   vid_a = vid-via
   vre_b = vre-vrb
   vie_b = vie-vib
   vrf_c = vrf-vrc
   vif_c = vif-vic
   qrad = ds13*vrd_a+ds16*vre_b+ds19*vrf_c
   qiad = ds13*vid_a+ds16*vie_b+ds19*vif_c
   qra = qrad+ds14*vra+ds17*vrb+ds20*vrc
   qia = qiad+ds14*via+ds17*vib+ds20*vic
   qrd = qrad+ds15*vrd+ds18*vre+ds21*vrf
   qid = qiad+ds15*vid+ds18*vie+ds21*vif
   qrbe = ds19*vrd_a+ds13*vre_b+ds16*vrf_c
   qibe = ds19*vid_a+ds13*vie_b+ds16*vif_c
   qrb = qrbe+ds20*vra+ds14*vrb+ds17*vrc
   qib = qibe+ds20*via+ds14*vib+ds17*vic
   qre = qrbe+ds21*vrd+ds15*vre+ds18*vrf
   qie = qibe+ds21*vid+ds15*vie+ds18*vif
   qrcf = ds16*vrd_a+ds19*vre_b+ds13*vrf_c
   qicf = ds16*vid_a+ds19*vie_b+ds13*vif_c
   qrc = qrcf+ds17*vra+ds20*vrb+ds14*vrc
   qic = qicf+ds17*via+ds20*vib+ds14*vic
   qrf = qrcf+ds18*vrd+ds21*vre+ds15*vrf
   qif = qicf+ds18*vid+ds21*vie+ds15*vif
   qr4a = qr4+qra
   qi4a = qi4+qia
   tr1 = qr1+c1*qr4a
   ti1 = qi1+c1*qi4a
   tr47ad = qr1+c2*qr4a
   ti47ad = qi1+c2*qi4a
   tr4d = tr47ad+qra
   ti4d = ti47ad+qia
   tr4_d = qr7+s1*qrd
   ti4_d = qi7+s1*qid
   tr4 = tr4d+tr4_d
   ti4 = ti4d+ti4_d
   trd = tr4d-tr4_d
   tid = ti4d-ti4_d
   tr7a = tr47ad+qr4
   ti7a = ti47ad+qi4
   tr7_a = s1*qr7-qrd
   ti7_a = s1*qi7-qid
   tr7 = tr7a+tr7_a
   ti7 = ti7a+ti7_a
   tra = tr7a-tr7_a
   tia = ti7a-ti7_a
   qr5b = qr5+qrb
   qi5b = qi5+qib
   tr2 = qr2+c1*qr5b
   ti2 = qi2+c1*qi5b
   tr58be = qr2+c2*qr5b
   ti58be = qi2+c2*qi5b
   tr5e = tr58be+qrb
   ti5e = ti58be+qib
   tr5_e = qr8+s1*qre
   ti5_e = qi8+s1*qie
   tr5 = tr5e+tr5_e
   ti5 = ti5e+ti5_e
   tre = tr5e-tr5_e
   tie = ti5e-ti5_e
   tr8b = tr58be+qr5
   ti8b = ti58be+qi5
   tr8_b = s1*qr8-qre
   ti8_b = s1*qi8-qie
   tr8 = tr8b+tr8_b
   ti8 = ti8b+ti8_b
   trb = tr8b-tr8_b
   tib = ti8b-ti8_b
   qr6c = qr6+qrc
   qi6c = qi6+qic
   tr3 = qr3+c1*qr6c
   ti3 = qi3+c1*qi6c
   tr69cf = qr3+c2*qr6c
   ti69cf = qi3+c2*qi6c
   tr6f = tr69cf+qrc
   ti6f = ti69cf+qic
   tr6_f = qr9+s1*qrf
   ti6_f = qi9+s1*qif
   tr6 = tr6f+tr6_f
   ti6 = ti6f+ti6_f
   trf = tr6f-tr6_f
   tif = ti6f-ti6_f
   tr9c = tr69cf+qr6
   ti9c = ti69cf+qi6
   tr9_c = s1*qr9-qrf
   ti9_c = s1*qi9-qif
   tr9 = tr9c+tr9_c
   ti9 = ti9c+ti9_c
   trc = tr9c-tr9_c
   tic = ti9c-ti9_c

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