## Why is the answer from the calculator and matlab different?

### soe (view profile)

on 11 Apr 2019
Latest activity Edited by John D'Errico

### John D'Errico (view profile)

on 15 Apr 2019
I tried to find the answer in the following:
a=mod(455^100,17)
7
and
the answer from the calculator is
1.
Could anyone help me explain the difference?

### John D'Errico (view profile)

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

Show 1 older comment
David Goodmanson

### David Goodmanson (view profile)

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

### Bjorn Gustavsson (view profile)

on 15 Apr 2019
Yeah, integer arithmetics (don't know how to best put this into English...) blows my mind out -> "that works??? Wow! How, why?"
John D'Errico

### John D'Errico (view profile)

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.)