telafft.c

用于Pattern Recongnition的一个图像处理程序－－Trace Transform。

       const double wa1[], const double wa2[])
{
    static const double taur=-0.5;
    static const double taui=0.866025403784439;
    int     i, k, ic;
    double  ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;

    for(k=0; k<l1; k++)
    {
	cr2=cc[(k+l1)*ido]+cc[(k+2*l1)*ido];
	ch[3*k*ido]=cc[k*ido]+cr2;
	ch[(3*k+2)*ido]=taui*(cc[(k+l1*2)*ido]-cc[(k+l1)*ido]);
	ch[ido-1+(3*k+1)*ido]=cc[k*ido]+taur*cr2;
    }
    if(ido==1)
	return;
    for(k=0; k<l1; k++)
    {
	for(i=2; i<ido; i+=2)
	{
	    ic=ido-i;
	    dr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+
		wa1[i-1]*cc[i+(k+l1)*ido];
	    di2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
	    dr3=wa2[i-2]*cc[i-1+(k+l1*2)*ido]+wa2[i-1]*cc[i+(k+l1*2)*ido];
	    di3=wa2[i-2]*cc[i+(k+l1*2)*ido]-wa2[i-1]*cc[i-1+(k+l1*2)*ido];
	    cr2=dr2+dr3;
	    ci2=di2+di3;
	    ch[i-1+3*k*ido]=cc[i-1+k*ido]+cr2;
	    ch[i+3*k*ido]=cc[i+k*ido]+ci2;
	    tr2=cc[i-1+k*ido]+taur*cr2;
	    ti2=cc[i+k*ido]+taur*ci2;
	    tr3=taui*(di2-di3);
	    ti3=taui*(dr3-dr2);
	    ch[i-1+(3*k+2)*ido]=tr2+tr3;
	    ch[ic-1+(3*k+1)*ido]=tr2-tr3;
	    ch[i+(3*k+2)*ido]=ti2+ti3;
	    ch[ic+(3*k+1)*ido]=ti3-ti2;
	}
    }
}				/*radf3*/


static void 
radb3(int ido, int l1, const double cc[], double ch[],
       const double wa1[], const double wa2[])
{
    static const double taur=-0.5;
    static const double taui=0.866025403784439;
    int     i, k, ic;
    double  ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;

    for(k=0; k<l1; k++)
    {
	tr2=2*cc[ido-1+(3*k+1)*ido];
	cr2=cc[3*k*ido]+taur*tr2;
	ch[k*ido]=cc[3*k*ido]+tr2;
	ci3=2*taui*cc[(3*k+2)*ido];
	ch[(k+l1)*ido]=cr2-ci3;
	ch[(k+2*l1)*ido]=cr2+ci3;
    }
    if(ido==1)
	return;
    for(k=0; k<l1; k++)
    {
	for(i=2; i<ido; i+=2)
	{
	    ic=ido-i;
	    tr2=cc[i-1+(3*k+2)*ido]+cc[ic-1+(3*k+1)*ido];
	    cr2=cc[i-1+3*k*ido]+taur*tr2;
	    ch[i-1+k*ido]=cc[i-1+3*k*ido]+tr2;
	    ti2=cc[i+(3*k+2)*ido]-cc[ic+(3*k+1)*ido];
	    ci2=cc[i+3*k*ido]+taur*ti2;
	    ch[i+k*ido]=cc[i+3*k*ido]+ti2;
	    cr3=taui*(cc[i-1+(3*k+2)*ido]-cc[ic-1+(3*k+1)*ido]);
	    ci3=taui*(cc[i+(3*k+2)*ido]+cc[ic+(3*k+1)*ido]);
	    dr2=cr2-ci3;
	    dr3=cr2+ci3;
	    di2=ci2+cr3;
	    di3=ci2-cr3;
	    ch[i-1+(k+l1)*ido]=wa1[i-2]*dr2-wa1[i-1]*di2;
	    ch[i+(k+l1)*ido]=wa1[i-2]*di2+wa1[i-1]*dr2;
	    ch[i-1+(k+2*l1)*ido]=wa2[i-2]*dr3-wa2[i-1]*di3;
	    ch[i+(k+2*l1)*ido]=wa2[i-2]*di3+wa2[i-1]*dr3;
	}
    }
}				/*radb3*/


