fast way to calculate a^(b!)mod(c)

Question

Here is what I've tried so far. I don't know what is wrong with this code, but it's giving the wrong answer for large cases: x,y,n > 10^5.

#include <stdio.h>
long long factorial(long long N)
{
  long long product = 1;
  for ( long long j=1;j<=N;j++)
    product *= j;
  return product;
}
long long power(long long x, unsigned long long y, long long p) 
{ 
    long long res = 1;      // Initialize result 

    x = x % p;  // Update x if it is more than or  
                // equal to p 

    while (y > 0) 
    { 
        // If y is odd, multiply x with result 
        if (y & 1) 
            res = (res*x) % p; 

        // y must be even now 
        y = y>>1; // y = y/2 
        x = (x*x) % p;   
    } 
    return res; 
} 
int main()  
{  
   long long A,B,C,test;
   scanf("%lld",&test);
   for (long long i = 0;i<test;i++)
   {
       scanf("%lld %lld %lld",&A,&B,&C);
        long long fa = factorial(B);
        long long res = power(A,fa,C);
        printf("%lld \n",res);
   }
   return 0;  
}

Any help or walkthrough will be appreciated.

Calculating `b!` that way won't help. This is related to modular arithmetic. So, all we would need in the end is some value between `0` and `c-1`. — nice_dev, Dec 17 '18 at 14:05
Have a read about how to do this in a memory efficient manner (so you won't experience overflow given your modulus isn't enormous): https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method — Matt Messersmith, Dec 17 '18 at 14:25
The largest factorial you can rely on fitting in a long long is 20! (2432902008176640000). — Andrew Morton, Dec 17 '18 at 14:45
You don't need to calculate the actual factorial. Since `a^b mod c = a^(b mod phi(c)) mod c`. I.e., calculate the factorial up to `mod phi(c)` where `phi` is [Euler's totient function](https://en.wikipedia.org/wiki/Euler%27s_totient_function). — Nico Schertler, Dec 17 '18 at 14:46
As for why it's giving the wrong answers: You say that n > 10^5, 10000! is around 10^456573. Compare that to the maximum value supported by `long long`. — Raymond Chen, Dec 17 '18 at 15:21
A duplicate maybe [Calculating (a^b)%MOD](https://stackoverflow.com/questions/11272437/calculating-abmod)? — gsamaras, Dec 18 '18 at 09:27

Spektre · Accepted Answer · 2018-12-18T09:38:43.007

well apart exploiting some knowledge from number theory like Nico Schertler suggests in his comment you can do this on small ints with naive approach too.

Your problem is that your sub result is not fitting into your variables:

a^(b!) mod c < c
b! > c

The factorial result is really huge and will not fit into your small int variables without the use of bigints. But you can change the computation slightly to this:

  10^(3!) mod c
= 10^(1*2*3) mod c
= (((10^1 mod c)^2 mod c)^3 mod c)

So the idea is instead of computing the factorial by multiplying i iterator of a loop use modpow on the a^... sub result. This way all the sub results still fits into your variables.

Here a small C++ code for this:

//---------------------------------------------------------------------------
typedef unsigned __int32 DWORD;
//---------------------------------------------------------------------------
DWORD mod(DWORD a,DWORD p)
    {
    DWORD bb;
    for (bb=p;(DWORD(a)>DWORD(bb))&&(!DWORD(bb&0x80000000));bb<<=1);
    for (;;)
        {
        if (DWORD(a)>=DWORD(bb)) a-=bb;
        if (bb==p) break;
        bb>>=1;
        }
    return a;
    }
//---------------------------------------------------------------------------
DWORD modadd(DWORD a,DWORD b,DWORD p)
    {
    DWORD d,cy;
    a=mod(a,p);
    b=mod(b,p);
    d=a+b;
    cy=((a>>1)+(b>>1)+(((a&1)+(b&1))>>1))&0x80000000;
    if (cy) d-=p;
    if (DWORD(d)>=DWORD(p)) d-=p;
    return d;
    }

//---------------------------------------------------------------------------
DWORD modsub(DWORD a,DWORD b,DWORD p)
    {
    DWORD d;
    a=mod(a,p);
    b=mod(b,p);
    d=a-b; if (DWORD(a)<DWORD(b)) d+=p;
    if (DWORD(d)>=DWORD(p)) d-=p;
    return d;
    }
//---------------------------------------------------------------------------
DWORD modmul(DWORD a,DWORD b,DWORD p)
    {
    int i;
    DWORD d;
    a=mod(a,p);
    for (d=0,i=0;i<32;i++)
        {
        if (DWORD(a&1))    d=modadd(d,b,p);
        a>>=1;
        b=modadd(b,b,p);
        }
    return d;
    }
//---------------------------------------------------------------------------
DWORD modpow(DWORD a,DWORD b,DWORD p)
    {
    int i;
    DWORD d=1;
    mod(a,p);
    for (i=0;i<32;i++)
        {
        d=modmul(d,d,p);
        if (DWORD(b&0x80000000)) d=modmul(d,a,p);
        b<<=1;
        }
    return d;
    }
//---------------------------------------------------------------------------
DWORD modpowfact(DWORD a,DWORD b,DWORD c) // returns a^(b!) mod c
    {
    DWORD i,y=mod(a,c);
    for (i=2;i<=b;i++)
     y=modpow(y,i,c);
    return y;
    }
//---------------------------------------------------------------------------

the code above returns these on my setup:

10^(0!) mod 1031 = 10
10^(1!) mod 1031 = 10
10^(2!) mod 1031 = 100
10^(3!) mod 1031 = 961
10^(4!) mod 1031 = 72
10^(5!) mod 1031 = 754

the modular arithmetic was taken from here:

Modular arithmetics and NTT (finite field DFT) optimizations

I used the slow non optimized ones so it can be compiled and run on any platform ...

score -1 · Answer 2 · answered Dec 17 '18 at 21:15

check this answer!

https://stackoverflow.com/a/26282858/7480909

I use this function to solve this problem

UVA 374 - Big Mod

http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=310
// a^b % T

// with Exponentiation by squaring (http://en.wikipedia.org/wiki/Exponentiation_by_squaring#Basic_method)

// if a very large use
// R=(unsigned long long)(R*a)%T;

int restOfPot(long long a,int b,int T)
{
    a%=T;
    long long R=1;

    while(b)
    {
        if(b&1) R=(R*a)%T;
        a=(a*a)%T;
        b>>=1;
    }

    return int(R);
}

This should be a comment...and it doesn't solve the problem. The problem is the OP is calling `factorial(B)` and then *passing that value around*. A good answer would show how to *not* do that. — Matt Messersmith, Dec 17 '18 at 21:46

fast way to calculate a^(b!)mod(c)

2 Answers2