最下位bit列から状態を復元するやつ(2)

XorShift128p_Inverse.cs

using System;
using System.Numerics;
using System.Linq;

using static System.Console;

namespace ConsoleAppCore
{
    class Program
    {
        static void Main(string[] args)
        {
            var matrix = new (uint s0, uint s1, uint s2, uint s3)[128];

            // 遷移行列Tに対し、各T^iの最下位bit生成に関係する行を取り出す。
            // 各基底に遷移関数を作用させたときの最下位bitを並べればよい。
            for (int k = 0; k < 32; k++)
            {
                var state = (1u << k, 0U, 0u, 0u);
                for (int i = 0; i < 128; i++) matrix[i].s0 |= ((state.GetRand() & 1) << k);
            }
            for (int k = 0; k < 32; k++)
            {
                var state = (0u, 1u << k, 0u, 0u);
                for (int i = 0; i < 128; i++) matrix[ i].s1 |= ((state.GetRand() & 1) << k);
            }
            for (int k = 0; k < 32; k++)
            {
                var state = (0u, 0U, 1u << k, 0u);
                for (int i = 0; i < 128; i++) matrix[i].s2 |= ((state.GetRand() & 1) << k);
            }
            for (int k = 0; k < 32; k++)
            {
                var state = (0u, 0U, 0u, 1u << k);
                for (int i = 0; i < 128; i++) matrix[i].s3 |= ((state.GetRand() & 1) << k);
            }

            // 逆行列を計算する。
            var inv = matrix.GetInv();

            // テスト
            {
                var state = (0xBEEFu, 0xDEADu, 0xFACEu, 0xBADu);

                // 右逆行列になっているか？
                {
                    var test = state.Products(matrix).Products(inv);
                    // stateが復元される。
                    WriteLine(test.ToSeed());
                    WriteLine();
                }

                // 左逆行列になっているか？
                {
                    var test = state.Products(inv).Products(matrix);
                    // stateが復元される。
                    WriteLine(test.ToSeed());
                    WriteLine();
                }

                // 実際に最下位bit列を生成する。
                {
                    var copy = state;

                    var result = (0u, 0u, 0u, 0u);
                    for (int i = 0; i < 32; i++)
                    {
                        result.Item1 |= (copy.GetRand() & 1) << i;
                    }
                    for (int i = 0; i < 32; i++)
                    {
                        result.Item2 |= (copy.GetRand() & 1) << i;
                    }
                    for (int i = 0; i < 32; i++)
                    {
                        result.Item3 |= (copy.GetRand() & 1) << i;
                    }
                    for (int i = 0; i < 32; i++)
                    {
                        result.Item4 |= (copy.GetRand() & 1) << i;
                    }
                    
                    // matrixが最下位bit列への写像になっていることの確認
                    WriteLine(result.ToSeed());
                    WriteLine(state.Products(matrix).ToSeed());
                    WriteLine();

                    // stateが復元される。
                    WriteLine(result.Products(inv).ToSeed());

                    WriteLine();
                }
            }

            //foreach (var row in inv) WriteLine(row.ToSeed());
        }
    }

