An Improved Method for Predicting Truncated Fibonacci LFSRs over Integer Residue Rings

// This code is used to recover the truncated ZUC's driving sequence.
// The library used is the NTL 11.4.3 library.
// The whole experimental process consists of 4 steps. Searching linear relations, Recovering the modulus, the coefficients and the initial states.

#include <NTL/vec_ZZ.h>
#include <NTL/mat_ZZ.h>
#include <NTL/mat_ZZ_p.h>
#include <NTL/ZZ_pX.h>
#include <NTL/ZZX.h>
#include <NTL/LLL.h>
#include <fstream>
#include <ctime>
#include <math.h>
#include <chrono>
#include <fstream>
NTL_CLIENT
using namespace std;

// Search polynomials that annihilate the sequence a
// The high-order scenario
// The modulus m is known
mat_ZZ Search_linear_relations_high_m(ZZ m, vec_ZZ y, int beta, int r, int t)
{
    auto start = std::chrono::high_resolution_clock::now();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = m;
    for (i = t; i < r + t; i++)
    {
        L[i][i] = power2_ZZ(beta);
        for (j = 0; j < t; j++)
        {
            L[i][j] = y[i + j - t] * power2_ZZ(beta);
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Choose 3 polynomials that annihilate the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t] / power2_ZZ(beta);
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    double time_in_seconds = duration.count();
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time_in_seconds << "s\n\n";

    return C;
}

// Search polynomials that annihilate the sequence a
// The high-order scenario
// The modulus m is unknown
mat_ZZ Search_linear_relations_high(vec_ZZ y, int alpha, int r, int t)
{
    auto start = std::chrono::high_resolution_clock::now();
    mat_ZZ L, C;
    int i, j;

    L.SetDims(r + t, r + t);
    C.SetDims(r, r);

    // Construct the lattice L
    for (i = 0; i < r + t; i++)
    {
        for (j = 0; j < r + t; j++)
        {
            L[i][j] = 0;
        }
    }
    for (i = 0; i < t; i++)
        L[i][i] = power2_ZZ(alpha);
    for (i = t; i < r + t; i++)
    {
        L[i][i] = 1;
        for (j = 0; j < t; j++)
        {
            L[i][j] = y[i + j - t];
        }
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Choose 3 polynomials that annihilate the sequence a
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            C[i][j] = L[i][j + t];
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    double time_in_seconds = duration.count();
    cout << "Step 1: time of searching polynomials that annihilate the sequence a: " << time_in_seconds << "s\n\n";

    return C;
}

// Recover the modulus by the resultant
ZZ Recover_the_modulus(mat_ZZ C, int n, int k, int r)
{
    auto start = std::chrono::high_resolution_clock::now();
    ZZ res, m_power, recovered_m;
    ZZX f, g;
    int i, j, flag;

    f.SetLength(r);
    g.SetLength(r);

    flag = 0;
    for (i = 0; i < r; i++)
        if (C[0][i] == 0)
            flag = flag + 1;
    if (flag == r)
        cout << "Fail to recover the modulus m!\n";

    flag = 2;
    m_power = 0;
    while (flag < r)
    {
        for (i = 0; i < r; i++)
            SetCoeff(f, i, C[flag - 1][i]);
        for (j = 0; j < flag - 1; j++)
        {
            for (i = 0; i < r; i++)
            {
                SetCoeff(g, i, C[j][i]);
            }
            res = resultant(f, g);
            m_power = GCD(m_power, res);
            //cout << m_power << "\n";
        }
        if (NumBits(m_power) == k * n)
            break;
        flag = flag + 1;
    }
    recovered_m = SqrRoot(SqrRoot(SqrRoot(SqrRoot(m_power))));

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    double time_in_seconds = duration.count();
    cout << "Step 2: time of recovering the modulus: " << time_in_seconds << "s\n\n";
    cout << "Number of the annihilating polynomials for recovering the modulus = " << flag << "\n";

    return recovered_m;
}

// Recover the coefficients by the Chinese Remainder Theorem
ZZ_pX Recover_the_coefficients(mat_ZZ C, ZZ m, int n, int r)
{
    auto start = std::chrono::high_resolution_clock::now();
    ZZ_pX f, g, h;
    mat_ZZ_p C_m;
    ZZ_p c;
    int i, j, mark, flag;
    c.init(m);

    C_m.SetDims(r, r);

    ofstream fout("C.txt");
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            fout << C[i][j] << " ";
        }
        fout << "\n";
    }
    fout.close();
    freopen("C.txt", "r", stdin);

    // Make the polynomials monic
    for (i = 0; i < r; i++)
    {
        for (j = 0; j < r; j++)
        {
            cin >> c;
            C_m[i][j] = c;
        }
        mark = r - 1;
        while (C_m[i][mark] == 0)
        {
            mark = mark - 1;
        }
        if (mark >= 0)
        {
            for (j = 0; j < r; j++)
            {
                C_m[i][j] = C_m[i][j] * inv(C_m[i][mark]);
            }
        }
    }

    f.SetLength(r);
    g.SetLength(r);

    for (i = 0; i < r; i++)
    {
        SetCoeff(f, i, C_m[0][i]);
        SetCoeff(g, i, C_m[1][i]);
    }

    h = GCD(f, g);
    flag = 2;

    while (deg(h) > n)
    {
        for (i = 0; i < r; i++)
            SetCoeff(f, i, C_m[flag][i]);
        h = GCD(f, h);
        flag = flag + 1;
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    double time_in_seconds = duration.count();
    cout << "Step 3: time of recovering the coefficients: " << time_in_seconds << "s\n\n";
    cout << "Number of the annihilating polynomials for recovering the coefficients = " << flag << "\n";

    return h;
}

// Recover the initial state by Kannan's embedding technique
// The high-order scenario
vec_ZZ Recover_the_initial_state(vec_ZZ y, ZZ_pX f, ZZ m, int n, int d, int beta)
{
    auto start = std::chrono::high_resolution_clock::now();
    vec_ZZ a, z;
    mat_ZZ Q, Q_power, L;
    ZZ b;
    int i, j, l;

    a.SetLength(n);
    z.SetLength(n);
    Q.SetDims(n, n);
    Q_power.SetDims(n, n);
    L.SetDims(d + 1, d + 1);

    Q_power = ident_mat_ZZ(n);

    Q[0][n - 1] = (-rep(f[0])) % m;
    for (i = 1; i < n; i++)
    {
        Q[i][i - 1] = 1;
        Q[i][n - 1] = (-rep(f[i])) % m;
    }

    // Construct the lattice L
    for (i = 1; i < n; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (l = 0; l < n; l++)
            {
                Q_power[j][l] = Q_power[j][l] % m;
            }
        }
    }

    for (i = n; i < d; i++)
    {
        Q_power = Q_power * Q;
        for (j = 0; j < n; j++)
        {
            for (l = 0; l < n; l++)
            {
                Q_power[j][l] = Q_power[j][l] % m;
            }
        }
        b = 0;
        for (j = 0; j < n; j++)
        {
            L[j + 1][i + 1] = Q_power[j][0];
            b += Q_power[j][0] * y[j];
        }
        L[0][i + 1] = (pow(2, beta) * (y[i] - b)) % m + pow(2, beta - 1);
    }

    L[0][0] = pow(2, beta - 1);
    for (i = 1; i <= n; i++)
    {
        L[0][i] = pow(2, beta - 1);
        L[i][i] = 1;
    }
    for (i = n + 1; i <= d; i++)
    {
        L[i][i] = m;
    }

    // Reduce the lattice L
    BKZ_FP(L, 0.99, 20);

    // Recover the initial state
    if (L[0][0] == -pow(2, beta - 1))
    {
        for (i = 0; i < n; i++)
        {
            z[i] = L[0][i + 1] + pow(2, beta - 1);
            a[i] = y[i] * pow(2, beta) + z[i];
        }
    }
    if (L[0][0] == pow(2, beta - 1))
    {
        for (i = 0; i < n; i++)
        {
            z[i] = pow(2, beta - 1) - L[0][i + 1];
            a[i] = y[i] * pow(2, beta) + z[i];
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration = end - start;
    double time_in_seconds = duration.count();
    cout << "Step 4: time of recovering the initial state: " << time_in_seconds << "s\n\n";

    return a;
}

void Experiment_2_1(vec_ZZ y, ZZ m, int n, int k)
{
    vec_ZZ recovered_a;
    mat_ZZ C;
    ZZ_pX f;
    int r, t, d, alpha, beta;
    int i;

    r = 71;
    t = 71;
    alpha = 17;
    beta = 14;
    d = 30;

    C = Search_linear_relations_high_m(m, y, beta, r, t);
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state(y, f, m, n, d, beta);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}

void Experiment_3_1(vec_ZZ y, int n, int k)
{
    vec_ZZ recovered_a;
    mat_ZZ C;
    ZZ m;
    ZZ_pX f;
    int r, t, d, alpha, beta;
    int i;

    r = 72;
    t = 70;
    alpha = 17;
    beta = 14;
    d = 30;

    C = Search_linear_relations_high(y, alpha, r, t);
    m = Recover_the_modulus(C, n, k, r);
    cout << "m = " << m << "\n\n\n";
    f = Recover_the_coefficients(C, m, n, r);
    cout << "f = " << f << "\n\n\n";
    recovered_a = Recover_the_initial_state(y, f, m, n, d, beta);
    cout << "Number of truncated digits for recovering the initial state = " << d << "\n";
    cout << "a = " << recovered_a << "\n\n\n\n\n";
}

void main()
{
    vec_ZZ y;
    ZZ m;
    int n, k, i, N;

    n = 16;
    k = 31;
    m = 2147483647;
    N = 141;

    y.SetLength(N);
    long Y[141] = { 312,6195,20905,36760,33275,128773,5120,82629,3434,13450,115969,78595,46080,21032,126472,116368,39147,118169,53608,16848,30654,130814,92427,66994,64940,95674,83420,12319,51686,121628,103101,4578,55561,114165,99105,66029,11183,33637,126822,45172,71211,75549,31249,100570,57872,40625,53029,98594,35382,114241,86331,100286,74451,71452,113487,22824,119007,42576,9420,18312,62834,3099,103626,124788,28318,119993,120470,1417,66987,121877,100551,36538,107257,95791,90747,11856,78683,27967,33439,112227,42492,76803,73986,127640,127988,113849,62073,121003,67257,62990,28618,14066,37196,88819,91188,119086,47509,14040,46365,117363,96302,61501,69537,19261,5302,29632,112053,60563,73212,123325,7485,52451,92665,108139,31895,6630,28136,97771,1874,10433,123318,98757,45897,34866,118234,75543,95739,41629,82990,25310,28330,21183,102442,39164,13744,125171,78275,42927,27280,113230,115192 };
    for (i = 0; i < N; i++)
        y[i] = ZZ(Y[i]);


    cout << "***Experiment 2.1***\n";
    Experiment_2_1(y, m, n, k);

    cout << "***Experiment 3.1***\n";
    Experiment_3_1(y, n, k);
}