/* csvorbis
 * Copyright (C) 2000 ymnk, JCraft,Inc.
 *  
 * Written by: 2000 ymnk<ymnk@jcraft.com>
 * Ported to C# from JOrbis by: Mark Crichton <crichton@gimp.org> 
 *   
 * Thanks go to the JOrbis team, for licencing the code under the
 * LGPL, making my job a lot easier.
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public License
 * as published by the Free Software Foundation; either version 2 of
 * the License, or (at your option) any later version.
   
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


using System;

namespace csvorbis 
{
	class Drft
	{
		int n;
		float[] trigcache;    
		int[] splitcache;

		internal void backward(float[] data)
		{
			if(n==1)return;
			drftb1(n,data,trigcache,trigcache,n,splitcache);
		}

		internal void init(int n)
		{
			this.n=n;
			trigcache=new float[3*n];
			splitcache=new int[32];
			fdrffti(n, trigcache, splitcache);
		}

		internal void clear()
		{
			if(trigcache!=null)trigcache=null;
			if(splitcache!=null)splitcache=null;
		}

		static int[] ntryh = { 4,2,3,5 };
		static float tpi = 6.28318530717958647692528676655900577f;
		static float hsqt2 = .70710678118654752440084436210485f;
		static float taui = .86602540378443864676372317075293618f;
		static float taur = -.5f;
		static float sqrt2 = 1.4142135623730950488016887242097f;

		static void drfti1(int n, float[] wa, int index, int[] ifac)
		{
			float arg,argh,argld,fi;
			int ntry=0,i,j=-1;
			int k1, l1, l2, ib;
			int ld, ii, ip, iis, nq, nr;
			int ido, ipm, nfm1;
			int nl=n;
			int nf=0;


			L101: j++;
			if (j < 4)
				ntry=ntryh[j];
			else
				ntry+=2;
							
			L104: nq=nl/ntry;
			nr=nl-ntry*nq;
			if(nr!=0) goto L101;
			
			nf++;
			ifac[nf+1]=ntry;
			nl=nq;
			if(ntry!=2) goto L107;
			if(nf==1) goto L107;
			
			for(i=1;i<nf;i++)
			{
				ib=nf-i+1;
				ifac[ib+1]=ifac[ib];
			}
			ifac[2] = 2;
			
			L107: if(nl!=1) goto L104;
			ifac[0]=n;
			ifac[1]=nf;
			argh=tpi/n;
			iis=0;
			nfm1=nf-1;
			l1=1;

			if(nfm1==0)return;

			for (k1=0;k1<nfm1;k1++)
			{
				ip=ifac[k1+2];
				ld=0;
				l2=l1*ip;
				ido=n/l2;
				ipm=ip-1;

				for (j=0;j<ipm;j++)
				{
					ld+=l1;
					i=iis;
					argld=(float)ld*argh;
					fi=0.0f;
					for (ii=2;ii<ido;ii+=2)
					{
						fi+=1.0f;
						arg=fi*argld;
						wa[index+i++]=(float)Math.Cos(arg);
						wa[index+i++]=(float)Math.Sin(arg);
					}
					iis+=ido;
				}
				l1=l2;
			}
		}

		static void fdrffti(int n, float[] wsave, int[] ifac)
		{
			if(n == 1) return;
			drfti1(n, wsave, n, ifac);
		}

		static void dradf2(int ido,int l1,float[] cc, float[] ch, float[] wa1, int index)
		{
			int i,k;
			float ti2,tr2;
			int t0,t1,t2,t3,t4,t5,t6;

			t1=0;
			t0=(t2=l1*ido);
			t3=ido<<1;
			for(k=0;k<l1;k++)
			{
				ch[t1<<1]=cc[t1]+cc[t2];
				ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
				t1+=ido;
				t2+=ido;
			}
    
			if(ido<2) return;

			if(ido==2) goto L105;
			
			t1=0;
			t2=t0;
			for(k=0;k<l1;k++)
			{
				t3=t2;
				t4=(t1<<1)+(ido<<1);
				t5=t1;
				t6=t1+t1;
				for(i=2;i<ido;i+=2)
				{
					t3+=2;
					t4-=2;
					t5+=2;
					t6+=2;
					tr2=wa1[index+i-2]*cc[t3-1]+wa1[index+i-1]*cc[t3];
					ti2=wa1[index+i-2]*cc[t3]-wa1[index+i-1]*cc[t3-1];
					ch[t6]=cc[t5]+ti2;
					ch[t4]=ti2-cc[t5];
					ch[t6-1]=cc[t5-1]+tr2;
					ch[t4-1]=cc[t5-1]-tr2;
				}
				t1+=ido;
				t2+=ido;
			}
			if(ido%2==1)return;
    
			L105: t3=(t2=(t1=ido)-1);
			t2+=t0;
			for(k=0;k<l1;k++)
			{
				ch[t1]=-cc[t2];
				ch[t1-1]=cc[t3];
				t1+=ido<<1;
				t2+=ido;
				t3+=ido;
			}
		}

		static void dradf4(int ido,int l1,float[] cc, float[] ch, 
			float[] wa1, int index1,
			float[] wa2, int index2,
			float[] wa3, int index3)
		{
			int i,k,t0,t1,t2,t3,t4,t5,t6;
			float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
			t0=l1*ido;
  
			t1=t0;
			t4=t1<<1;
			t2=t1+(t1<<1);
			t3=0;

			for(k=0;k<l1;k++)
			{
				tr1=cc[t1]+cc[t2];
				tr2=cc[t3]+cc[t4];

				ch[t5=t3<<2]=tr1+tr2;
				ch[(ido<<2)+t5-1]=tr2-tr1;
				ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
				ch[t5]=cc[t2]-cc[t1];
      
				t1+=ido;
				t2+=ido;
				t3+=ido;
				t4+=ido;
			}
			if(ido<2)return;

			if(ido==2) goto L105;
			
			t1=0;
			for(k=0;k<l1;k++)
			{
				t2=t1;
				t4=t1<<2;
				t5=(t6=ido<<1)+t4;
				for(i=2;i<ido;i+=2)
				{
					t3=(t2+=2);
					t4+=2;
					t5-=2;

					t3+=t0;
					cr2=wa1[index1+i-2]*cc[t3-1]+wa1[index1+i-1]*cc[t3];
					ci2=wa1[index1+i-2]*cc[t3]-wa1[index1+i-1]*cc[t3-1];
					t3+=t0;
					cr3=wa2[index2+i-2]*cc[t3-1]+wa2[index2+i-1]*cc[t3];
					ci3=wa2[index2+i-2]*cc[t3]-wa2[index2+i-1]*cc[t3-1];
					t3+=t0;
					cr4=wa3[index3+i-2]*cc[t3-1]+wa3[index3+i-1]*cc[t3];
					ci4=wa3[index3+i-2]*cc[t3]-wa3[index3+i-1]*cc[t3-1];

					tr1=cr2+cr4;
					tr4=cr4-cr2;
					ti1=ci2+ci4;
					ti4=ci2-ci4;

					ti2=cc[t2]+ci3;
					ti3=cc[t2]-ci3;
					tr2=cc[t2-1]+cr3;
					tr3=cc[t2-1]-cr3;
	  
					ch[t4-1]=tr1+tr2;
					ch[t4]=ti1+ti2;

					ch[t5-1]=tr3-ti4;
					ch[t5]=tr4-ti3;

					ch[t4+t6-1]=ti4+tr3;
					ch[t4+t6]=tr4+ti3;

					ch[t5+t6-1]=tr2-tr1;
					ch[t5+t6]=ti1-ti2;
				}
				t1+=ido;
			}
			if((ido&1)!=0)return;

			L105: t2=(t1=t0+ido-1)+(t0<<1);
			t3=ido<<2;
			t4=ido;
			t5=ido<<1;
			t6=ido;

			for(k=0;k<l1;k++)
			{
				ti1=-hsqt2*(cc[t1]+cc[t2]);
				tr1=hsqt2*(cc[t1]-cc[t2]);

				ch[t4-1]=tr1+cc[t6-1];
				ch[t4+t5-1]=cc[t6-1]-tr1;

				ch[t4]=ti1-cc[t1+t0];
				ch[t4+t5]=ti1+cc[t1+t0];
      
				t1+=ido;
				t2+=ido;
				t4+=t3;
				t6+=ido;
			}
		}

		static void dradfg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
			float[] c2, float[] ch, float[] ch2, float[] wa, int index)
		{
			int idij,ipph,i,j,k,l,ic,ik,iis;
			int t0,t1,t2=0,t3,t4,t5,t6,t7,t8,t9,t10;
			float dc2,ai1,ai2,ar1,ar2,ds2;
			int nbd;
			float dcp=0,arg,dsp=0,ar1h,ar2h;
			int idp2,ipp2;
  
			arg=tpi/(float)ip;
			dcp=(float)Math.Cos(arg);
			dsp=(float)Math.Sin(arg);
			ipph=(ip+1)>>1;
			ipp2=ip;
			idp2=ido;
			nbd=(ido-1)>>1;
			t0=l1*ido;
			t10=ip*ido;

			if(ido==1) goto L119;
			for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];

			t1=0;
			for(j=1;j<ip;j++)
			{
				t1+=t0;
				t2=t1;
				for(k=0;k<l1;k++)
				{
					ch[t2]=c1[t2];
					t2+=ido;
				}
			}

			iis=-ido;
			t1=0;
			if(nbd>l1)
			{
				for(j=1;j<ip;j++)
				{
					t1+=t0;
					iis+=ido;
					t2= -ido+t1;
					for(k=0;k<l1;k++)
					{
						idij=iis-1;
						t2+=ido;
						t3=t2;
						for(i=2;i<ido;i+=2)
						{
							idij+=2;
							t3+=2;
							ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
							ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
						}
					}
				}
			}
			else
			{

				for(j=1;j<ip;j++)
				{
					iis+=ido;
					idij=iis-1;
					t1+=t0;
					t2=t1;
					for(i=2;i<ido;i+=2)
					{
						idij+=2;
						t2+=2;
						t3=t2;
						for(k=0;k<l1;k++)
						{
							ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
							ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
							t3+=ido;
						}
					}
				}
			}

			t1=0;
			t2=ipp2*t0;
			if(nbd<l1)
			{
				for(j=1;j<ipph;j++)
				{
					t1+=t0;
					t2-=t0;
					t3=t1;
					t4=t2;
					for(i=2;i<ido;i+=2)
					{
						t3+=2;
						t4+=2;
						t5=t3-ido;
						t6=t4-ido;
						for(k=0;k<l1;k++)
						{
							t5+=ido;
							t6+=ido;
							c1[t5-1]=ch[t5-1]+ch[t6-1];
							c1[t6-1]=ch[t5]-ch[t6];
							c1[t5]=ch[t5]+ch[t6];
							c1[t6]=ch[t6-1]-ch[t5-1];
						}
					}
				}
			}
			else
			{
				for(j=1;j<ipph;j++)
				{
					t1+=t0;
					t2-=t0;
					t3=t1;
					t4=t2;
					for(k=0;k<l1;k++)
					{
						t5=t3;
						t6=t4;
						for(i=2;i<ido;i+=2)
						{
							t5+=2;
							t6+=2;
							c1[t5-1]=ch[t5-1]+ch[t6-1];
							c1[t6-1]=ch[t5]-ch[t6];
							c1[t5]=ch[t5]+ch[t6];
							c1[t6]=ch[t6-1]-ch[t5-1];
						}
						t3+=ido;
						t4+=ido;
					}
				}
			}

			L119:		for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];

			t1=0;
			t2=ipp2*idl1;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1-ido;
				t4=t2-ido;
				for(k=0;k<l1;k++)
				{
					t3+=ido;
					t4+=ido;
					c1[t3]=ch[t3]+ch[t4];
					c1[t4]=ch[t4]-ch[t3];
				}
			}

			ar1=1.0f;
			ai1=0.0f;
			t1=0;
			t2=ipp2*idl1;
			t3=(ip-1)*idl1;
			for(l=1;l<ipph;l++)
			{
				t1+=idl1;
				t2-=idl1;
				ar1h=dcp*ar1-dsp*ai1;
				ai1=dcp*ai1+dsp*ar1;
				ar1=ar1h;
				t4=t1;
				t5=t2;
				t6=t3;
				t7=idl1;

				for(ik=0;ik<idl1;ik++)
				{
					ch2[t4++]=c2[ik]+ar1*c2[t7++];
					ch2[t5++]=ai1*c2[t6++];
				}

				dc2=ar1;
				ds2=ai1;
				ar2=ar1;
				ai2=ai1;

				t4=idl1;
				t5=(ipp2-1)*idl1;
				for(j=2;j<ipph;j++)
				{
					t4+=idl1;
					t5-=idl1;

					ar2h=dc2*ar2-ds2*ai2;
					ai2=dc2*ai2+ds2*ar2;
					ar2=ar2h;

					t6=t1;
					t7=t2;
					t8=t4;
					t9=t5;
					for(ik=0;ik<idl1;ik++)
					{
						ch2[t6]+=ar2*c2[t8++];
						t6++;
						ch2[t7]+=ai2*c2[t9++];
						t7++;
					}
				}
			}
			t1=0;
			for(j=1;j<ipph;j++)
			{
				t1+=idl1;
				t2=t1;
				for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
			}

			if(ido<l1) goto L132;
			
			t1=0;
			t2=0;
			for(k=0;k<l1;k++)
			{
				t3=t1;
				t4=t2;
				for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
				t1+=ido;
				t2+=t10;
			}
			
			goto L135;

			L132:			for(i=0;i<ido;i++)
							{
								t1=i;
								t2=i;
								for(k=0;k<l1;k++)
								{
									cc[t2]=ch[t1];
									t1+=ido;
									t2+=t10;
								}
							}
								
			L135:	t1=0;
			t2=ido<<1;
			t3=0;
			t4=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t2;
				t3+=t0;
				t4-=t0;

				t5=t1;
				t6=t3;
				t7=t4;

				for(k=0;k<l1;k++)
				{
					cc[t5-1]=ch[t6];
					cc[t5]=ch[t7];
					t5+=t10;
					t6+=ido;
					t7+=ido;
				}
			}

			if(ido==1)return;
			if(nbd<l1) goto L141;
								
			t1=-ido;
			t3=0;
			t4=0;
			t5=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t2;
				t3+=t2;
				t4+=t0;
				t5-=t0;
				t6=t1;
				t7=t3;
				t8=t4;
				t9=t5;
				for(k=0;k<l1;k++)
				{
					for(i=2;i<ido;i+=2)
					{
						ic=idp2-i;
						cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
						cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
						cc[i+t7]=ch[i+t8]+ch[i+t9];
						cc[ic+t6]=ch[i+t9]-ch[i+t8];
					}
					t6+=t10;
					t7+=t10;
					t8+=ido;
					t9+=ido;
				}
			}
			return;
								
			L141:	t1=-ido;
			t3=0;
			t4=0;
			t5=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t2;
				t3+=t2;
				t4+=t0;
				t5-=t0;
				for(i=2;i<ido;i+=2)
				{
					t6=idp2+t1-i;
					t7=i+t3;
					t8=i+t4;
					t9=i+t5;
					for(k=0;k<l1;k++)
					{
						cc[t7-1]=ch[t8-1]+ch[t9-1];
						cc[t6-1]=ch[t8-1]-ch[t9-1];
						cc[t7]=ch[t8]+ch[t9];
						cc[t6]=ch[t9]-ch[t8];
						t6+=t10;
						t7+=t10;
						t8+=ido;
						t9+=ido;
					}
				}
			}
		}

		static void drftf1(int n,float[] c, float[] ch, float[] wa, int[] ifac)
		{
			int i,k1,l1,l2;
			int na,kh,nf;
			int ip,iw,ido,idl1,ix2,ix3;

			nf=ifac[1];
			na=1;
			l2=n;
			iw=n;

			for(k1=0;k1<nf;k1++)
			{
				kh=nf-k1;
				ip=ifac[kh+1];
				l1=l2/ip;
				ido=n/l2;
				idl1=ido*l1;
				iw-=(ip-1)*ido;
				na=1-na;

				if(ip!=4) goto L102;

				ix2=iw+ido;
				ix3=ix2+ido;
				if(na!=0)
					dradf4(ido,l1,ch,c,wa,iw-1,wa,ix2-1,wa,ix3-1);
				else
					dradf4(ido,l1,c,ch,wa,iw-1,wa,ix2-1,wa,ix3-1);
				goto L110;

			L102: if(ip!=2) goto L104;
				if(na!=0) goto L103;
				
				dradf2(ido,l1,c,ch,wa, iw-1);
				goto L110;

			L103: dradf2(ido,l1,ch,c,wa, iw-1);
				goto L110;

			L104: if(ido==1)na=1-na;
				if(na!=0) goto L109;
				dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,iw-1);
				na=1;
				goto L110;

			L109: dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,iw-1);
				na=0;
				l2=l1;
				break;

			L110: l2=l1;
			}

			if(na==1)return;
			for(i=0;i<n;i++)c[i]=ch[i];
		}

		static void dradb2(int ido,int l1,float[] cc,float[] ch,float[] wa1, int index)
		{
			int i,k,t0,t1,t2,t3,t4,t5,t6;
			float ti2,tr2;

			t0=l1*ido;
    
			t1=0;
			t2=0;
			t3=(ido<<1)-1;
			for(k=0;k<l1;k++)
			{
				ch[t1]=cc[t2]+cc[t3+t2];
				ch[t1+t0]=cc[t2]-cc[t3+t2];
				t2=(t1+=ido)<<1;
			}

			if(ido<2)return;
			if(ido==2) goto L105;
			
			t1=0;
			t2=0;
			for(k=0;k<l1;k++)
			{
				t3=t1;
				t5=(t4=t2)+(ido<<1);
				t6=t0+t1;
				for(i=2;i<ido;i+=2)
				{
					t3+=2;
					t4+=2;
					t5-=2;
					t6+=2;
					ch[t3-1]=cc[t4-1]+cc[t5-1];
					tr2=cc[t4-1]-cc[t5-1];
					ch[t3]=cc[t4]-cc[t5];
					ti2=cc[t4]+cc[t5];
					ch[t6-1]=wa1[index+i-2]*tr2-wa1[index+i-1]*ti2;
					ch[t6]=wa1[index+i-2]*ti2+wa1[index+i-1]*tr2;
				}
				t2=(t1+=ido)<<1;
			}
			if((ido%2)==1)return;
			

			L105: t1=ido-1;
			t2=ido-1;
			for(k=0;k<l1;k++)
			{
				ch[t1]=cc[t2]+cc[t2];
				ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
				t1+=ido;
				t2+=ido<<1;
			}
		}

		static void dradb3(int ido,int l1,float[] cc,float[] ch,
			float[] wa1, int index1,
			float[] wa2, int index2)
		{
			int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
			float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
			t0=l1*ido;

			t1=0;
			t2=t0<<1;
			t3=ido<<1;
			t4=ido+(ido<<1);
			t5=0;
			for(k=0;k<l1;k++)
			{
				tr2=cc[t3-1]+cc[t3-1];
				cr2=cc[t5]+(taur*tr2);
				ch[t1]=cc[t5]+tr2;
				ci3=taui*(cc[t3]+cc[t3]);
				ch[t1+t0]=cr2-ci3;
				ch[t1+t2]=cr2+ci3;
				t1+=ido;
				t3+=t4;
				t5+=t4;
			}

			if(ido==1)return;

			t1=0;
			t3=ido<<1;
			for(k=0;k<l1;k++)
			{
				t7=t1+(t1<<1);
				t6=(t5=t7+t3);
				t8=t1;
				t10=(t9=t1+t0)+t0;

				for(i=2;i<ido;i+=2)
				{
					t5+=2;
					t6-=2;
					t7+=2;
					t8+=2;
					t9+=2;
					t10+=2;
					tr2=cc[t5-1]+cc[t6-1];
					cr2=cc[t7-1]+(taur*tr2);
					ch[t8-1]=cc[t7-1]+tr2;
					ti2=cc[t5]-cc[t6];
					ci2=cc[t7]+(taur*ti2);
					ch[t8]=cc[t7]+ti2;
					cr3=taui*(cc[t5-1]-cc[t6-1]);
					ci3=taui*(cc[t5]+cc[t6]);
					dr2=cr2-ci3;
					dr3=cr2+ci3;
					di2=ci2+cr3;
					di3=ci2-cr3;
					ch[t9-1]=wa1[index1+i-2]*dr2-wa1[index1+i-1]*di2;
					ch[t9]=wa1[index1+i-2]*di2+wa1[index1+i-1]*dr2;
					ch[t10-1]=wa2[index2+i-2]*dr3-wa2[index2+i-1]*di3;
					ch[t10]=wa2[index2+i-2]*di3+wa2[index2+i-1]*dr3;
				}
				t1+=ido;
			}
		}

		static void dradb4(int ido,int l1,float[] cc,float[] ch,
			float[] wa1, int index1,
			float[] wa2, int index2,
			float[] wa3, int index3)
		{
			int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
			float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
			t0=l1*ido;
  
			t1=0;
			t2=ido<<2;
			t3=0;
			t6=ido<<1;
			for(k=0;k<l1;k++)
			{
				t4=t3+t6;
				t5=t1;
				tr3=cc[t4-1]+cc[t4-1];
				tr4=cc[t4]+cc[t4]; 
				tr1=cc[t3]-cc[(t4+=t6)-1];
				tr2=cc[t3]+cc[t4-1];
				ch[t5]=tr2+tr3;
				ch[t5+=t0]=tr1-tr4;
				ch[t5+=t0]=tr2-tr3;
				ch[t5+=t0]=tr1+tr4;
				t1+=ido;
				t3+=t2;
			}

			if(ido<2)return;
			if(ido==2) goto L105;
			
			t1=0;
			for(k=0;k<l1;k++)
			{
				t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
				t7=t1;
				for(i=2;i<ido;i+=2)
				{
					t2+=2;
					t3+=2;
					t4-=2;
					t5-=2;
					t7+=2;
					ti1=cc[t2]+cc[t5];
					ti2=cc[t2]-cc[t5];
					ti3=cc[t3]-cc[t4];
					tr4=cc[t3]+cc[t4];
					tr1=cc[t2-1]-cc[t5-1];
					tr2=cc[t2-1]+cc[t5-1];
					ti4=cc[t3-1]-cc[t4-1];
					tr3=cc[t3-1]+cc[t4-1];
					ch[t7-1]=tr2+tr3;
					cr3=tr2-tr3;
					ch[t7]=ti2+ti3;
					ci3=ti2-ti3;
					cr2=tr1-tr4;
					cr4=tr1+tr4;
					ci2=ti1+ti4;
					ci4=ti1-ti4;

					ch[(t8=t7+t0)-1]=wa1[index1+i-2]*cr2-wa1[index1+i-1]*ci2;
					ch[t8]=wa1[index1+i-2]*ci2+wa1[index1+i-1]*cr2;
					ch[(t8+=t0)-1]=wa2[index2+i-2]*cr3-wa2[index2+i-1]*ci3;
					ch[t8]=wa2[index2+i-2]*ci3+wa2[index2+i-1]*cr3;
					ch[(t8+=t0)-1]=wa3[index3+i-2]*cr4-wa3[index3+i-1]*ci4;
					ch[t8]=wa3[index3+i-2]*ci4+wa3[index3+i-1]*cr4;
				}
				t1+=ido;
			}
			if(ido%2 == 1)return;
			

			L105: t1=ido;
			t2=ido<<2;
			t3=ido-1;
			t4=ido+(ido<<1);
			for(k=0;k<l1;k++)
			{
				t5=t3;
				ti1=cc[t1]+cc[t4];
				ti2=cc[t4]-cc[t1];
				tr1=cc[t1-1]-cc[t4-1];
				tr2=cc[t1-1]+cc[t4-1];
				ch[t5]=tr2+tr2;
				ch[t5+=t0]=sqrt2*(tr1-ti1);
				ch[t5+=t0]=ti2+ti2;
				ch[t5+=t0]=-sqrt2*(tr1+ti1);
      
				t3+=ido;
				t1+=t2;
				t4+=t2;
			}
		}

		static void dradbg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
			float[] c2,float[] ch,float[] ch2,float[] wa, int index )
		{

			int idij,ipph=0,i,j,k,l,ik,iis,t0=0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10=0,
				t11,t12;
			float dc2,ai1,ai2,ar1,ar2,ds2;
			int nbd=0;
			float dcp=0,arg,dsp=0,ar1h,ar2h;
			int ipp2=0;

			t10=ip*ido;
			t0=l1*ido;
			arg=tpi/(float)ip;
			dcp=(float)Math.Cos(arg);
			dsp=(float)Math.Sin(arg);
			nbd=(int)((uint)(ido-1)>>1);
			ipp2=ip;
			ipph=(int)((uint)(ip+1)>>1);
			if(ido<l1) goto L103;
										
			t1=0;
			t2=0;
			for(k=0;k<l1;k++)
			{
				t3=t1;
				t4=t2;
				for(i=0;i<ido;i++)
				{
					ch[t3]=cc[t4];
					t3++;
					t4++;
				}
				t1+=ido;
				t2+=t10;
			}
			goto L106;

			L103: t1=0;
			for(i=0;i<ido;i++)
			{
				t2=t1;
				t3=t1;
				for(k=0;k<l1;k++)
				{
					ch[t2]=cc[t3];
					t2+=ido;
					t3+=t10;
				}
				t1++;
			}
			L106: t1=0;
			t2=ipp2*t0;
			t7=(t5=ido<<1);
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;
				t6=t5;
				for(k=0;k<l1;k++)
				{
					ch[t3]=cc[t6-1]+cc[t6-1];
					ch[t4]=cc[t6]+cc[t6];
					t3+=ido;
					t4+=ido;
					t6+=t10;
				}
				t5+=t7;
			}
			if (ido == 1) goto L116;
			if(nbd<l1) goto L112;
										
			t1=0;
			t2=ipp2*t0;
			t7=0;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;

				t7+=(ido<<1);
				t8=t7;
				for(k=0;k<l1;k++)
				{
					t5=t3;
					t6=t4;
					t9=t8;
					t11=t8;
					for(i=2;i<ido;i+=2)
					{
						t5+=2;
						t6+=2;
						t9+=2;
						t11-=2;
						ch[t5-1]=cc[t9-1]+cc[t11-1];
						ch[t6-1]=cc[t9-1]-cc[t11-1];
						ch[t5]=cc[t9]-cc[t11];
						ch[t6]=cc[t9]+cc[t11];
					}
					t3+=ido;
					t4+=ido;
					t8+=t10;
				}
			}
			goto L116;
							
			L112: t1=0;
			t2=ipp2*t0;
			t7=0;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;
				t7+=(ido<<1);
				t8=t7;
				t9=t7;
				for(i=2;i<ido;i+=2)
				{
					t3+=2;
					t4+=2;
					t8+=2;
					t9-=2;
					t5=t3;
					t6=t4;
					t11=t8;
					t12=t9;
					for(k=0;k<l1;k++)
					{
						ch[t5-1]=cc[t11-1]+cc[t12-1];
						ch[t6-1]=cc[t11-1]-cc[t12-1];
						ch[t5]=cc[t11]-cc[t12];
						ch[t6]=cc[t11]+cc[t12];
						t5+=ido;
						t6+=ido;
						t11+=t10;
						t12+=t10;
					}
				}
			}
			L116: ar1=1.0f;
			ai1=0.0f;
			t1=0;
			t9=(t2=ipp2*idl1);
			t3=(ip-1)*idl1;
			for(l=1;l<ipph;l++)
			{
				t1+=idl1;
				t2-=idl1;
	  
				ar1h=dcp*ar1-dsp*ai1;
				ai1=dcp*ai1+dsp*ar1;
				ar1=ar1h;
				t4=t1;
				t5=t2;
				t6=0;
				t7=idl1;
				t8=t3;
				for(ik=0;ik<idl1;ik++)
				{
					c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
					c2[t5++]=ai1*ch2[t8++];
				}
				dc2=ar1;
				ds2=ai1;
				ar2=ar1;
				ai2=ai1;

				t6=idl1;
				t7=t9-idl1;
				for(j=2;j<ipph;j++)
				{
					t6+=idl1;
					t7-=idl1;
					ar2h=dc2*ar2-ds2*ai2;
					ai2=dc2*ai2+ds2*ar2;
					ar2=ar2h;
					t4=t1;
					t5=t2;
					t11=t6;
					t12=t7;
					for(ik=0;ik<idl1;ik++)
					{
						c2[t4]+=ar2*ch2[t11++];
						t4++;
						c2[t5]+=ai2*ch2[t12++];
						t5++;
					}
				}
			}

			t1=0;
			for(j=1;j<ipph;j++)
			{
				t1+=idl1;
				t2=t1;
				for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
			}

			t1=0;
			t2=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;
				for(k=0;k<l1;k++)
				{
					ch[t3]=c1[t3]-c1[t4];
					ch[t4]=c1[t3]+c1[t4];
					t3+=ido;
					t4+=ido;
				}
			}

			if(ido==1) goto L132;
			if(nbd<l1) goto L128;
										
			t1=0;
			t2=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;
				for(k=0;k<l1;k++)
				{
					t5=t3;
					t6=t4;
					for(i=2;i<ido;i+=2)
					{
						t5+=2;
						t6+=2;
						ch[t5-1]=c1[t5-1]-c1[t6];
						ch[t6-1]=c1[t5-1]+c1[t6];
						ch[t5]=c1[t5]+c1[t6-1];
						ch[t6]=c1[t5]-c1[t6-1];
					}
					t3+=ido;
					t4+=ido;
				}
			}
			goto L132;
			L128: t1=0;
			t2=ipp2*t0;
			for(j=1;j<ipph;j++)
			{
				t1+=t0;
				t2-=t0;
				t3=t1;
				t4=t2;
				for(i=2;i<ido;i+=2)
				{
					t3+=2;
					t4+=2;
					t5=t3;
					t6=t4;
					for(k=0;k<l1;k++)
					{
						ch[t5-1]=c1[t5-1]-c1[t6];
						ch[t6-1]=c1[t5-1]+c1[t6];
						ch[t5]=c1[t5]+c1[t6-1];
						ch[t6]=c1[t5]-c1[t6-1];
						t5+=ido;
						t6+=ido;
					}
				}
			}
				
			L132: if(ido==1)return;

			for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];

			t1=0;
			for(j=1;j<ip;j++)
			{
				t2=(t1+=t0);
				for(k=0;k<l1;k++)
				{
					c1[t2]=ch[t2];
					t2+=ido;
				}
			}

			if(nbd>l1) goto L139;
										
			iis= -ido-1;
			t1=0;
			for(j=1;j<ip;j++)
			{
				iis+=ido;
				t1+=t0;
				idij=iis;
				t2=t1;
				for(i=2;i<ido;i+=2)
				{
					t2+=2;
					idij+=2;
					t3=t2;
					for(k=0;k<l1;k++)
					{
						c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
						c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
						t3+=ido;
					}
				}
			}
			return;

			L139: iis= -ido-1;
			t1=0;
			for(j=1;j<ip;j++)
			{
				iis+=ido;
				t1+=t0;
				t2=t1;
				for(k=0;k<l1;k++)
				{
					idij=iis;
					t3=t2;
					for(i=2;i<ido;i+=2)
					{
						idij+=2;
						t3+=2;
						c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
						c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
					}
					t2+=ido;
				}
			}
		}

		static void drftb1(int n, float[] c, float[] ch, float[] wa, int index, int[] ifac)
		{
			int i,k1,l1,l2=0;
			int na;
			int nf,ip=0,iw,ix2,ix3,ido=0,idl1=0;

			nf=ifac[1];
			na=0;
			l1=1;
			iw=1;

			for(k1=0;k1<nf;k1++)
			{
				ip=ifac[k1 + 2];
				l2=ip*l1;
				ido=n/l2;
				idl1=ido*l1;
				if(ip!=4) goto L103;
				ix2=iw+ido;
				ix3=ix2+ido;

				if(na!=0)
					dradb4(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
				else
					dradb4(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
				na=1-na;
				goto L115;							
			L103: if(ip!=2) goto L106;
											
				if(na!=0)
					dradb2(ido,l1,ch,c,wa,index+iw-1);
				else
					dradb2(ido,l1,c,ch,wa,index+iw-1);
				na=1-na;
				goto L115;

			L106: if(ip!=3) goto L109;
											
				ix2=iw+ido;
				if(na!=0)
					dradb3(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1);
				else
					dradb3(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1);
				na=1-na;
				goto L115;
										
			L109: if(na!=0)
					  dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,index+iw-1);
				  else
					  dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,index+iw-1);
				if(ido==1)na=1-na;
				l1=l2;
				iw+=(ip-1)*ido;

			L115: l1=l2;
				iw+=(ip-1)*ido;
			}
			if(na==0)return;
			for(i=0;i<n;i++)c[i]=ch[i];
		}
	}
}