题意

定义特殊的加法 R(x3,x3)=P(x1,x2)+Q(x1,x2)R(x_3,x_3) = P(x_1,x_2) + Q(x_1,x_2)R(x3,x3)=P(x1,x2)+Q(x1,x2)

x3=k2x1x2(<mtext> </mtext>mod<mtext> </mtext>p)x_3 = k^2-x_1-x_2(\bmod p)x3=k2x1x2(modp)

y3=k(x1x3)y1(<mtext> </mtext>mod<mtext> </mtext>p)y_3 = k(x_1-x_3) - y_1 (\bmod p)y3=k(x1x3)y1(modp)

其中,若PQP\neq QP=Qk=y2y1x2x1(<mtext> </mtext>mod<mtext> </mtext>p)k = \frac{y_2-y_1}{x_2-x_1} (\bmod p)k=x2x1y2y1(modp),否则k=3x12+12y1(<mtext> </mtext>mod<mtext> </mtext>p)k = \frac{3x_1^2+1}{2y_1}(\bmod p)k=2y13x12+1(modp)

现在给P0(x,y)P_0(x,y)P0(x,y)nnn

nnnP0P_0P0相加的和P0+P0+...+P0P_0+P_0+...+P_0P0+P0+...+P0,其中相加的加法用上述定义的加法

方法

模拟/实现

因为这里是特殊的加法,所以乘法也就不再是数乘,nP0(x,y)(nx,ny)nP_0(x,y) \neq (nx,ny) nP0(x,y)=(nx,ny)

把题目给的递推式文字转换为代码, 我们可以通过for循环,相加点实现计算,

代码

/**
 * struct Point {
 *	int x;
 *	int y;
 * };
 */

class Solution {
public:
    const int mod = 1000000007;
    long long add(long long v1,long long v2){ // 带模运算的加法
        return ((v1+v2)%mod+mod)%mod; 
    }
    long long mypow(long long v,long long p){ // 带模运算的数值快速幂
        long long r = 1;
        while(p){
            if(p%2)(r*=v)%=mod;
            (v*=v)%=mod;
            p/=2;
        }
        return r;
    }
    Point add(Point P1,Point P2){ // 按照上面数学公式直接翻译的点加法
        long long k ;
        if(P1.x == P2.x && P1.y == P2.y){
            // 注意overflow
            k = add(P1.x*(long long)P1.x%mod*3,1)*mypow(2*P1.y%mod,mod-2)%mod; // $k = \frac{3x_1^2+1}{2y_1}
        }else{
            k = add(P2.y,-P1.y)*mypow(add(P2.x,- P1.x),mod-2)%mod; // 若$P\neq Q$ 则 $k = \frac{y_2-y_1}{x_2-x_1} (\bmod p)$
        }
        int x3 = add(k*k,-P1.x-P2.x); // $x_3 = k^2-x_1-x_2(\bmod p)$
        int y3 = add(k*add(P1.x,-x3),-P1.y); // $y_3 = k(x_1-x_3) - y_1 (\bmod p)$
        return Point{x3,y3};
    }
    /**
     * 
     * @param P Point类 
     * @param n int整型 
     * @return Point类
     */
    Point NTimesPoint(Point P, int n) {
        Point r = P;
        for(int i = 1;i<n;i++){ // n个点逐个相加
            r = add(r,P);
        }
        return r;
    }
};

复杂度分析

时间复杂度: 因为nnn的范围达到了10910^9109,并不能在时间复杂度内完成, O(n)O(n)O(n)

空间复杂度: 只用记录当前加和的值O(1)O(1)O(1)

快速幂

对于普通的加法,同一个数多次相加,我们可以用乘法一次完成计算

对于普通的乘法,同一个数多次相乘,我们可以用快速幂O(log(n))O(log(n))O(log(n))次完成计算

因为这里是特殊的加法,我们不具备一个特殊的乘法来一次完成计算

所以套用普通乘法的思路,用快速幂的形式去计算,也就是如下表格的拆分,每个数表示成其二进制形式的整幂次和,然后依次计算20P,21P,22P,23P,<mtext> </mtext>,2iP2^0P,2^1P,2^2P,2^3P,\cdots,2^iP20P,21P,22P,23P,,2iP ,在对应位满足时,将当前计算的内容加进结果里,这样只需要O(log(n))O(log(n))O(log(n))的时间复杂度即可计算出nPnPnP的值

n 1 2 3 4 5 6 7 8 10 ...
P 2P 2P+P 4P 4P+P 4P+2P 4P+2P+P 8P 8P+P ...