static void 
radf4(int ido, int l1, const double cc[], double ch[],
       const double wa1[], const double wa2[], const double wa3[])
{
    static const double hsqt2=0.7071067811865475;
    int     i, k, ic;
    double  ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
            tr3, tr4;

    for(k=0; k<l1; k++)
    {
	tr1=cc[(k+l1)*ido]+cc[(k+3*l1)*ido];
	tr2=cc[k*ido]+cc[(k+2*l1)*ido];
	ch[4*k*ido]=tr1+tr2;
	ch[ido-1+(4*k+3)*ido]=tr2-tr1;
	ch[ido-1+(4*k+1)*ido]=cc[k*ido]-cc[(k+2*l1)*ido];
	ch[(4*k+2)*ido]=cc[(k+3*l1)*ido]-cc[(k+l1)*ido];
    }
    if(ido<2)
	return;
    if(ido !=2)
    {
	for(k=0; k<l1; k++)
	{
	    for(i=2; i<ido; i+=2)
	    {
		ic=ido-i;
		cr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
		ci2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
		cr3=wa2[i-2]*cc[i-1+(k+2*l1)*ido]+wa2[i-1]*cc[i+(k+2*l1)*ido];
		ci3=wa2[i-2]*cc[i+(k+2*l1)*ido]-wa2[i-1]*cc[i-1+(k+2*l1)*ido];
		cr4=wa3[i-2]*cc[i-1+(k+3*l1)*ido]+wa3[i-1]*cc[i+(k+3*l1)*ido];
		ci4=wa3[i-2]*cc[i+(k+3*l1)*ido]-wa3[i-1]*cc[i-1+(k+3*l1)*ido];
		tr1=cr2+cr4;
		tr4=cr4-cr2;
		ti1=ci2+ci4;
		ti4=ci2-ci4;
		ti2=cc[i+k*ido]+ci3;
		ti3=cc[i+k*ido]-ci3;
		tr2=cc[i-1+k*ido]+cr3;
		tr3=cc[i-1+k*ido]-cr3;
		ch[i-1+4*k*ido]=tr1+tr2;
		ch[ic-1+(4*k+3)*ido]=tr2-tr1;
		ch[i+4*k*ido]=ti1+ti2;
		ch[ic+(4*k+3)*ido]=ti1-ti2;
		ch[i-1+(4*k+2)*ido]=ti4+tr3;
		ch[ic-1+(4*k+1)*ido]=tr3-ti4;
		ch[i+(4*k+2)*ido]=tr4+ti3;
		ch[ic+(4*k+1)*ido]=tr4-ti3;
	    }
	}
	if(ido%2==1)
	    return;
    }
    for(k=0; k<l1; k++)
    {
	ti1=-hsqt2*(cc[ido-1+(k+l1)*ido]+cc[ido-1+(k+3*l1)*ido]);
	tr1=hsqt2*(cc[ido-1+(k+l1)*ido]-cc[ido-1+(k+3*l1)*ido]);
	ch[ido-1+4*k*ido]=tr1+cc[ido-1+k*ido];
	ch[ido-1+(4*k+2)*ido]=cc[ido-1+k*ido]-tr1;
	ch[(4*k+1)*ido]=ti1-cc[ido-1+(k+2*l1)*ido];
	ch[(4*k+3)*ido]=ti1+cc[ido-1+(k+2*l1)*ido];
    }
}				/*radf4*/