    public static class XorShift128p
    {
        public static uint GetRand(this ref (uint s1, uint s2, uint s3, uint s4) state)
        {
            var t1 = state.s1 ^ (state.s1 << 11);
            var t2 = state.s4 ^ (state.s4 >> 19) ^ t1 ^ (t1 >> 8);

            state = (state.s2, state.s3, state.s4, t2);

            return t2;
        }
        public static string ToSeed(this (uint s1, uint s2, uint s3, uint s4) state)
            => $"({state.s1:X8}, {state.s2:X8}, {state.s3:X8}, {state.s4:X8})";
    }
    public static class MatCalc
    {
        // 掃き出し法にて逆行列を求める。
        public static (uint s0, uint s1, uint s2, uint s3)[] GetInv(this (uint s0, uint s1, uint s2, uint s3)[] mat)
        {
            mat = mat.ToArray();

            // 初期状態として単位行列を作る。
            var inv = new (uint s0, uint s1, uint s2, uint s3)[128];
            for (int i = 0; i < 32; i++)
            {
                inv[i].s0 = (1U << i);
                inv[i + 32].s1 = (1U << i);
                inv[i + 64].s2 = (1U << i);
                inv[i + 96].s3 = (1U << i);
            }

            for (int i = 0; i < 32; i++)
            {
                // i行目が1である列を探す。
                var r = mat.GetPoppedRow(i);

                // 見つけた列をi列目と交換。
                if (r != i)
                {
                    mat.SwapRows(i, r);
                    inv.SwapRows(i, r);
                }

                // i列目以外のi行目が全て0になるように引いて(足して)いく。
                // (2元体なので加算=減算)。
                for (int j = 0; j < 128; j++)
                {
                    if (j == i) continue;

                    if ((mat[j].s0 & (1UL << i)) != 0)
                    {
                        mat[j].s0 ^= mat[i].s0;
                        mat[j].s1 ^= mat[i].s1;
                        mat[j].s2 ^= mat[i].s2;
                        mat[j].s3 ^= mat[i].s3;

                        inv[j].s0 ^= inv[i].s0;
                        inv[j].s1 ^= inv[i].s1;
                        inv[j].s2 ^= inv[i].s2;
                        inv[j].s3 ^= inv[i].s3;
                    }
                }
            }

            for (int i = 0; i < 32; i++)
            {
                var k = i + 32;
                var r = mat.GetPoppedRow(k);
                if (r != k)
                {
                    mat.SwapRows(k, r);
                    inv.SwapRows(k, r);
                }

                for (int j = 0; j < 128; j++)
                {
                    if (j == k) continue;

                    if ((mat[j].s1 & (1U << i)) != 0)
                    {
                        mat[j].s0 ^= mat[k].s0;
                        mat[j].s1 ^= mat[k].s1;
                        mat[j].s2 ^= mat[k].s2;
                        mat[j].s3 ^= mat[k].s3;

                        inv[j].s0 ^= inv[k].s0;
                        inv[j].s1 ^= inv[k].s1;
                        inv[j].s2 ^= inv[k].s2;
                        inv[j].s3 ^= inv[k].s3;
                    }
                }
            }

            for (int i = 0; i < 32; i++)
            {
                var k = i + 64;
                var r = mat.GetPoppedRow(k);
                if (r != k)
                {
                    mat.SwapRows(k, r);
                    inv.SwapRows(k, r);
                }

                for (int j = 0; j < 128; j++)
                {
                    if (j == k) continue;

                    if ((mat[j].s2 & (1U << i)) != 0)
                    {
                        mat[j].s0 ^= mat[k].s0;
                        mat[j].s1 ^= mat[k].s1;
                        mat[j].s2 ^= mat[k].s2;
                        mat[j].s3 ^= mat[k].s3;

                        inv[j].s0 ^= inv[k].s0;
                        inv[j].s1 ^= inv[k].s1;
                        inv[j].s2 ^= inv[k].s2;
                        inv[j].s3 ^= inv[k].s3;
                    }
                }
            }

            for (int i = 0; i < 32; i++)
            {
                var k = i + 96;
                var r = mat.GetPoppedRow(k);
                if (r != k)
                {
                    mat.SwapRows(k, r);
                    inv.SwapRows(k, r);
                }

                for (int j = 0; j < 128; j++)
                {
                    if (j == k) continue;

                    if ((mat[j].s3 & (1U << i)) != 0)
                    {
                        mat[j].s0 ^= mat[k].s0;
                        mat[j].s1 ^= mat[k].s1;
                        mat[j].s2 ^= mat[k].s2;
                        mat[j].s3 ^= mat[k].s3;

                        inv[j].s0 ^= inv[k].s0;
                        inv[j].s1 ^= inv[k].s1;
                        inv[j].s2 ^= inv[k].s2;
                        inv[j].s3 ^= inv[k].s3;
                    }
                }
            }
            
            return inv;
        }

        private static int GetPoppedRow(this (uint s0, uint s1, uint s2, uint s3)[] mat, int i)
        {
            for(int k=i; k<128; k++)
            {
                if (i < 32)
                { if ((mat[k].s0 & (1U << i)) != 0) return k; }
                else if(i < 64)
                { if ((mat[k].s1 & (1U << (i - 32))) != 0) return k; }
                else if (i < 96)
                { if ((mat[k].s2 & (1U << (i - 64))) != 0) return k; }
                else
                { if ((mat[k].s3 & (1U << (i - 96))) != 0) return k; }
            }

            // ここには来ないはず。
            throw new Exception();
        }

