libsnark：Lagrange interpolation at Roots of unity

libsnark\depends\libfqfft\libfqfft\evaluation_domain\domains\basic_radix2_domain_aux.tcc中的_basic_radix2_evaluate_all_lagrange_polynomials(…)函数的计算依据为《Behavior of the Lagrange Interpolants in the Roots of unity》论文中如下图的公式：

根据论文公式，令$w=e(\frac{1}{m})=e^{j\ast \frac{2\ast \pi }{m}}=cos(\frac{2\ast \pi }{m})+j\ast sin(\frac{2\ast \pi }{m})$
$$\Downarrow$$
$$L_{m-1}(f,z)=\frac{z^m-1}{m}\ast \sum _{k=0}^{n-1}\frac{w^k\ast f(w^k)}{z-w^k}$$
令$z=t$，可以推导出：
$u[k]=\frac{t^m-1}{m}\ast w^k\ast \frac{1}{t-w^k}$
$S[k]=f(w^k)$
使得 $L_{m-1}(f,t)=u\cdot S$(点乘)

/**
 * Compute the m Lagrange coefficients, relative to the set S={omega^{0},...,omega^{m-1}}, at the field element t.
 */
template<typename FieldT>
std::vector<FieldT> _basic_radix2_evaluate_all_lagrange_polynomials(const size_t m, const FieldT &t)
{
    if (m == 1)
    {
        return std::vector<FieldT>(1, FieldT::one());
    }

    if (m != (1u << libff::log2(m))) throw DomainSizeException("expected m == (1u << log2(m))");

    const FieldT omega = libff::get_root_of_unity<FieldT>(m);

    std::vector<FieldT> u(m, FieldT::zero());

    /*
     If t equals one of the roots of unity in S={omega^{0},...,omega^{m-1}}
     then output 1 at the right place, and 0 elsewhere
     */

    if ((t^m) == (FieldT::one()))
    {
        FieldT omega_i = FieldT::one();
        for (size_t i = 0; i < m; ++i)
        {
            if (omega_i == t) // i.e., t equals omega^i
            {
                u[i] = FieldT::one();
                return u;
            }

            omega_i *= omega;
        }
    }

    /*
     Otherwise, if t does not equal any of the roots of unity in S,
     then compute each L_{i,S}(t) as Z_{S}(t) * v_i / (t-\omega^i)
     where:
     - Z_{S}(t) = \prod_{j} (t-\omega^j) = (t^m-1), and
     - v_{i} = 1 / \prod_{j \neq i} (\omega^i-\omega^j).
     Below we use the fact that v_{0} = 1/m and v_{i+1} = \omega * v_{i}. 
     // 计算依据见《Behavior of the Lagrange Interpolants in the Roots of unity.pdf》中的Ln-1(f, z)的系数
     // S=[f(1),f(w),f(w^2),....,f(w^m-1)]
     */

    const FieldT Z = (t^m)-FieldT::one();
    FieldT l = Z * FieldT(m).inverse();
    FieldT r = FieldT::one();
    for (size_t i = 0; i < m; ++i)
    {
        u[i] = l * (t - r).inverse();
        l *= omega;
        r *= omega;
    }

    return u;
}