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;
    }
}

2010/09/25

Toom-Cook 2, Random BigIntegers, bit length BigInteger

In the previous, the first blog, Toom Cook 2, TC2, was to fast. A lot of 1's were mutiplied with a lot of 1's.
1111 1111 * 1111 1111 =  1111 1110 0000 0001
If the numbers are cut into pieces, limbs, a lot of multiplications by one and zero are done.
To get realistic results, random BigIntegers have to be used.
Other gadgets in the code below are the bitlength of a BigInteger, and the bitlength of a byte.


//---------RaNDom----------bitLength--------------TC2---------------------
using System;
using System.Numerics;
using System.Diagnostics;

class TC2_test
{
    private static Stopwatch sw = new Stopwatch();

    static void Main()
    {
        Console.WindowHeight = Console.LargestWindowHeight;
        test_RND();
        test_MTP();

        Console.ReadLine();
    }

    private static void test_RND()
    {
        BigInteger U, V;
        for (int i = 1; i < 200; i += 20)
        {
            U = RND(i); V = RND(i);
            Console.WriteLine(i + "  " + bitLength(U) + "  " + U);
            Console.WriteLine(i + "  " + bitLength(V) + "  " + V);
        }
    }

    private static void test_MTP()
    {
        BigInteger U, V;
        for (int i = 0; i < 10; i++)
        {
            U = RND(100000);
            V = RND(100000);
            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();
        }
    }

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

    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 int seed;
    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 = BigInteger.Abs(U).ToByteArray();
        int i = bytes.Length - 1;
        return (i << 3) + bitLength(bytes[i]);
    }

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

}
//---------RaNDom----------bitLength--------------TC2---------------------

2010/09/21

BigInteger Long Multiplication and Division

It's 2010, I used VB.Net 2005 Express and the Java BigInteger Library (vjslib.dll),
and wrote the classical multiplication, Knuth section 4.3.1, and division algorithms.
It's still 2010, downloaded C# 2010 Express.
Wrote them again, and found out it's useless, C# has them.
Next step: Karatsuba Multiplication (Toom Cook 2).
Still have to test it, but the blog is started.
TC2 is about 30 times faster in the example code below,
squaring seems to be 10% faster.
2E85TTGXQJHZ

//---------------------------------------------------------------
using System;
using System.Numerics;
using System.Diagnostics;
namespace TC2Blog
{
    class Program
    {
        static void Main()
        {
            Stopwatch sw = new Stopwatch();
            BigInteger U = (BigInteger.One << 1000000) - 1;
            BigInteger V = -U;
            sw.Start();
            BigInteger P1 = MTP(U, V);
            sw.Stop();
            Console.WriteLine(sw.ElapsedMilliseconds);
            sw.Restart();
            BigInteger P2 = U * V;
            sw.Stop();
            Console.WriteLine(sw.ElapsedMilliseconds);
            sw.Restart();
            BigInteger P3 = BigInteger.Pow(U, 2);
            sw.Stop();
            Console.WriteLine(sw.ElapsedMilliseconds);
            if (P1.CompareTo(P2) != 0)
                Console.WriteLine("MTP WRONG");
            if (BigInteger.Abs(P1).CompareTo(P3) != 0)
                Console.WriteLine("SQUARE WRONG");
            Console.WriteLine("Ready");
            Console.ReadLine();
        }

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

        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;
        }
    }
}
//---------------------------------------------------------------