        private static void SwapRows(this (uint s0, uint s1, uint s2, uint s3)[] mat, int r1, int r2)
        {
            var temp = mat[r1];
            mat[r1] = mat[r2];
            mat[r2] = temp;
        }

        // ベクトルに行列を作用させる。
        public static (uint s0, uint s1, uint s2, uint s3) Products(this (uint s0, uint s1, uint s2, uint s3) state, (uint s0, uint s1, uint s2, uint s3)[] matrix)
        {
            var result = (0U, 0U, 0u, 0u);
            for (int i = 0; i < 32; i++)
            {
                var (s0, s1, s2, s3) = matrix[i];
                // 論理積を取り、立っているビットの偶奇を見ることで畳み込み計算ができる。
                // convolution; 畳み込み
                var conv = (uint)((s0 & state.s0).PopCount() 
                    + (s1 & state.s1).PopCount()
                    + (s2 & state.s2).PopCount()
                + (s3 & state.s3).PopCount()) & 1;

                result.Item1 |= (conv << i);
            }
            for (int i = 0; i < 32; i++)
            {
                var (s0, s1, s2, s3) = matrix[i + 32];
                var conv = (uint)((s0 & state.s0).PopCount()
                    + (s1 & state.s1).PopCount()
                    + (s2 & state.s2).PopCount()
                + (s3 & state.s3).PopCount()) & 1;

                result.Item2 |= (conv << i);
            }
            for (int i = 0; i < 32; i++)
            {
                var (s0, s1, s2, s3) = matrix[i + 64];
                var conv = (uint)((s0 & state.s0).PopCount()
                    + (s1 & state.s1).PopCount()
                    + (s2 & state.s2).PopCount()
                + (s3 & state.s3).PopCount()) & 1;

                result.Item3 |= (conv << i);
            }
            for (int i = 0; i < 32; i++)
            {
                var (s0, s1, s2, s3) = matrix[i + 96];
                var conv = (uint)((s0 & state.s0).PopCount()
                    + (s1 & state.s1).PopCount()
                    + (s2 & state.s2).PopCount()
                + (s3 & state.s3).PopCount()) & 1;

                result.Item4 |= (conv << i);
            }

            return result;
        }
    }

    public static class MyBitOps
    {
        public static int PopCount(this ulong n)
            => BitOperations.PopCount(n);
        public static int PopCount(this uint n)
            => BitOperations.PopCount(n);
    }
}

元の状態に復元する逆行列(s0,s1,s2,s3の順に並んでいる)

