Anonymous

/********************************************************************************                                                                              **   (C) 1997-2004 by Ernst W. Mayer.                                           **                                                                              **  This program is free software; you can redistribute it and/or modify it     **  under the terms of the GNU 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 General Public License for   **  more details.                                                               **                                                                              **  You should have received a copy of the GNU General Public License along     **  with this program; see the file GPL.txt.  If not, you may view one at       **  http://www.fsf.org/licenses/licenses.html, or obtain one by writing to the  **  Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA     **  02111-1307, USA.                                                            **                                                                              ********************************************************************************/#include "Mlucas.h"/***************/void radix16_dif_pass1(double a[], int n){/*!...Acronym: DIF = Decimation In Frequency!!...Subroutine to perform an initial radix-16 complex DIF FFT pass on the data in the length-N real vector A.!!   The data are stored in a 1-D zero-offset array, with 2^PAD_BITS 8-byte padding elements inserted!   between every block of 2^DAT_BITS contiguous data. The array padding is to prevent data being accessed!   in strides that are large powers of two and thus to minimize cache thrashing!   (in cache-based microprocessor architectures) or bank conflicts (in supercomputers.)!!   See the documentation in radix16_dif_pass for further details on storage and indexing.*/	int j,j1;	static int n16,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15, first_entry=TRUE;	static double c = 0.9238795325112867561281831, s = 0.3826834323650897717284599;	/* exp[i*(twopi/16)]	*/	double rt,it	,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32;	if(!first_entry && (n >> 4) != n16)	/* New runlength?	*/	{	  first_entry=TRUE;	}/*...initialize things upon first entry	*/	if(first_entry)	{	  first_entry=FALSE;	  n16=n/16;	  p1 = n16 + ( (n16 >> DAT_BITS) << PAD_BITS );	  p2 = p1 +p1;	  p3 = p2 +p1;	  p4 = p3 +p1;	  p5 = p4 +p1;	  p6 = p5 +p1;	  p7 = p6 +p1;	  p8 = p7 +p1;	  p9 = p8 +p1;	  p10= p9 +p1;	  p11= p10+p1;	  p12= p11+p1;	  p13= p12+p1;	  p14= p13+p1;	  p15= p14+p1;	}/*...The radix-16 pass is here.	*/      for(j=0; j < n16; j += 2)      {	  j1 = j  + ( (j  >> DAT_BITS) << PAD_BITS );/*       gather the needed data (16 64-bit complex, i.e. 16 64-bit reals) and do the first set of four length-4 transforms.	*//*...Block 1:	*/	    t1 =a[j1    ];	t2 =a[j1+1    ];	    rt =a[j1+p8 ];	it =a[j1+1+p8 ];	    t3 =t1 -rt;  	t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p4 ];	t6 =a[j1+1+p4 ];	    rt =a[j1+p12];	it =a[j1+1+p12];	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 +t8;	t3 =t3 -t8;			t8 =t4 -rt;	t4 =t4 +rt;/*...Block 2:	*/	    t9 =a[j1+p2 ];	t10=a[j1+1+p2 ];	    rt =a[j1+p10];	it =a[j1+1+p10];	    t11=t9 -rt;  	t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ];	t14=a[j1+1+p6 ];	    rt =a[j1+p14];	it =a[j1+1+p14];	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11+t16;	t11=t11-t16;			t16=t12-rt;	t12=t12+rt;/*...Block 3:	*/	    t17=a[j1+p1 ];	t18=a[j1+1+p1 ];	    rt =a[j1+p9 ];	it =a[j1+1+p9 ];	    t19=t17-rt;  	t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p5 ];	t22=a[j1+1+p5 ];	    rt =a[j1+p13];	it =a[j1+1+p13];	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19+t24;	t19=t19-t24;			t24=t20-rt;	t20=t20+rt;/*...Block 4:	*/	    t25=a[j1+p3 ];	t26=a[j1+1+p3 ];	    rt =a[j1+p11];	it =a[j1+1+p11];	    t27=t25-rt;  	t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p7 ];	t30=a[j1+1+p7 ];	    rt =a[j1+p15];	it =a[j1+1+p15];	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27+t32;	t27=t27-t32;			t32=t28-rt;	t28=t28+rt;/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(i* 1*twopi/16) =       ( c, s), exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 3*twopi/16) =       ( s, c) (for inputs to transform block 2)  	1, exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 4*twopi/16) =       ( 0, 1), exp(i* 6*twopi/16) = ISRT2*(-1, 1) (for inputs to transform block 3)  	1, exp(i* 3*twopi/16) =       ( s, c), exp(i* 6*twopi/16) = ISRT2*(-1, 1), exp(i* 9*twopi/16) =       (-c,-s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;	    a[j1    ]=t1+t17;	a[j1+1    ]=t2+t18;	    a[j1+p1 ]=t1-t17;	a[j1+1+p1 ]=t2-t18;	    a[j1+p2 ]=t9 -t26;	a[j1+1+p2 ]=t10+t25;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p3 ]=t9 +t26;	a[j1+1+p3 ]=t10-t25;/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 +t14;	t5 =t5 -t14;		/* twiddle mpy by E^4 = I	*/			t14=t6 -rt;	t6 =t6 +rt;	    rt =(t21-t22)*ISRT2;t22=(t21+t22)*ISRT2;	t21=rt;	/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t21,22,29,30	*/	    rt =(t30+t29)*ISRT2;it =(t30-t29)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;	    a[j1+p4 ]=t5+t21;	a[j1+1+p4 ]=t6+t22;	    a[j1+p5 ]=t5-t21;	a[j1+1+p5 ]=t6-t22;	    a[j1+p6 ]=t13-t30;	a[j1+1+p6 ]=t14+t29;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p7 ]=t13+t30;	a[j1+1+p7 ]=t14-t29;/*...Block 2: t3,11,19,27	*/	    rt =(t11-t12)*ISRT2;it =(t11+t12)*ISRT2;		/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t11,12	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c - t20*s;	t20=t20*c + t19*s;	t19=rt;	/* twiddle mpy by E^1	*/	    rt =t27*s - t28*c;	it =t28*s + t27*c;		/* twiddle mpy by E^3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;	    a[j1+p8 ]=t3+t19;	a[j1+1+p8 ]=t4+t20;	    a[j1+p9 ]=t3-t19;	a[j1+1+p9 ]=t4-t20;	    a[j1+p10]=t11-t28;	a[j1+1+p10]=t12+t27;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p11]=t11+t28;	a[j1+1+p11]=t12-t27;/*...Block 4: t7,15,23,31	*/	    rt =(t16+t15)*ISRT2;it =(t16-t15)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s - t24*c;	t24=t24*s + t23*c;	t23=rt;	/* twiddle mpy by E^3	*/	    rt =t31*c - t32*s;	it =t32*c + t31*s;		/* twiddle mpy by E^1 = -E^9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;/* Note: t23+rt = t23*(s+1)	*/	    a[j1+p12]=t7+t23;	a[j1+1+p12]=t8+t24;	    a[j1+p13]=t7-t23;	a[j1+1+p13]=t8-t24;	    a[j1+p14]=t15-t32;	a[j1+1+p14]=t16+t31;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p15]=t15+t32;	a[j1+1+p15]=t16-t31;	}}/***************/int radix16_ditN_cy_dif1(double a[], int n, int nwt, int nwt_bits, double wt0[], double wt1[], int si[], int nrt, int nrt_bits, struct complex rt0[],                         struct complex rt1[], double base[], double baseinv[], int iter, double *fracmax, int p){/*!...Acronym: DWT = Discrete Weighted Transform, DIT = Decimation In Time, DIF = Decimation In Frequency!!...Performs a final radix-16 complex DIT pass, an inverse DWT weighting, a carry propagation,!   a forward DWT weighting, and an initial radix-16 complex DIF pass on the data in the length-N real vector A.!!   Data enter and are returned in the A-array.!!   See the documentation in mers_mod_square and radix16_dif_pass for further details on the array!   storage scheme, and radix16_ditN_cy_dif1 for details on the reduced-length weights array scheme.*/	int n16,ithread,i,j,j1,iroot,root_incr,k1,k2,k,khi,l,outer;	static int *_i, *bjmodn0 = 0x0, *bjmodn1 = 0x0, *bjmodn2 = 0x0, *bjmodn3 = 0x0, *bjmodn4 = 0x0, *bjmodn5 = 0x0, *bjmodn6 = 0x0, *bjmodn7 = 0x0, *bjmodn8 = 0x0, *bjmodn9 = 0x0, *bjmodnA = 0x0, *bjmodnB = 0x0, *bjmodnC = 0x0, *bjmodnD = 0x0, *bjmodnE = 0x0, *bjmodnF = 0x0	, *bjmodnini = 0x0, *jstart = 0x0, *jhi = 0x0, *_col = 0x0, *_co2 = 0x0, *_co3 = 0x0;	static int psave=0,bw,sw,nm1,nrtm1,pini,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15;	const  double one_half[3] = {1.0, 0.5, 0.25};	/* Needed for small-weights-tables scheme */	static double c = 0.9238795325112867561281831, s = 0.3826834323650897717284599, radix_inv,n2inv	, *cy0 = 0x0, *cy1 = 0x0, *cy2 = 0x0, *cy3 = 0x0, *cy4 = 0x0, *cy5 = 0x0, *cy6 = 0x0, *cy7 = 0x0, *cy8 = 0x0, *cy9 = 0x0, *cyA = 0x0, *cyB = 0x0, *cyC = 0x0, *cyD = 0x0, *cyE = 0x0, *cyF = 0x0	, *_maxerr = 0x0;	double rt,it#if !USE_NEW_DIT		,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15		,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15#endif	,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32	,a1p0r,a1p1r,a1p2r,a1p3r,a1p4r,a1p5r,a1p6r,a1p7r,a1p8r,a1p9r,a1pAr,a1pBr,a1pCr,a1pDr,a1pEr,a1pFr	,a1p0i,a1p1i,a1p2i,a1p3i,a1p4i,a1p5i,a1p6i,a1p7i,a1p8i,a1p9i,a1pAi,a1pBi,a1pCi,a1pDi,a1pEi,a1pFi	,temp,frac,scale;	double maxerr = 0.0;#if PFETCH	double *add0, *add;#endif	int err;	FILE *f1,*f2;	static int first_entry=TRUE;/*...stuff for the reduced-length DWT weights array is here:	*/	int n_div_nwt;	int col,co2,co3,m,m2,n_minus_sil,n_minus_silp1,sinwt,sinwtm1;	double wt,wtinv,wtl,wtlp1,wtn,wtnm1,wtA,wtB,wtC;#ifdef CY_DBG	int ierr, npad;	static double*b = 0x0;#endif/*...change n16 and n_div_wt to non-static to work around a gcc compiler bug. */	n16   = n/16;	n_div_nwt = n16 >> nwt_bits;	if(NTHREADS > 1)	{		assert(n16      %NTHREADS == 0);		assert(n_div_nwt%NTHREADS == 0);	}	if((n_div_nwt << nwt_bits) != n16)	{	  if(INTERACT)printf("FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  f1 = fopen(   OFILE,"a");	  f2 = fopen(STATFILE,"a");	  fprintf(f1,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fprintf(f2,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fclose(f1);	  fclose(f2);	  err=1;	  return(err);	}	if(p != psave)	{	  first_entry=TRUE;	}/*...initialize things upon first entry	*/	if(first_entry)	{#ifdef CY_DBG	npad = n + ( (n >> DAT_BITS) << PAD_BITS );	/* length of padded data array.	*/	printf("CY_DBG: allocating scratch B-array of %d doubles...\n", npad);	b = (double *)malloc(npad*sizeof(double));	if(!b){ printf("FATAL: unable to allocate array B.\n"); exit(EXIT_FAILURE); }#endif	  psave = p;	  first_entry=FALSE;	  radix_inv = (double)((long double)1.0L/(long double)16.0L);	  n2inv     = (double)((long double)1.0L/(long double)(n/2));	  bw    = p%n;		/* Number of bigwords in the Crandall/Fagin mixed-radix representation = (Mersenne exponent) mod (vector length).	*/	  sw    = n - bw;	/* Number of smallwords.	*/	  nm1   = n-1;	  nrtm1 = nrt - 1;/*   constant index offsets for load/stores are here.	*/	  pini = n16/NTHREADS;	  pini += ( (pini >> DAT_BITS) << PAD_BITS );	  p1 = n16 + ( (n16 >> DAT_BITS) << PAD_BITS );	  p2 = p1 +p1;	  p3 = p2 +p1;	  p4 = p3 +p1;	  p5 = p4 +p1;	  p6 = p5 +p1;	  p7 = p6 +p1;	  p8 = p7 +p1;	  p9 = p8 +p1;	  p10= p9 +p1;	  p11= p10+p1;	  p12= p11+p1;	  p13= p12+p1;	  p14= p13+p1;	  p15= p14+p1;	  if(cy0)	/* If it's a new exponent of a range test, need to deallocate these. */	  {	    free((void *)_i     ); _i      = 0x0;	    free((void *)bjmodn0); bjmodn0 = 0x0;	    free((void *)bjmodn1); bjmodn1 = 0x0;	    free((void *)bjmodn2); bjmodn2 = 0x0;	    free((void *)bjmodn3); bjmodn3 = 0x0;	    free((void *)bjmodn4); bjmodn4 = 0x0;	    free((void *)bjmodn5); bjmodn5 = 0x0;	    free((void *)bjmodn6); bjmodn6 = 0x0;	    free((void *)bjmodn7); bjmodn7 = 0x0;	    free((void *)bjmodn8); bjmodn8 = 0x0;	    free((void *)bjmodn9); bjmodn9 = 0x0;	    free((void *)bjmodnA); bjmodnA = 0x0;	    free((void *)bjmodnB); bjmodnB = 0x0;	    free((void *)bjmodnC); bjmodnC = 0x0;	    free((void *)bjmodnD); bjmodnD = 0x0;	    free((void *)bjmodnE); bjmodnE = 0x0;	    free((void *)bjmodnF); bjmodnF = 0x0;	    free((void *)cy0    ); cy0     = 0x0;	    free((void *)cy1    ); cy1     = 0x0;	    free((void *)cy2    ); cy2     = 0x0;	    free((void *)cy3    ); cy3     = 0x0;	    free((void *)cy4    ); cy4     = 0x0;	    free((void *)cy5    ); cy5     = 0x0;	    free((void *)cy6    ); cy6     = 0x0;	    free((void *)cy7    ); cy7     = 0x0;	    free((void *)cy8    ); cy8     = 0x0;	    free((void *)cy9    ); cy9     = 0x0;	    free((void *)cyA    ); cyA     = 0x0;	    free((void *)cyB    ); cyB     = 0x0;	    free((void *)cyC    ); cyC     = 0x0;	    free((void *)cyD    ); cyD     = 0x0;	    free((void *)cyE    ); cyE     = 0x0;	    free((void *)cyF    ); cyF     = 0x0;	    free((void *)jstart ); jstart  = 0x0;	    free((void *)jhi    ); jhi     = 0x0;	    free((void *)_maxerr); _maxerr = 0x0;	    free((void *)_col   ); _col    = 0x0;	    free((void *)_co2   ); _co2    = 0x0;	    free((void *)_co3   ); _co3    = 0x0;	    free((void *)bjmodnini); bjmodnini = 0x0;	  }#ifdef CY_DBG	printf("CY_DBG: Allocating arrays...\n");#endif	  _i      = (int *)malloc(NTHREADS*sizeof(int)); if(!_i     ){ printf("FATAL: unable to allocate array _i      in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn0 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn0){ printf("FATAL: unable to allocate array bjmodn0 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn1 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn1){ printf("FATAL: unable to allocate array bjmodn1 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn2 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn2){ printf("FATAL: unable to allocate array bjmodn2 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn3 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn3){ printf("FATAL: unable to allocate array bjmodn3 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn4 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn4){ printf("FATAL: unable to allocate array bjmodn4 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn5 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn5){ printf("FATAL: unable to allocate array bjmodn5 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn6 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn6){ printf("FATAL: unable to allocate array bjmodn6 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn7 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn7){ printf("FATAL: unable to allocate array bjmodn7 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn8 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn8){ printf("FATAL: unable to allocate array bjmodn8 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn9 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn9){ printf("FATAL: unable to allocate array bjmodn9 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnA = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnA){ printf("FATAL: unable to allocate array bjmodnA in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnB = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnB){ printf("FATAL: unable to allocate array bjmodnB in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnC = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnC){ printf("FATAL: unable to allocate array bjmodnC in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnD = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnD){ printf("FATAL: unable to allocate array bjmodnD in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnE = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnE){ printf("FATAL: unable to allocate array bjmodnE in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnF = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnF){ printf("FATAL: unable to allocate array bjmodnF in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  jstart  = (int *)malloc(NTHREADS*sizeof(int)); if(!jstart ){ printf("FATAL: unable to allocate array jstart  in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  jhi     = (int *)malloc(NTHREADS*sizeof(int)); if(!jhi    ){ printf("FATAL: unable to allocate array jhi     in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _col    = (int *)malloc(NTHREADS*sizeof(int)); if(!_col   ){ printf("FATAL: unable to allocate array _col    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _co2    = (int *)malloc(NTHREADS*sizeof(int)); if(!_co2   ){ printf("FATAL: unable to allocate array _co2    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _co3    = (int *)malloc(NTHREADS*sizeof(int)); if(!_co3   ){ printf("FATAL: unable to allocate array _co3    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy0     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn0){ printf("FATAL: unable to allocate array bjmodn0 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy1     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn1){ printf("FATAL: unable to allocate array bjmodn1 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy2     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn2){ printf("FATAL: unable to allocate array bjmodn2 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy3     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn3){ printf("FATAL: unable to allocate array bjmodn3 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy4     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn4){ printf("FATAL: unable to allocate array bjmodn4 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy5     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn5){ printf("FATAL: unable to allocate array bjmodn5 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy6     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn6){ printf("FATAL: unable to allocate array bjmodn6 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy7     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn7){ printf("FATAL: unable to allocate array bjmodn7 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy8     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn8){ printf("FATAL: unable to allocate array bjmodn8 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy9     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn9){ printf("FATAL: unable to allocate array bjmodn9 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyA     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnA){ printf("FATAL: unable to allocate array bjmodnA in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyB     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnB){ printf("FATAL: unable to allocate array bjmodnB in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyC     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnC){ printf("FATAL: unable to allocate array bjmodnC in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyD     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnD){ printf("FATAL: unable to allocate array bjmodnD in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyE     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnE){ printf("FATAL: unable to allocate array bjmodnE in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyF     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnF){ printf("FATAL: unable to allocate array bjmodnF in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _maxerr = (double *)malloc(NTHREADS*sizeof(double)); if(!_maxerr){ printf("FATAL: unable to allocate array _maxerr in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnini = (int *)malloc((NTHREADS + 1)*sizeof(int)); if(!bjmodnini){ printf("FATAL: unable to allocate array bjmodnini in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnini[0] = 0;	  bjmodnini[1] = 0;#ifdef CY_DBG	printf("CY_DBG: Init bjmodnini...\n");#endif	  for(j=0; j < n16/NTHREADS; j++)	  {	    bjmodnini[1] -= sw; bjmodnini[1] = bjmodnini[1] + ( (-(int)((uint32)bjmodnini[1] >> 31)) & n);	  }	  if(NTHREADS > 1)	  {	    for(ithread = 2; ithread <= NTHREADS; ithread++)	    {	      bjmodnini[ithread] = bjmodnini[ithread-1] + bjmodnini[1] - n; bjmodnini[ithread] = bjmodnini[ithread] + ( (-(int)((uint32)bjmodnini[ithread] >> 31)) & n);            }	  }	}/*...The radix-16 final DIT pass is here.	*/#ifdef CY_DBG	printf("CY_DBG: Iter = %d; Copying A-array into B...\n", iter);	for(i=0; i < n; i++)	{		j = i + ( (i >> DAT_BITS) << PAD_BITS );		b[j] = a[j];	}	printf("CY_DBG: Calling radix16_ditN_cy_dif1_nochk...\n");	ierr = radix16_ditN_cy_dif1_nochk(b,n,nwt,nwt_bits,wt0,wt1,si,nrt,nrt_bits,rt0,rt1,base,baseinv,iter,p);	printf("IERR = %d\n", ierr);#endif    for(ithread = 0; ithread < NTHREADS; ithread++)    {	/* init carries	*/	cy0[ithread] = 0;	cy1[ithread] = 0;	cy2[ithread] = 0;	cy3[ithread] = 0;	cy4[ithread] = 0;	cy5[ithread] = 0;	cy6[ithread] = 0;	cy7[ithread] = 0;	cy8[ithread] = 0;	cy9[ithread] = 0;	cyA[ithread] = 0;	cyB[ithread] = 0;	cyC[ithread] = 0;	cyD[ithread] = 0;	cyE[ithread] = 0;	cyF[ithread] = 0;    }    cy0[0] = -2;	*fracmax=0;	/* init max. fractional error	*/	iroot = 0;	/* init sincos array index	*/	root_incr = 1;	/* init sincos array index increment (set = 1 for normal carry pass, = 0 for wrapper pass)	*/	scale = n2inv;	/* init inverse-weight scale factor  (set = 2/n for normal carry pass, = 1 for wrapper pass)	*/    for(ithread = 0; ithread < NTHREADS; ithread++)    {	_maxerr[ithread] = 0.0;    }for(outer=0; outer <= 1; outer++){	_i[0] = 1;		/* Pointer to the BASE and BASEINV arrays. lowest-order digit is always a bigword (_i[0] = 1).	*/	if(NTHREADS > 1)	{	  for(ithread = 1; ithread < NTHREADS; ithread++)	  {	    _i[ithread] = ((uint32)(sw - bjmodnini[ithread]) >> 31);          }	}	// Moved this inside the outer-loop, so on cleanup pass can use it to reset _col,_co2,_co3 starting values,	// then simply overwrite it with 1 prior to starting the k-loop.	khi = n_div_nwt/NTHREADS;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	bjmodn0[ithread] = bjmodnini[ithread];	bjmodn1[ithread] = bjmodn0[ithread] + bjmodnini[NTHREADS] - n; bjmodn1[ithread] = bjmodn1[ithread] + ( (-(int)((uint32)bjmodn1[ithread] >> 31)) & n);	bjmodn2[ithread] = bjmodn1[ithread] + bjmodnini[NTHREADS] - n; bjmodn2[ithread] = bjmodn2[ithread] + ( (-(int)((uint32)bjmodn2[ithread] >> 31)) & n);	bjmodn3[ithread] = bjmodn2[ithread] + bjmodnini[NTHREADS] - n; bjmodn3[ithread] = bjmodn3[ithread] + ( (-(int)((uint32)bjmodn3[ithread] >> 31)) & n);	bjmodn4[ithread] = bjmodn3[ithread] + bjmodnini[NTHREADS] - n; bjmodn4[ithread] = bjmodn4[ithread] + ( (-(int)((uint32)bjmodn4[ithread] >> 31)) & n);	bjmodn5[ithread] = bjmodn4[ithread] + bjmodnini[NTHREADS] - n; bjmodn5[ithread] = bjmodn5[ithread] + ( (-(int)((uint32)bjmodn5[ithread] >> 31)) & n);	bjmodn6[ithread] = bjmodn5[ithread] + bjmodnini[NTHREADS] - n; bjmodn6[ithread] = bjmodn6[ithread] + ( (-(int)((uint32)bjmodn6[ithread] >> 31)) & n);	bjmodn7[ithread] = bjmodn6[ithread] + bjmodnini[NTHREADS] - n; bjmodn7[ithread] = bjmodn7[ithread] + ( (-(int)((uint32)bjmodn7[ithread] >> 31)) & n);	bjmodn8[ithread] = bjmodn7[ithread] + bjmodnini[NTHREADS] - n; bjmodn8[ithread] = bjmodn8[ithread] + ( (-(int)((uint32)bjmodn8[ithread] >> 31)) & n);	bjmodn9[ithread] = bjmodn8[ithread] + bjmodnini[NTHREADS] - n; bjmodn9[ithread] = bjmodn9[ithread] + ( (-(int)((uint32)bjmodn9[ithread] >> 31)) & n);	bjmodnA[ithread] = bjmodn9[ithread] + bjmodnini[NTHREADS] - n; bjmodnA[ithread] = bjmodnA[ithread] + ( (-(int)((uint32)bjmodnA[ithread] >> 31)) & n);	bjmodnB[ithread] = bjmodnA[ithread] + bjmodnini[NTHREADS] - n; bjmodnB[ithread] = bjmodnB[ithread] + ( (-(int)((uint32)bjmodnB[ithread] >> 31)) & n);	bjmodnC[ithread] = bjmodnB[ithread] + bjmodnini[NTHREADS] - n; bjmodnC[ithread] = bjmodnC[ithread] + ( (-(int)((uint32)bjmodnC[ithread] >> 31)) & n);	bjmodnD[ithread] = bjmodnC[ithread] + bjmodnini[NTHREADS] - n; bjmodnD[ithread] = bjmodnD[ithread] + ( (-(int)((uint32)bjmodnD[ithread] >> 31)) & n);	bjmodnE[ithread] = bjmodnD[ithread] + bjmodnini[NTHREADS] - n; bjmodnE[ithread] = bjmodnE[ithread] + ( (-(int)((uint32)bjmodnE[ithread] >> 31)) & n);	bjmodnF[ithread] = bjmodnE[ithread] + bjmodnini[NTHREADS] - n; bjmodnF[ithread] = bjmodnF[ithread] + ( (-(int)((uint32)bjmodnF[ithread] >> 31)) & n);	jstart[ithread] = ithread*n16/NTHREADS;	if(root_incr==0)		jhi   [ithread] = jstart[ithread] + 3;		/* Cleanup loop assumes carryins propagate at most 4 words up. */	else		jhi   [ithread] = jstart[ithread] + nwt-1;	_col[ithread] = ithread*(khi*16);			/* col gets incremented by RADIX_VEC[0] on every pass through the k-loop */	_co2[ithread] = (n>>nwt_bits)-1+16 - _col[ithread];	/* co2 gets decremented by RADIX_VEC[0] on every pass through the k-loop */	_co3[ithread] = _co2[ithread]-16;			/* At the start of each new j-loop, co3=co2-RADIX_VEC[0]	*/    }    // Move this cleanup-pass-specific khi setting here, since need regular-pass khi value for above inits:    if(root_incr==0) khi = 1;#if defined(USE_THREADS) && defined(_OPENMP)/*printf("Using %d threads in carry step\n", NTHREADS);*/  omp_set_num_threads(NTHREADS);  // Needed to remove the prefetch-address vars add0 & add for this to compile properly:  #pragma omp parallel private(i,j,j1,k,l,col,co2,co3,m,m2,n_minus_sil,n_minus_silp1,sinwt,sinwtm1,wtl,wtlp1,wtn,wtnm1,wt,wtinv,wtA,wtB,wtC,maxerr,rt,it,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32,a1p0r,a1p1r,a1p2r,a1p3r,a1p4r,a1p5r,a1p6r,a1p7r,a1p8r,a1p9r,a1pAr,a1pBr,a1pCr,a1pDr,a1pEr,a1pFr,a1p0i,a1p1i,a1p2i,a1p3i,a1p4i,a1p5i,a1p6i,a1p7i,a1p8i,a1p9i,a1pAi,a1pBi,a1pCi,a1pDi,a1pEi,a1pFi,temp,frac) default(shared)  {    #pragma omp for schedule(dynamic) nowait#endif    for(ithread = 0; ithread < NTHREADS; ithread++)    {	i = _i[ithread];	/* Pointer to the BASE and BASEINV arrays.	*//*printf("KHI = %d\n",khi);*/	for(k=1; k <= khi; k++)	/* Do n/(radix(1)*nwt) outer loop executions...	*/	{/*printf("JLO = %d, JHI = %d\n",jstart[ithread],jhi[ithread]);*/	  for(j=jstart[ithread]; j < jhi[ithread]; j += 2)	/* Each inner loop execution processes (radix(1)*nwt) array data.	*/	  {	    j1 = j + ( (j >> DAT_BITS) << PAD_BITS );	/* padded-array fetch index is here	*//*	here are the needed sincos data - these will be processed below in bit-reversed order.	*/#if USE_NEW_DIT/*...Block 1:	*/	    t1 =a[j1    ];	t2 =a[j1+1    ];	    rt =a[j1+p1 ];	it =a[j1+1+p1 ];	    t3 =t1 -rt;		t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ];	t6 =a[j1+1+p2 ];	    rt =a[j1+p3 ];	it =a[j1+1+p3 ];	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ];	t10=a[j1+1+p4 ];	    rt =a[j1+p5 ];	it =a[j1+1+p5 ];	    t11=t9 -rt;		t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ];	t14=a[j1+1+p6 ];	    rt =a[j1+p7 ];	it =a[j1+1+p7 ];	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ];	t18=a[j1+1+p8 ];	    rt =a[j1+p9 ];	it =a[j1+1+p9 ];	    t19=t17-rt;		t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10];	t22=a[j1+1+p10];	    rt =a[j1+p11];	it =a[j1+1+p11];	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12];	t26=a[j1+1+p12];	    rt =a[j1+p13];	it =a[j1+1+p13];	    t27=t25-rt;		t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14];	t30=a[j1+1+p14];	    rt =a[j1+p15];	it =a[j1+1+p15];	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#else	    l = iroot;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l += iroot;			/* 2*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c1=t1*t3-t2*t4;	s1=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 4*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c2=t1*t3-t2*t4;	s2=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 8*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c4=t1*t3-t2*t4;	s4=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+(iroot << 2)+iroot;		/* 13*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c8=t1*t3-t2*t4;	s8=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c13=t1*t3-t2*t4;	s13=t1*t4+t2*t3;	      t1=c1*c4; t2=c1*s4; t3=s1*c4; t4=s1*s4;	      c3=t1+t4; s3=t2-t3; c5=t1-t4; s5=t2+t3;	      t1=c1*c8; t2=c1*s8; t3=s1*c8; t4=s1*s8;	      c7=t1+t4; s7=t2-t3; c9=t1-t4; s9=t2+t3;	      t1=c1*c13; t2=c1*s13; t3=s1*c13; t4=s1*s13;	      c12=t1+t4; s12=t2-t3; c14=t1-t4; s14=t2+t3;	      t1=c2*c8; t2=c2*s8; t3=s2*c8; t4=s2*s8;	      c6=t1+t4; s6=t2-t3; c10=t1-t4; s10=t2+t3;	      t1=c2*c13; t2=c2*s13; t3=s2*c13; t4=s2*s13;	      c11=t1+t4; s11=t2-t3; c15=t1-t4; s15=t2+t3;/*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*//*...Block 1:	*/	    t1 =a[j1    ];			t2 =a[j1+1    ];	    rt =a[j1+p1 ]*c8 +a[j1+1+p1 ]*s8;	it =a[j1+1+p1 ]*c8 -a[j1+p1 ]*s8;	    t3 =t1 -rt;  	t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ]*c4 +a[j1+1+p2 ]*s4;	t6 =a[j1+1+p2 ]*c4 -a[j1+p2 ]*s4;	    rt =a[j1+p3 ]*c12+a[j1+1+p3 ]*s12;	it =a[j1+1+p3 ]*c12-a[j1+p3 ]*s12;	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ]*c2 +a[j1+1+p4 ]*s2 ;	t10=a[j1+1+p4 ]*c2 -a[j1+p4 ]*s2;	    rt =a[j1+p5 ]*c10+a[j1+1+p5 ]*s10;	it =a[j1+1+p5 ]*c10-a[j1+p5 ]*s10;	    t11=t9 -rt;  	t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ]*c6 +a[j1+1+p6 ]*s6;	t14=a[j1+1+p6 ]*c6 -a[j1+p6 ]*s6;	    rt =a[j1+p7 ]*c14+a[j1+1+p7 ]*s14;	it =a[j1+1+p7 ]*c14-a[j1+p7 ]*s14;	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ]*c1 +a[j1+1+p8 ]*s1;	t18=a[j1+1+p8 ]*c1 -a[j1+p8 ]*s1;	    rt =a[j1+p9 ]*c9 +a[j1+1+p9 ]*s9;	it =a[j1+1+p9 ]*c9 -a[j1+p9 ]*s9;	    t19=t17-rt;  	t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10]*c5 +a[j1+1+p10]*s5 ;	t22=a[j1+1+p10]*c5 -a[j1+p10]*s5;	    rt =a[j1+p11]*c13+a[j1+1+p11]*s13;	it =a[j1+1+p11]*c13-a[j1+p11]*s13;	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12]*c3 +a[j1+1+p12]*s3 ;	t26=a[j1+1+p12]*c3 -a[j1+p12]*s3;	    rt =a[j1+p13]*c11+a[j1+1+p13]*s11;	it =a[j1+1+p13]*c11-a[j1+p13]*s11;	    t27=t25-rt;  	t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14]*c7 +a[j1+1+p14]*s7 ;	t30=a[j1+1+p14]*c7 -a[j1+p14]*s7;	    rt =a[j1+p15]*c15+a[j1+1+p15]*s15;	it =a[j1+1+p15]*c15-a[j1+p15]*s15;	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(-i* 1*twopi/16) =       ( c,-s), exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 3*twopi/16) =       ( s,-c) (for inputs to transform block 2)  	1, exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 4*twopi/16) =       ( 0,-1), exp(-i* 6*twopi/16) = ISRT2*(-1,-1) (for inputs to transform block 3)  	1, exp(-i* 3*twopi/16) =       ( s,-c), exp(-i* 6*twopi/16) = ISRT2*(-1,-1), exp(-i* 9*twopi/16) =       (-c, s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;	    a1p0r =t1+t17;	a1p0i =t2+t18;	    a1p8r =t1-t17;	a1p8i =t2-t18;	    a1p4r =t9 +t26;	a1p4i =t10-t25;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pCr=t9 -t26;	a1pCi=t10+t25;/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 -t14;	t5 =t5 +t14;		/* twiddle mpy by E^4 =-I	*/			t14=t6 +rt;	t6 =t6 -rt;	    rt =(t22+t21)*ISRT2;t22=(t22-t21)*ISRT2;	t21=rt;	/* twiddle mpy by E^-2	*/	    rt =(t29-t30)*ISRT2;it =(t29+t30)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^-6 by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;	    a1p2r =t5+t21;	a1p2i =t6+t22;	    a1pAr=t5-t21;	a1pAi=t6-t22;	    a1p6r =t13+t30;	a1p6i =t14-t29;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pEr=t13-t30;	a1pEi=t14+t29;/*...Block 2: t3,11,19,27	*/	    rt =(t12+t11)*ISRT2;it =(t12-t11)*ISRT2;		/* twiddle mpy by E^-2	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c + t20*s;	t20=t20*c - t19*s;	t19=rt;	/* twiddle mpy by E^-1	*/	    rt =t27*s + t28*c;	it =t28*s - t27*c;		/* twiddle mpy by E^-3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;	    a1p1r =t3+t19;	a1p1i =t4+t20;	    a1p9r =t3-t19;	a1p9i =t4-t20;	    a1p5r =t11+t28;	a1p5i =t12-t27;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pDr=t11-t28;	a1pDi=t12+t27;/*...Block 4: t7,15,23,31	*/	    rt =(t15-t16)*ISRT2;it =(t15+t16)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s + t24*c;	t24=t24*s - t23*c;	t23=rt;	/* twiddle mpy by E^-3	*/	    rt =t31*c + t32*s;	it =t32*c - t31*s;		/* twiddle mpy by E^-1 = -E^-9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;	    a1p3r =t7+t23;	a1p3i =t8+t24;	    a1pBr=t7-t23;	a1pBi=t8-t24;	    a1p7r =t15+t32;	a1p7i =t16-t31;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pFr=t15-t32;	a1pFi=t16+t31;/*...and combine those to complete the radix-16 transform and do the carries. Since the outputs would     normally be getting dispatched to 16 separate blocks of the A-array, we need 16 separate carries.	*/	    l= j & (nwt-1);			/* We want (S*J mod N) - SI(L) for all 16 carries, so precompute	*/	    n_minus_sil   = n-si[l  ];		/* N - SI(L) and for each J, find N - (B*J mod N) - SI(L)		*/	    n_minus_silp1 = n-si[l+1];		/* For the inverse weight, want (S*(N - J) mod N) - SI(NWT - L) =	*/	    sinwt   = si[nwt-l  ];		/*	= N - (S*J mod N) - SI(NWT - L) = (B*J mod N) - SI(NWT - L).	*/	    sinwtm1 = si[nwt-l-1];	    wtl     =wt0[    l  ];	    wtn     =wt0[nwt-l  ]*scale;	/* Include 1/(n/2) scale factor of inverse transform here...	*/	    wtlp1   =wt0[    l+1];	    wtnm1   =wt0[nwt-l-1]*scale;	/* ...and here.	*//*...set0 is slightly different from others:	*/	col = _col[ithread];	co2 = _co2[ithread];	co3 = _co3[ithread];	   cmplx_carry_norm_pow2_errcheck0(a1p0r,a1p0i,cy0[ithread],bjmodn0[ithread]);	    cmplx_carry_norm_pow2_errcheck(a1p1r,a1p1i,cy1[ithread],bjmodn1[ithread],0x1);	    cmplx_carry_norm_pow2_errcheck(a1p2r,a1p2i,cy2[ithread],bjmodn2[ithread],0x2);	    cmplx_carry_norm_pow2_errcheck(a1p3r,a1p3i,cy3[ithread],bjmodn3[ithread],0x3);	    cmplx_carry_norm_pow2_errcheck(a1p4r,a1p4i,cy4[ithread],bjmodn4[ithread],0x4);	    cmplx_carry_norm_pow2_errcheck(a1p5r,a1p5i,cy5[ithread],bjmodn5[ithread],0x5);	    cmplx_carry_norm_pow2_errcheck(a1p6r,a1p6i,cy6[ithread],bjmodn6[ithread],0x6);	    cmplx_carry_norm_pow2_errcheck(a1p7r,a1p7i,cy7[ithread],bjmodn7[ithread],0x7);	    cmplx_carry_norm_pow2_errcheck(a1p8r,a1p8i,cy8[ithread],bjmodn8[ithread],0x8);	    cmplx_carry_norm_pow2_errcheck(a1p9r,a1p9i,cy9[ithread],bjmodn9[ithread],0x9);	    cmplx_carry_norm_pow2_errcheck(a1pAr,a1pAi,cyA[ithread],bjmodnA[ithread],0xA);	    cmplx_carry_norm_pow2_errcheck(a1pBr,a1pBi,cyB[ithread],bjmodnB[ithread],0xB);	    cmplx_carry_norm_pow2_errcheck(a1pCr,a1pCi,cyC[ithread],bjmodnC[ithread],0xC);	    cmplx_carry_norm_pow2_errcheck(a1pDr,a1pDi,cyD[ithread],bjmodnD[ithread],0xD);	    cmplx_carry_norm_pow2_errcheck(a1pEr,a1pEi,cyE[ithread],bjmodnE[ithread],0xE);	    cmplx_carry_norm_pow2_errcheck(a1pFr,a1pFi,cyF[ithread],bjmodnF[ithread],0xF);	    i =((uint32)(sw - bjmodn0[ithread]) >> 31);	/* get ready for the next set...	*/	    _co2[ithread] = _co3[ithread];	/* For all data but the first set in each j-block, co2=co3. Thus, after the first block of data is done						  (and only then: for all subsequent blocks it's superfluous), this assignment decrements co2 by radix(1).	*//*...The radix-16 DIF pass is here:	*//*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*/#if PFETCHadd0 = &a[j1];prefetch_p_doubles(add0);#endif/*...Block 1:	*/	    t3 =a1p0r -a1p8r;  	t1 =a1p0r +a1p8r;	    t4 =a1p0i -a1p8i;	t2 =a1p0i +a1p8i;	    t7 =a1p4r -a1pCr;  t5 =a1p4r +a1pCr;	    t8 =a1p4i -a1pCi;  t6 =a1p4i +a1pCi;#if PFETCHadd = add0+p1;prefetch_p_doubles(add);#endif	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 +t8;	t3 =t3 -t8;			t8 =t4 -rt;	t4 =t4 +rt;#if PFETCHadd = add0+p2;prefetch_p_doubles(add);#endif/*...Block 2:	*/	    t11=a1p2r -a1pAr;  t9 =a1p2r +a1pAr;	    t12=a1p2i -a1pAi;	t10=a1p2i +a1pAi;	    t15=a1p6r -a1pEr;  t13=a1p6r +a1pEr;	    t16=a1p6i -a1pEi;	t14=a1p6i +a1pEi;#if PFETCHadd = add0+p3;prefetch_p_doubles(add);#endif	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11+t16;	t11=t11-t16;			t16=t12-rt;	t12=t12+rt;#if PFETCHadd = add0+p4;prefetch_p_doubles(add);#endif/*...Block 3:	*/	    t19=a1p1r -a1p9r;  	t17=a1p1r +a1p9r;	    t20=a1p1i -a1p9i;	t18=a1p1i +a1p9i;	    t23=a1p5r -a1pDr;  t21=a1p5r +a1pDr;	    t24=a1p5i -a1pDi;	t22=a1p5i +a1pDi;#if PFETCHadd = add0+p5;prefetch_p_doubles(add);#endif	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19+t24;	t19=t19-t24;			t24=t20-rt;	t20=t20+rt;#if PFETCHadd = add0+p6;prefetch_p_doubles(add);#endif/*...Block 4:	*/	    t27=a1p3r -a1pBr;  t25=a1p3r +a1pBr;	    t28=a1p3i -a1pBi;	t26=a1p3i +a1pBi;	    t31=a1p7r -a1pFr;  t29=a1p7r +a1pFr;	    t32=a1p7i -a1pFi;	t30=a1p7i +a1pFi;#if PFETCHadd = add0+p7;prefetch_p_doubles(add);#endif	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27+t32;	t27=t27-t32;			t32=t28-rt;	t28=t28+rt;#if PFETCHadd = add0+p8;prefetch_p_doubles(add);#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(i* 1*twopi/16) =       ( c, s), exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 3*twopi/16) =       ( s, c) (for inputs to transform block 2)  	1, exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 4*twopi/16) =       ( 0, 1), exp(i* 6*twopi/16) = ISRT2*(-1, 1) (for inputs to transform block 3)  	1, exp(i* 3*twopi/16) =       ( s, c), exp(i* 6*twopi/16) = ISRT2*(-1, 1), exp(i* 9*twopi/16) =       (-c,-s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;#if PFETCHadd = add0+p9;prefetch_p_doubles(add);#endif	    a[j1    ]=t1+t17;	a[j1+1    ]=t2+t18;	    a[j1+p1 ]=t1-t17;	a[j1+1+p1 ]=t2-t18;	    a[j1+p2 ]=t9 -t26;	a[j1+1+p2 ]=t10+t25;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p3 ]=t9 +t26;	a[j1+1+p3 ]=t10-t25;#if PFETCHadd = add0+p10;prefetch_p_doubles(add);#endif/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 +t14;	t5 =t5 -t14;		/* twiddle mpy by E^4 = I	*/			t14=t6 -rt;	t6 =t6 +rt;	    rt =(t21-t22)*ISRT2;t22=(t21+t22)*ISRT2;	t21=rt;	/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t21,22,29,30	*/	    rt =(t30+t29)*ISRT2;it =(t30-t29)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;#if PFETCHadd = add0+p11;prefetch_p_doubles(add);#endif	    a[j1+p4 ]=t5+t21;	a[j1+1+p4 ]=t6+t22;	    a[j1+p5 ]=t5-t21;	a[j1+1+p5 ]=t6-t22;	    a[j1+p6 ]=t13-t30;	a[j1+1+p6 ]=t14+t29;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p7 ]=t13+t30;	a[j1+1+p7 ]=t14-t29;#if PFETCHadd = add0+p12;prefetch_p_doubles(add);#endif/*...Block 2: t3,11,19,27	*/	    rt =(t11-t12)*ISRT2;it =(t11+t12)*ISRT2;		/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t11,12	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c - t20*s;	t20=t20*c + t19*s;	t19=rt;	/* twiddle mpy by E^1	*/	    rt =t27*s - t28*c;	it =t28*s + t27*c;		/* twiddle mpy by E^3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;#if PFETCHadd = add0+p13;prefetch_p_doubles(add);#endif	    a[j1+p8 ]=t3+t19;	a[j1+1+p8 ]=t4+t20;	    a[j1+p9 ]=t3-t19;	a[j1+1+p9 ]=t4-t20;	    a[j1+p10]=t11-t28;	a[j1+1+p10]=t12+t27;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p11]=t11+t28;	a[j1+1+p11]=t12-t27;#if PFETCHadd = add0+p14;prefetch_p_doubles(add);#endif/*...Block 4: t7,15,23,31	*/	    rt =(t16+t15)*ISRT2;it =(t16-t15)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s - t24*c;	t24=t24*s + t23*c;	t23=rt;	/* twiddle mpy by E^3	*/	    rt =t31*c - t32*s;	it =t32*c + t31*s;		/* twiddle mpy by E^1 = -E^9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;			/* Note: t23+rt = t23*(s+1)	*/#if PFETCHadd = add0+p15;prefetch_p_doubles(add);#endif	    a[j1+p12]=t7+t23;	a[j1+1+p12]=t8+t24;	    a[j1+p13]=t7-t23;	a[j1+1+p13]=t8-t24;	    a[j1+p14]=t15-t32;	a[j1+1+p14]=t16+t31;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p15]=t15+t32;	a[j1+1+p15]=t16-t31;#ifdef CY_DBG//#if 0	if(a[j1      ] != b[j1      ]){ printf("j1   = %d :p0 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn0 = %d\n", j1  , a[j1      ], b[j1      ], ithread, i, l, col, co2, co3, bjmodn0[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1    ] != b[j1+1    ]){ printf("j1+1 = %d :p0 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn0 = %d\n", j1+1, a[j1+1    ], b[j1+1    ], ithread, i, l, col, co2, co3, bjmodn0[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p1 ] != b[j1  +p1 ]){ printf("j1   = %d :p1 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn1 = %d\n", j1  , a[j1  +p1 ], b[j1  +p1 ], ithread, i, l, col, co2, co3, bjmodn1[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p1 ] != b[j1+1+p1 ]){ printf("j1+1 = %d :p1 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn1 = %d\n", j1+1, a[j1+1+p1 ], b[j1+1+p1 ], ithread, i, l, col, co2, co3, bjmodn1[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p2 ] != b[j1  +p2 ]){ printf("j1   = %d :p2 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn2 = %d\n", j1  , a[j1  +p2 ], b[j1  +p2 ], ithread, i, l, col, co2, co3, bjmodn2[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p2 ] != b[j1+1+p2 ]){ printf("j1+1 = %d :p2 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn2 = %d\n", j1+1, a[j1+1+p2 ], b[j1+1+p2 ], ithread, i, l, col, co2, co3, bjmodn2[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p3 ] != b[j1  +p3 ]){ printf("j1   = %d :p3 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn3 = %d\n", j1  , a[j1  +p3 ], b[j1  +p3 ], ithread, i, l, col, co2, co3, bjmodn3[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p3 ] != b[j1+1+p3 ]){ printf("j1+1 = %d :p3 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn3 = %d\n", j1+1, a[j1+1+p3 ], b[j1+1+p3 ], ithread, i, l, col, co2, co3, bjmodn3[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p4 ] != b[j1  +p4 ]){ printf("j1   = %d :p4 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn4 = %d\n", j1  , a[j1  +p4 ], b[j1  +p4 ], ithread, i, l, col, co2, co3, bjmodn4[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p4 ] != b[j1+1+p4 ]){ printf("j1+1 = %d :p4 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn4 = %d\n", j1+1, a[j1+1+p4 ], b[j1+1+p4 ], ithread, i, l, col, co2, co3, bjmodn4[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p5 ] != b[j1  +p5 ]){ printf("j1   = %d :p5 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn5 = %d\n", j1  , a[j1  +p5 ], b[j1  +p5 ], ithread, i, l, col, co2, co3, bjmodn5[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p5 ] != b[j1+1+p5 ]){ printf("j1+1 = %d :p5 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn5 = %d\n", j1+1, a[j1+1+p5 ], b[j1+1+p5 ], ithread, i, l, col, co2, co3, bjmodn5[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p6 ] != b[j1  +p6 ]){ printf("j1   = %d :p6 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn6 = %d\n", j1  , a[j1  +p6 ], b[j1  +p6 ], ithread, i, l, col, co2, co3, bjmodn6[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p6 ] != b[j1+1+p6 ]){ printf("j1+1 = %d :p6 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn6 = %d\n", j1+1, a[j1+1+p6 ], b[j1+1+p6 ], ithread, i, l, col, co2, co3, bjmodn6[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p7 ] != b[j1  +p7 ]){ printf("j1   = %d :p7 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn7 = %d\n", j1  , a[j1  +p7 ], b[j1  +p7 ], ithread, i, l, col, co2, co3, bjmodn7[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p7 ] != b[j1+1+p7 ]){ printf("j1+1 = %d :p7 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn7 = %d\n", j1+1, a[j1+1+p7 ], b[j1+1+p7 ], ithread, i, l, col, co2, co3, bjmodn7[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p8 ] != b[j1  +p8 ]){ printf("j1   = %d :p8 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn8 = %d\n", j1  , a[j1  +p8 ], b[j1  +p8 ], ithread, i, l, col, co2, co3, bjmodn8[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p8 ] != b[j1+1+p8 ]){ printf("j1+1 = %d :p8 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn8 = %d\n", j1+1, a[j1+1+p8 ], b[j1+1+p8 ], ithread, i, l, col, co2, co3, bjmodn8[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p9 ] != b[j1  +p9 ]){ printf("j1   = %d :p9 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn9 = %d\n", j1  , a[j1  +p9 ], b[j1  +p9 ], ithread, i, l, col, co2, co3, bjmodn9[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p9 ] != b[j1+1+p9 ]){ printf("j1+1 = %d :p9 : A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodn9 = %d\n", j1+1, a[j1+1+p9 ], b[j1+1+p9 ], ithread, i, l, col, co2, co3, bjmodn9[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p10] != b[j1  +p10]){ printf("j1   = %d :p10: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnA = %d\n", j1  , a[j1  +p10], b[j1  +p10], ithread, i, l, col, co2, co3, bjmodnA[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p10] != b[j1+1+p10]){ printf("j1+1 = %d :p10: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnA = %d\n", j1+1, a[j1+1+p10], b[j1+1+p10], ithread, i, l, col, co2, co3, bjmodnA[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p11] != b[j1  +p11]){ printf("j1   = %d :p11: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnB = %d\n", j1  , a[j1  +p11], b[j1  +p11], ithread, i, l, col, co2, co3, bjmodnB[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p11] != b[j1+1+p11]){ printf("j1+1 = %d :p11: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnB = %d\n", j1+1, a[j1+1+p11], b[j1+1+p11], ithread, i, l, col, co2, co3, bjmodnB[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p12] != b[j1  +p12]){ printf("j1   = %d :p12: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnC = %d\n", j1  , a[j1  +p12], b[j1  +p12], ithread, i, l, col, co2, co3, bjmodnC[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p12] != b[j1+1+p12]){ printf("j1+1 = %d :p12: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnC = %d\n", j1+1, a[j1+1+p12], b[j1+1+p12], ithread, i, l, col, co2, co3, bjmodnC[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p13] != b[j1  +p13]){ printf("j1   = %d :p13: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnD = %d\n", j1  , a[j1  +p13], b[j1  +p13], ithread, i, l, col, co2, co3, bjmodnD[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p13] != b[j1+1+p13]){ printf("j1+1 = %d :p13: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnD = %d\n", j1+1, a[j1+1+p13], b[j1+1+p13], ithread, i, l, col, co2, co3, bjmodnD[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p14] != b[j1  +p14]){ printf("j1   = %d :p14: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnE = %d\n", j1  , a[j1  +p14], b[j1  +p14], ithread, i, l, col, co2, co3, bjmodnE[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p14] != b[j1+1+p14]){ printf("j1+1 = %d :p14: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnE = %d\n", j1+1, a[j1+1+p14], b[j1+1+p14], ithread, i, l, col, co2, co3, bjmodnE[ithread]); exit(EXIT_FAILURE); }	if(a[j1  +p15] != b[j1  +p15]){ printf("j1   = %d :p15: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnF = %d\n", j1  , a[j1  +p15], b[j1  +p15], ithread, i, l, col, co2, co3, bjmodnF[ithread]); exit(EXIT_FAILURE); }	if(a[j1+1+p15] != b[j1+1+p15]){ printf("j1+1 = %d :p15: A = %20.15f, B = %20.15f, thread = %d, i = %d, l = %d, col/co2/co3 = %d, %d, %d, bjmodnF = %d\n", j1+1, a[j1+1+p15], b[j1+1+p15], ithread, i, l, col, co2, co3, bjmodnF[ithread]); exit(EXIT_FAILURE); }#endif	  iroot += root_incr;		/* increment sincos index.	*/	  }	  jstart[ithread] += nwt;	  jhi   [ithread] += nwt;	  _col[ithread] += 16;	  _co3[ithread] -= 16;	  // Since will lose separate maxerr values when threads are merged, save them after each pass.	  if(_maxerr[ithread] < maxerr)		  _maxerr[ithread] = maxerr;	}    }	/******* END OF PARALLEL FOR-LOOP ********/#if defined(USE_THREADS) && defined(_OPENMP)  }#endif	if(root_incr==0)break;/*   Wraparound carry cleanup loop is here:     The cleanup carries from the end of each length-N/16 block into the begining of the next     can all be neatly processed as follows:     (1) Invert the radix-16 forward DIF FFT of the first block of 16 complex elements in A and unweight;     (2) Propagate cleanup carries among the real and imaginary parts of the 16 outputs of (1);     (3) Reweight and perform a radix-16 forward DIF FFT on the result of (2);     (4) If any of the exit carries from (2) are nonzero, advance to the next 16 elements and repeat (1-4).*/	t1  = cyF[NTHREADS - 1];    if(NTHREADS > 1)    {      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyF[ithread] = cyF[ithread-1]; }; cyF[0] = cyE[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyE[ithread] = cyE[ithread-1]; }; cyE[0] = cyD[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyD[ithread] = cyD[ithread-1]; }; cyD[0] = cyC[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyC[ithread] = cyC[ithread-1]; }; cyC[0] = cyB[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyB[ithread] = cyB[ithread-1]; }; cyB[0] = cyA[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyA[ithread] = cyA[ithread-1]; }; cyA[0] = cy9[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy9[ithread] = cy9[ithread-1]; }; cy9[0] = cy8[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy8[ithread] = cy8[ithread-1]; }; cy8[0] = cy7[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy7[ithread] = cy7[ithread-1]; }; cy7[0] = cy6[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy6[ithread] = cy6[ithread-1]; }; cy6[0] = cy5[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy5[ithread] = cy5[ithread-1]; }; cy5[0] = cy4[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy4[ithread] = cy4[ithread-1]; }; cy4[0] = cy3[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy3[ithread] = cy3[ithread-1]; }; cy3[0] = cy2[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy2[ithread] = cy2[ithread-1]; }; cy2[0] = cy1[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy1[ithread] = cy1[ithread-1]; }; cy1[0] = cy0[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy0[ithread] = cy0[ithread-1]; }; cy0[0] = t1;    }    else    {	cyF[0]= cyE[0];	cyE[0]= cyD[0];	cyD[0]= cyC[0];	cyC[0]= cyB[0];	cyB[0]= cyA[0];	cyA[0]= cy9[0];	cy9[0]= cy8[0];	cy8[0]= cy7[0];	cy7[0]= cy6[0];	cy6[0]= cy5[0];	cy5[0]= cy4[0];	cy4[0]= cy3[0];	cy3[0]= cy2[0];	cy2[0]= cy1[0];	cy1[0]= cy0[0];	cy0[0]= t1;    }	iroot = 0;	root_incr = 0;	scale = 1;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	for(j = ithread*pini; j <= ithread*pini + 3; j++)	{	  a[j    ] *= radix_inv;	  a[j+p1 ] *= radix_inv;	  a[j+p2 ] *= radix_inv;	  a[j+p3 ] *= radix_inv;	  a[j+p4 ] *= radix_inv;	  a[j+p5 ] *= radix_inv;	  a[j+p6 ] *= radix_inv;	  a[j+p7 ] *= radix_inv;	  a[j+p8 ] *= radix_inv;	  a[j+p9 ] *= radix_inv;	  a[j+p10] *= radix_inv;	  a[j+p11] *= radix_inv;	  a[j+p12] *= radix_inv;	  a[j+p13] *= radix_inv;	  a[j+p14] *= radix_inv;	  a[j+p15] *= radix_inv;	}    }}/***********************************///exit(EXIT_FAILURE);/***********************************//*printf("cy0[0] = %20.10f\n", cy0[0]);printf("cy1[0] = %20.10f\n", cy1[0]);printf("cy2[0] = %20.10f\n", cy2[0]);printf("cy3[0] = %20.10f\n", cy3[0]);printf("cy4[0] = %20.10f\n", cy4[0]);printf("cy5[0] = %20.10f\n", cy5[0]);printf("cy6[0] = %20.10f\n", cy6[0]);printf("cy7[0] = %20.10f\n", cy7[0]);printf("cy8[0] = %20.10f\n", cy8[0]);printf("cy9[0] = %20.10f\n", cy9[0]);printf("cyA[0] = %20.10f\n", cyA[0]);printf("cyB[0] = %20.10f\n", cyB[0]);printf("cyC[0] = %20.10f\n", cyC[0]);printf("cyD[0] = %20.10f\n", cyD[0]);printf("cyE[0] = %20.10f\n", cyE[0]);printf("cyF[0] = %20.10f\n", cyF[0]);*/    t1 = 0;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	t1 += fabs(cy0[ithread])+fabs(cy1[ithread])+fabs(cy2[ithread])+fabs(cy3[ithread])+fabs(cy4[ithread])+fabs(cy5[ithread])+fabs(cy6[ithread])+fabs(cy7[ithread])		+fabs(cy8[ithread])+fabs(cy9[ithread])+fabs(cyA[ithread])+fabs(cyB[ithread])+fabs(cyC[ithread])+fabs(cyD[ithread])+fabs(cyE[ithread])+fabs(cyF[ithread]);	if(*fracmax < _maxerr[ithread])		*fracmax = _maxerr[ithread];    }	if(t1 != 0.0)	{	    if(INTERACT)printf("FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    f1 = fopen(   OFILE,"a");	    f2 = fopen(STATFILE,"a");	    fprintf(f1,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    fprintf(f2,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    fclose(f1);	    fclose(f2);	    err=1;	    return(err);	}	/*printf("ifracmax = %20llu\n",ifracmax);*/	return(0);}/**************/int radix16_ditN_cy_dif1_nochk(double a[], int n, int nwt, int nwt_bits, double wt0[], double wt1[], int si[], int nrt, int nrt_bits, struct complex rt0[],                         struct complex rt1[], double base[], double baseinv[], int iter                 , int p){/*!...Acronym: DWT = Discrete Weighted Transform, DIT = Decimation In Time, DIF = Decimation In Frequency!!...Performs a final radix-16 complex DIT pass, an inverse DWT weighting, a carry propagation,!   a forward DWT weighting, and an initial radix-16 complex DIF pass on the data in the length-N real vector A.!!   Data enter and are returned in the A-array.!!   See the documentation in mers_mod_square and radix16_dif_pass for further details on the array!   storage scheme, and radix16_ditN_cy_dif1 for details on the reduced-length weights array scheme.*/#if defined(USE_THREADS) && defined(_OPENMP)	// **************** MULTI-THREADED VERSION ********************	int n16,ithread,i,j,j1,iroot,root_incr,k1,k2,k,khi,l,outer;	static int *_i, *bjmodn0 = 0x0, *bjmodn1 = 0x0, *bjmodn2 = 0x0, *bjmodn3 = 0x0, *bjmodn4 = 0x0, *bjmodn5 = 0x0, *bjmodn6 = 0x0, *bjmodn7 = 0x0, *bjmodn8 = 0x0, *bjmodn9 = 0x0, *bjmodnA = 0x0, *bjmodnB = 0x0, *bjmodnC = 0x0, *bjmodnD = 0x0, *bjmodnE = 0x0, *bjmodnF = 0x0	, *bjmodnini = 0x0, *jstart = 0x0, *jhi = 0x0, *_col = 0x0, *_co2 = 0x0, *_co3 = 0x0;	static int psave=0,bw,sw,nm1,nrtm1,pini,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15;	const  double one_half[3] = {1.0, 0.5, 0.25};	/* Needed for small-weights-tables scheme */	static double c = 0.9238795325112867561281831, s = 0.3826834323650897717284599, radix_inv,n2inv	, *cy0 = 0x0, *cy1 = 0x0, *cy2 = 0x0, *cy3 = 0x0, *cy4 = 0x0, *cy5 = 0x0, *cy6 = 0x0, *cy7 = 0x0, *cy8 = 0x0, *cy9 = 0x0, *cyA = 0x0, *cyB = 0x0, *cyC = 0x0, *cyD = 0x0, *cyE = 0x0, *cyF = 0x0;	double rt,it#if !USE_NEW_DIT		,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15		,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15#endif	,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32	,a1p0r,a1p1r,a1p2r,a1p3r,a1p4r,a1p5r,a1p6r,a1p7r,a1p8r,a1p9r,a1pAr,a1pBr,a1pCr,a1pDr,a1pEr,a1pFr	,a1p0i,a1p1i,a1p2i,a1p3i,a1p4i,a1p5i,a1p6i,a1p7i,a1p8i,a1p9i,a1pAi,a1pBi,a1pCi,a1pDi,a1pEi,a1pFi	,temp,scale;#if PFETCH	double *add0, *add;#endif	int err;	FILE *f1,*f2;	static int first_entry=TRUE;/*...stuff for the reduced-length DWT weights array is here:	*/	int n_div_nwt;	int col,co2,co3,m,m2,n_minus_sil,n_minus_silp1,sinwt,sinwtm1;	double wt,wtinv,wtl,wtlp1,wtn,wtnm1,wtA,wtB,wtC;#ifdef CY_DBG	int ierr, npad;	static double*b = 0x0;#endif/*...change n16 and n_div_wt to non-static to work around a gcc compiler bug. */	n16   = n/16;	n_div_nwt = n16 >> nwt_bits;	if(NTHREADS > 1)	{		assert(n16      %NTHREADS == 0);		assert(n_div_nwt%NTHREADS == 0);	}	if((n_div_nwt << nwt_bits) != n16)	{	  if(INTERACT)printf("FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  f1 = fopen(   OFILE,"a");	  f2 = fopen(STATFILE,"a");	  fprintf(f1,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fprintf(f2,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fclose(f1);	  fclose(f2);	  err=1;	  return(err);	}	if(p != psave)	{	  first_entry=TRUE;	}/*...initialize things upon first entry	*/	if(first_entry)	{#ifdef CY_DBG	npad = n + ( (n >> DAT_BITS) << PAD_BITS );	/* length of padded data array.	*/	printf("CY_DBG: allocating scratch B-array of %d doubles...\n", npad);	b = (double *)malloc(npad*sizeof(double));	if(!b){ printf("FATAL: unable to allocate array B.\n"); exit(EXIT_FAILURE); }#endif	  psave = p;	  first_entry=FALSE;	  radix_inv = (double)((long double)1.0L/(long double)16.0L);	  n2inv     = (double)((long double)1.0L/(long double)(n/2));	  bw    = p%n;		/* Number of bigwords in the Crandall/Fagin mixed-radix representation = (Mersenne exponent) mod (vector length).	*/	  sw    = n - bw;	/* Number of smallwords.	*/	  nm1   = n-1;	  nrtm1 = nrt - 1;/*   constant index offsets for load/stores are here.	*/	  pini = n16/NTHREADS;	  pini += ( (pini >> DAT_BITS) << PAD_BITS );	  p1 = n16 + ( (n16 >> DAT_BITS) << PAD_BITS );	  p2 = p1 +p1;	  p3 = p2 +p1;	  p4 = p3 +p1;	  p5 = p4 +p1;	  p6 = p5 +p1;	  p7 = p6 +p1;	  p8 = p7 +p1;	  p9 = p8 +p1;	  p10= p9 +p1;	  p11= p10+p1;	  p12= p11+p1;	  p13= p12+p1;	  p14= p13+p1;	  p15= p14+p1;	  if(cy0)	/* If it's a new exponent of a range test, need to deallocate these. */	  {	    free((void *)_i     ); _i      = 0x0;	    free((void *)bjmodn0); bjmodn0 = 0x0;	    free((void *)bjmodn1); bjmodn1 = 0x0;	    free((void *)bjmodn2); bjmodn2 = 0x0;	    free((void *)bjmodn3); bjmodn3 = 0x0;	    free((void *)bjmodn4); bjmodn4 = 0x0;	    free((void *)bjmodn5); bjmodn5 = 0x0;	    free((void *)bjmodn6); bjmodn6 = 0x0;	    free((void *)bjmodn7); bjmodn7 = 0x0;	    free((void *)bjmodn8); bjmodn8 = 0x0;	    free((void *)bjmodn9); bjmodn9 = 0x0;	    free((void *)bjmodnA); bjmodnA = 0x0;	    free((void *)bjmodnB); bjmodnB = 0x0;	    free((void *)bjmodnC); bjmodnC = 0x0;	    free((void *)bjmodnD); bjmodnD = 0x0;	    free((void *)bjmodnE); bjmodnE = 0x0;	    free((void *)bjmodnF); bjmodnF = 0x0;	    free((void *)cy0    ); cy0     = 0x0;	    free((void *)cy1    ); cy1     = 0x0;	    free((void *)cy2    ); cy2     = 0x0;	    free((void *)cy3    ); cy3     = 0x0;	    free((void *)cy4    ); cy4     = 0x0;	    free((void *)cy5    ); cy5     = 0x0;	    free((void *)cy6    ); cy6     = 0x0;	    free((void *)cy7    ); cy7     = 0x0;	    free((void *)cy8    ); cy8     = 0x0;	    free((void *)cy9    ); cy9     = 0x0;	    free((void *)cyA    ); cyA     = 0x0;	    free((void *)cyB    ); cyB     = 0x0;	    free((void *)cyC    ); cyC     = 0x0;	    free((void *)cyD    ); cyD     = 0x0;	    free((void *)cyE    ); cyE     = 0x0;	    free((void *)cyF    ); cyF     = 0x0;	    free((void *)jstart ); jstart  = 0x0;	    free((void *)jhi    ); jhi     = 0x0;	    free((void *)_col   ); _col    = 0x0;	    free((void *)_co2   ); _co2    = 0x0;	    free((void *)_co3   ); _co3    = 0x0;	    free((void *)bjmodnini); bjmodnini = 0x0;	  }#ifdef CY_DBG	printf("CY_DBG: Allocating arrays...\n");#endif	  _i      = (int *)malloc(NTHREADS*sizeof(int)); if(!_i     ){ printf("FATAL: unable to allocate array _i      in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn0 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn0){ printf("FATAL: unable to allocate array bjmodn0 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn1 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn1){ printf("FATAL: unable to allocate array bjmodn1 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn2 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn2){ printf("FATAL: unable to allocate array bjmodn2 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn3 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn3){ printf("FATAL: unable to allocate array bjmodn3 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn4 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn4){ printf("FATAL: unable to allocate array bjmodn4 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn5 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn5){ printf("FATAL: unable to allocate array bjmodn5 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn6 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn6){ printf("FATAL: unable to allocate array bjmodn6 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn7 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn7){ printf("FATAL: unable to allocate array bjmodn7 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn8 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn8){ printf("FATAL: unable to allocate array bjmodn8 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodn9 = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodn9){ printf("FATAL: unable to allocate array bjmodn9 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnA = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnA){ printf("FATAL: unable to allocate array bjmodnA in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnB = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnB){ printf("FATAL: unable to allocate array bjmodnB in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnC = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnC){ printf("FATAL: unable to allocate array bjmodnC in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnD = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnD){ printf("FATAL: unable to allocate array bjmodnD in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnE = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnE){ printf("FATAL: unable to allocate array bjmodnE in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnF = (int *)malloc(NTHREADS*sizeof(int)); if(!bjmodnF){ printf("FATAL: unable to allocate array bjmodnF in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  jstart  = (int *)malloc(NTHREADS*sizeof(int)); if(!jstart ){ printf("FATAL: unable to allocate array jstart  in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  jhi     = (int *)malloc(NTHREADS*sizeof(int)); if(!jhi    ){ printf("FATAL: unable to allocate array jhi     in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _col    = (int *)malloc(NTHREADS*sizeof(int)); if(!_col   ){ printf("FATAL: unable to allocate array _col    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _co2    = (int *)malloc(NTHREADS*sizeof(int)); if(!_co2   ){ printf("FATAL: unable to allocate array _co2    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  _co3    = (int *)malloc(NTHREADS*sizeof(int)); if(!_co3   ){ printf("FATAL: unable to allocate array _co3    in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy0     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn0){ printf("FATAL: unable to allocate array bjmodn0 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy1     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn1){ printf("FATAL: unable to allocate array bjmodn1 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy2     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn2){ printf("FATAL: unable to allocate array bjmodn2 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy3     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn3){ printf("FATAL: unable to allocate array bjmodn3 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy4     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn4){ printf("FATAL: unable to allocate array bjmodn4 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy5     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn5){ printf("FATAL: unable to allocate array bjmodn5 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy6     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn6){ printf("FATAL: unable to allocate array bjmodn6 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy7     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn7){ printf("FATAL: unable to allocate array bjmodn7 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy8     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn8){ printf("FATAL: unable to allocate array bjmodn8 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cy9     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodn9){ printf("FATAL: unable to allocate array bjmodn9 in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyA     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnA){ printf("FATAL: unable to allocate array bjmodnA in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyB     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnB){ printf("FATAL: unable to allocate array bjmodnB in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyC     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnC){ printf("FATAL: unable to allocate array bjmodnC in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyD     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnD){ printf("FATAL: unable to allocate array bjmodnD in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyE     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnE){ printf("FATAL: unable to allocate array bjmodnE in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  cyF     = (double *)malloc(NTHREADS*sizeof(double)); if(!bjmodnF){ printf("FATAL: unable to allocate array bjmodnF in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnini = (int *)malloc((NTHREADS + 1)*sizeof(int)); if(!bjmodnini){ printf("FATAL: unable to allocate array bjmodnini in radix16_ditN_cy_dif1.\n"); exit(EXIT_FAILURE); }	  bjmodnini[0] = 0;	  bjmodnini[1] = 0;#ifdef CY_DBG	printf("CY_DBG: Init bjmodnini...\n");#endif	  for(j=0; j < n16/NTHREADS; j++)	  {	    bjmodnini[1] -= sw; bjmodnini[1] = bjmodnini[1] + ( (-(int)((uint32)bjmodnini[1] >> 31)) & n);	  }	  if(NTHREADS > 1)	  {	    for(ithread = 2; ithread <= NTHREADS; ithread++)	    {	      bjmodnini[ithread] = bjmodnini[ithread-1] + bjmodnini[1] - n; bjmodnini[ithread] = bjmodnini[ithread] + ( (-(int)((uint32)bjmodnini[ithread] >> 31)) & n);            }	  }	}/*...The radix-16 final DIT pass is here.	*/#ifdef CY_DBG	printf("CY_DBG: Iter = %d; Copying A-array into B...\n", iter);	for(i=0; i < n; i++)	{		j = i + ( (i >> DAT_BITS) << PAD_BITS );		b[j] = a[j];	}	printf("CY_DBG: Calling radix16_ditN_cy_dif1_nochk...\n");	ierr = radix16_ditN_cy_dif1_nochk(b,n,nwt,nwt_bits,wt0,wt1,si,nrt,nrt_bits,rt0,rt1,base,baseinv,iter,p);	printf("IERR = %d\n", ierr);#endif    for(ithread = 0; ithread < NTHREADS; ithread++)    {	/* init carries	*/	cy0[ithread] = 0;	cy1[ithread] = 0;	cy2[ithread] = 0;	cy3[ithread] = 0;	cy4[ithread] = 0;	cy5[ithread] = 0;	cy6[ithread] = 0;	cy7[ithread] = 0;	cy8[ithread] = 0;	cy9[ithread] = 0;	cyA[ithread] = 0;	cyB[ithread] = 0;	cyC[ithread] = 0;	cyD[ithread] = 0;	cyE[ithread] = 0;	cyF[ithread] = 0;    }    cy0[0] = -2;	iroot = 0;	/* init sincos array index	*/	root_incr = 1;	/* init sincos array index increment (set = 1 for normal carry pass, = 0 for wrapper pass)	*/	scale = n2inv;	/* init inverse-weight scale factor  (set = 2/n for normal carry pass, = 1 for wrapper pass)	*/for(outer=0; outer <= 1; outer++){	_i[0] = 1;		/* Pointer to the BASE and BASEINV arrays. lowest-order digit is always a bigword (_i[0] = 1).	*/	if(NTHREADS > 1)	{	  for(ithread = 1; ithread < NTHREADS; ithread++)	  {	    _i[ithread] = ((uint32)(sw - bjmodnini[ithread]) >> 31);          }	}	// Moved this inside the outer-loop, so on cleanup pass can use it to reset _col,_co2,_co3 starting values,	// then simply overwrite it with 1 prior to starting the k-loop.	khi = n_div_nwt/NTHREADS;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	bjmodn0[ithread] = bjmodnini[ithread];	bjmodn1[ithread] = bjmodn0[ithread] + bjmodnini[NTHREADS] - n; bjmodn1[ithread] = bjmodn1[ithread] + ( (-(int)((uint32)bjmodn1[ithread] >> 31)) & n);	bjmodn2[ithread] = bjmodn1[ithread] + bjmodnini[NTHREADS] - n; bjmodn2[ithread] = bjmodn2[ithread] + ( (-(int)((uint32)bjmodn2[ithread] >> 31)) & n);	bjmodn3[ithread] = bjmodn2[ithread] + bjmodnini[NTHREADS] - n; bjmodn3[ithread] = bjmodn3[ithread] + ( (-(int)((uint32)bjmodn3[ithread] >> 31)) & n);	bjmodn4[ithread] = bjmodn3[ithread] + bjmodnini[NTHREADS] - n; bjmodn4[ithread] = bjmodn4[ithread] + ( (-(int)((uint32)bjmodn4[ithread] >> 31)) & n);	bjmodn5[ithread] = bjmodn4[ithread] + bjmodnini[NTHREADS] - n; bjmodn5[ithread] = bjmodn5[ithread] + ( (-(int)((uint32)bjmodn5[ithread] >> 31)) & n);	bjmodn6[ithread] = bjmodn5[ithread] + bjmodnini[NTHREADS] - n; bjmodn6[ithread] = bjmodn6[ithread] + ( (-(int)((uint32)bjmodn6[ithread] >> 31)) & n);	bjmodn7[ithread] = bjmodn6[ithread] + bjmodnini[NTHREADS] - n; bjmodn7[ithread] = bjmodn7[ithread] + ( (-(int)((uint32)bjmodn7[ithread] >> 31)) & n);	bjmodn8[ithread] = bjmodn7[ithread] + bjmodnini[NTHREADS] - n; bjmodn8[ithread] = bjmodn8[ithread] + ( (-(int)((uint32)bjmodn8[ithread] >> 31)) & n);	bjmodn9[ithread] = bjmodn8[ithread] + bjmodnini[NTHREADS] - n; bjmodn9[ithread] = bjmodn9[ithread] + ( (-(int)((uint32)bjmodn9[ithread] >> 31)) & n);	bjmodnA[ithread] = bjmodn9[ithread] + bjmodnini[NTHREADS] - n; bjmodnA[ithread] = bjmodnA[ithread] + ( (-(int)((uint32)bjmodnA[ithread] >> 31)) & n);	bjmodnB[ithread] = bjmodnA[ithread] + bjmodnini[NTHREADS] - n; bjmodnB[ithread] = bjmodnB[ithread] + ( (-(int)((uint32)bjmodnB[ithread] >> 31)) & n);	bjmodnC[ithread] = bjmodnB[ithread] + bjmodnini[NTHREADS] - n; bjmodnC[ithread] = bjmodnC[ithread] + ( (-(int)((uint32)bjmodnC[ithread] >> 31)) & n);	bjmodnD[ithread] = bjmodnC[ithread] + bjmodnini[NTHREADS] - n; bjmodnD[ithread] = bjmodnD[ithread] + ( (-(int)((uint32)bjmodnD[ithread] >> 31)) & n);	bjmodnE[ithread] = bjmodnD[ithread] + bjmodnini[NTHREADS] - n; bjmodnE[ithread] = bjmodnE[ithread] + ( (-(int)((uint32)bjmodnE[ithread] >> 31)) & n);	bjmodnF[ithread] = bjmodnE[ithread] + bjmodnini[NTHREADS] - n; bjmodnF[ithread] = bjmodnF[ithread] + ( (-(int)((uint32)bjmodnF[ithread] >> 31)) & n);	jstart[ithread] = ithread*n16/NTHREADS;	if(root_incr==0)		jhi   [ithread] = jstart[ithread] + 3;		/* Cleanup loop assumes carryins propagate at most 4 words up. */	else		jhi   [ithread] = jstart[ithread] + nwt-1;	_col[ithread] = ithread*(khi*16);			/* col gets incremented by RADIX_VEC[0] on every pass through the k-loop */	_co2[ithread] = (n>>nwt_bits)-1+16 - _col[ithread];	/* co2 gets decremented by RADIX_VEC[0] on every pass through the k-loop */	_co3[ithread] = _co2[ithread]-16;			/* At the start of each new j-loop, co3=co2-RADIX_VEC[0]	*/    }    // Move this cleanup-pass-specific khi setting here, since need regular-pass khi value for above inits:    if(root_incr==0) khi = 1;/*printf("Using %d threads in carry step\n", NTHREADS);*/  omp_set_num_threads(NTHREADS);  // Needed to remove the prefetch-address vars add0 & add for this to compile properly:  #pragma omp parallel private(i,j,j1,k,l,col,co2,co3,m,m2,n_minus_sil,n_minus_silp1,sinwt,sinwtm1,wtl,wtlp1,wtn,wtnm1,wt,wtinv,wtA,wtB,wtC,rt,it,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32,a1p0r,a1p1r,a1p2r,a1p3r,a1p4r,a1p5r,a1p6r,a1p7r,a1p8r,a1p9r,a1pAr,a1pBr,a1pCr,a1pDr,a1pEr,a1pFr,a1p0i,a1p1i,a1p2i,a1p3i,a1p4i,a1p5i,a1p6i,a1p7i,a1p8i,a1p9i,a1pAi,a1pBi,a1pCi,a1pDi,a1pEi,a1pFi,temp) default(shared)  {    #pragma omp for schedule(dynamic) nowait    for(ithread = 0; ithread < NTHREADS; ithread++)    {	i = _i[ithread];	/* Pointer to the BASE and BASEINV arrays.	*//*printf("KHI = %d\n",khi);*/	for(k=1; k <= khi; k++)	/* Do n/(radix(1)*nwt) outer loop executions...	*/	{/*printf("JLO = %d, JHI = %d\n",jstart[ithread],jhi[ithread]);*/	  for(j=jstart[ithread]; j < jhi[ithread]; j += 2)	/* Each inner loop execution processes (radix(1)*nwt) array data.	*/	  {	    j1 = j + ( (j >> DAT_BITS) << PAD_BITS );	/* padded-array fetch index is here	*//*	here are the needed sincos data - these will be processed below in bit-reversed order.	*/#if USE_NEW_DIT/*...Block 1:	*/	    t1 =a[j1    ];	t2 =a[j1+1    ];	    rt =a[j1+p1 ];	it =a[j1+1+p1 ];	    t3 =t1 -rt;		t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ];	t6 =a[j1+1+p2 ];	    rt =a[j1+p3 ];	it =a[j1+1+p3 ];	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ];	t10=a[j1+1+p4 ];	    rt =a[j1+p5 ];	it =a[j1+1+p5 ];	    t11=t9 -rt;		t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ];	t14=a[j1+1+p6 ];	    rt =a[j1+p7 ];	it =a[j1+1+p7 ];	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ];	t18=a[j1+1+p8 ];	    rt =a[j1+p9 ];	it =a[j1+1+p9 ];	    t19=t17-rt;		t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10];	t22=a[j1+1+p10];	    rt =a[j1+p11];	it =a[j1+1+p11];	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12];	t26=a[j1+1+p12];	    rt =a[j1+p13];	it =a[j1+1+p13];	    t27=t25-rt;		t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14];	t30=a[j1+1+p14];	    rt =a[j1+p15];	it =a[j1+1+p15];	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#else	    l = iroot;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l += iroot;			/* 2*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c1=t1*t3-t2*t4;	s1=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 4*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c2=t1*t3-t2*t4;	s2=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 8*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c4=t1*t3-t2*t4;	s4=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+(iroot << 2)+iroot;		/* 13*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c8=t1*t3-t2*t4;	s8=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c13=t1*t3-t2*t4;	s13=t1*t4+t2*t3;	      t1=c1*c4; t2=c1*s4; t3=s1*c4; t4=s1*s4;	      c3=t1+t4; s3=t2-t3; c5=t1-t4; s5=t2+t3;	      t1=c1*c8; t2=c1*s8; t3=s1*c8; t4=s1*s8;	      c7=t1+t4; s7=t2-t3; c9=t1-t4; s9=t2+t3;	      t1=c1*c13; t2=c1*s13; t3=s1*c13; t4=s1*s13;	      c12=t1+t4; s12=t2-t3; c14=t1-t4; s14=t2+t3;	      t1=c2*c8; t2=c2*s8; t3=s2*c8; t4=s2*s8;	      c6=t1+t4; s6=t2-t3; c10=t1-t4; s10=t2+t3;	      t1=c2*c13; t2=c2*s13; t3=s2*c13; t4=s2*s13;	      c11=t1+t4; s11=t2-t3; c15=t1-t4; s15=t2+t3;/*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*//*...Block 1:	*/	    t1 =a[j1    ];			t2 =a[j1+1    ];	    rt =a[j1+p1 ]*c8 +a[j1+1+p1 ]*s8;	it =a[j1+1+p1 ]*c8 -a[j1+p1 ]*s8;	    t3 =t1 -rt;  	t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ]*c4 +a[j1+1+p2 ]*s4;	t6 =a[j1+1+p2 ]*c4 -a[j1+p2 ]*s4;	    rt =a[j1+p3 ]*c12+a[j1+1+p3 ]*s12;	it =a[j1+1+p3 ]*c12-a[j1+p3 ]*s12;	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ]*c2 +a[j1+1+p4 ]*s2 ;	t10=a[j1+1+p4 ]*c2 -a[j1+p4 ]*s2;	    rt =a[j1+p5 ]*c10+a[j1+1+p5 ]*s10;	it =a[j1+1+p5 ]*c10-a[j1+p5 ]*s10;	    t11=t9 -rt;  	t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ]*c6 +a[j1+1+p6 ]*s6;	t14=a[j1+1+p6 ]*c6 -a[j1+p6 ]*s6;	    rt =a[j1+p7 ]*c14+a[j1+1+p7 ]*s14;	it =a[j1+1+p7 ]*c14-a[j1+p7 ]*s14;	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ]*c1 +a[j1+1+p8 ]*s1;	t18=a[j1+1+p8 ]*c1 -a[j1+p8 ]*s1;	    rt =a[j1+p9 ]*c9 +a[j1+1+p9 ]*s9;	it =a[j1+1+p9 ]*c9 -a[j1+p9 ]*s9;	    t19=t17-rt;  	t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10]*c5 +a[j1+1+p10]*s5 ;	t22=a[j1+1+p10]*c5 -a[j1+p10]*s5;	    rt =a[j1+p11]*c13+a[j1+1+p11]*s13;	it =a[j1+1+p11]*c13-a[j1+p11]*s13;	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12]*c3 +a[j1+1+p12]*s3 ;	t26=a[j1+1+p12]*c3 -a[j1+p12]*s3;	    rt =a[j1+p13]*c11+a[j1+1+p13]*s11;	it =a[j1+1+p13]*c11-a[j1+p13]*s11;	    t27=t25-rt;  	t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14]*c7 +a[j1+1+p14]*s7 ;	t30=a[j1+1+p14]*c7 -a[j1+p14]*s7;	    rt =a[j1+p15]*c15+a[j1+1+p15]*s15;	it =a[j1+1+p15]*c15-a[j1+p15]*s15;	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(-i* 1*twopi/16) =       ( c,-s), exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 3*twopi/16) =       ( s,-c) (for inputs to transform block 2)  	1, exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 4*twopi/16) =       ( 0,-1), exp(-i* 6*twopi/16) = ISRT2*(-1,-1) (for inputs to transform block 3)  	1, exp(-i* 3*twopi/16) =       ( s,-c), exp(-i* 6*twopi/16) = ISRT2*(-1,-1), exp(-i* 9*twopi/16) =       (-c, s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;	    a1p0r =t1+t17;	a1p0i =t2+t18;	    a1p8r =t1-t17;	a1p8i =t2-t18;	    a1p4r =t9 +t26;	a1p4i =t10-t25;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pCr=t9 -t26;	a1pCi=t10+t25;/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 -t14;	t5 =t5 +t14;		/* twiddle mpy by E^4 =-I	*/			t14=t6 +rt;	t6 =t6 -rt;	    rt =(t22+t21)*ISRT2;t22=(t22-t21)*ISRT2;	t21=rt;	/* twiddle mpy by E^-2	*/	    rt =(t29-t30)*ISRT2;it =(t29+t30)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^-6 by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;	    a1p2r =t5+t21;	a1p2i =t6+t22;	    a1pAr=t5-t21;	a1pAi=t6-t22;	    a1p6r =t13+t30;	a1p6i =t14-t29;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pEr=t13-t30;	a1pEi=t14+t29;/*...Block 2: t3,11,19,27	*/	    rt =(t12+t11)*ISRT2;it =(t12-t11)*ISRT2;		/* twiddle mpy by E^-2	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c + t20*s;	t20=t20*c - t19*s;	t19=rt;	/* twiddle mpy by E^-1	*/	    rt =t27*s + t28*c;	it =t28*s - t27*c;		/* twiddle mpy by E^-3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;	    a1p1r =t3+t19;	a1p1i =t4+t20;	    a1p9r =t3-t19;	a1p9i =t4-t20;	    a1p5r =t11+t28;	a1p5i =t12-t27;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pDr=t11-t28;	a1pDi=t12+t27;/*...Block 4: t7,15,23,31	*/	    rt =(t15-t16)*ISRT2;it =(t15+t16)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s + t24*c;	t24=t24*s - t23*c;	t23=rt;	/* twiddle mpy by E^-3	*/	    rt =t31*c + t32*s;	it =t32*c - t31*s;		/* twiddle mpy by E^-1 = -E^-9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;	    a1p3r =t7+t23;	a1p3i =t8+t24;	    a1pBr=t7-t23;	a1pBi=t8-t24;	    a1p7r =t15+t32;	a1p7i =t16-t31;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pFr=t15-t32;	a1pFi=t16+t31;/*...and combine those to complete the radix-16 transform and do the carries. Since the outputs would     normally be getting dispatched to 16 separate blocks of the A-array, we need 16 separate carries.	*/	    l= j & (nwt-1);			/* We want (S*J mod N) - SI(L) for all 16 carries, so precompute	*/	    n_minus_sil   = n-si[l  ];		/* N - SI(L) and for each J, find N - (B*J mod N) - SI(L)		*/	    n_minus_silp1 = n-si[l+1];		/* For the inverse weight, want (S*(N - J) mod N) - SI(NWT - L) =	*/	    sinwt   = si[nwt-l  ];		/*	= N - (S*J mod N) - SI(NWT - L) = (B*J mod N) - SI(NWT - L).	*/	    sinwtm1 = si[nwt-l-1];	    wtl     =wt0[    l  ];	    wtn     =wt0[nwt-l  ]*scale;	/* Include 1/(n/2) scale factor of inverse transform here...	*/	    wtlp1   =wt0[    l+1];	    wtnm1   =wt0[nwt-l-1]*scale;	/* ...and here.	*//*...set0 is slightly different from others:	*/	col = _col[ithread];	co2 = _co2[ithread];	co3 = _co3[ithread];	   cmplx_carry_norm_pow2_nocheck0(a1p0r,a1p0i,cy0[ithread],bjmodn0[ithread]);	    cmplx_carry_norm_pow2_nocheck(a1p1r,a1p1i,cy1[ithread],bjmodn1[ithread],0x1);	    cmplx_carry_norm_pow2_nocheck(a1p2r,a1p2i,cy2[ithread],bjmodn2[ithread],0x2);	    cmplx_carry_norm_pow2_nocheck(a1p3r,a1p3i,cy3[ithread],bjmodn3[ithread],0x3);	    cmplx_carry_norm_pow2_nocheck(a1p4r,a1p4i,cy4[ithread],bjmodn4[ithread],0x4);	    cmplx_carry_norm_pow2_nocheck(a1p5r,a1p5i,cy5[ithread],bjmodn5[ithread],0x5);	    cmplx_carry_norm_pow2_nocheck(a1p6r,a1p6i,cy6[ithread],bjmodn6[ithread],0x6);	    cmplx_carry_norm_pow2_nocheck(a1p7r,a1p7i,cy7[ithread],bjmodn7[ithread],0x7);	    cmplx_carry_norm_pow2_nocheck(a1p8r,a1p8i,cy8[ithread],bjmodn8[ithread],0x8);	    cmplx_carry_norm_pow2_nocheck(a1p9r,a1p9i,cy9[ithread],bjmodn9[ithread],0x9);	    cmplx_carry_norm_pow2_nocheck(a1pAr,a1pAi,cyA[ithread],bjmodnA[ithread],0xA);	    cmplx_carry_norm_pow2_nocheck(a1pBr,a1pBi,cyB[ithread],bjmodnB[ithread],0xB);	    cmplx_carry_norm_pow2_nocheck(a1pCr,a1pCi,cyC[ithread],bjmodnC[ithread],0xC);	    cmplx_carry_norm_pow2_nocheck(a1pDr,a1pDi,cyD[ithread],bjmodnD[ithread],0xD);	    cmplx_carry_norm_pow2_nocheck(a1pEr,a1pEi,cyE[ithread],bjmodnE[ithread],0xE);	    cmplx_carry_norm_pow2_nocheck(a1pFr,a1pFi,cyF[ithread],bjmodnF[ithread],0xF);	    i =((uint32)(sw - bjmodn0[ithread]) >> 31);	/* get ready for the next set...	*/	    _co2[ithread] = _co3[ithread];	/* For all data but the first set in each j-block, co2=co3. Thus, after the first block of data is done						  (and only then: for all subsequent blocks it's superfluous), this assignment decrements co2 by radix(1).	*//*...The radix-16 DIF pass is here:	*//*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*/#if PFETCHadd0 = &a[j1];prefetch_p_doubles(add0);#endif/*...Block 1:	*/	    t3 =a1p0r -a1p8r;  	t1 =a1p0r +a1p8r;	    t4 =a1p0i -a1p8i;	t2 =a1p0i +a1p8i;	    t7 =a1p4r -a1pCr;  t5 =a1p4r +a1pCr;	    t8 =a1p4i -a1pCi;  t6 =a1p4i +a1pCi;#if PFETCHadd = add0+p1;prefetch_p_doubles(add);#endif	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 +t8;	t3 =t3 -t8;			t8 =t4 -rt;	t4 =t4 +rt;#if PFETCHadd = add0+p2;prefetch_p_doubles(add);#endif/*...Block 2:	*/	    t11=a1p2r -a1pAr;  t9 =a1p2r +a1pAr;	    t12=a1p2i -a1pAi;	t10=a1p2i +a1pAi;	    t15=a1p6r -a1pEr;  t13=a1p6r +a1pEr;	    t16=a1p6i -a1pEi;	t14=a1p6i +a1pEi;#if PFETCHadd = add0+p3;prefetch_p_doubles(add);#endif	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11+t16;	t11=t11-t16;			t16=t12-rt;	t12=t12+rt;#if PFETCHadd = add0+p4;prefetch_p_doubles(add);#endif/*...Block 3:	*/	    t19=a1p1r -a1p9r;  	t17=a1p1r +a1p9r;	    t20=a1p1i -a1p9i;	t18=a1p1i +a1p9i;	    t23=a1p5r -a1pDr;  t21=a1p5r +a1pDr;	    t24=a1p5i -a1pDi;	t22=a1p5i +a1pDi;#if PFETCHadd = add0+p5;prefetch_p_doubles(add);#endif	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19+t24;	t19=t19-t24;			t24=t20-rt;	t20=t20+rt;#if PFETCHadd = add0+p6;prefetch_p_doubles(add);#endif/*...Block 4:	*/	    t27=a1p3r -a1pBr;  t25=a1p3r +a1pBr;	    t28=a1p3i -a1pBi;	t26=a1p3i +a1pBi;	    t31=a1p7r -a1pFr;  t29=a1p7r +a1pFr;	    t32=a1p7i -a1pFi;	t30=a1p7i +a1pFi;#if PFETCHadd = add0+p7;prefetch_p_doubles(add);#endif	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27+t32;	t27=t27-t32;			t32=t28-rt;	t28=t28+rt;#if PFETCHadd = add0+p8;prefetch_p_doubles(add);#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(i* 1*twopi/16) =       ( c, s), exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 3*twopi/16) =       ( s, c) (for inputs to transform block 2)  	1, exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 4*twopi/16) =       ( 0, 1), exp(i* 6*twopi/16) = ISRT2*(-1, 1) (for inputs to transform block 3)  	1, exp(i* 3*twopi/16) =       ( s, c), exp(i* 6*twopi/16) = ISRT2*(-1, 1), exp(i* 9*twopi/16) =       (-c,-s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;#if PFETCHadd = add0+p9;prefetch_p_doubles(add);#endif	    a[j1    ]=t1+t17;	a[j1+1    ]=t2+t18;	    a[j1+p1 ]=t1-t17;	a[j1+1+p1 ]=t2-t18;	    a[j1+p2 ]=t9 -t26;	a[j1+1+p2 ]=t10+t25;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p3 ]=t9 +t26;	a[j1+1+p3 ]=t10-t25;#if PFETCHadd = add0+p10;prefetch_p_doubles(add);#endif/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 +t14;	t5 =t5 -t14;		/* twiddle mpy by E^4 = I	*/			t14=t6 -rt;	t6 =t6 +rt;	    rt =(t21-t22)*ISRT2;t22=(t21+t22)*ISRT2;	t21=rt;	/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t21,22,29,30	*/	    rt =(t30+t29)*ISRT2;it =(t30-t29)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;#if PFETCHadd = add0+p11;prefetch_p_doubles(add);#endif	    a[j1+p4 ]=t5+t21;	a[j1+1+p4 ]=t6+t22;	    a[j1+p5 ]=t5-t21;	a[j1+1+p5 ]=t6-t22;	    a[j1+p6 ]=t13-t30;	a[j1+1+p6 ]=t14+t29;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p7 ]=t13+t30;	a[j1+1+p7 ]=t14-t29;#if PFETCHadd = add0+p12;prefetch_p_doubles(add);#endif/*...Block 2: t3,11,19,27	*/	    rt =(t11-t12)*ISRT2;it =(t11+t12)*ISRT2;		/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t11,12	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c - t20*s;	t20=t20*c + t19*s;	t19=rt;	/* twiddle mpy by E^1	*/	    rt =t27*s - t28*c;	it =t28*s + t27*c;		/* twiddle mpy by E^3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;#if PFETCHadd = add0+p13;prefetch_p_doubles(add);#endif	    a[j1+p8 ]=t3+t19;	a[j1+1+p8 ]=t4+t20;	    a[j1+p9 ]=t3-t19;	a[j1+1+p9 ]=t4-t20;	    a[j1+p10]=t11-t28;	a[j1+1+p10]=t12+t27;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p11]=t11+t28;	a[j1+1+p11]=t12-t27;#if PFETCHadd = add0+p14;prefetch_p_doubles(add);#endif/*...Block 4: t7,15,23,31	*/	    rt =(t16+t15)*ISRT2;it =(t16-t15)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s - t24*c;	t24=t24*s + t23*c;	t23=rt;	/* twiddle mpy by E^3	*/	    rt =t31*c - t32*s;	it =t32*c + t31*s;		/* twiddle mpy by E^1 = -E^9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;			/* Note: t23+rt = t23*(s+1)	*/#if PFETCHadd = add0+p15;prefetch_p_doubles(add);#endif	    a[j1+p12]=t7+t23;	a[j1+1+p12]=t8+t24;	    a[j1+p13]=t7-t23;	a[j1+1+p13]=t8-t24;	    a[j1+p14]=t15-t32;	a[j1+1+p14]=t16+t31;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p15]=t15+t32;	a[j1+1+p15]=t16-t31;	  iroot += root_incr;		/* increment sincos index.	*/	  }	  jstart[ithread] += nwt;	  jhi   [ithread] += nwt;	  _col[ithread] += 16;	  _co3[ithread] -= 16;	}    }	/******* END OF PARALLEL FOR-LOOP ********/  }	if(root_incr==0)break;/*   Wraparound carry cleanup loop is here:     The cleanup carries from the end of each length-N/16 block into the begining of the next     can all be neatly processed as follows:     (1) Invert the radix-16 forward DIF FFT of the first block of 16 complex elements in A and unweight;     (2) Propagate cleanup carries among the real and imaginary parts of the 16 outputs of (1);     (3) Reweight and perform a radix-16 forward DIF FFT on the result of (2);     (4) If any of the exit carries from (2) are nonzero, advance to the next 16 elements and repeat (1-4).*/	t1  = cyF[NTHREADS - 1];    if(NTHREADS > 1)    {      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyF[ithread] = cyF[ithread-1]; }; cyF[0] = cyE[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyE[ithread] = cyE[ithread-1]; }; cyE[0] = cyD[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyD[ithread] = cyD[ithread-1]; }; cyD[0] = cyC[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyC[ithread] = cyC[ithread-1]; }; cyC[0] = cyB[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyB[ithread] = cyB[ithread-1]; }; cyB[0] = cyA[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cyA[ithread] = cyA[ithread-1]; }; cyA[0] = cy9[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy9[ithread] = cy9[ithread-1]; }; cy9[0] = cy8[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy8[ithread] = cy8[ithread-1]; }; cy8[0] = cy7[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy7[ithread] = cy7[ithread-1]; }; cy7[0] = cy6[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy6[ithread] = cy6[ithread-1]; }; cy6[0] = cy5[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy5[ithread] = cy5[ithread-1]; }; cy5[0] = cy4[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy4[ithread] = cy4[ithread-1]; }; cy4[0] = cy3[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy3[ithread] = cy3[ithread-1]; }; cy3[0] = cy2[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy2[ithread] = cy2[ithread-1]; }; cy2[0] = cy1[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy1[ithread] = cy1[ithread-1]; }; cy1[0] = cy0[NTHREADS - 1];      for(ithread = NTHREADS - 1; ithread > 0; ithread--){ cy0[ithread] = cy0[ithread-1]; }; cy0[0] = t1;    }    else    {	cyF[0]= cyE[0];	cyE[0]= cyD[0];	cyD[0]= cyC[0];	cyC[0]= cyB[0];	cyB[0]= cyA[0];	cyA[0]= cy9[0];	cy9[0]= cy8[0];	cy8[0]= cy7[0];	cy7[0]= cy6[0];	cy6[0]= cy5[0];	cy5[0]= cy4[0];	cy4[0]= cy3[0];	cy3[0]= cy2[0];	cy2[0]= cy1[0];	cy1[0]= cy0[0];	cy0[0]= t1;    }	iroot = 0;	root_incr = 0;	scale = 1;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	for(j = ithread*pini; j <= ithread*pini + 3; j++)	{	  a[j    ] *= radix_inv;	  a[j+p1 ] *= radix_inv;	  a[j+p2 ] *= radix_inv;	  a[j+p3 ] *= radix_inv;	  a[j+p4 ] *= radix_inv;	  a[j+p5 ] *= radix_inv;	  a[j+p6 ] *= radix_inv;	  a[j+p7 ] *= radix_inv;	  a[j+p8 ] *= radix_inv;	  a[j+p9 ] *= radix_inv;	  a[j+p10] *= radix_inv;	  a[j+p11] *= radix_inv;	  a[j+p12] *= radix_inv;	  a[j+p13] *= radix_inv;	  a[j+p14] *= radix_inv;	  a[j+p15] *= radix_inv;	}    }}    t1 = 0;    for(ithread = 0; ithread < NTHREADS; ithread++)    {	t1 += fabs(cy0[ithread])+fabs(cy1[ithread])+fabs(cy2[ithread])+fabs(cy3[ithread])+fabs(cy4[ithread])+fabs(cy5[ithread])+fabs(cy6[ithread])+fabs(cy7[ithread])		+fabs(cy8[ithread])+fabs(cy9[ithread])+fabs(cyA[ithread])+fabs(cyB[ithread])+fabs(cyC[ithread])+fabs(cyD[ithread])+fabs(cyE[ithread])+fabs(cyF[ithread]);    }	if(t1 != 0.0)	{	    if(INTERACT)printf("FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    f1 = fopen(   OFILE,"a");	    f2 = fopen(STATFILE,"a");	    fprintf(f1,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    fprintf(f2,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1.\n",iter);	    fclose(f1);	    fclose(f2);	    err=1;	    return(err);	}	return(0);#else	// **************** SINGLE-THREADED VERSION ********************	int n16,bjmodn0,bjmodn1,bjmodn2,bjmodn3,bjmodn4,bjmodn5,bjmodn6,bjmodn7,bjmodn8,bjmodn9,bjmodnA	,bjmodnB,bjmodnC,bjmodnD,bjmodnE,bjmodnF,i,j,j1,jstart,jhi,iroot,root_incr,k1,k2,k,khi,l,outer;	static int psave=0,bw,sw,bjmodnini,nm1,nrtm1,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15;	const double one_half[3] = {1.0, 0.5, 0.25};	/* Needed for small-weights-tables scheme */	static double c = 0.9238795325112867561281831, s = 0.3826834323650897717284599, radix_inv,n2inv;	double rt,it#if !USE_NEW_DIT		,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15		,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15#endif	,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32	,a1p0r,a1p1r,a1p2r,a1p3r,a1p4r,a1p5r,a1p6r,a1p7r,a1p8r,a1p9r,a1pAr,a1pBr,a1pCr,a1pDr,a1pEr,a1pFr	,a1p0i,a1p1i,a1p2i,a1p3i,a1p4i,a1p5i,a1p6i,a1p7i,a1p8i,a1p9i,a1pAi,a1pBi,a1pCi,a1pDi,a1pEi,a1pFi		,cy0,cy1,cy2,cy3,cy4,cy5,cy6,cy7,cy8,cy9,cyA,cyB,cyC,cyD,cyE,cyF,temp,scale;#if PFETCH	double *add0, *add;#endif	int err;	FILE *f1,*f2;	static int first_entry=TRUE;/*...stuff for the reduced-length DWT weights array is here:	*/	int n_div_nwt;	int col,co2,co3,m,m2,n_minus_sil,n_minus_silp1,sinwt,sinwtm1;	double wt,wtinv,wtl,wtlp1,wtn,wtnm1,wtA,wtB,wtC;/*...change n16 and n_div_wt to non-static to work around a gcc compiler bug. */	n16   = n/16;	n_div_nwt = n16 >> nwt_bits;	if((n_div_nwt << nwt_bits) != n16)	{	  if(INTERACT)printf("FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  f1 = fopen(   OFILE,"a");	  f2 = fopen(STATFILE,"a");	  fprintf(f1,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fprintf(f2,"FATAL: iter = %10d; NWT_BITS does not divide N/16 in radix16_ditN_cy_dif1.\n",iter);	  fclose(f1);	  fclose(f2);	  err=1;	  return(err);	}	if(p != psave)	{	  first_entry=TRUE;	}/*...initialize things upon first entry	*/	if(first_entry)	{	  psave = p;	  first_entry=FALSE;	  radix_inv = (double)((long double)1.0L/(long double)16.0L);	  n2inv     = (double)((long double)1.0L/(long double)(n/2));	  bw    = p%n;		/* Number of bigwords in the Crandall/Fagin mixed-radix representation = (Mersenne exponent) mod (vector length).	*/	  sw    = n - bw;	/* Number of smallwords.	*/	  nm1   = n-1;	  nrtm1 = nrt - 1;/*   constant index offsets for load/stores are here.	*/	  p1 = n16 + ( (n16 >> DAT_BITS) << PAD_BITS );	  p2 = p1 +p1;	  p3 = p2 +p1;	  p4 = p3 +p1;	  p5 = p4 +p1;	  p6 = p5 +p1;	  p7 = p6 +p1;	  p8 = p7 +p1;	  p9 = p8 +p1;	  p10= p9 +p1;	  p11= p10+p1;	  p12= p11+p1;	  p13= p12+p1;	  p14= p13+p1;	  p15= p14+p1;	  bjmodnini=0;	  for(j=0; j < n16; j++)	  {	    bjmodnini -= sw; bjmodnini = bjmodnini + ( (-(int)((uint32)bjmodnini >> 31)) & n);	  }	}/*...The radix-16 final DIT pass is here.	*/	cy0 =-2;		/* init carries	*/	cy1 = 0;	cy2 = 0;	cy3 = 0;	cy4 = 0;	cy5 = 0;	cy6 = 0;	cy7 = 0;	cy8 = 0;	cy9 = 0;	cyA= 0;	cyB= 0;	cyC= 0;	cyD= 0;	cyE= 0;	cyF= 0;	iroot = 0;	/* init sincos array index	*/	root_incr = 1;	/* init sincos array index increment (set = 1 for normal carry pass, = 0 for wrapper pass)	*/	scale = n2inv;	/* init inverse-weight scale factor  (set = 2/n for normal carry pass, = 1 for wrapper pass)	*/	jstart = 0;	jhi = jstart+nwt-1;	khi = n_div_nwt;for(outer=0; outer <= 1; outer++){	i = 1;		/* Pointer to the BASE and BASEINV arrays. lowest-order digit is always a bigword (i = 1).	*/	bjmodn0 = 0;	bjmodn1 = bjmodnini;	bjmodn2 = bjmodn1 +bjmodnini-n; bjmodn2 = bjmodn2 + ( (-(int)((uint32)bjmodn2 >> 31)) & n);	bjmodn3 = bjmodn2 +bjmodnini-n; bjmodn3 = bjmodn3 + ( (-(int)((uint32)bjmodn3 >> 31)) & n);	bjmodn4 = bjmodn3 +bjmodnini-n; bjmodn4 = bjmodn4 + ( (-(int)((uint32)bjmodn4 >> 31)) & n);	bjmodn5 = bjmodn4 +bjmodnini-n; bjmodn5 = bjmodn5 + ( (-(int)((uint32)bjmodn5 >> 31)) & n);	bjmodn6 = bjmodn5 +bjmodnini-n; bjmodn6 = bjmodn6 + ( (-(int)((uint32)bjmodn6 >> 31)) & n);	bjmodn7 = bjmodn6 +bjmodnini-n; bjmodn7 = bjmodn7 + ( (-(int)((uint32)bjmodn7 >> 31)) & n);	bjmodn8 = bjmodn7 +bjmodnini-n; bjmodn8 = bjmodn8 + ( (-(int)((uint32)bjmodn8 >> 31)) & n);	bjmodn9 = bjmodn8 +bjmodnini-n; bjmodn9 = bjmodn9 + ( (-(int)((uint32)bjmodn9 >> 31)) & n);	bjmodnA = bjmodn9 +bjmodnini-n; bjmodnA = bjmodnA + ( (-(int)((uint32)bjmodnA >> 31)) & n);	bjmodnB = bjmodnA +bjmodnini-n; bjmodnB = bjmodnB + ( (-(int)((uint32)bjmodnB >> 31)) & n);	bjmodnC = bjmodnB +bjmodnini-n; bjmodnC = bjmodnC + ( (-(int)((uint32)bjmodnC >> 31)) & n);	bjmodnD = bjmodnC +bjmodnini-n; bjmodnD = bjmodnD + ( (-(int)((uint32)bjmodnD >> 31)) & n);	bjmodnE = bjmodnD +bjmodnini-n; bjmodnE = bjmodnE + ( (-(int)((uint32)bjmodnE >> 31)) & n);	bjmodnF = bjmodnE +bjmodnini-n; bjmodnF = bjmodnF + ( (-(int)((uint32)bjmodnF >> 31)) & n);	col=0;	co2=(n >> nwt_bits)-1+16;	co3=co2-16;		/* At the start of each new j-loop, co3=co2-radix(1)	*/	for(k=1; k <= khi; k++)	/* Do n/(radix(1)*nwt) outer loop executions...	*/	{	  for(j=jstart; j<jhi; j += 2)	/* Each inner loop execution processes (radix(1)*nwt) array data.	*/	  {	    j1 = j + ( (j >> DAT_BITS) << PAD_BITS );	/* padded-array fetch index is here	*//*	here are the needed sincos data - these will be processed below in bit-reversed order.	*/#if USE_NEW_DIT/*...Block 1:	*/	    t1 =a[j1    ];	t2 =a[j1+1    ];	    rt =a[j1+p1 ];	it =a[j1+1+p1 ];	    t3 =t1 -rt;		t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ];	t6 =a[j1+1+p2 ];	    rt =a[j1+p3 ];	it =a[j1+1+p3 ];	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ];	t10=a[j1+1+p4 ];	    rt =a[j1+p5 ];	it =a[j1+1+p5 ];	    t11=t9 -rt;		t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ];	t14=a[j1+1+p6 ];	    rt =a[j1+p7 ];	it =a[j1+1+p7 ];	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ];	t18=a[j1+1+p8 ];	    rt =a[j1+p9 ];	it =a[j1+1+p9 ];	    t19=t17-rt;		t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10];	t22=a[j1+1+p10];	    rt =a[j1+p11];	it =a[j1+1+p11];	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12];	t26=a[j1+1+p12];	    rt =a[j1+p13];	it =a[j1+1+p13];	    t27=t25-rt;		t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14];	t30=a[j1+1+p14];	    rt =a[j1+p15];	it =a[j1+1+p15];	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#else	    l = iroot;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l += iroot;				/* 2*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c1=t1*t3-t2*t4;	s1=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 4*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c2=t1*t3-t2*t4;	s2=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+l;				/* 8*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c4=t1*t3-t2*t4;	s4=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      l = l+(iroot << 2)+iroot;		/* 13*iroot	*/	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c8=t1*t3-t2*t4;	s8=t1*t4+t2*t3;	      k1=(l & nrtm1);	      k2=(l >> nrt_bits);	      t1=rt0[k1].re;	t2=rt0[k1].im;	      t3=rt1[k2].re;	t4=rt1[k2].im;	      c13=t1*t3-t2*t4;	s13=t1*t4+t2*t3;	      t1=c1*c4; t2=c1*s4; t3=s1*c4; t4=s1*s4;	      c3=t1+t4; s3=t2-t3; c5=t1-t4; s5=t2+t3;	      t1=c1*c8; t2=c1*s8; t3=s1*c8; t4=s1*s8;	      c7=t1+t4; s7=t2-t3; c9=t1-t4; s9=t2+t3;	      t1=c1*c13; t2=c1*s13; t3=s1*c13; t4=s1*s13;	      c12=t1+t4; s12=t2-t3; c14=t1-t4; s14=t2+t3;	      t1=c2*c8; t2=c2*s8; t3=s2*c8; t4=s2*s8;	      c6=t1+t4; s6=t2-t3; c10=t1-t4; s10=t2+t3;	      t1=c2*c13; t2=c2*s13; t3=s2*c13; t4=s2*s13;	      c11=t1+t4; s11=t2-t3; c15=t1-t4; s15=t2+t3;/*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*//*...Block 1:	*/	    t1 =a[j1    ];			t2 =a[j1+1    ];	    rt =a[j1+p1 ]*c8 +a[j1+1+p1 ]*s8;	it =a[j1+1+p1 ]*c8 -a[j1+p1 ]*s8;	    t3 =t1 -rt;  	t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ]*c4 +a[j1+1+p2 ]*s4;	t6 =a[j1+1+p2 ]*c4 -a[j1+p2 ]*s4;	    rt =a[j1+p3 ]*c12+a[j1+1+p3 ]*s12;	it =a[j1+1+p3 ]*c12-a[j1+p3 ]*s12;	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ]*c2 +a[j1+1+p4 ]*s2 ;	t10=a[j1+1+p4 ]*c2 -a[j1+p4 ]*s2;	    rt =a[j1+p5 ]*c10+a[j1+1+p5 ]*s10;	it =a[j1+1+p5 ]*c10-a[j1+p5 ]*s10;	    t11=t9 -rt;  	t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ]*c6 +a[j1+1+p6 ]*s6;	t14=a[j1+1+p6 ]*c6 -a[j1+p6 ]*s6;	    rt =a[j1+p7 ]*c14+a[j1+1+p7 ]*s14;	it =a[j1+1+p7 ]*c14-a[j1+p7 ]*s14;	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ]*c1 +a[j1+1+p8 ]*s1;	t18=a[j1+1+p8 ]*c1 -a[j1+p8 ]*s1;	    rt =a[j1+p9 ]*c9 +a[j1+1+p9 ]*s9;	it =a[j1+1+p9 ]*c9 -a[j1+p9 ]*s9;	    t19=t17-rt;  	t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10]*c5 +a[j1+1+p10]*s5 ;	t22=a[j1+1+p10]*c5 -a[j1+p10]*s5;	    rt =a[j1+p11]*c13+a[j1+1+p11]*s13;	it =a[j1+1+p11]*c13-a[j1+p11]*s13;	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12]*c3 +a[j1+1+p12]*s3 ;	t26=a[j1+1+p12]*c3 -a[j1+p12]*s3;	    rt =a[j1+p13]*c11+a[j1+1+p13]*s11;	it =a[j1+1+p13]*c11-a[j1+p13]*s11;	    t27=t25-rt;  	t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14]*c7 +a[j1+1+p14]*s7 ;	t30=a[j1+1+p14]*c7 -a[j1+p14]*s7;	    rt =a[j1+p15]*c15+a[j1+1+p15]*s15;	it =a[j1+1+p15]*c15-a[j1+p15]*s15;	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(-i* 1*twopi/16) =       ( c,-s), exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 3*twopi/16) =       ( s,-c) (for inputs to transform block 2)  	1, exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 4*twopi/16) =       ( 0,-1), exp(-i* 6*twopi/16) = ISRT2*(-1,-1) (for inputs to transform block 3)  	1, exp(-i* 3*twopi/16) =       ( s,-c), exp(-i* 6*twopi/16) = ISRT2*(-1,-1), exp(-i* 9*twopi/16) =       (-c, s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;	    a1p0r =t1+t17;	a1p0i =t2+t18;	    a1p8r =t1-t17;	a1p8i =t2-t18;	    a1p4r =t9 +t26;	a1p4i =t10-t25;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pCr=t9 -t26;	a1pCi=t10+t25;/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 -t14;	t5 =t5 +t14;		/* twiddle mpy by E^4 =-I	*/			t14=t6 +rt;	t6 =t6 -rt;	    rt =(t22+t21)*ISRT2;t22=(t22-t21)*ISRT2;	t21=rt;	/* twiddle mpy by E^-2	*/	    rt =(t29-t30)*ISRT2;it =(t29+t30)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^-6 by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;	    a1p2r =t5+t21;	a1p2i =t6+t22;	    a1pAr=t5-t21;	a1pAi=t6-t22;	    a1p6r =t13+t30;	a1p6i =t14-t29;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pEr=t13-t30;	a1pEi=t14+t29;/*...Block 2: t3,11,19,27	*/	    rt =(t12+t11)*ISRT2;it =(t12-t11)*ISRT2;		/* twiddle mpy by E^-2	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c + t20*s;	t20=t20*c - t19*s;	t19=rt;	/* twiddle mpy by E^-1	*/	    rt =t27*s + t28*c;	it =t28*s - t27*c;		/* twiddle mpy by E^-3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;	    a1p1r =t3+t19;	a1p1i =t4+t20;	    a1p9r =t3-t19;	a1p9i =t4-t20;	    a1p5r =t11+t28;	a1p5i =t12-t27;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pDr=t11-t28;	a1pDi=t12+t27;/*...Block 4: t7,15,23,31	*/	    rt =(t15-t16)*ISRT2;it =(t15+t16)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s + t24*c;	t24=t24*s - t23*c;	t23=rt;	/* twiddle mpy by E^-3	*/	    rt =t31*c + t32*s;	it =t32*c - t31*s;		/* twiddle mpy by E^-1 = -E^-9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;	    a1p3r =t7+t23;	a1p3i =t8+t24;	    a1pBr=t7-t23;	a1pBi=t8-t24;	    a1p7r =t15+t32;	a1p7i =t16-t31;	/* mpy by E^-4 = -I is inlined here...	*/	    a1pFr=t15-t32;	a1pFi=t16+t31;/*...and combine those to complete the radix-16 transform and do the carries. Since the outputs would     normally be getting dispatched to 16 separate blocks of the A-array, we need 16 separate carries.	*/	    l= j & (nwt-1);			/* We want (S*J mod N) - SI(L) for all 16 carries, so precompute	*/	    n_minus_sil   = n-si[l  ];		/* N - SI(L) and for each J, find N - (B*J mod N) - SI(L)		*/	    n_minus_silp1 = n-si[l+1];		/* For the inverse weight, want (S*(N - J) mod N) - SI(NWT - L) =	*/	    sinwt   = si[nwt-l  ];		/*	= N - (S*J mod N) - SI(NWT - L) = (B*J mod N) - SI(NWT - L).	*/	    sinwtm1 = si[nwt-l-1];	    wtl     =wt0[    l  ];	    wtn     =wt0[nwt-l  ]*scale;	/* Include 1/(n/2) scale factor of inverse transform here...	*/	    wtlp1   =wt0[    l+1];	    wtnm1   =wt0[nwt-l-1]*scale;	/* ...and here.	*//*...set0 is slightly different from others:	*/	   cmplx_carry_norm_pow2_nocheck0(a1p0r,a1p0i,cy0,bjmodn0);	    cmplx_carry_norm_pow2_nocheck(a1p1r,a1p1i,cy1,bjmodn1,0x1);	    cmplx_carry_norm_pow2_nocheck(a1p2r,a1p2i,cy2,bjmodn2,0x2);	    cmplx_carry_norm_pow2_nocheck(a1p3r,a1p3i,cy3,bjmodn3,0x3);	    cmplx_carry_norm_pow2_nocheck(a1p4r,a1p4i,cy4,bjmodn4,0x4);	    cmplx_carry_norm_pow2_nocheck(a1p5r,a1p5i,cy5,bjmodn5,0x5);	    cmplx_carry_norm_pow2_nocheck(a1p6r,a1p6i,cy6,bjmodn6,0x6);	    cmplx_carry_norm_pow2_nocheck(a1p7r,a1p7i,cy7,bjmodn7,0x7);	    cmplx_carry_norm_pow2_nocheck(a1p8r,a1p8i,cy8,bjmodn8,0x8);	    cmplx_carry_norm_pow2_nocheck(a1p9r,a1p9i,cy9,bjmodn9,0x9);	    cmplx_carry_norm_pow2_nocheck(a1pAr,a1pAi,cyA,bjmodnA,0xA);	    cmplx_carry_norm_pow2_nocheck(a1pBr,a1pBi,cyB,bjmodnB,0xB);	    cmplx_carry_norm_pow2_nocheck(a1pCr,a1pCi,cyC,bjmodnC,0xC);	    cmplx_carry_norm_pow2_nocheck(a1pDr,a1pDi,cyD,bjmodnD,0xD);	    cmplx_carry_norm_pow2_nocheck(a1pEr,a1pEi,cyE,bjmodnE,0xE);	    cmplx_carry_norm_pow2_nocheck(a1pFr,a1pFi,cyF,bjmodnF,0xF);	    i =((uint32)(sw - bjmodn0) >> 31);	/* get ready for the next set...	*/	    co2=co3;	/* For all data but the first set in each j-block, co2=co3. Thus, after the first block of data is done			   (and only then: for all subsequent blocks it's superfluous), this assignment decrements co2 by radix(1).	*//*...The radix-16 DIF pass is here:	*//*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*/#if PFETCHadd0 = &a[j1];prefetch_p_doubles(add0);#endif/*...Block 1:	*/	    t3 =a1p0r -a1p8r;  	t1 =a1p0r +a1p8r;	    t4 =a1p0i -a1p8i;	t2 =a1p0i +a1p8i;	    t7 =a1p4r -a1pCr;  t5 =a1p4r +a1pCr;	    t8 =a1p4i -a1pCi;  t6 =a1p4i +a1pCi;#if PFETCHadd = add0+p1;prefetch_p_doubles(add);#endif	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 +t8;	t3 =t3 -t8;			t8 =t4 -rt;	t4 =t4 +rt;#if PFETCHadd = add0+p2;prefetch_p_doubles(add);#endif/*...Block 2:	*/	    t11=a1p2r -a1pAr;  t9 =a1p2r +a1pAr;	    t12=a1p2i -a1pAi;	t10=a1p2i +a1pAi;	    t15=a1p6r -a1pEr;  t13=a1p6r +a1pEr;	    t16=a1p6i -a1pEi;	t14=a1p6i +a1pEi;#if PFETCHadd = add0+p3;prefetch_p_doubles(add);#endif	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11+t16;	t11=t11-t16;			t16=t12-rt;	t12=t12+rt;#if PFETCHadd = add0+p4;prefetch_p_doubles(add);#endif/*...Block 3:	*/	    t19=a1p1r -a1p9r;  	t17=a1p1r +a1p9r;	    t20=a1p1i -a1p9i;	t18=a1p1i +a1p9i;	    t23=a1p5r -a1pDr;  t21=a1p5r +a1pDr;	    t24=a1p5i -a1pDi;	t22=a1p5i +a1pDi;#if PFETCHadd = add0+p5;prefetch_p_doubles(add);#endif	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19+t24;	t19=t19-t24;			t24=t20-rt;	t20=t20+rt;#if PFETCHadd = add0+p6;prefetch_p_doubles(add);#endif/*...Block 4:	*/	    t27=a1p3r -a1pBr;  t25=a1p3r +a1pBr;	    t28=a1p3i -a1pBi;	t26=a1p3i +a1pBi;	    t31=a1p7r -a1pFr;  t29=a1p7r +a1pFr;	    t32=a1p7i -a1pFi;	t30=a1p7i +a1pFi;#if PFETCHadd = add0+p7;prefetch_p_doubles(add);#endif	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27+t32;	t27=t27-t32;			t32=t28-rt;	t28=t28+rt;#if PFETCHadd = add0+p8;prefetch_p_doubles(add);#endif/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(i* 1*twopi/16) =       ( c, s), exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 3*twopi/16) =       ( s, c) (for inputs to transform block 2)  	1, exp(i* 2*twopi/16) = ISRT2*( 1, 1), exp(i* 4*twopi/16) =       ( 0, 1), exp(i* 6*twopi/16) = ISRT2*(-1, 1) (for inputs to transform block 3)  	1, exp(i* 3*twopi/16) =       ( s, c), exp(i* 6*twopi/16) = ISRT2*(-1, 1), exp(i* 9*twopi/16) =       (-c,-s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;#if PFETCHadd = add0+p9;prefetch_p_doubles(add);#endif	    a[j1    ]=t1+t17;	a[j1+1    ]=t2+t18;	    a[j1+p1 ]=t1-t17;	a[j1+1+p1 ]=t2-t18;	    a[j1+p2 ]=t9 -t26;	a[j1+1+p2 ]=t10+t25;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p3 ]=t9 +t26;	a[j1+1+p3 ]=t10-t25;#if PFETCHadd = add0+p10;prefetch_p_doubles(add);#endif/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 +t14;	t5 =t5 -t14;		/* twiddle mpy by E^4 = I	*/			t14=t6 -rt;	t6 =t6 +rt;	    rt =(t21-t22)*ISRT2;t22=(t21+t22)*ISRT2;	t21=rt;	/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t21,22,29,30	*/	    rt =(t30+t29)*ISRT2;it =(t30-t29)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;#if PFETCHadd = add0+p11;prefetch_p_doubles(add);#endif	    a[j1+p4 ]=t5+t21;	a[j1+1+p4 ]=t6+t22;	    a[j1+p5 ]=t5-t21;	a[j1+1+p5 ]=t6-t22;	    a[j1+p6 ]=t13-t30;	a[j1+1+p6 ]=t14+t29;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p7 ]=t13+t30;	a[j1+1+p7 ]=t14-t29;#if PFETCHadd = add0+p12;prefetch_p_doubles(add);#endif/*...Block 2: t3,11,19,27	*/	    rt =(t11-t12)*ISRT2;it =(t11+t12)*ISRT2;		/* twiddle mpy by E^2; could also do mpy by ISRT2 directly on t11,12	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c - t20*s;	t20=t20*c + t19*s;	t19=rt;	/* twiddle mpy by E^1	*/	    rt =t27*s - t28*c;	it =t28*s + t27*c;		/* twiddle mpy by E^3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;#if PFETCHadd = add0+p13;prefetch_p_doubles(add);#endif	    a[j1+p8 ]=t3+t19;	a[j1+1+p8 ]=t4+t20;	    a[j1+p9 ]=t3-t19;	a[j1+1+p9 ]=t4-t20;	    a[j1+p10]=t11-t28;	a[j1+1+p10]=t12+t27;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p11]=t11+t28;	a[j1+1+p11]=t12-t27;#if PFETCHadd = add0+p14;prefetch_p_doubles(add);#endif/*...Block 4: t7,15,23,31	*/	    rt =(t16+t15)*ISRT2;it =(t16-t15)*ISRT2;		/* twiddle mpy by -E^6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s - t24*c;	t24=t24*s + t23*c;	t23=rt;	/* twiddle mpy by E^3	*/	    rt =t31*c - t32*s;	it =t32*c + t31*s;		/* twiddle mpy by E^1 = -E^9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;			/* Note: t23+rt = t23*(s+1)	*/#if PFETCHadd = add0+p15;prefetch_p_doubles(add);#endif	    a[j1+p12]=t7+t23;	a[j1+1+p12]=t8+t24;	    a[j1+p13]=t7-t23;	a[j1+1+p13]=t8-t24;	    a[j1+p14]=t15-t32;	a[j1+1+p14]=t16+t31;	/* mpy by E^4=i is inlined here...	*/	    a[j1+p15]=t15+t32;	a[j1+1+p15]=t16-t31;	  iroot += root_incr;		/* increment sincos index.	*/	  }	  jstart += nwt;	  jhi    += nwt;	  col += 16;	  co3 -= 16;	}	if(root_incr==0)break;/*   Wraparound carry cleanup loop is here:     The cleanup carries from the end of each length-N/16 block into the begining of the next     can all be neatly processed as follows:     (1) Invert the radix-16 forward DIF FFT of the first block of 16 complex elements in A and unweight;     (2) Propagate cleanup carries among the real and imaginary parts of the 16 outputs of (1);     (3) Reweight and perform a radix-16 forward DIF FFT on the result of (2);     (4) If any of the exit carries from (2) are nonzero, advance to the next 16 elements and repeat (1-4).*/	t1  = cyF;	cyF= cyE;	cyE= cyD;	cyD= cyC;	cyC= cyB;	cyB= cyA;	cyA= cy9;	cy9 = cy8;	cy8 = cy7;	cy7 = cy6;	cy6 = cy5;	cy5 = cy4;	cy4 = cy3;	cy3 = cy2;	cy2 = cy1;	cy1 = cy0;	cy0 = t1;	iroot = 0;	root_incr = 0;	scale = 1;	jstart = 0;	jhi = 3;	khi = 1;	for(j=0; j<=jhi; j++)	{	  a[j    ] *= radix_inv;	  a[j+p1 ] *= radix_inv;	  a[j+p2 ] *= radix_inv;	  a[j+p3 ] *= radix_inv;	  a[j+p4 ] *= radix_inv;	  a[j+p5 ] *= radix_inv;	  a[j+p6 ] *= radix_inv;	  a[j+p7 ] *= radix_inv;	  a[j+p8 ] *= radix_inv;	  a[j+p9 ] *= radix_inv;	  a[j+p10] *= radix_inv;	  a[j+p11] *= radix_inv;	  a[j+p12] *= radix_inv;	  a[j+p13] *= radix_inv;	  a[j+p14] *= radix_inv;	  a[j+p15] *= radix_inv;	}}	if(fabs(cy0)+fabs(cy1)+fabs(cy2)+fabs(cy3)+fabs(cy4)+fabs(cy5)+fabs(cy6)+fabs(cy7)		+fabs(cy8)+fabs(cy9)+fabs(cyA)+fabs(cyB)+fabs(cyC)+fabs(cyD)+fabs(cyE)+fabs(cyF) != 0.0)	{	    if(INTERACT)printf("FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1_nochk.\n",iter);	    f1 = fopen(   OFILE,"a");	    f2 = fopen(STATFILE,"a");	    fprintf(f1,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1_nochk.\n",iter);	    fprintf(f2,"FATAL: iter = %10d; nonzero exit carry in radix16_ditN_cy_dif1_nochk.\n",iter);	    fclose(f1);	    fclose(f2);	    err=1;	    return(err);	}	return(0);#endif	// **************** END MULTITHREAD #IFDEF ********************}/**************/void radix16_dit_pass1(double a[], int n){/*!...Acronym: DIT = Decimation In Time!!...Subroutine to perform an initial radix-16 complex inverse DIT FFT pass on the data in the length-N real vector A.!!   See the documentation in radix16_dif_pass for further details.*/	int j,j1;	static int n16,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15, first_entry=TRUE;	static double c = 0.9238795325112867561281831, s = 0.3826834323650897717284599;	/* exp[i*(twopi/16)]	*/	double rt,it	,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32;	if(!first_entry && (n >> 4) != n16)	/* New runlength?	*/	{	  first_entry=TRUE;	}/*...initialize things upon first entry	*/	if(first_entry)	{	  first_entry=FALSE;	  n16=n/16;	  p1 = n16 + ( (n16 >> DAT_BITS) << PAD_BITS );	  p2 = p1 +p1;	  p3 = p2 +p1;	  p4 = p3 +p1;	  p5 = p4 +p1;	  p6 = p5 +p1;	  p7 = p6 +p1;	  p8 = p7 +p1;	  p9 = p8 +p1;	  p10= p9 +p1;	  p11= p10+p1;	  p12= p11+p1;	  p13= p12+p1;	  p14= p13+p1;	  p15= p14+p1;	}/*...The radix-16 pass is here.	*/      for(j=0; j < n16; j += 2)      {	  j1 = j  + ( (j  >> DAT_BITS) << PAD_BITS );/*       gather the needed data (16 64-bit complex, i.e. 32 64-bit reals) and do the first set of four length-4 transforms.	*//*...Block 1:	*/	    t1 =a[j1    ];	t2 =a[j1+1    ];	    rt =a[j1+p1 ];	it =a[j1+1+p1 ];	    t3 =t1 -rt;  	t1 =t1 +rt;	    t4 =t2 -it;		t2 =t2 +it;	    t5 =a[j1+p2 ];	t6 =a[j1+1+p2 ];	    rt =a[j1+p3 ];	it =a[j1+1+p3 ];	    t7 =t5 -rt;  	t5 =t5 +rt;	    t8 =t6 -it;  	t6 =t6 +it;	    rt =t5;	t5 =t1 -rt;	t1 =t1 +rt;	    it =t6;	t6 =t2 -it;	t2 =t2 +it;	    rt =t7;	t7 =t3 -t8;	t3 =t3 +t8;			t8 =t4 +rt;	t4 =t4 -rt;/*...Block 2:	*/	    t9 =a[j1+p4 ];	t10=a[j1+1+p4 ];	    rt =a[j1+p5 ];	it =a[j1+1+p5 ];	    t11=t9 -rt;  	t9 =t9 +rt;	    t12=t10-it;		t10=t10+it;	    t13=a[j1+p6 ];	t14=a[j1+1+p6 ];	    rt =a[j1+p7 ];	it =a[j1+1+p7 ];	    t15=t13-rt;  	t13=t13+rt;	    t16=t14-it;		t14=t14+it;	    rt =t13;	t13=t9 -rt;	t9 =t9 +rt;	    it =t14;	t14=t10-it;	t10=t10+it;	    rt =t15;	t15=t11-t16;	t11=t11+t16;			t16=t12+rt;	t12=t12-rt;/*...Block 3:	*/	    t17=a[j1+p8 ];	t18=a[j1+1+p8 ];	    rt =a[j1+p9 ];	it =a[j1+1+p9 ];	    t19=t17-rt;  	t17=t17+rt;	    t20=t18-it;		t18=t18+it;	    t21=a[j1+p10];	t22=a[j1+1+p10];	    rt =a[j1+p11];	it =a[j1+1+p11];	    t23=t21-rt;  	t21=t21+rt;	    t24=t22-it;		t22=t22+it;	    rt =t21;	t21=t17-rt;	t17=t17+rt;	    it =t22;	t22=t18-it;	t18=t18+it;	    rt =t23;	t23=t19-t24;	t19=t19+t24;			t24=t20+rt;	t20=t20-rt;/*...Block 4:	*/	    t25=a[j1+p12];	t26=a[j1+1+p12];	    rt =a[j1+p13];	it =a[j1+1+p13];	    t27=t25-rt;  	t25=t25+rt;	    t28=t26-it;		t26=t26+it;	    t29=a[j1+p14];	t30=a[j1+1+p14];	    rt =a[j1+p15];	it =a[j1+1+p15];	    t31=t29-rt;  	t29=t29+rt;	    t32=t30-it;		t30=t30+it;	    rt =t29;	t29=t25-rt;	t25=t25+rt;	    it =t30;	t30=t26-it;	t26=t26+it;	    rt =t31;	t31=t27-t32;	t27=t27+t32;			t32=t28+rt;	t28=t28-rt;/*...and now do four more radix-4 transforms, including the internal twiddle factors:  	1, exp(-i* 1*twopi/16) =       ( c,-s), exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 3*twopi/16) =       ( s,-c) (for inputs to transform block 2)  	1, exp(-i* 2*twopi/16) = ISRT2*( 1,-1), exp(-i* 4*twopi/16) =       ( 0,-1), exp(-i* 6*twopi/16) = ISRT2*(-1,-1) (for inputs to transform block 3)  	1, exp(-i* 3*twopi/16) =       ( s,-c), exp(-i* 6*twopi/16) = ISRT2*(-1,-1), exp(-i* 9*twopi/16) =       (-c, s) (for inputs to transform block 4).    (This is 4 real*complex and 4 complex*complex multiplies (= 24 FMUL), compared to 6 and 4, respectively (= 28 FMUL) for my old scheme.)     I.e. do similar as above, except inputs a(j1  +p0:15:1) are replaced by t0:30:2,                                             a(j1+1+p0:15:1) are replaced by t1:31:2, and v.v. for outputs,     and only the last 3 inputs to each of the radix-4 transforms 2 through 4 are multiplied by non-unity twiddles.	*//*...Block 1: t1,9,17,25	*/	    rt =t9;	t9 =t1 -rt;	t1 =t1 +rt;	    it =t10;	t10=t2 -it;	t2 =t2 +it;	    rt =t25;	t25=t17-rt;	t17=t17+rt;	    it =t26;	t26=t18-it;	t18=t18+it;	    a[j1    ]=t1+t17;	a[j1+1    ]=t2+t18;	    a[j1+p8 ]=t1-t17;	a[j1+1+p8 ]=t2-t18;	    a[j1+p4 ]=t9 +t26;	a[j1+1+p4 ]=t10-t25;	/* mpy by E^-4 = -I is inlined here...	*/	    a[j1+p12]=t9 -t26;	a[j1+1+p12]=t10+t25;/*...Block 3: t5,13,21,29	*/	    rt =t13;	t13=t5 -t14;	t5 =t5 +t14;		/* twiddle mpy by E^4 =-I	*/			t14=t6 +rt;	t6 =t6 -rt;	    rt =(t22+t21)*ISRT2;t22=(t22-t21)*ISRT2;	t21=rt;	/* twiddle mpy by E^-2	*/	    rt =(t29-t30)*ISRT2;it =(t29+t30)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t29=t21+rt;		t21=t21-rt;			/* ...and get E^-6 by flipping signs here.	*/	    t30=t22+it;		t22=t22-it;	    a[j1+p2 ]=t5+t21;	a[j1+1+p2 ]=t6+t22;	    a[j1+p10]=t5-t21;	a[j1+1+p10]=t6-t22;	    a[j1+p6 ]=t13+t30;	a[j1+1+p6 ]=t14-t29;	/* mpy by E^-4 = -I is inlined here...	*/	    a[j1+p14]=t13-t30;	a[j1+1+p14]=t14+t29;/*...Block 2: t3,11,19,27	*/	    rt =(t12+t11)*ISRT2;it =(t12-t11)*ISRT2;		/* twiddle mpy by E^-2	*/	    t11=t3 -rt;		t3 =t3 +rt;	    t12=t4 -it;		t4 =t4 +it;	    rt =t19*c + t20*s;	t20=t20*c - t19*s;	t19=rt;	/* twiddle mpy by E^-1	*/	    rt =t27*s + t28*c;	it =t28*s - t27*c;		/* twiddle mpy by E^-3	*/	    t27=t19-rt;		t19=t19+rt;	    t28=t20-it;		t20=t20+it;	    a[j1+p1 ]=t3+t19;	a[j1+1+p1 ]=t4+t20;	    a[j1+p9 ]=t3-t19;	a[j1+1+p9 ]=t4-t20;	    a[j1+p5 ]=t11+t28;	a[j1+1+p5 ]=t12-t27;	/* mpy by E^-4 = -I is inlined here...	*/	    a[j1+p13]=t11-t28;	a[j1+1+p13]=t12+t27;/*...Block 4: t7,15,23,31	*/	    rt =(t15-t16)*ISRT2;it =(t15+t16)*ISRT2;		/* twiddle mpy by E^2 = -E^-6 is here...	*/	    t15=t7 +rt;		t7 =t7 -rt;			/* ...and get E^6=(i-1)/sqrt by flipping signs here.	*/	    t16=t8 +it;		t8 =t8 -it;	    rt =t23*s + t24*c;	t24=t24*s - t23*c;	t23=rt;	/* twiddle mpy by E^-3	*/	    rt =t31*c + t32*s;	it =t32*c - t31*s;		/* twiddle mpy by E^-1 = -E^-9...	*/	    t31=t23+rt;		t23=t23-rt;			/* ...and get E^9 by flipping signs here.	*/	    t32=t24+it;		t24=t24-it;	    a[j1+p3 ]=t7+t23;	a[j1+1+p3 ]=t8+t24;	    a[j1+p11]=t7-t23;	a[j1+1+p11]=t8-t24;	    a[j1+p7 ]=t15+t32;	a[j1+1+p7 ]=t16-t31;	/* mpy by E^-4 = -I is inlined here...	*/	    a[j1+p15]=t15-t32;	a[j1+1+p15]=t16+t31;	}}