Asked by soe
on 11 Apr 2019

I tried to find the answer in the following:

a=mod(455^100,17)

The answer from MATLAB is

7

and

the answer from the calculator is

1.

Could anyone help me explain the difference?

Answer by John D'Errico
on 11 Apr 2019

Because 455^100 is a computation done in DOUBLE PRECISION. Assuch, it goes well beyond the limits of what a double can represent exactly.

455^100

ans =

6.3262e+265

But we can do it using extended computations, using the symbolic toolbox, or my own VPI toolbox.

sym(455)^100

ans =

63261526279114425171518999214678518568130878514925691067972899819122977043674229777870415964659361018260513068441824908812079611837251889975323134728654421081265983719189871544107740294078906928547031851680435583183069234007069037151183010792010463774204254150390625

MATLAB computations using a double are NOT carried out to full precision. If they were, you would wait forever to do some simple computations. So a double can represent integers as large as 2^53-1 exactly. (Actually, it gets 2^53 correct also.) But beyond that point, that number will not be represented correctly, exactly, as an integer. And of course 455^100 is well beyond that point, almost to the point where it will overflow into an inf.

The symbolic toolbox has a powermod function in it though, if you have that TB. (VPI also supplies a powermod utility.)

powermod(sym(455),100,17)

ans =

1

David Goodmanson
on 14 Apr 2019

Hi soe,

If you are ok with sticking to the limits of double precision, you can go with

(a*m + b)^N (mod m) = b^N (mod m)

since expanding out the left hand side gives a bunch of coefficients times positive powers of m, all of which go away mod m, leaving only the term b^N.

% find mod(455^100,17)

mod(455,17)

ans = 13 % = b

% find mod(13^100,17), still too big

% find mod((13^10)^10,17)

mod(13^10,17) % this works within the limits of double precision

ans = 16

mod(16^10,17) % this also works within the limits of double precision

ans = 1

Same basic idea as Bjorn's.

Bjorn Gustavsson
on 15 Apr 2019

John D'Errico
on 15 Apr 2019

The powermod trick works because computation of a modulus follows a few simple rules. The most important is:

mod(a*b,p) = mod(mod(a,p) * mod(b,p), p)

So you can compute a modulus of a power easily. For example, suppose we wanted to compute mod(455^100,17)? There is an easy way, and a not so easy way. Lets do it the hard way, first.

a = 455;

p = 17;

d = 100;

M = 1;

for i = 1:100

M = mod(M*a,p);

end

M

M =

1

And that indeed is the modulus mod(455^100,17).

Easy, peasy. But things would have been more nasty, had the exponent d been something like 1e15. Thus what is mod(455^1e15,17)?

So how might I have done this differently? We first compute the binary expansion of 100.

find(flip(dec2bin(100)- '0')) - 1

ans =

2 5 6

That is, 100 = 2^2 + 2^5 + 2^6 = 4 + 32 + 64.

What does that tell us about 455^100?

455^100 = 455^4 * 455^32 * 455^64

Does that help us in any way? Actually, it does. In fact, it allows us to compute the desired power trivially and on paper. That is because 455^4 is not too large of a number that it overflows integer computations with a double. So we see this:

mod(455^4,17)

ans =

1

Is that nice? It is if we recognize that

455^32 = (455^4)^8

So if mod(455^4,17) == 1, then mod(455^32,17) = mod(1^8,17) = 1

Likewise, 455^64 = (455^32)^2.

So trivially, we find that

mod(455,^100,17) = mod(1*1*1,17) = 1.

Now, could I have done this differently, but just as easily, IF we understand how the binomial expansion applies to mod computations. Thus:

mod(a + k*p)^n,p) = mod(a^n,p)

What does that tell us?

mod(455,17)

ans =

13

So 455 is congruent to 13, modulo 17. Therefore, 455 is also congruent to -4, modulo 17. So the binomial trick I showed above tells us that

mod(455^100,17) = mod((-4)^100,17) = mod(4^100,17)

The negative goes away there since the exponent is even. Next, we might notice that

mod(4^10,17)

ans =

16

But 16 is also minus 1, mod 17. That quickly brings us to the conclusion that

mod(455^100,17) = mod((-1)^10,17) = mod(1,17) = 1

Again, no large scale computation was ever needed.

That was an easy case. For larger exponents, or where we need to do more than computer powers of the number 1, we do it in very much the same way. Consider d=1e15? As long as we can compute a binary expansion of the exponent, there is no problem at all.

binexp = find(flip(dec2bin(1e15)- '0')) - 1

binexp =

15 17 18 22 23 26 29 31 33 34 35 36 37 38 40 42 43 47 48 49

That is,

1e15 = 2^15 + 2^17 + 2^18 + ... + 2^48 + 2^49

sum(2.^binexp) == 1e15

ans =

logical

1

Now, we compute 455^1e15 as just a product of only 20 terms, done in a loop.

numel(binexp)

ans =

20

Is this an important thing? Well, actually, it can be terribly useful. There are a variety of tests to identify the primality of some massively HUGE numbers.

For a simple one, consider this number:

X = sym('9999999999999999999999999999999999999999594780488384396490825872803373354740335455085706895782052009');

Is X a prime number? Well, as it turns out, it is, and isprime will say so.

isprime(X)

ans =

logical

1

We can also use some other tests, like the Fermat test, which uses powermod.

X

X =

9999999999999999999999999999999999999999594780488384396490825872803373354740335455085706895782052009

>> powermod(2,X-1,X)

ans =

1

If the above test results in 1, then this is evidence for the primality of X.

However, the Fermat test is not a 100% test for primality, and there are many things you can read about what are called Fermat liars or Fermat pseudo-primes.

In the end, why do we care about huge prime numbers? Well, they form the foundation of many encryption schemes.

(I think this has diverged a bit from the original question, but it happens to touch on a topic of personal interest to me, that of Quasi-Woodall primes. So I could go on here for hours.)

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.