C++值元编程

——永远不要在OJ上使用值元编程，过于简单的没有优势，能有优势的编译错误。

背景

2019年10月，我在学习算法。有一道作业题，输入规模很小，可以用打表法解决。具体方案有以下三种：

1. 运行时预处理，生成所需的表格，根据输入直接找到对应项，稍加处理后输出；

2. 一个程序生成表格，作为提交程序的一部分，后续与方法1相同，这样就省去了运行时计算的步骤；

3. 以上两种方法结合，编译期计算表格，运行时直接查询，即元编程（metaprogramming）。

做题当然是用方法1或2，但是元编程已经埋下了种子。时隔大半年，我来补上这个坑。

题目

北京大学OpenJudge 百练4119 复杂的整数划分问题

描述

将正整数 (n) 表示成一系列正整数之和，(n = n_1 + n_2 + ... + n_k)，其中 (n_1 geq n_2 geq ... geq n_k geq 1)，(k geq 1)。正整数 (n) 的这种表示称为正整数 (n) 的划分。

输入

标准的输入包含若干组测试数据。每组测试数据是一行输入数据，包括两个整数 (N) 和 (K)。（ (0 le N leq 50)，(0 le K leq N) ）

输出

对于每组测试数据，输出以下三行数据：

第一行: (N) 划分成 (K) 个正整数之和的划分数目

第二行: (N) 划分成若干个不同正整数之和的划分数目

第三行: (N) 划分成若干个奇正整数之和的划分数目

样例输入

5 2

样例输出

2
3
3

提示

第一行: 4+1，3+2

第二行: 5，4+1，3+2

第三行: 5，1+1+3，1+1+1+1+1+1

解答

标准的动态规划题。用dp[c][i][j]表示把i分成c个正整数之和的方法数，其中每个数都不超过j。

第一行。初始化：由 (i leq j) 是否成立决定dp[1][i][j]的值，当 (i leq j) 时为1，划分为 (i = i)，否则无法划分，值为0。

递推：为了求dp[c][i][j]，对 (i = i_1 + i_2 + ... + i_c)，(i_1 geq i_2 geq ... geq i_c) 中的最大数 (i_1) 分类讨论，最小为 (1)，最大不超过 (i - 1)，因为 (c geq 2)，同时不超过 (j)，因为定义。最大数为 (n) 时，对于把 (i - n) 分成 (c - 1) 个数，每个数不超过 (n) 的划分，追加上 (n) 可得 (i) 的一个划分。(n) 只有这些取值，没有漏；对于不同的 (n)，由于最大数不一样，两个划分也不一样，没有多。故递推式为：

[dp[c][i][j] = sum_{n=1}^{min{i-1,j}}dp[c-1][i-n][n] ]

dp[K][N][N]即为所求ans1[K][N]。

第二行。可以把递推式中的dp[c - 1][i - n][n]修改为dp[c - 1][i - n][n - 1]后重新计算。由于只需一个与c无关的结果，可以省去c这一维度，相应地改变递推顺序，每轮累加。

另一种方法是利用已经计算好的ans1数组。设 (i = i_1 + i_2 + ... + + i_{c-1} + i_c)，其中 (i_1 ge i_2 ge ... ge i_{c+1} ge i_c ge 0)，则 (i_1 - left( c-1 ight) geq i_2 - left( c-2 ight) geq ... geq i_{c-1} - 1 geq i_c ge 0)，且 (left( i_1 - left( c-1 ight) ight) + left( i_2 - left( c-2 ight) ight) + ... + left( i_{c-1} - 1 ight) + left( i_c ight) = i - frac {c left( c-1 ight)} {2})，故把i划分成c个不同正整数之和的划分数目等于ans[c][i - c * (c - 1) / 2]，遍历c累加即得结果。

第三行。想法与第二行相似，也是找一个对应，此处从略。另外，数学上可以证明，第二行和第三行的结果一定是一样的。

#include <iostream>
#include <algorithm>

constexpr int max = 50;
int dp[max + 1][max + 1][max + 1] = { 0 };
int ans1[max + 1][max + 1] = { 0 };
int ans2[max + 1] = { 0 };
int ans3[max + 1] = { 0 };

