Pagina's

2010/09/26

Toom Cook 3

using System;
using System.Numerics;
using System.Diagnostics;

class Program
{
    private static Stopwatch sw = new Stopwatch();
    private static ushort[] qr3;
    private static int seed;

    static void Main()
    {
        create_qr3();
        test_MTP();
        Console.ReadLine();
    }

    private static void test_MTP()
    {
        BigInteger U, V;
        for (int i = 0; i < 10; i++)
        {
            seed += i;
            U = RND(200000);
            V = -RND(200000);
            sw.Start();
            BigInteger P1 = MTP(U, V);
            sw.Stop();
            Console.WriteLine("MTP " + sw.ElapsedMilliseconds);
            sw.Restart();
            BigInteger P2 = U * V;
            sw.Stop();
            Console.WriteLine("U*V " + sw.ElapsedMilliseconds);
            if (P1 != P2) Console.WriteLine("MTP WRONG");
            sw.Reset();
            Console.WriteLine(BitLength(P1));
        }
    }

    private static BigInteger D3(BigInteger U)
    {
        ushort t = 0;
        byte[] bytes = (U.Sign * U).ToByteArray();
        for (int i = bytes.Length - 1; i >= 0; i--)
        {
            t <<= 8;
            t = qr3[t | bytes[i]];
            bytes[i] = (byte)(t >> 8);
        }
        return U.Sign * new BigInteger(bytes);
    }

    private static void create_qr3()
    {
        int i, j = 0, k;
        qr3 = new ushort[byte.MaxValue + (1 << 9) + 1];
        for (i = 0; i <= byte.MaxValue + (1 << 9); i += 3)
        {
            for (k = 0; k < 3; k++) qr3[i + k] = (ushort)(j | k);
            j += 1 << 8;
        }
    }

    private static BigInteger MTP(BigInteger U, BigInteger V)
    {
        int n = BigInteger.Max(U, V).ToByteArray().Length << 3;
        if (n > 10000) return TC3(U, V, n);
        if (n > 3200) return TC2(U, V, n);
        return U * V;
    }

    private static BigInteger TC3(BigInteger U2, BigInteger V2, int n)
    {
        n = n / 3;// (int)((long)(n) * 0x55555556 >> 32);

        BigInteger Mask = (BigInteger.One << n) - 1;
        BigInteger U0 = U2 & Mask; U2 >>= n;
        BigInteger U1 = U2 & Mask; U2 >>= n;
        BigInteger V0 = V2 & Mask; V2 >>= n;
        BigInteger V1 = V2 & Mask; V2 >>= n;
        BigInteger W0 = MTP(U0, V0);
        BigInteger W4 = MTP(U2, V2);
        BigInteger P3 = MTP((U2 << 2) + (U1 << 1) + U0, (V2 << 2) + (V1 << 1) + V0);
        U2 += U0;
        V2 += V0;
        BigInteger P2 = MTP(U2 - U1, V2 - V1);
        BigInteger P1 = MTP(U2 + U1, V2 + V1);
        BigInteger W2 = ((P1 + P2) >> 1) - (W0 + W4);
        BigInteger W3 = W0 - P1;
        W3 = D3((W3 + P3 - P2 >> 1) + W3) - (W4 << 1);
        BigInteger W1 = P1 - (W4 + W3 + W2 + W0);
        return (((((((W4 << n) + W3) << n) + W2) << n) + W1) << n) + W0;
    }

    private static BigInteger TC2(BigInteger U1, BigInteger V1, int n)
    {
        n = n >> 1;
        BigInteger Mask = (BigInteger.One << n) - 1;
        BigInteger U0 = U1 & Mask; U1 >>= n;
        BigInteger V0 = V1 & Mask; V1 >>= n;
        BigInteger P0 = MTP(U0, V0);
        BigInteger P2 = MTP(U1, V1);
        BigInteger P1 = MTP(U0 + U1, V0 + V1) - (P0 + P2);
        return (((P2 << n) + P1) << n) + P0;
    }

    private static BigInteger RND(int n)
    {
        if (n == 0) return 0;
        if (seed == int.MaxValue) seed = 0; else seed++;
        Random rand = new Random(seed);
        byte[] bytes = new byte[((n - 1) >> 3) + 2];
        rand.NextBytes(bytes);
        bytes[bytes.Length - 1] = 0;
        n = ((bytes.Length - 1) << 3) - n;
        bytes[bytes.Length - 2] >>= n;
        bytes[bytes.Length - 2] |= (byte)(128 >> n);
        return new BigInteger(bytes);
    }

    private static int BitLength(BigInteger U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return (i << 3) + bitLength(bytes[i]);
    }

    private static int bitLength(byte b)
    {
        return b < 8 ? b < 2 ? b < 1 ? 0 : 1 : b < 4 ? 2 : 3 : b < 32 ? b < 16 ? 4 : 5 : b < 64 ? 6 : 7;
    }
}

No comments:

Post a Comment