(FFD3C800, 3F985D65, 10046D8B, 10000000)
(5F51051F, 06CABE94, AACC6042, 82007048)
(FEF760C5, C79606F2, 39E96831, ABC0A1B7)
(67EFDA35, 2058504E, E8053930, C1D46F6B)
(C09D20C3, 8782C17D, 01CB0EF3, 565D52F6)
(0EBA7E48, 71B58959, 2A7CEE82, 0FFA8B16)
(2D3DDA6E, 610B409F, 48BBA9BB, 125C8CE1)
(BBC7C1ED, 2E0C9610, FC86D8C7, DDE78D3D)
(3AFE1972, 6DA74F7D, E3CCBDDE, 8B3314CA)
(F45C13B9, A1EB726C, 81C66763, BDDA8C0A)
(BEEBB8CC, B7D11308, A31605AD, 03ADA194)
(D93E7D15, C89FA989, 7D96B677, 4278D782)
(4C5F7375, A8DA2F0C, EF27BDF9, 6BF414C5)
(1E573796, 023AC64F, E4506698, 0467B5A6)
(63A20074, E9266084, 2DDE4D2B, BB0CE524)
(6AFE96B1, DC1A087D, DA982FFE, A64A1052)
(8AC49C3E, 5003FF9E, 208251A7, D29FB3D2)
(9FA6B5E7, 547608A9, A34AF89A, 38C41BAE)
(587A6F81, B555CE2D, 28AD4D9F, 0A3F8E98)
(E79EC22E, 148E707A, 2E641ABB, 63666299)
(087F4DF4, C49DB2C4, 7A41E2E4, 41FBCF50)
(49812A4F, 4007853C, 192A5644, C94C42DC)
(1D2CA3BA, 0386D219, 08FB6FBA, 500926A0)
(04729B90, F285B351, 4AD1D26B, CD12D1D5)
(8EE577FB, 650C7D8B, E064E2B9, 64042EDA)
(2C1E9EB7, CB804BD8, E85AE01A, DBF515F6)
(CD36BCC5, B8B1C3FF, C10D0FA0, 83968EB4)
(1B70F4F7, DA2CCBC0, 791632D7, 179FF8DF)
(A9005B4C, 492A87A1, 2F4A97DF, D487D776)
(68960B80, 7DE23C7F, CBDDBCE8, 1E7856FB)
(EB94BAF6, 1755992C, 92B44E1D, F1CCCCBC)
(3772C09B, 1D919F41, 269F4C6E, 5EA699B4)
(FFA79000, 7F30BACB, 2008DB16, 20000000)
(439E8A3F, F410AB77, 55DE1837, 0400E090)
(00D2418B, 76A9DBBA, 739408D0, 5781436F)
(32E3346B, B93576C3, D04CAAD3, 83A8DED6)
(813A4186, 0F0582FB, 03961DE7, ACBAA5EC)
(1D74FC90, E36B12B2, 54F9DD04, 1FF5162C)
(5A7BB4DC, C216813E, 91775376, 24B919C2)
(8AB303DB, A59CFA7E, F94B693D, BBCF1A7A)
(88C0B2E5, 22CB48A5, C7DFA30F, 16662994)
(1584A773, BA533286, 03CA1674, 7BB51814)
(7DD77198, 6FA22611, 462C0B5B, 075B4329)
(B27CFA2A, 913F5313, FB2D6CEF, 84F1AF04)
(98BEE6EA, 51B45E18, DE4F7BF3, D7E8298B)
(3CAE6F2C, 04758C9E, C8A0CD30, 08CF6B4D)
(3A7880E9, 2BC91757, 5BFA42E4, 7619CA49)
(28C1AD63, 41B1C6A5, B576874E, 4C9420A4)
(E8B5B87D, 59822962, 41427BFD, A53F67A5)
(3F4D6BCE, A8EC1153, 4695F134, 7188375D)
(B0F4DF02, 6AAB9C5A, 515A9B3F, 147F1D30)
(CF3D845C, 291CE0F5, 5CC83576, C6CCC532)
(10FE9BE8, 893B6588, F483C5C9, 83F79EA0)
(6E3ED49F, 798ADC27, 3212743B, 929885B9)
(3A594774, 070DA432, 11F6DF74, A0124D40)
(F5D9B721, 1C8EB0FD, 95E57C64, 9A25A3AB)
(1DCAEFF6, CA18FB17, C0C9C572, C8085DB5)
(A501BD6F, 6E8541EF, D0F31886, B7EA2BEC)
(6751F98B, 88E651A0, 825CC7F2, 072D1D68)
(36E1E9EE, B4599780, F22C65AF, 2F3FF1BE)
(AF3C3699, 6BD0D91C, 5ED3F70D, A90FAEED)
(D12C1700, FBC478FE, 97BB79D0, 3CF0ADF7)
(2A15F5ED, D72EE406, 252E4489, E3999978)
(6EE58136, 3B233E82, 4D3E98DC, BD4D3368)
(FF4F2000, FE617597, 4011B62C, 40000000)
(873D147E, E82156EE, ABBC306F, 0801C120)
(01A48316, ED53B774, E72811A0, AF0286DE)
(98FAE8D7, 8BEF3BD9, A0DF8D14, 0751BDAC)
(FF48030D, E78ED3A8, 076AE37D, 59754BD9)
(3AE9F920, C6D62564, A9F3BA09, 3FEA2C58)
(B4F769B8, 842D027C, 22EEA6ED, 49723385)
(E85A87B7, B2BC22A2, F2D00AC8, 779E34F4)
(118165CA, 4596914B, 8FBF461E, 2CCC5329)
(2B094EE6, 74A6650C, 07942CE9, F76A3028)
(FBAEE330, DF444C22, 8C5816B6, 0EB68652)
(99C57455, DBFB7078, F61C016C, 09E35E08)
(CC414DD5, 5AED6A6E, BCD82F55, AFD05316)
(795CDE58, 08EB193C, 91419A60, 119ED69B)
(74F101D2, 57922EAE, B7F485C8, EC339492)
(51835AC6, 83638D4A, 6AED0E9C, 99284149)
(2C57F0FB, 4A81849A, 82C22F49, 4A7ECF4B)
(7E9AD79C, 51D822A6, 8D2BE269, E3106EBA)
(61E9BE04, D55738B5, A2B5367E, 28FE3A60)
(634788B9, ABBC17B4, B9D6B25F, 8D998A65)
(DCC1B7D1, EBF31D4F, E9415320, 07EF3D40)
(2141293F, 0A906E11, 646230C5, 25310B73)
(898E0EE9, F79E9E3B, 23AB665B, 40249A81)
(168FEE43, C098B7A4, 2B8C207B, 344B4756)
(C6A95FED, 6DB42071, 81D55256, 9010BB6A)
(B73FFADF, 248F5580, A1A0E9BF, 6FD457D8)
(CEA3F316, 11CCA340, 04B98FE5, 0E5A3AD1)
(6DC3D3DC, 68B32F00, E458CB5F, 5E7FE37D)
(A344ED33, 2E246466, BDE136A9, 521F5DDB)
(A2582E00, F788F1FD, 2F76F3A1, 79E15BEF)
(A9176BDB, 57D81E53, 4A1A51A0, C73332F1)
(20F7826D, 8FC3AB5B, 9A3BE90B, 7A9A66D1)
(FE9E4000, FCC2EB2F, 80236C59, 80000000)
(0E7A28FC, D042ADDD, 577860DF, 10038241)
(FE75862D, 2322B8B7, CE16FBF2, 5E050DBC)
(31F5D1AE, 17DE77B3, 41BF1A29, 0EA37B59)
(FE90061A, CF1DA751, 0ED5C6FB, B2EA97B2)
(75D3F240, 8DAC4AC8, 53E77413, 7FD458B1)
(69EED370, 085A04F9, 45DD4DDB, 92E4670A)
(D0B50F6E, 65784545, E5A01591, EF3C69E9)
(2302CB94, 8B2D2296, 1F7E8C3C, 5998A653)
(AB2E1DCD, 10C91C47, 0F6E8161, EED46051)
(F75DC660, BE889845, 18B02D6D, 1D6D0CA5)
(338AE8AA, B7F6E0F1, EC3802D9, 13C6BC11)
(65BE1BAB, 4C5F0282, 79F68619, 5FA0A62C)
(F2B9BCB0, 11D63278, 228334C0, 233DAD37)
(14DE83A5, 56A18B03, 6FAFD323, D8672924)
(5E3A358D, FF42CCCB, D59CC58A, 32508293)
(58AFE1F6, 95030934, 05845E92, 94FD9E97)
(00092F39, 5A359313, 1A111C61, C620DD74)
(C3D37C08, AAAE716A, 456A6CFD, 51FC74C1)
(3BB39173, AEFDF937, 73EBBC0C, 1B3314CA)
(B9836FA2, D7E63A9F, D282A641, 0FDE7A81)
(4282527E, 1520DC22, C8C4618A, 4A6216E6)
(131C1DD2, EF3D3C77, 4756CCB7, 80493502)
(2D1FDC86, 81316F48, 571840F7, 68968EAC)
(706E3FDB, 22ED96BC, 03EC7C1F, 202176D4)
(6E7FF5BE, 491EAB01, 4341D37E, DFA8AFB1)
(9D47E62C, 23994681, 09731FCA, 1CB475A2)
(DB87A7B8, D1665E00, C8B196BE, BCFFC6FB)
(4689DA66, 5C48C8CD, 7BC26D52, A43EBBB7)
(44B05C00, EF11E3FB, 5EEDE743, F3C2B7DE)
(AF1257B7, 5635EAF8, 94727BF3, 8E6665E3)
(41EF04DA, 1F8756B6, 3477D217, F534CDA3)