/*  ssqcpinv.h by K.Tsuru  */
// since ver. 2.18
/***************************************************************
SN library
Template class of square coupling method for inverse R.
This is a special case of binary splitting method.
It evaluates power series
S_L = A_0/R + A_1 / R^2 + A_2 / R^3 + ...
It can be used to evaluate a polynomial of SLong or SInteger
 e.g. radix conversion of decimal.

algorithm:

(0.A[0]A[1]A[2]...)_R =
                        A[k]
S_L = \sum_{k=0}^{L-1}--------
                       R^(k+1)

                        A[2k]        A[2k+1]
 = \sum_{k=0}^{L/2-1}{---------- + ----------- }
                       R^(2k+1)     R^(2k+2)

                        A[2k] R + A[2k+1]
 = \sum_{k=0}^{L/2-1} -----------------------
                           (R^2)^(k+1)

Then repeating while L < 1
      L/2 <-- L
      R <-- R*R
      A[k] <-- A[2k]*R+A[2k+1]
yields
S_L = A[0]/R
***************************************************************/

#ifndef SQUARE_COUPLING_INV_TEMPLATE_H
#define SQUARE_COUPLING_INV_TEMPLATE_H

template <class SqCType> class SquareCouplingInverseR {
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);
    }
    SquareCouplingInverseR() : ZERO(0.0), isOk(false), L(0), userL(0) {}
    SquareCouplingInverseR(SqCType *a, const SqCType& r, ulong size)
        : ZERO(0.0), isOk(false), A(a), R(r) {
        setL(size);
    }
    virtual ~SquareCouplingInverseR(){}
    void putTogether();
    bool isOK() const { return isOk; }

    // sum is obtained S_L = A[0]/R
    // It returns the value of A[0].
    SqCType getA() {
        if(!isOk) putTogether();
        return A[0];
    }
    // It returns the value of R.
    SqCType getR() {
        if(!isOk) putTogether();
        return R;
    }
};
/***************************************
* Implementation of class SquareCouplingInverseR
****************************************/
template <class SqCType> void SquareCouplingInverseR<SqCType> ::freeHalf(ulong from, ulong to) {
    for (ulong i = from; i < to; i++) A[i].SizeZero();
}
template <class SqCType> void SquareCouplingInverseR<SqCType>::abort() {
    R.SetError(R.SYNTAX_ERR, "class SquareCouplingInverseR", 1602);
}

template <class SqCType> void SquareCouplingInverseR<SqCType>::putTogether() {
    if (L == 0) abort();

    // 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] * R + A[j+1];

    // here j == uj-2
    A[j/2] = element(j) * R + element(j+1);
    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) * R + element(i+1);
#endif

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

    while (n > 1uL) {
        for (ulong i = 0; i < n ; i += 2) A[i/2] = A[i] * R + A[i+1];
        n /= 2;
        freeHalf(n, 2 * n - 1);
        R *= R;
    }
    isOk = true;
}

#endif // SQUARE_COUPLING_INV_TEMPLATE_H