int main()
{
    int num, k;
    for (int i = 1; i <= max; ++i)
        for (int j = 1; j <= max; ++j)
            dp[1][i][j] = i <= j;
    for (int cnt = 2; cnt <= max; ++cnt)
        for (int i = 1; i <= max; ++i)
            for (int j = 1; j <= max; ++j)
            {
                auto min = std::min(i - 1, j);
                for (int n = 1; n <= min; ++n)
                    dp[cnt][i][j] += dp[cnt - 1][i - n][n];
            }
    for (int cnt = 1; cnt <= max; ++cnt)
        for (int i = 1; i <= max; ++i)
            ans1[cnt][i] = dp[cnt][i][i];
    for (int i = 1; i <= max; ++i)
        for (int cnt = 1; cnt <= i; ++cnt)
        {
            int j = i - cnt * (cnt - 1) / 2;
            if (j <= 0)
                break;
            ans2[i] += ans1[cnt][j];
        }
    for (int i = 1; i <= max; ++i)
        for (int cnt = 1; cnt <= i; ++cnt)
        {
            int j = i + cnt;
            if (j % 2)
                continue;
            j /= 2;
            ans3[i] += ans1[cnt][j];
        }
    
    while (std::cin >> num)
    {
        std::cin >> k;
        std::cout << ans1[k][num] << std::endl;
        std::cout << ans2[num] << std::endl;
        std::cout << ans3[num] << std::endl;
    }
}

值元编程基础

元编程是指计算机程序能把其他程序作为它们的数据的编程技术。在目前的C++中，元编程体现为用代码生成代码，包括宏与模板。当我们使用了std::vector<int>中的任何一个名字时，std::vector类模板就用模板参数int, std::allocator<int>实例化为std::vector<int, std::allocator<int>>模板类，这是一种元编程，不过我们通常不这么讲。

狭义的C++模板元编程（template metaprogramming，TMP）包括值元编程、类型元编程，以及两者的相交。本文讨论的是值元编程，即为编译期值编程。

在C++中有两套工具可用于值元编程：模板和constexpr。C++模板是图灵完全的，这是模板被引入C++以后才被发现的，并不是C++模板的初衷，因此用模板做计算在C++中算不上一等用法，导致其语法比较冗长复杂。constexpr的初衷是提供纯正的编译期常量，后来才取消对计算的限制，但不能保证计算一定在编译期完成。总之，这两套工具都不完美，所以本文都会涉及。

严格来说，constexpr不符合上述对元编程的定义，但它确实可以提供运行时程序需要的数据，所以也归入元编程的类别。

constexpr式值元编程

从constexpr开始讲，是因为它与我们在C++中惯用的编程范式——过程式范式是一致的。

constexpr关键字在C++11中被引入。当时，constexpr函数中只能包含一条求值语句，就是return语句，返回值可以用于初始化constexpr变量，作模板参数等用途。如果需要分支语句，用三目运算符?:；如果需要循环语句，用函数递归实现。比如，计算阶乘：

constexpr int factorial(int n)
{
    return n <= 1 ? 1 : (n * factorial(n - 1));
}

对于编译期常量i，factorial(i)产生编译期常量；对于运行时值j，factorial(j)产生运行时值，也就是说，constexpr可以视为对既有函数的附加修饰。

然而，多数函数不止有一句return语句，constexpr对函数体的限制使它很难用于中等复杂的计算任务，为此C++14放宽了限制，允许定义局部变量，允许if-else、switch-case、while、for等控制流。factorial函数可以改写为：

constexpr int factorial(int n)
{
    int result = 1;
    for (; n > 1; --n)
        result *= n;
    return result;
}

也许你会觉得factorial函数的递归版本比循环版本易懂，那是因为你学习递归时接触的第一个例子就是它。对于C++开发者来说，大多数情况下首选的还是循环。

计算单个constexpr值用C++14就足够了，但是传递数组需要C++17，因为std::array的operator[]从C++17开始才是constexpr的。

整数划分问题的constexpr元编程实现需要C++17标准：