static void 
radb4(int ido, int l1, const double cc[], double ch[],
       const double wa1[], const double wa2[], const double wa3[])
{
    static const double sqrt2=1.414213562373095;
    int     i, k, ic;
    double  ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
            tr3, tr4;

    for(k=0; k<l1; k++)
    {
	tr1=cc[4*k*ido]-cc[ido-1+(4*k+3)*ido];
	tr2=cc[4*k*ido]+cc[ido-1+(4*k+3)*ido];
	tr3=cc[ido-1+(4*k+1)*ido]+cc[ido-1+(4*k+1)*ido];
	tr4=cc[(4*k+2)*ido]+cc[(4*k+2)*ido];
	ch[k*ido]=tr2+tr3;
	ch[(k+l1)*ido]=tr1-tr4;
	ch[(k+2*l1)*ido]=tr2-tr3;
	ch[(k+3*l1)*ido]=tr1+tr4;
    }
    if(ido<2)
	return;
    if(ido !=2)
    {
	for(k=0; k<l1;++k)
	{
	    for(i=2; i<ido; i+=2)
	    {
		ic=ido-i;
		ti1=cc[i+4*k*ido]+cc[ic+(4*k+3)*ido];
		ti2=cc[i+4*k*ido]-cc[ic+(4*k+3)*ido];
		ti3=cc[i+(4*k+2)*ido]-cc[ic+(4*k+1)*ido];
		tr4=cc[i+(4*k+2)*ido]+cc[ic+(4*k+1)*ido];
		tr1=cc[i-1+4*k*ido]-cc[ic-1+(4*k+3)*ido];
		tr2=cc[i-1+4*k*ido]+cc[ic-1+(4*k+3)*ido];
		ti4=cc[i-1+(4*k+2)*ido]-cc[ic-1+(4*k+1)*ido];
		tr3=cc[i-1+(4*k+2)*ido]+cc[ic-1+(4*k+1)*ido];
		ch[i-1+k*ido]=tr2+tr3;
		cr3=tr2-tr3;
		ch[i+k*ido]=ti2+ti3;
		ci3=ti2-ti3;
		cr2=tr1-tr4;
		cr4=tr1+tr4;
		ci2=ti1+ti4;
		ci4=ti1-ti4;
		ch[i-1+(k+l1)*ido]=wa1[i-2]*cr2-wa1[i-1]*ci2;
		ch[i+(k+l1)*ido]=wa1[i-2]*ci2+wa1[i-1]*cr2;
		ch[i-1+(k+2*l1)*ido]=wa2[i-2]*cr3-wa2[i-1]*ci3;
		ch[i+(k+2*l1)*ido]=wa2[i-2]*ci3+wa2[i-1]*cr3;
		ch[i-1+(k+3*l1)*ido]=wa3[i-2]*cr4-wa3[i-1]*ci4;
		ch[i+(k+3*l1)*ido]=wa3[i-2]*ci4+wa3[i-1]*cr4;
	    }
	}
	if(ido%2==1)
	    return;
    }
    for(k=0; k<l1; k++)
    {
	ti1=cc[(4*k+1)*ido]+cc[(4*k+3)*ido];
	ti2=cc[(4*k+3)*ido]-cc[(4*k+1)*ido];
	tr1=cc[ido-1+4*k*ido]-cc[ido-1+(4*k+2)*ido];
	tr2=cc[ido-1+4*k*ido]+cc[ido-1+(4*k+2)*ido];
	ch[ido-1+k*ido]=tr2+tr2;
	ch[ido-1+(k+l1)*ido]=sqrt2*(tr1-ti1);
	ch[ido-1+(k+2*l1)*ido]=ti2+ti2;
	ch[ido-1+(k+3*l1)*ido]=-sqrt2*(tr1+ti1);
    }
}				/*radb4*/


