题解 | #牛牛的函数2#

思路：

题目的主要信息：

定义函数 $f(x) = (a+1)*x^a+(a+2)*x^{a+1}+...+b*x^{b-1}$
已知 $a$ 与 $b$ ，求 $f(n)\%10000000033$

方法一：暴力解法（超时）
具体做法：
写一个循环算幂的函数，然后遍历 $a$ 到 $b$ ，将计算的幂结果按照公式相乘再相加并取模。
因为 $0<a<=b-1<=10^{20}$ ，超出了long long的表示范围，因此该方法就算不超时也会有部分过大的案例无法通过。

class Solution {
public:
    long long cal(int x, int y){ //计算x的y次幂
        long long res = 1;
        for(int i = 0; i < y; i++)
            res *= x;
        return res;
    }
    long long solve(string sa, string sb, int n) {
        long long a = atoll(sa.c_str()); //转换long long型，不完全
        long long b = atoll(sb.c_str());
        long long res = 0;
        for(int i = a; i < b; i++)
            res = (res + (i + 1) * cal(n, i)) % 10000000033; //公式叠加取模
        return res;
    }
};

复杂度分析：

时间复杂度： $O(b^2-a^2)$ ，内外循环加起来一共运算次数是 $a+a+1+a+2+...+b-1+b=(a+b)*(b-a+1)/2$
空间复杂度： $O(1)$ ，常数空间

方法二：数列求和+快速幂+快速乘法+逆元
具体做法：
这道题的数据很大，我们不仅要考虑超时的问题，还应该考虑运算时爆long long的问题，因此我们需要多次优化：

对于string类型表示数大于long long表示范围，我们采用__int128_t来表示数字
对于相加运算次数，我们可以用数列求和解决
对于幂运算次数，我们可以用快速幂解决
对于运算时爆long long的问题，我们可以用快速乘法，以及费马小定理的逆元解决

数列求和。 $f(x) = (a+1)*x^a+(a+2)*x^{a+1}+...+b*x^{b-1}$ 这是一个差比数列，我们 $xf(x) = (a+1)*x^{a+1}+(a+2)*x^{a+2}+...+b*x^{b}$ ，令第一个式子减去第二个式子得： $(1-x)f(x)=(a+1)*x^a+x^{a+1}+x^{a+2}+...+x^{b-1}-b*x^b$ ，经过整理求和后我们有： $f(x)=(b*x^{b+1}-(b+1)*x^b+(a+1)*x^a-a*x^{a+1})/(x-1)^2$ ，我们的任务就从循环相加变成了计算上述式子。
快速幂。如果我们要计算 $5^{10}$ ，常规的算法是 $5*5=25$ ，然后再 $25*5=125$ ，如此往下，一共是 $9$ 次运算，即 $n-1$ 次。但是我们可以考虑这样： $5*5=25$ (二次）、 $25*25=625$ （四次）、 $625*625=...$ （八次），这是一个二分的思维，运算次数缩减到了 $log_2n$ 次，公式如下：
快速乘法。直接计算 $x^a$ 会超出long long的表示范围，因此我们可以考虑用加法来代替乘法，并在这其中取模。就比如 $a*b=a*(b_1+b_2+b_3+...)$ ，其中 $b_i$ 是数字 $b$ 的二进制各位，假设 $a=5$ ， $b=110101$ ，我们有 $a*b=a*(100000*1+10000*1+1000*0+100*1+10*0+1*1)$ ，如下表所示可以换成加法运算并在加法中取模：
逆元。为了计算 $(b*x^{b+1}-(b+1)*x^b+(a+1)*x^a-a*x^{a+1})/(x-1)^2$ ，我们用到了一个数学定理（费马小定理）：如果 $p$ 是一个质数，而整数 $x$ 不是 $p$ 的倍数，则有 $x^{p-1} \% p \equiv 1$ ，我们取的模10000000033正式一个质数。上述公式两边同时除 $x^2$ 得 $x^{p-3} \% p \equiv 1/x^2 \% p$ ，如此我们有 $y/x^2 \% mod \equiv y \% mod * 1/x^2 \% mod$
$\equiv y \% mod * x^{mod-3} \% mod$ ，其中 $x$ 与 $y$ 对应上述计算公式中的： $y=b*x^{b+1}-(b+1)*x^b+(a+1)*x^a-a*x^{a+1}$ , $x=(x-1)$ ，这样我们就可以用逆元的方式计算上述公式。

需要注意当 $n=1$ 时， $f(1)=(a+b)*(b-a+1)/2$ ，要单独讨论，这个也可以用逆元来解， $(a+b)*(b-a+1)/2 \equiv (a+b)*(b-a+1) \% mod * (2)^{mod-2} \% mod$

class Solution {
public:
    const __int128_t mod = 10000000033;
    __int128_t transform(string s){ //字符串转换成数字，超出long long范围使用__int128_t
        __int128_t res = 0, x = 1;
        for(int i = s.length() - 1; i >= 0; i--){
            res += x * (s[i] - '0');
            x *= 10;
        }
        return res;
    }
    __int128_t fast(__int128_t x, __int128_t y){ //快速乘法
        __int128_t res = 0;
        x %= mod;
        y %= mod;
        while(y){
            if(y & 1){
                res += x; //加法代替乘法，防止爆long long
                if(res >= mod)
                    res -= mod;
            }
            y = y >> 1;
            x = x << 1;
            if(x >= mod)
                x -= mod;
        }
        return res;
    }
    __int128_t Pow(__int128_t x, __int128_t y){ //快速幂
        __int128_t res = 1;
        while(y){
            if(y & 1) //可以再往上乘一个
                res = fast(res, x);
            x = fast(x, x); //叠加
            y = y >> 1; //减少乘次数
        }
        return res;
    }
    long long solve(string sa, string sb, int n) {
        if(n == 0)
            return 0;
        __int128_t a = transform(sa); //转换
        __int128_t b = transform(sb);
        if(n == 1)
            return fast(fast(a + b + 1, b - a), Pow(2, mod - 2));
        //根据公式推算
        __int128_t n_pow_a = Pow(n, a);
        __int128_t n_pow_b = Pow(n, b);
        __int128_t temp1 = (fast(fast(n_pow_a, n), a) - fast(n_pow_a, a + 1) + mod) % mod;
        __int128_t temp2 = (fast(fast(n_pow_b, n), b) - fast(n_pow_b, b + 1) + mod) % mod;
        return fast((temp2 - temp1 + mod) % mod, Pow(n - 1, mod - 3)); //逆元
    }
};

复杂度分析：

时间复杂度： $O(L+log_2(log_2n))$ ，其中 $L$ 是两个字符串长度的较长者，为字符串转换数字的时间复杂度，这里的 $n=max(b,a,mod-2)$ ，不管是快速幂还是快速乘法都是 $O(log_2n)$ ，这里嵌套使用
空间复杂度： $O(1)$ ，常数个变量，无额外空间使用