integer overflow in matrix exponentiation

Question

I have been trying to produce several sequences with require matrix exponentiation. My code produces correct values for small cases.I have provided below 2 different matrix which I have been using for 2 different sequences.

Matrix1

[2  1 -2 -1]
[1  0  0  0]
[0  1  0  0]
[0  0  1  0]

Martrix2

[4 1]
[-4 0]

I have tried exponentiation before with non-negative values in the matrix,and was getting the desired result.

To accommodate large value of result I am using MOD=1000000007.

As far as I have figured out, the problem arises when one of the values in the matrix is reduced by MOD and others don't,that in turn makes the final result wrong.Also negative values are also getting large so I am also not sure what to do with them, and my final answer is negative for many values.

I would appriciate any link/resource/help on this issue.
Here is the code, it generates Fibonacci product.

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <vector>
#include <string>
#include <cstring>
#include <cmath>
#include <list>
#include <map>
#include<conio.h>
#include <cassert>
#define pb push_back
#define mp make_pair
#define all(x) x.begin(),x.end()
#define rep(i,n) for(int i=0;i<n;i++)
#define fill(x,y) memset(x,y,sizeof(x))
#define Size(x) (int(x.size()))
#define FOR(k,a,b) for(typeof(a) k=(a); k < (b); ++k)
typedef long long ll;
using namespace std;
#define dim 4


struct matrix
{
    int a[dim][dim];
};

#define MOD 1000000007

matrix mul(matrix x, matrix y)
{
    matrix res;
    FOR(a, 0, dim) FOR(b, 0, dim) res.a[a][b] = 0;
    FOR(a, 0, dim) FOR(b, 0, dim) FOR(c, 0, dim)
    {
        res.a[a][c] += (x.a[a][b] * ll(y.a[b][c])) % MOD;
        res.a[a][c] %= MOD;
    }
    return res;
}

matrix power(matrix m, int n)
{
    if (n == 1) return m;
    matrix u = mul(m, m);
    u = power(u, n / 2);
    if (n & 1) u = mul(u, m);
    return u;
}

matrix M, C, D, RP, A;

int main()
{
    FOR(a, 0, dim) FOR(b, 0, dim) M.a[a][b] = 0;
    M.a[0][0] = 2;
    M.a[1][0] = 1;
    M.a[2][0] = -2;
    M.a[3][0] = -1;
    M.a[0][1] = M.a[1][2] = M.a[2][3] = 1;
    C.a[0][0] = 96;
    C.a[0][1] = 44;
    C.a[0][2] = 18;
    C.a[0][3] = 5;
    int nt;
    scanf("%d", & nt);
    while (nt--)
    {
        int n;
        scanf("%d", & n);
        if (n <= 5)
        {
            if (n == 5) printf("96\n");
            if (n == 4) printf("44\n");
            if (n == 2) printf("5\n");
            if (n == 3) printf("18\n");
            if (n == 1) printf("2\n");
            continue;
        }
        if (n > 5)
        {
            int rs = n - 5;
            RP = power(M, rs);
            A = mul(C, RP);
            printf("%d\n", A.a[0][0]);
        }
    }
    getch();
    return 0;
}