static void 
radf5(int ido, int l1, const double cc[], double ch[],
       const double wa1[], const double wa2[], const double wa3[], const double wa4[])
{
    static const double tr11=0.309016994374947;
    static const double ti11=0.951056516295154;
    static const double tr12=-0.809016994374947;
    static const double ti12=0.587785252292473;
    int     i, k, ic;
    double  ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
            dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;

    for(k=0; k<l1; k++)
    {
	cr2=cc[(k+4*l1)*ido]+cc[(k+l1)*ido];
	ci5=cc[(k+4*l1)*ido]-cc[(k+l1)*ido];
	cr3=cc[(k+3*l1)*ido]+cc[(k+2*l1)*ido];
	ci4=cc[(k+3*l1)*ido]-cc[(k+2*l1)*ido];
	ch[5*k*ido]=cc[k*ido]+cr2+cr3;
	ch[ido-1+(5*k+1)*ido]=cc[k*ido]+tr11*cr2+tr12*cr3;
	ch[(5*k+2)*ido]=ti11*ci5+ti12*ci4;
	ch[ido-1+(5*k+3)*ido]=cc[k*ido]+tr12*cr2+tr11*cr3;
	ch[(5*k+4)*ido]=ti12*ci5-ti11*ci4;
    }
    if(ido==1)
	return;
    for(k=0; k<l1;++k)
    {
	for(i=2; i<ido; i+=2)
	{
	    ic=ido-i;
	    dr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
	    di2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
	    dr3=wa2[i-2]*cc[i-1+(k+2*l1)*ido]+wa2[i-1]*cc[i+(k+2*l1)*ido];
	    di3=wa2[i-2]*cc[i+(k+2*l1)*ido]-wa2[i-1]*cc[i-1+(k+2*l1)*ido];
	    dr4=wa3[i-2]*cc[i-1+(k+3*l1)*ido]+wa3[i-1]*cc[i+(k+3*l1)*ido];
	    di4=wa3[i-2]*cc[i+(k+3*l1)*ido]-wa3[i-1]*cc[i-1+(k+3*l1)*ido];
	    dr5=wa4[i-2]*cc[i-1+(k+4*l1)*ido]+wa4[i-1]*cc[i+(k+4*l1)*ido];
	    di5=wa4[i-2]*cc[i+(k+4*l1)*ido]-wa4[i-1]*cc[i-1+(k+4*l1)*ido];
	    cr2=dr2+dr5;
	    ci5=dr5-dr2;
	    cr5=di2-di5;
	    ci2=di2+di5;
	    cr3=dr3+dr4;
	    ci4=dr4-dr3;
	    cr4=di3-di4;
	    ci3=di3+di4;
	    ch[i-1+5*k*ido]=cc[i-1+k*ido]+cr2+cr3;
	    ch[i+5*k*ido]=cc[i+k*ido]+ci2+ci3;
	    tr2=cc[i-1+k*ido]+tr11*cr2+tr12*cr3;
	    ti2=cc[i+k*ido]+tr11*ci2+tr12*ci3;
	    tr3=cc[i-1+k*ido]+tr12*cr2+tr11*cr3;
	    ti3=cc[i+k*ido]+tr12*ci2+tr11*ci3;
	    tr5=ti11*cr5+ti12*cr4;
	    ti5=ti11*ci5+ti12*ci4;
	    tr4=ti12*cr5-ti11*cr4;
	    ti4=ti12*ci5-ti11*ci4;
	    ch[i-1+(5*k+2)*ido]=tr2+tr5;
	    ch[ic-1+(5*k+1)*ido]=tr2-tr5;
	    ch[i+(5*k+2)*ido]=ti2+ti5;
	    ch[ic+(5*k+1)*ido]=ti5-ti2;
	    ch[i-1+(5*k+4)*ido]=tr3+tr4;
	    ch[ic-1+(5*k+3)*ido]=tr3-tr4;
	    ch[i+(5*k+4)*ido]=ti3+ti4;
	    ch[ic+(5*k+3)*ido]=ti4-ti3;
	}
    }
}				/*radf5*/