#include <iostream>
#include <utility>
#include <array>

constexpr int MAX = 50;

constexpr auto calculate_ans1()
{
    std::array<std::array<std::array<int, MAX + 1>, MAX + 1>, MAX + 1> dp{};
    std::array<std::array<int, MAX + 1>, MAX + 1> ans1{};
    constexpr int max = MAX;
    for (int i = 1; i <= max; ++i)
        for (int j = 1; j <= max; ++j)
            dp[1][i][j] = i <= j;
    for (int cnt = 2; cnt <= max; ++cnt)
        for (int i = 1; i <= max; ++i)
            for (int j = 1; j <= max; ++j)
            {
                auto min = std::min(i - 1, j);
                for (int n = 1; n <= min; ++n)
                    dp[cnt][i][j] += dp[cnt - 1][i - n][n];
            }
    for (int cnt = 1; cnt <= max; ++cnt)
        for (int i = 1; i <= max; ++i)
            ans1[cnt][i] = dp[cnt][i][i];
    return ans1;
}

constexpr auto calculate_ans2()
{
    constexpr auto ans1 = calculate_ans1();
    std::array<int, MAX + 1> ans2{};
    constexpr int max = MAX;
    for (int i = 1; i <= max; ++i)
        for (int cnt = 1; cnt <= i; ++cnt)
        {
            int j = i - cnt * (cnt - 1) / 2;
            if (j <= 0)
                break;
            ans2[i] += ans1[cnt][j];
        }
    return ans2;
}

int main()
{
    constexpr auto ans1 = calculate_ans1();
    constexpr auto ans2 = calculate_ans2();

    for (int cnt = 1; cnt <= 10; ++cnt)
    {
        for (int i = 1; i <= 10; ++i)
            std::cout << ans1[cnt][i] << ' ';+
        std::cout << std::endl;
    }
    std::cout << std::endl;
    for (int i = 1; i <= 50; ++i)
        std::cout << ans2[i] << ' ';
    std::cout << std::endl;

    int num, k;
    while (std::cin >> num)
    {
        std::cin >> k;
        std::cout << ans1[k][num] << std::endl;
        std::cout << ans2[num] << std::endl;
        std::cout << ans2[num] << std::endl;
    }
}

模板式值元编程

模板式与C++11中的constexpr式类似，必须把循环化为递归。事实上C++模板是一门函数式编程语言，对值元编程和类型元编程都是如此。

程序控制流有三种基本结构：顺序、分支与循环。

顺序

在函数式编程中，数据都是不可变的，函数总是接受若干参数，返回若干结果，参数和结果是不同的变量；修改原来的变量是不允许的。对于C++模板这门语言，函数是类模板，也称“元函数”（metafunction）；参数是模板参数；运算结果是模板类中定义的静态编译期常量（在C++11以前，常用enum来定义；C++11开始用constexpr）。

比如，对于参数 (x)，计算 (x + 1) 和 (x ^ 2) 的元函数：

template<int X>
struct PlusOne
{
    static constexpr int value = X + 1;
};

template<int X>
struct Square
{
    static constexpr int value = X * X;
};

这里假定运算数的类型为int。从C++17开始，可以用auto声明非类型模板参数。

顺序结构，是对数据依次进行多个操作，可以用函数嵌套来实现：

std::cout << PlusOne<1>::value << std::endl;
std::cout << Square<2>::value << std::endl;
std::cout << Square<PlusOne<3>::value>::value << std::endl;
std::cout << PlusOne<Square<4>::value>::value << std::endl;

或者借助constexpr函数，回归熟悉的过程式范式：

template<int X>
struct SquareAndIncrease
{
    static constexpr int calculate()
    {
        int x = X;
        x = x * x;
        x = x + 1;
        return x;
    }
    static constexpr int value = calculate();
};

void f()
{
    std::cout << SquareAndIncrease<5>::value << std::endl;
}

过程式方法同样可以用于分支和循环结构，以下省略；函数式方法可以相似地用于值元编程与类型元编程，所以我更青睐（主要还是逼格更高）。

分支

