模乘逆

Modular multiplicative inverse

我正在计算 ((A^B)/C)%M,但是当 A、B、C、M 的数量很大时,我的代码不起作用。当 A、B、C、D 较小时,此代码给出正确答案 int

这里有什么问题?

这里C和M互质

示例输入 2 3 4 5 示例输出 2

这些输入的代码失败 969109092 60139073 122541116 75884463

C程序

#include <stdio.h>

int d,x,y;

模指数 (A^B)%M

int power(int A, int B, int M)
{
    long long int result=1;
    while(B>0)
    {
        if(B % 2 ==1)
        {
            result=(result * A)%M;
        }
        A=(A*A)%M;
        B=B/2;
    }
    return result;
}

模乘逆

void extendedEuclid(int A, int B)
{
    if(B == 0)
    {
        d = A;
        x = 1;
        y = 0;
    }
    else
    {
        extendedEuclid(B,A%B);
        int temp = x;
        x = y;
        y = temp - (A/B)*y;
    }
}

int modInv(int A, int M)
{
    extendedEuclid(A,M);
    return (x%M+M)%M;
}

main()

int main()
{
    int A,B,C,M;
    scanf("%d %d %d %d",&A,&B,&C,&M);
    int inv = modInv(C,M)%M;
    printf("%d\n",inv);
    long long int p = (power(A,B,M))%M;
    printf("%d\n",p);
    long long int ans = (p * inv)%M;
    //printf("%d",((modInv(C,M)*(power(A,B,M))))%M);
    printf("%lld",ans);
    return 0;
}

int 的值可能不够大,请尝试使用 long 或 double。

小心因为power returns an int not long long int

代码至少存在以下问题:

intA*A 中溢出。代码需要使用更广泛的数学计算乘积 A*A。这就是为什么代码适用于小值而不适用于大值。

// A=(A*A)%M;
A = ((long long)A*A) % M;
// or 
A = (1LL*A*A) % M;

打印说明符错误。这意味着编译器警告未完全启用。节省时间,全部启用。

long long int p = (power(A,B,M))%M;
// printf("%d\n",p);
printf("%lld\n",p);

代码有负值错误。不要修补 int 漏洞,而是使用 unsigned 类型。

unsigned power(unsigned A, unsigned B, unsigned M) {
  unsigned long long result = 1;
  ...

power(A,0,1) 中的极端案例失败。当 M==1.

时,result 应为 0
// long long int result=1;
long long int result=1%M;

带有评论中指出的修复的测试版本:

#include <stdio.h>

int d,x,y;

int power(int A, int B, int M)
{
    long long int result=1;
    long long int S = A;            /* fix */
    while(B>0)
    {
        if(B % 2 ==1)
        {
            result=(result * S)%M;  /* fix */
        }
        S=(S*S)%M;                  /* fix */
        B=B/2;
    }
    return (int)result;
}

void extendedEuclid(int A, int B)
{
int temp;                           /* C */
    if(B == 0)
    {
        d = A;
        x = 1;
        y = 0;
    }
    else
    {
        extendedEuclid(B,A%B);
        temp = x;
        x = y;
        y = temp - (A/B)*y;
    }
}

int modInv(int A, int M)
{
    extendedEuclid(A,M);
/*  x = x%M;                        ** not needed */
    if (x < 0)                      /* fix */
        x += M;                     /* fix */
    return (x);                     /* fix */
}

int main()
{
    int A,B,C,M;                    /* C */
    int inv, p, ans;                /* C */
    A = 969109092;                  /* 2^2 × 3^2 ×7 × 1249 × 3079 */
    B =  60139073;                  /* 60139073 */
    C = 122541116;                  /* 2^2 × 1621 × 18899 */
    M =  75884463;                  /* 3^2 × 8431607 */

    inv = modInv(C,M)%M;            /* 15543920 */
    printf("%d\n",inv);
    p = power(A,B,M)%M;             /*  6704397 */
    printf("%d\n",p);
    ans = (unsigned)(((unsigned long long)p * inv)%M);  /* fix 22271562 */
    printf("%d\n",ans);
    return 0;
}

你可以试试mod_inv C函数:

// return a modular multiplicative inverse of n with respect to the modulus.
// return 0 if the linear congruence has no solutions.
unsigned mod_inv(unsigned ra, unsigned rb) {
    unsigned rc, sa = 1, sb = 0, sc, i = 0;
    if (rb > 1) do {
            rc = ra % rb;
            sc = sa - (ra / rb) * sb;
            sa = sb, sb = sc;
            ra = rb, rb = rc;
        } while (++i, rc);
    sa *= (i *= ra == 1) != 0;
    sa += (i & 1) * sb;
    return sa;
}

这基本上是标准算法,为了避免溢出,符号存储在 d 变量中,您可以使用结构来执行此操作。此外,当 n = 1mod = 0 时,输出为 0,而不是1,我认为我们没有多少计算可以执行 modulo 0.

The modular multiplicative inverse of an integer N modulo m is an integer n such as the inverse of N modulo m equals n, if a modular inverse exists then it is unique. To calculate the value of the modulo inverse, use the extended euclidean algorithm which finds solutions to the Bezout identity.

用法示例:

#include <assert.h>

int main(void) {
    unsigned n, mod, res;

    n = 52, mod = 107;
    res = mod_inv(n, mod);
    assert(res == 35); //  35 is a solution of the linear congruence.

    n = 66, mod = 123;
    res = mod_inv(n, mod);
    assert(res == 0); // 66 does note have an inverse modulo 123.
}

/*
    n =     7 and mod = 45    then res = 13    so 1 == (   13 * 7    ) % 45   
    n =    52 and mod = 107   then res = 35    so 1 == (   35 * 52   ) % 107  
    n =   213 and mod = 155   then res = 147   so 1 == (  147 * 213  ) % 155  
    n =   392 and mod = 45    then res = 38    so 1 == (   38 * 392  ) % 45   
    n =   687 and mod = 662   then res = 53    so 1 == (   53 * 687  ) % 662  
    n =   451 and mod = 799   then res = 512   so 1 == (  512 * 451  ) % 799  
    n =  1630 and mod = 259   then res = 167   so 1 == (  167 * 1630 ) % 259  
    n =  4277 and mod = 4722  then res = 191   so 1 == (  191 * 4277 ) % 4722
*/

Source