51nod - 1436 方程的解数

请计算这个方程组有多少合法的整数解，答案比较大，对m取余后输出。

对于样例，有三组解{1, 1}, {3, 1}, {1, 3}。

 

Input

单组测试数据。
第一行包含四个整数 n, k, l, m (2 ≤ n ≤ 10^18, 0 ≤ k ≤ 10^18, 0 ≤ l ≤ 64, 1 ≤ m ≤ 10^9 + 7)。

Output

对于每一组数据输出答案占一行。

Input示例

2 1 2 10

Output示例

3

思路：

此题由于方程中仅含位操作，所以首先可以考虑将k以二进制的形式分开来看每一位。
假设当前考虑的是k的第i位，那么这一位将由且仅由n个ai的第i位决定，其中，n个数的总情况数为2^n。

若k的第i位为0，那么a1, a2...an中，在第i位上必然不存在相邻的1，设所有满足不存在相邻1的情况数为x；

那么满足使k的第i位为1的情况数就是2^n - x，我们设y = 2^n - x。

于是若k在二进制的0 ~ L位中，存在p个1和q个0，ans即为 y^p * x^q。

对于解x的值，当n = 1，x = 2；当n = 2，x = 3；当n更大的时候，可以考虑第n位的两种情况，若第n位为0那么前n-1位只要满足“不存在相邻1”即可；若第n位为1，那么第n-1位必须为0，而前n-2位满足“不存在相邻1”即可。

于是得到转移方程：dp[n] = dp[n - 1] + dp[n - 2]，但是n比较大，无法O(n)递推，但是很庆幸这个转移方程就是Fibonacci数列，很容易想到使用矩阵快速幂进行加速求得某项的值，于是此题得解。

以上内容写的比较啰嗦...

但是友情提示一句，需要特判考虑无解的情况：）

 

#include <iostream>  
#include <cstring>  
#include <algorithm>  
#include <cstdio>  
using namespace std;  

typedef unsigned long long ull;
typedef long long ll;

struct Node  
{  
    ll val[3][3];  
};  

int m;  
ll n, k, mod;  

inline Node mul(Node a, Node b)  
{  
    Node result;  
    for (int i = 1; i <= 2; i++)  
	{
        for (int j = 1; j <= 2;j++)  
        {  
            result.val[i][j] = 0;  
            for (int k = 1; k <= 2; k++)  
			{
                 result.val[i][j] = (result.val[i][j] + a.val[i][k] * b.val[k][j] % mod) % mod;             
			}
        }  
	}

    return result;  
}

inline Node del(Node tmp, ll p)  
{  
    Node result;  
    result.val[1][1] = 1;
	result.val[2][2] = 1; 
	result.val[1][2] = 0;
	result.val[2][1] = 0;  
    while (p)  
    {  
        if (p & 1) 
		{
			result = mul(result, tmp);  
		}
        tmp = mul(tmp, tmp);  
        p >>= 1;  
    }  

    return result;  
}  

inline ll prim(ll a, ll b)  
{  
    ll result = 1;  
    while (b)  
    {  
        if (b & 1) 
		{
			result = (result * a) % mod;  
		}
        a = (a * a) % mod;  
        b >>= 1;  
    }  

    return result;  
}  

int main()  
{  
    Node tmp;  
	cin >> n >> k >> m >> mod;
    ull t = 1ULL << m;  
    if (mod == 1 || (k >= t && m != 64))  
    {  
		cout << 0 << endl;
		return 0;  
    }  
    tmp.val[1][1] = 1;
	tmp.val[1][2] = 1;
	tmp.val[2][1] = 1; 
	tmp.val[2][2] = 0;  
    tmp = del(tmp, n);  
    ll x = (tmp.val[1][1] + tmp.val[2][1]) % mod;  
    ll y = (prim(2, n) % mod - x % mod + mod) % mod;  
    int num1 = 0, num0 = 0;  
    while (k)  
    {  
        if (k & 1) 
		{
			num1++;  
		}
        k >>= 1;  
    }  
    num0 = m - num1;  
    ll result = (prim(x, num0) % mod * prim(y, num1) % mod) % mod;  
    cout << result << endl;  

	return 0;  
}