C++模板元编程实现分支的方式是模板特化与模板参数匹配，用一个额外的带默认值的bool类型模板参数作匹配规则，特化false或true的情形，另一种情形留给主模板。

比如，计算 (x) 的绝对值：

template<int X, bool Pos = (X > 0)>
struct AbsoluteHelper
{
    static constexpr int value = X;
};

template<int X>
struct AbsoluteHelper<X, false>
{
    static constexpr int value = -X;
};

如果你怕用户瞎写模板参数，可以再包装一层：

template<int X>
struct Absolute : AbsoluteHelper<X> { };

void g()
{
    std::cout << Absolute<6>::value << std::endl;
    std::cout << Absolute<-7>::value << std::endl;
}

标准库提供了std::conditional及其辅助类型std::conditional_t用于模板分支：

template<bool B, class T, class F>
struct conditional;

定义了成员类型type，当B == true时为T，否则为F。

模板匹配实际上是在处理switch-case的分支，bool只是其中一种简单情况。对于对应关系不太规则的分支语句，可以用一个constexpr函数把参数映射到一个整数或枚举上：

enum class Port_t
{
    PortB, PortC, PortD, PortError,
};

constexpr Port_t portMap(int pin)
{
    Port_t result = Port_t::PortError;
    if (pin < 0)
        ;
    else if (pin < 8)
        result = Port_t::PortD;
    else if (pin < 14)
        result = Port_t::PortB;
    else if (pin < 20)
        result = Port_t::PortC;
    return result;
}

template<int Pin, Port_t Port = portMap(Pin)>
struct PinOperation;

template<int Pin>
struct PinOperation<Pin, Port_t::PortB> { /* ... */ };

template<int Pin>
struct PinOperation<Pin, Port_t::PortC> { /* ... */ };

template<int Pin>
struct PinOperation<Pin, Port_t::PortD> { /* ... */ };

如果同一个模板有两个参数分别处理两种分支（这已经从分支上升到模式匹配了），或同时处理分支和循环的特化，总之有两个或以上维度的特化，需要注意两个维度的特化是否会同时满足，如果有这样的情形但没有提供两参数都特化的模板特化，编译会出错。见problem2::Accumulator，它不需要提供两个参数同时特化的版本。

循环

如前所述，循环要化为递归，循环的开始与结束是递归的起始与终点或两者对调，递归终点的模板需要特化。比如，还是计算阶乘：

template<int N>
struct Factorial
{
    static constexpr int value = N * Factorial<N - 1>::value;
};

template<>
struct Factorial<0>
{
    static constexpr int value = 1;
};

或许阶乘的递归定义很大程度上来源于数学，那就再看一个平方和的例子：

template<int N>
struct SquareSum
{
    static constexpr int value = SquareSum<N - 1>::value + N * N;
};

template<>
struct SquareSum<0>
{
    static constexpr int value = 0;
};

（(1^2 + 2^2 + cdots + n^2 = frac {n left( n + 1 ight) left( 2n + 1 ight)} {6})）

好吧，还是挺数学的，去下面看实例感觉一下吧，那里还有break——哦不，被我放到思考题中去了。

加群是交换群，求和顺序不影响结果，上面这样的顺序写起来方便。有些运算符不满足交换律，需要逆转顺序。还以平方和为例：

template<int N, int Cur = 0>
struct SquareSumR
{
    static constexpr int value = Cur * Cur + SquareSumR<N, Cur + 1>::value;
};

template<int N>
struct SquareSumR<N, N>
{
    static constexpr int value = N * N;
};

递归

递归在过程式中是一种高级的结构，它可以直接转化为函数式的递归，后面会提到两者的异同。

比如，计算平方根，这个例子来源于C++ Templates: The Complete Guide 2e：

// primary template for main recursive step
template<int N, int LO = 1, int HI = N>
struct Sqrt {
    // compute the midpoint, rounded up
    static constexpr auto mid = (LO + HI + 1) / 2;
    // search a not too large value in a halved interval
    using SubT = std::conditional_t<(N < mid * mid),
                                   Sqrt<N, LO, mid - 1>,
                                   Sqrt<N, mid, HI>>;
    static constexpr auto value = SubT::value;
};
// partial specialization for end of recursion criterion
template<int N, int S>
struct Sqrt<N, S, S> {
    static constexpr auto value = S;
};

