/*  ssqcoupl.h by K.Tsuru  */
// since ver. 2.18
/***************************************************************
SN library
Template class of square coupling method.

It evaluates power series
S_L = A_0 + A_1 * R + A_2 * R^2 + ...
It can be used to evaluate a polynomial of SLong or SInteger
 e.g. radix conversion.

This is a generalization of divide radix conversion 
of "SLong SInteger::DConvToDec() const;" in "simdcvdc.cpp".
***************************************************************/

#ifndef SQUARE_COUPLING_TEMPLATE_H
#define SQUARE_COUPLING_TEMPLATE_H

template <class SqCType> class SquareCoupling {
private:
    const SqCType ZERO;
    void abort();   // abort by syntax error
    void freeHalf(ulong from, ulong to);    // free memory of disused elements.
    bool isOk;
    SqCType *A, R;
    ulong L, userL; // L/2 < userL =<  L
    void setL(ulong s) {
        L = userL = s;
        if ( userL & (userL - 1) ) L = ceilpow2(userL); // L != 2^n It changes into 2^n.
        SCalcInfo::upToTerm = L;
    }
    SqCType element(ulong k) {
        return (k < userL) ? A[k] : ZERO;
    }
public:
    // A function sets coefficients A[].
    // caution : The content of a[] is destroied.
    void setAR(SqCType *a, const SqCType& r, ulong size) {
        isOk = false;   // reset for second usage
        A = a;  R = r;  setL(size);
    }
    SquareCoupling() : ZERO(0.0), isOk(false), L(0), userL(0) {}
    SquareCoupling(SqCType *a, const SqCType& r, ulong size)
        : ZERO(0.0), isOk(false), A(a), R(r) {
        setL(size);
    }
    virtual ~SquareCoupling(){}
    void putTogether();
    bool isOK() const { return isOk; }
    SqCType getValue(); //It returns the value of A[0].
};
/***************************************
* Implementation of class SquareCoupling
****************************************/
template <class SqCType> void SquareCoupling<SqCType> ::freeHalf(ulong from, ulong to) {
    for (ulong i = from; i < to; i++) A[i].SizeZero();
}
template <class SqCType> void SquareCoupling<SqCType>::abort() {
    R.SetError(R.SYNTAX_ERR, "class SquareCoupling", 1601);
}
template <class SqCType> void SquareCoupling<SqCType>::putTogether() {
    if (L == 0) abort();
    SqCType Q(R);

    // First step "element()" is used because A[i] (userL =< i < L) are undecided.
#if 1
    ulong j, uj = userL;
    if (userL & 1) uj++;
    for ( j = 0; j < uj - 2; j += 2) A[j/2] = A[j+1] * Q + A[j];

    // here j == uj-2
    A[j/2] = element(j+1) * Q + element(j);
    j += 2;
    for ( ; j< L ; j += 2) A[j/2].SetZero();
#else
    // simple version
    for (ulong i = 0; i < L; i += 2) A[i/2] = element(i+1) * Q + element(i);
#endif

    ulong n = L / 2;
    freeHalf(n, userL);
    Q *= Q;

    while (1) {
        for (ulong i = 0; i < n ; i += 2) A[i/2] = A[i+1] * Q + A[i];
        n /= 2;
        freeHalf(n, 2 * n - 1);
        if (n == 1) break;
        Q *= Q;
    }
    isOk = true;
}
template <class SqCType> SqCType SquareCoupling<SqCType>::getValue() {
    if (!isOk) putTogether();
    return A[0];
}

#endif // SQUARE_COUPLING_TEMPLATE_H