图解:

假设我们要计算7P

- r(结果) p(幂次) v(基数) 对应
初始化 1 7 vvv 未计算的状态
幂次为7,模2余1,乘上基数 1v=v1 \cdot v =v1v=v 7/2=3 (v)2=v2(v)^2=v^2(v)2=v2 +1P
幂次为3,模2余1,乘上基数 vv2=v3v \cdot v^2 = v^3vv2=v3 3/2=1 (v2)2=v4(v^2)^2=v^4(v2)2=v4 +2P
幂次为1,模2余1,乘上基数 v3v4=v3v^3 \cdot v^4 = v^3v3v4=v3 1/2=0 (v4)2=v8(v^4)^2=v^8(v4)2=v8 +4P
幂次为0 结束 - - - -

代码

/**
 * struct Point {
 *	int x;
 *	int y;
 * };
 */

class Solution {
public:
    const int mod = 1000000007;
    long long add(long long v1,long long v2){ // 带模运算的加法
        // Non-negative add result
        return ((v1+v2)%mod+mod)%mod;
    }
    long long mypow(long long v,long long p){ // 带模运算的快速幂
        long long r = 1;
        while(p){ // 见上面图解的快速幂
            if(p%2)(r*=v)%=mod;
            (v*=v)%=mod;
            p/=2;
        }
        return r;
    }
    Point add(Point P1,Point P2){ // 按照上面数学公式直接翻译的点加法
        long long k ;
        if(P1.x == P2.x && P1.y == P2.y){
            // 注意overflow
            k = add(P1.x*(long long)P1.x%mod*3,1)*mypow(2*P1.y%mod,mod-2)%mod; // $k = \frac{3x_1^2+1}{2y_1}(\bmod p)$
        }else{
            k = add(P2.y,-P1.y)*mypow(add(P2.x,- P1.x),mod-2)%mod; // 若$P\neq Q$ 则 $k = \frac{y_2-y_1}{x_2-x_1} (\bmod p)$
        }
        int x3 = add(k*k,-P1.x-P2.x); // $x_3 = k^2-x_1-x_2(\bmod p)$
        int y3 = add(k*add(P1.x,-x3),-P1.y); // $y_3 = k(x_1-x_3) - y_1 (\bmod p)$
        return Point{x3,y3};
    }
    /**
     * 
     * @param P Point类 
     * @param n int整型 
     * @return Point类
     */
    Point NTimesPoint(Point P, int n) {
        Point r = P;
        n-=1;
        while(n){ // 见上面图解的快速幂
            if(n%2)r=add(r,P);
            P=add(P,P);
            n/=2;
        }
        return r;
    }
};

需要注意的点

  1. PointPointPointx,yx,yx,y字段是int类型,做乘法可能溢出,所以需要转换成longlonglong longlonglong参与计算,计算过程中注意保持加法结果的点是取模后的非负值

  2. 因为是特殊的加法,初始结果不能取(0,0)(0,0)(0,0) 因为(0,0)+PP(0,0) + P \neq P(0,0)+P=P, 所以这里采取nP=P+(n1)PnP=P+(n-1)PnP=P+(n1)P初始结果是PPP,剩余n1n-1n1PPP来计算

复杂度分析

时间复杂度: 上面分析了,我们对幂次的值按二进制分解,分解的个数就是需要计算的轮次数,O(log(n))O(log(n))O(log(n))

空间复杂度: 因为空间上只需要记录当前的计算结果 和 2i2^i2iPPP的和,所以是常数级别的空间复杂度O(1)O(1)O(1)

知识点

乘法模逆元

题目中要求分数的模其意义是,如果

ab=1(<mtext> </mtext>mod<mtext> </mtext>p)a \cdot b = 1 (\bmod p)ab=1(modp)

那么bbbaaa互为模逆元

有费马小定理ap1=1(<mtext> </mtext>mod<mtext> </mtext>p)a^{p-1}=1(\bmod p)ap1=1(modp)

所以 aap2=1(<mtext> </mtext>mod<mtext> </mtext>p)a\cdot a^{p-2} = 1 (\bmod p)aap2=1(modp)

所以我们可以通过ap2(<mtext> </mtext>mod<mtext> </mtext>p)a^{p-2} (\bmod p)ap2(modp) 来计算aaa的模逆元

快速幂

不论是处理特殊加法还是这里计算模逆元,都会用到快速幂,快速幂的核心代码

        while(power){
            if(power%2) result = mul(result,base);
            base = mul(base,base);
            power/=2;
        }