这个递归很容易化为循环，有助于你对循环化递归的理解。

存储

实际应用中我们可能不需要把所有计算出来的值存储起来，但在打表的题目中需要。存储一系列数据需要用循环，循环的实现方式依然是递归。比如，存储阶乘（Factorial类模板见上）：

template<int N>
inline void storeFactorial(int* dst)
{
    storeFactorial<N - 1>(dst);
    dst[N] = Factorial<N>::value;
}

template<>
inline void storeFactorial<-1>(int* dst)
{
    ;
}

void h()
{
    constexpr int MAX = 10;
    int factorial[MAX + 1];
    storeFactorial<MAX>(factorial);
    for (int i = 0; i <= MAX; ++i)
        std::cout << factorial[i] << ' ';
    std::cout << std::endl;
}

多维数组同理，例子见下方。注意，函数模板不能偏特化，但有静态方法的类模板可以，这个静态方法就充当原来的模板函数。

虽然我们是对数组中的元素挨个赋值的，但编译器的生成代码不会这么做，即使不能优化成所有数据一起用memcpy，至少能做到一段一段拷贝。

类内定义的函数隐式成为inline，手动写上inline没有语法上的意义，但是对于一些编译器，写上以后函数被内联的可能性更高，所以写inline是一个好习惯。

解答

#include <iostream>
#include <algorithm>

constexpr int MAX = 50;

namespace problem1
{

template<int Count, int Num, int Max>
struct Partition;

template<int Count, int Num, int Loop>
struct Accumulator
{
    static constexpr int value = Accumulator<Count, Num, Loop - 1>::value + Partition<Count, Num - Loop, Loop>::value;
};

template<int Count, int Num>
struct Accumulator<Count, Num, 0>
{
    static constexpr int value = 0;
};

template<int Count, int Num, int Max = Num>
struct Partition
{
    static constexpr int value = Accumulator<Count - 1, Num, std::min(Num - 1, Max)>::value;
};

template<int Num, int Max>
struct Partition<1, Num, Max>
{
    static constexpr int value = Num <= Max;
};

template<int Count, int Num>
struct Store
{
    static inline void store(int* dst)
    {
        Store<Count, Num - 1>::store(dst);
        dst[Num] = Partition<Count, Num>::value;
    }
};

template<int Count>
struct Store<Count, 0>
{
    static inline void store(int* dst)
    {
        ;
    }
};

template<int Count>
inline void store(int (*dst)[MAX + 1])
{
    store<Count - 1>(dst);
    Store<Count, MAX>::store(dst[Count]);
}

template<>
inline void store<0>(int (*dst)[MAX + 1])
{
    ;
}

inline void store(int(*dst)[MAX + 1])
{
    store<MAX>(dst);
}

}

namespace problem2
{

template<int Num, int Count = Num, int Helper = Num - Count * (Count - 1) / 2, bool Valid = (Helper > 0)>
struct Accumulator
{
    static constexpr int value = Accumulator<Num, Count - 1>::value + problem1::Partition<Count, Helper>::value;
};

template<int Num, int Count, int Helper>
struct Accumulator<Num, Count, Helper, false>
{
    static constexpr int value = Accumulator<Num, Count - 1>::value;
};

template<int Num, int Helper, bool Valid>
struct Accumulator<Num, 0, Helper, Valid>
{
    static constexpr int value = 0;
};

template<int Num>
inline void store(int* dst)
{
    store<Num - 1>(dst);
    dst[Num] = Accumulator<Num>::value;
}

template<>
inline void store<0>(int* dst)
{
    ;
}

inline void store(int* dst)
{
    store<MAX>(dst);
}

}

int ans1[MAX + 1][MAX + 1];
int ans2[MAX + 1];