static void 
radb5(int ido, int l1, const double cc[], double ch[],
       const double wa1[], const double wa2[], const double wa3[], const double wa4[])
{
    static const double tr11=0.309016994374947;
    static const double ti11=0.951056516295154;
    static const double tr12=-0.809016994374947;
    static const double ti12=0.587785252292473;
    int     i, k, ic;
    double  ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
            ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;

    for(k=0; k<l1; k++)
    {
	ti5=2*cc[(5*k+2)*ido];
	ti4=2*cc[(5*k+4)*ido];
	tr2=2*cc[ido-1+(5*k+1)*ido];
	tr3=2*cc[ido-1+(5*k+3)*ido];
	ch[k*ido]=cc[5*k*ido]+tr2+tr3;
	cr2=cc[5*k*ido]+tr11*tr2+tr12*tr3;
	cr3=cc[5*k*ido]+tr12*tr2+tr11*tr3;
	ci5=ti11*ti5+ti12*ti4;
	ci4=ti12*ti5-ti11*ti4;
	ch[(k+l1)*ido]=cr2-ci5;
	ch[(k+2*l1)*ido]=cr3-ci4;
	ch[(k+3*l1)*ido]=cr3+ci4;
	ch[(k+4*l1)*ido]=cr2+ci5;
    }
    if(ido==1)
	return;
    for(k=0; k<l1;++k)
    {
	for(i=2; i<ido; i+=2)
	{
	    ic=ido-i;
	    ti5=cc[i+(5*k+2)*ido]+cc[ic+(5*k+1)*ido];
	    ti2=cc[i+(5*k+2)*ido]-cc[ic+(5*k+1)*ido];
	    ti4=cc[i+(5*k+4)*ido]+cc[ic+(5*k+3)*ido];
	    ti3=cc[i+(5*k+4)*ido]-cc[ic+(5*k+3)*ido];
	    tr5=cc[i-1+(5*k+2)*ido]-cc[ic-1+(5*k+1)*ido];
	    tr2=cc[i-1+(5*k+2)*ido]+cc[ic-1+(5*k+1)*ido];
	    tr4=cc[i-1+(5*k+4)*ido]-cc[ic-1+(5*k+3)*ido];
	    tr3=cc[i-1+(5*k+4)*ido]+cc[ic-1+(5*k+3)*ido];
	    ch[i-1+k*ido]=cc[i-1+5*k*ido]+tr2+tr3;
	    ch[i+k*ido]=cc[i+5*k*ido]+ti2+ti3;
	    cr2=cc[i-1+5*k*ido]+tr11*tr2+tr12*tr3;

	    ci2=cc[i+5*k*ido]+tr11*ti2+tr12*ti3;
	    cr3=cc[i-1+5*k*ido]+tr12*tr2+tr11*tr3;

	    ci3=cc[i+5*k*ido]+tr12*ti2+tr11*ti3;
	    cr5=ti11*tr5+ti12*tr4;
	    ci5=ti11*ti5+ti12*ti4;
	    cr4=ti12*tr5-ti11*tr4;
	    ci4=ti12*ti5-ti11*ti4;
	    dr3=cr3-ci4;
	    dr4=cr3+ci4;
	    di3=ci3+cr4;
	    di4=ci3-cr4;
	    dr5=cr2+ci5;
	    dr2=cr2-ci5;
	    di5=ci2-cr5;
	    di2=ci2+cr5;
	    ch[i-1+(k+l1)*ido]=wa1[i-2]*dr2-wa1[i-1]*di2;
	    ch[i+(k+l1)*ido]=wa1[i-2]*di2+wa1[i-1]*dr2;
	    ch[i-1+(k+2*l1)*ido]=wa2[i-2]*dr3-wa2[i-1]*di3;
	    ch[i+(k+2*l1)*ido]=wa2[i-2]*di3+wa2[i-1]*dr3;
	    ch[i-1+(k+3*l1)*ido]=wa3[i-2]*dr4-wa3[i-1]*di4;
	    ch[i+(k+3*l1)*ido]=wa3[i-2]*di4+wa3[i-1]*dr4;
	    ch[i-1+(k+4*l1)*ido]=wa4[i-2]*dr5-wa4[i-1]*di5;
	    ch[i+(k+4*l1)*ido]=wa4[i-2]*di5+wa4[i-1]*dr5;
	}
    }
}				/*radb5*/


