Why is the answer from the calculator and matlab different?

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?

[SWERA IRFAN]

WHY MATLAB ANSWER IS DIFFERENT FROM CALCULATOR?

Answer 1

[John D'Errico]

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

[Bjorn Gustavsson]

If I recall correctly it follows from mod's compatibility with multiplication that you can calculate mod of all powers of 455 to 455^n something like this:

M = 1;
for i1 = 1:n
  M(i1+1) = mod(M(i1)*455,17);
end

If that is the case then in addition to John's answer it means (possibly) that your calculator parses expressions differently than matlab does. Your calculator might have a grammar that implements this as the expression for "mod of a power of an integer" while matab might parse the input argument as an expression to be calculated sent to the function mod.

[David Goodmanson]

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]

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]

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.

https://en.wikipedia.org/wiki/Fermat_primality_test

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