int main()
{
    problem1::store(ans1);
    problem2::store(ans2);
    int num, k;
    while (std::cin >> num)
    {
        std::cin >> k;
        std::cout << ans1[k][num] << std::endl;
        std::cout << ans2[num] << std::endl;
        std::cout << ans2[num] << std::endl;
    }
}

请对照运行时版本自行理解。

讨论

constexpr

constexpr不保证计算在编译期完成，大部分编译器在Debug模式下把所有可以推迟的constexpr计算都推迟到运行时完成。但constexpr可以作为一个强有力的优化提示，原本在最高优化等级都不会编译期计算的代码，在有了constexpr后编译器会尽力帮你计算。如果编译器实在做不到，根据你是否强制编译期求值，编译器会给出错误或推迟到运行时计算。在不同的编译器中，这类行为的表现是不同的——众所周知MSVC对constexpr的支持不好。

目前（C++17）没有任何方法可以检查一个表达式是否是编译期求值的，但是有方法可以让编译器对于非编译期求值表达式给出一个错误，把期望constexpr的表达式放入模板参数或static_assert表达式都是可行的：如果编译期求值，则编译通过；否则编译错误。

（C++20：consteval、is_constant_evaluated）

模板

如果我们把Sqrt中的递归替换为如下语句：

static constexpr auto value = (N < mid * mid) ? Sqrt<N, LO, mid - 1>::value
                                              : Sqrt<N, mid, HI>::value;

显然计算结果是相同的，看上去还更简洁。但是问题在于，编译器会把Sqrt<N, LO, mid - 1>和Sqrt<N, mid, HI>两个类都实例化出来，尽管只有一个模板类的value会被使用到。这些类模板实例继续导致其他实例产生，最终将产生 (O left( n log n ight)) 个实例。相比之下，把两个类型名字传给std::conditional并不会导致类模板被实例化，std::conditional只是定义一个类型别名，对该类型求::value才会实例化它，一共产生 (O left( log n ight)) 个实例。

还有一个很常见的工具是变参模板，我没有介绍是因为暂时没有用到，而且我怕写出非多项式复杂度的元程序。如果我还有机会写一篇类型元编程的话，肯定会包含在其中的。

函数式

循环的一次迭代往往需要上一次迭代的结果，对应地在递归中就是函数对一个参数的结果依赖于对其他 (n) 个参数的结果。有些问题用递归解决比较直观，但是如果 (n geq 2)，计算过程就会指数爆炸，比如：

int fibonacci(int n)
{
    if (n <= 2)
        return 1;
    else
        return fibonacci(n - 2) + fibonacci(n - 1);
}

计算fibonacci(30)已经需要一点点时间了，而计算fibonacci(46)（4字节带符号整型能容纳的最大斐波那契数）就很慢了。把这种递归转化为循环，就是设计一个动态规划算法的过程。然而函数式中的递归与过程式中的循环可能有相同的渐近复杂度：

template<int N>
struct Fibonacci
{
    static constexpr int value = Fibonacci<N - 2>::value + Fibonacci<N - 1>::value;
};

template<>
struct Fibonacci<1>
{
    static constexpr int value = 1;
};

template<>
struct Fibonacci<2>
{
    static constexpr int value = 1;
};

因为只有Fibonacci<1>到Fibonacci<46>这46个类模板被实例化，是 (O left( n ight)) 复杂度的。

在题目中，由于表中的所有数据都有可能用到，并且运行时不能执行计算，所以要把所有数据都计算出来。实际问题中可能只需要其中一个值，比如我现在就想知道不同整数的划分问题对 (50) 的答案是多少，就写：

std::cout << problem2::Accumulator<50>::value << std::endl;

那么problem1::Partition的Count参数就不会超过10，不信的话你可以加一句static_assert。实例化的模板数量一共只有2000多个，而在完整的问题中这个数量要翻100倍不止。这种性质称为惰性求值，即用到了才求值。惰性求值是必需的，总不能穷尽模板参数的所有可能组合一一实例化出来吧？

函数式编程语言可以在运行时实现这些特性。

性能

我愧对这个小标题，因为C++值元编程根本没有性能，时间和空间都是。类型元编程也许是必需，至于值元编程，emm，做点简单的计算就可以了，这整篇文章都是反面教材。