static void 
radfg(int ido, int ip, int l1, int idl1,
       double cc[], double c1[], double c2[], double ch[], double ch2[], const double wa[])
{
    static const double twopi=6.28318530717959;
    int     idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
    double  dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;

    arg=twopi / ip;
    dcp=cos(arg);
    dsp=sin(arg);
    ipph=(ip+1)/ 2;
    nbd=(ido-1)/ 2;
    if(ido !=1)
    {
	for(ik=0; ik<idl1; ik++)
	    ch2[ik]=c2[ik];
	for(j=1; j<ip; j++)
	    for(k=0; k<l1; k++)
		ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido];
	if(nbd<=l1)
	{
	    is=-ido;
	    for(j=1; j<ip; j++)
	    {
		is+=ido;
		idij=is-1;
		for(i=2; i<ido; i+=2)
		{
		    idij+=2;
		    for(k=0; k<l1; k++)
		    {
			ch[i-1+(k+j*l1)*ido]=
			    wa[idij-1]*c1[i-1+(k+j*l1)*ido]+wa[idij]*c1[i+(k+j*l1)*ido];
			ch[i+(k+j*l1)*ido]=
			    wa[idij-1]*c1[i+(k+j*l1)*ido]-wa[idij]*c1[i-1+(k+j*l1)*ido];
		    }
		}
	    }
	}
	else
	{
	    is=-ido;
	    for(j=1; j<ip; j++)
	    {
		is+=ido;
		for(k=0; k<l1; k++)
		{
		    idij=is-1;
		    for(i=2; i<ido; i+=2)
		    {
			idij+=2;
			ch[i-1+(k+j*l1)*ido]=
			    wa[idij-1]*c1[i-1+(k+j*l1)*ido]+wa[idij]*c1[i+(k+j*l1)*ido];
			ch[i+(k+j*l1)*ido]=
			    wa[idij-1]*c1[i+(k+j*l1)*ido]-wa[idij]*c1[i-1+(k+j*l1)*ido];
		    }
		}
	    }
	}
	if(nbd>=l1)
	{
	    for(j=1; j<ipph; j++)
	    {
		jc=ip-j;
		for(k=0; k<l1; k++)
		{
		    for(i=2; i<ido; i+=2)
		    {
			c1[i-1+(k+j*l1)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
			c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
			c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
			c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
		    }
		}
	    }
	}
	else
	{
	    for(j=1; j<ipph; j++)
	    {
		jc=ip-j;
		for(i=2; i<ido; i+=2)
		{
		    for(k=0; k<l1; k++)
		    {
			c1[i-1+(k+j*l1)*ido]=
			    ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
			c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
			c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
			c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
		    }
		}
	    }
	}
    }
    else
    {				/*now ido==1*/
	for(ik=0; ik<idl1; ik++)
	    c2[ik]=ch2[ik];
    }
    for(j=1; j<ipph; j++)
    {
	jc=ip-j;
	for(k=0; k<l1; k++)
	{
	    c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido]+ch[(k+jc*l1)*ido];
	    c1[(k+jc*l1)*ido]=ch[(k+jc*l1)*ido]-ch[(k+j*l1)*ido];
	}
    }

    ar1=1;
    ai1=0;
    for(l=1; l<ipph; l++)
    {
	lc=ip-l;
	ar1h=dcp*ar1-dsp*ai1;
	ai1=dcp*ai1+dsp*ar1;
	ar1=ar1h;
	for(ik=0; ik<idl1; ik++)
	{
	    ch2[ik+l*idl1]=c2[ik]+ar1*c2[ik+idl1];
	    ch2[ik+lc*idl1]=ai1*c2[ik+(ip-1)*idl1];
	}
	dc2=ar1;
	ds2=ai1;
	ar2=ar1;
	ai2=ai1;
	for(j=2; j<ipph; j++)
	{
	    jc=ip-j;
	    ar2h=dc2*ar2-ds2*ai2;
	    ai2=dc2*ai2+ds2*ar2;
	    ar2=ar2h;
	    for(ik=0; ik<idl1; ik++)
	    {
		ch2[ik+l*idl1]+=ar2*c2[ik+j*idl1];
		ch2[ik+lc*idl1]+=ai2*c2[ik+jc*idl1];
	    }
	}
    }
    for(j=1; j<ipph; j++)
	for(ik=0; ik<idl1; ik++)
	    ch2[ik]+=c2[ik+j*idl1];

    if(ido>=l1)
    {
	for(k=0; k<l1; k++)
	{
	    for(i=0; i<ido; i++)
	    {
		cc[i+k*ip*ido]=ch[i+k*ido];
	    }
	}
    }
    else
    {
	for(i=0; i<ido; i++)
	{
	    for(k=0; k<l1; k++)
	    {
		cc[i+k*ip*ido]=ch[i+k*ido];
	    }
	}
    }
    for(j=1; j<ipph; j++)
    {
	jc=ip-j;
	j2=2*j;
	for(k=0; k<l1; k++)
	{
	    cc[ido-1+(j2-1+k*ip)*ido]=
		ch[(k+j*l1)*ido];
	    cc[(j2+k*ip)*ido]=