Using modulo to reduce your numbers to avoid overflow in a matrix exponentiation is just mathematically wrong. You can't get meaningful results that way. — nneonneo, Sep 08 '12 at 05:02
@nneonneo: You're wrong -- look up "finite field arithmetic". If `MOD` is prime and only non-negative numbers in the range `0 .. MOD-1` are involved in a series of additions and multiplications, the answer will be "correct" (i.e. the same as doing the operations to infinite precision, and then reducing the result modulo `MOD`). The problem here is that negative numbers are involved -- TTBOMK, these can't be used directly with FFA. — j_random_hacker, Sep 08 '12 at 07:21
@j_random_hacker: I understand finite field arithmetic. However, note that in modular arithmetic mod *any positive integer*, any sequences of additions and multiplications yields a correct result modulo your modulus. Nothing here prevents negative integers from being used, and the result is *still correct*. But, depending on your interpretation of the field's equivalence classes, they might come out 'positive' -- they're still equivalent. You don't need to arbitrarily restrict the input ranges -- 9999 * 9999 == 9 * 9 == 1 mod 10, for example. — nneonneo, Sep 08 '12 at 07:36
@nneonneo: Interesting, thanks. In that case, what is the design flaw as you see it? It seems a reasonable approach to use even with negative numbers, provided we remember that any given cell potentially "needs" an arbitrary number of `MOD`s subtracted from it. — j_random_hacker, Sep 08 '12 at 07:44
The reason why I said it was mathematically wrong is that when the numbers overflow and are clipped by modulo, you lose information about the original numbers permanently. Recording the "number of mods" is effectively just saying we'll use more bits to store the number -- there's no reason to use mod at all in that case! OP indicated that he used them 'to accommodate large value of result' -- MOD does exactly the opposite, in that it absolutely destroys the magnitude of the result beyond the value of N. — nneonneo, Sep 08 '12 at 07:47
As an example, if the 'correct' output value is 23384026197294446691258957323460528314494920687616, and I tell you it is 422922539 mod 1000000007, you can't even tell me how many digits the original value contained -- that information is permanently lost. If you say "let's just store the number of mods", you have to store the value 23384026033606264456015106131354785395011, which is *still* too big to fit in a 32- (or even 64-bit) value. At this point, you either use floats (and sacrifice precision) or go bigint (and sacrifice memory). — nneonneo, Sep 08 '12 at 07:50
@nneonneo: Sorry about my manners... :-/ I was assuming it is known that each final answer will be < `MOD`, and the problem was large *intermediate* results. (This is a common problem when calculating e.g. numbers of combinations in the naive way.) I'm not suggesting storing "the number of MODs". But I see now that even if we could likewise bound the answers from below when negative inputs are present ("no final answer below `-MOD`") we still need to "count" 0 or 1 MODs per output number. — j_random_hacker, Sep 08 '12 at 08:44
Here is a clue for most of the above commenters: When you see “1000000007” and some form of “modulo” in a question, it is a contest or homework problem, and the modulo is part of the specification. If the questioner has not stated this, they are seeking an undeserved advantage. No direct answer should be given. Some teaching answers might be offered. And the question is answered numerous times before on Stack Overflow. — Eric Postpischil, Sep 08 '12 at 08:44
possible duplicate of [I want to generate the nth term of the sequence 1,3,8,22,60 ,164 in Order(1) or order of (nlogn)](http://stackoverflow.com/questions/11301992/i-want-to-generate-the-nth-term-of-the-sequence-1-3-8-22-60-164-in-order1-or) — Eric Postpischil, Sep 08 '12 at 08:50
Cute. Yes, I recall the 1000000007 or similarly large primes appearing in several ACM-style problems in the past. — nneonneo, Sep 08 '12 at 10:55
@nneonneo: this is a contest problem,i agree that information would be lost if we take mod, but consider the following case at n=32 RP.a[0][0]=6566290 RP.a[1][0]=-2136000 ,now A.a[0][0]=96 A.a[1][0]=44 so 96*6566260=630360960 44*(-2136000)=-9398400 (with en the range and less than mod) at n=33 RP.a[0][0]=10996580 RP.a[0][1]=-3598250 , now 10996580*44=1055671680(greater than mod so it gets reduced),but abs(RP.a[1][0]*44) is less than MOD , so finally i get a huge negative answer.What approach should be used to contain the result?trying other approaches too based on the above hints/comments. — k53sc, Sep 08 '12 at 12:41
@EricPostpischil: Nobody is "seeking an undeserved advantage"or anything mate.If you know any resource for ex. a 200 page book on the topic i would be glad to through it. — k53sc, Sep 08 '12 at 12:44
@EricPostpischil:its is not "possible duplicate of I want to generate the nth term of the sequence 1,3,8,22,60 ,164 in Order(1) or order of (nlogn) " the sequence generated here is 5,18,44,96,195 — k53sc, Sep 08 '12 at 13:41
@k53sc: You asked for help on a contest problem without disclosing it is a contest problem. Answers to contest problems should be your own work. This question is a duplicate because changing the constants in a contest problem does not change the underlying problem. Searching Stack Overflow for “1000000007” shows dozens of similar questions. — Eric Postpischil, Sep 08 '12 at 15:35

integer overflow in matrix exponentiation

0 Answers0