思考题2用GCC编译，大概需要10分钟；用MSVC编译，出现我闻所未闻的错误：

因为编译器是32位的，4GB内存用完了就爆了。

停机问题

一个很有趣的问题是编译器对于死循环的行为。根据图灵停机问题，编译器无法判断它要编译的元程序是否包含死循环，那么它在遇到死循环时会怎样表现呢？当然不能跟着元程序一起死循环，constexpr的循环次数与模板的嵌套深度都是有限制的。在GCC中，可以用-fconstexpr-depth、-fconstexpr-loop-limit和-ftemplate-depth等命令行参数来控制。

思考题

1. problem2::Accumulator从Count == 0到Count == Num都要实例化，但其实只需实例化到 (O left( sqrt{n} ight)) 就可以了，试改写之。

2. 洛谷 NOIp2016提高组D2T1 组合数问题，用元编程实现。

   • 只需完成 (n leq 100, m leq 100) 的任务点；

   • 使用64位编译器（指编译器本身而非目标代码），给编译器亿点点时间；

   • 不要去网站上提交，我已经试过了，编译错误。

   • 测试数据下载。

题目描述

组合数 (inom {n} {m}) 表示的是从 (n) 个物品中选出 (m) 个物品的方法数。举个例子，从 (left( 1, 2, 3 ight)) 三个物品中选择两个物品可以有 (left( 1, 2 ight), left( 1, 3 ight), left( 2, 3 ight)) 这三种选择方法。根据组合数的定义，我们可以给出计算组合数 (inom {n} {m}) 的一般公式

[inom {n} {m} = frac {n!} {m! left( n-m ight) !} \,, ]

其中 (n! = 1 imes 2 imes cdots imes n)；特别地，定义 (0! = 1)。

小葱想知道如果给定 (n)，(m) 和 (k)，对于所有的 (0 leq i leq n, 0 leq j leq min left( i, m ight)) 有多少对 (left( i, j ight)) 满足 (k mid inom {i} {j})。

输入格式

第一行有个两个整数 (t, k)，其中 (t) 代表该测试点总共有多少组测试数据，(k) 的意义见问题描述。

接下来 (t) 行每行两个整数 (n, m)，其中 (n, m) 的意义见问题描述。

输出格式

共 (t) 行，每行一个整数代表所有的 (0 leq i leq n, 0 leq j leq min left( i, m ight)) 有多少对 (left( i, j ight)) 满足 (k mid inom {i} {j})。

输入输出样例

【输入#1】

1 2
3 3

【输出#1】

1

【输入#2】

2 5
4 5
6 7

【输出#2】

0 7

说明/提示

【样例1说明】

在所有可能的情况中，只有 (inom {2} {1} = 2) 一种情况是 (2) 的倍数。

【子任务】

测试点	(n)	(m)	(k)	(t)
1	(leq 3)	$ leq 3$	(= 2)	$ = 1$
2			(= 3)	(leq 10^4)
3	(leq 7)	$ leq 7$	(= 4)	$ = 1$
4			(= 5)	(leq 10^4)
5	(leq 10)	$ leq 10$	(= 6)	$ = 1$
6			(= 7)	(leq 10^4)
7	(leq 20)	$ leq 100$	(= 8)	$ = 1$
8			(= 9)	(leq 10^4)
9	(leq 25)	$ leq 2000$	(=10)	$ = 1$
10			(=11)	(leq 10^4)
11	(leq 60)	$ leq 20$	(=12)	$ = 1$
12			(=13)	(leq 10^4)
13	(leq 100)	$ leq 25$	(=14)	$ = 1$
14			(=15)	(leq 10^4)
15		$ leq 60$	(=16)	$ = 1$
16			(=17)	(leq 10^4)
17	(leq 2000)	$ leq 100$	(=18)	$ = 1$
18			(=19)	(leq 10^4)
19		$ leq 2000$	(=20)	$ = 1$
20			(=21)	(leq 10^4)

• 对于全部的测试点，保证 (0 leq n, m leq 2 imes 10^3, 1 leq t leq 10^4)。