#define SVAR	cdftci_build_svar(ICA, ICB, JCA, JCB, ITA, ITB, JTA, JTB, \
				NBasis, NAlpha, NBeta, NB2)
#define FW_GRAD(C,T,I,J,OI,OJ)	\
	do { \
	tmp1 = *((C) + (I)); \
	tmp2 = *((T) + (J)); \
	*(jA + en * (OJ) + (OI)) = 0.0; \
	*((C) + (I)) = tmp1 + 2.0 * epsilon; \
	*((T) + (J)) = tmp2 + 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-1.0 / 12.0) * (-1.0 / 12.0) * SVAR; \
	*((T) + (J)) = tmp2 + 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-1.0 / 12.0) * (2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-1.0 / 12.0) * (-2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-1.0 / 12.0) * (1.0 / 12.0) * SVAR; \
	*((C) + (I)) = tmp1 + 1.0 * epsilon; \
	*((T) + (J)) = tmp2 + 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (2.0 / 3.0) * (-1.0 / 12.0) * SVAR; \
	*((T) + (J)) = tmp2 + 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (2.0 / 3.0) * (2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (2.0 / 3.0) * (-2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (2.0 / 3.0) * (1.0 / 12.0) * SVAR; \
	*((C) + (I)) = tmp1 - 1.0 * epsilon; \
	*((T) + (J)) = tmp2 + 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-2.0 / 3.0) * (-1.0 / 12.0) * SVAR; \
	*((T) + (J)) = tmp2 + 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-2.0 / 3.0) * (2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-2.0 / 3.0) * (-2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (-2.0 / 3.0) * (1.0 / 12.0) * SVAR; \
	*((C) + (I)) = tmp1 - 2.0 * epsilon; \
	*((T) + (J)) = tmp2 + 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (1.0 / 12.0) * (-1.0 / 12.0) * SVAR; \
	*((T) + (J)) = tmp2 + 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (1.0 / 12.0) * (2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 1.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (1.0 / 12.0) * (-2.0 / 3.0) * SVAR; \
	*((T) + (J)) = tmp2 - 2.0 * epsilon; \
	*(jA + en * (OJ) + (OI)) += (1.0 / 12.0) * (1.0 / 12.0) * SVAR; \
	*((C) + (I)) = tmp1; \
	*((T) + (J)) = tmp2; \
	} while(0)
    for(int ii = 0; ii < NBasis * NAlpha; ++ii) {
	oi = ii;
	for(int jj = 0; jj < NBasis * NAlpha; ++jj) {
	    oj = jj;
	    FW_GRAD(ICA, ITA, ii, jj, oi, oj);
	}
	for(int jj = 0; jj < NBasis * NBeta; ++jj) {
	    oj = jj + NBasis * NAlpha;
	    FW_GRAD(ICA, ITB, ii, jj, oi, oj);
	}
    }
    for(int ii = 0; ii < NBasis * NBeta; ++ii) {
	oi = ii + NBasis * NAlpha;
	for(int jj = 0; jj < NBasis * NAlpha; ++jj) {
	    oj = jj;
	    FW_GRAD(ICB, ITA, ii, jj, oi, oj);
	}
	for(int jj = 0; jj < NBasis * NBeta; ++jj) {
	    oj = jj + NBasis * NAlpha;
	    FW_GRAD(ICB, ITB, ii, jj, oi, oj);
	}
    }
#undef FW_GRAD
#undef SVAR