Pagina's

2017/09/25

Direct Search Factorization


// MathWorld: DSF is the simplest (and most simple-minded) prime factorization algorithm.
// It consists of searching for factors of a number by systemetically performing trial divisions,
// usually using a sequence of increasing numbers. Multiples of small primes are commonly excluded
// to reduce the number of trial divisors, but just including them is sometimes faster than the time
// required to exclude them. DSF is very inefficient, and can be used only with fairly small numbers.

using System;                             // i7-4790@3.6GHz
using System.Collections.Generic;
class Program
{
    static void Main()
    {
        findPF(~0uL);                     // warm up
        time(4294967291);                 // largest prime < 2^32          50 µs
        time(4294967291 * 4294967291uL);  //                       6.5 s  
        time(3183958073 * 5793651691uL);  // 2 primes              4.8 s
        time(~0uL - 58);                  // largest prime < 2^64  6.5 s
        time(~0uL);                       // (2^64)-1                     100 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n == 1) return pf;
        ulong d = 1; uint b = 0, rt = sqrt(n);
        if (n <= ~0u) { findPF((uint)n, rt, 1, 0, pf); return pf; }
        for (; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > rt) { pf.Add(n); return pf; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return pf;
                    rt = sqrt(n);
                    if (n <= ~0u) { findPF((uint)n, rt, (uint)d, b, pf); return pf; }
                }
            }
            while (d + 30 <= rt && n % (d + 06) > 0 && n % (d + 10) > 0
                                && n % (d + 12) > 0 && n % (d + 16) > 0
                                && n % (d + 18) > 0 && n % (d + 22) > 0
                                && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static void findPF(uint n, uint rt, uint d, uint b, List<ulong> pf)
    {
        for (; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > rt) { pf.Add(n); return; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return;
                    rt = (uint)Math.Sqrt(n);
                }
            }
            while (d + 30 <= rt && n % (d + 06) > 0 && n % (d + 10) > 0
                                && n % (d + 12) > 0 && n % (d + 16) > 0
                                && n % (d + 18) > 0 && n % (d + 22) > 0
                                && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        ulong n1 = n0 - 1 >> 24; int s = 25;
        if (n1 > 255) { n1 >>= 4; s = 29; }
        if (n1 > 15) { n1 >>= 2; s |= 2; }
        if (n1 > 3) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}



///////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                               //
//                                  Eight threads, four cores.                                   //          
//                                                                                               //
///////////////////////////////////////////////////////////////////////////////////////////////////

using System;
using System.Collections.Generic;
using System.Threading.Tasks;
class Program
{
    static void Main()
    {
        findPF(~0uL);                     // warm up
        time(4294967291);                 // largest prime < 2^32          90 µs
        time(4294967291 * 4294967291uL);  //                              110 µs 
        time(3183958073 * 5793651691uL);  // 2 primes              0.9 s
        time(~0uL - 58);                  // largest prime < 2^64  1.1 s
        time(~0uL);                       // (2^64)-1                      80 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n > 1) { recFindF(n, pf); pf.Sort(); } return pf;
    }

    static void recFindF(ulong n, List<ulong> pf)
    {
        ulong f = findF(n); if (f == n) pf.Add(n);
        else { recFindF(f, pf); recFindF(n / f, pf); }
    }

    static ulong findF(ulong n)
    {
        uint r = sqrt(n); if (n % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            ulong d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                       k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31;
            ulong d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90;
            while (d0 <= r)
            {
                if (n % d0 == 0 || n % d1 == 0 || n % d2 == 0 || n % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n = n % d0 == 0 ? d0 : n % d1 == 0 ? d1 : n % d2 == 0 ? d2 : d3;
                        return;
                    }
                d0 += 120; d1 += 120; d2 += 120; d3 += 120;
            }
        }); return n;
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        ulong n1 = n0 - 1 >> 24; int s = 25;
        if (n1 > 255) { n1 >>= 4; s = 29; }
        if (n1 > 15) { n1 >>= 2; s |= 2; }
        if (n1 > 3) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}


///////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                               //
//                                               ?                                               //
//                                                                                               //
///////////////////////////////////////////////////////////////////////////////////////////////////

using System;
using System.Collections.Generic;
using System.Threading.Tasks;
class Program
{
    static void Main()
    {
        sqrt(1uL << 52);
        ulong[] n = new ulong[] { ~0u - 4, 0, ~0uL - 172, ~0uL - 58, ~0uL };
        n[1] = n[0] * n[0];
        findPF(n[0]); time(n[0]);  // largest prime < 2^32          50 µs
        findPF(n[1]); time(n[1]);  //                               80 µs 
        findPF(n[2]); time(n[2]);  // 2 primes              0.78 s
        findPF(n[3]); time(n[3]);  // largest prime < 2^64  1.05 s
        findPF(n[4]); time(n[4]);  // (2^64)-1                      70 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n > 1)
            if (n <= ~0u) findPF((uint)n, pf);
            else { recFindF(n, pf, 1); pf.Sort(); }
        return pf;
    }

    static void findPF(uint n, List<ulong> pf)
    {
        for (uint d = 1, r = (uint)Math.Sqrt(n), b = 0; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > r) { pf.Add(n); return; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return;
                    r = (uint)Math.Sqrt(n);
                }
            }
            while (d + 30 <= r && n % (d + 06) > 0 && n % (d + 10) > 0
                               && n % (d + 12) > 0 && n % (d + 16) > 0
                               && n % (d + 18) > 0 && n % (d + 22) > 0
                               && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static void recFindF(ulong n, List<ulong> pf, int i)
    {
        ulong f = findF(n);
        if (f == n) { while (i > 0) { pf.Add(n); i--; } }
        else if (n / f == f) { recFindF(f, pf, i << 1); }
        else { recFindF(f, pf, i); recFindF(n / f, pf, i); }
    }

    static ulong findF(ulong n0)
    {
        if (n0 <= ~0u) return findF((uint)n0);
        uint r = sqrt(n0); if (n0 % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            ulong d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                       k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31,
                  d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90, n1 = n0;
            for (; d0 <= r; d0 += 120, d1 += 120, d2 += 120, d3 += 120)
                if (n1 % d0 == 0 || n1 % d1 == 0 || n1 % d2 == 0 || n1 % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n0 = n0 % d0 == 0 ? d0 : n0 % d1 == 0 ? d1 : n0 % d2 == 0 ? d2 : d3;
                        return;
                    }
        }); return n0;
    }

    static uint findF(uint n0)
    {
        uint r = (uint)Math.Sqrt(n0); if (n0 % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            uint d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                      k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31,
                 d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90, n1 = n0;
            for (; d0 <= r; d0 += 120, d1 += 120, d2 += 120, d3 += 120)
                if (n1 % d0 == 0 || n1 % d1 == 0 || n1 % d2 == 0 || n1 % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n0 = n0 % d0 == 0 ? d0 : n0 % d1 == 0 ? d1 : n0 % d2 == 0 ? d2 : d3;
                        return;
                    }
        }); return n0;
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        uint n1 = (uint)(n0 >> 48); int s = 25;
        if (n1 > 255u) { n1 >>= 8; s = 29; }
        if (n1 > 15u) { n1 >>= 4; s |= 2; }
        if (n1 > 3u) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}