Euler0549_simple.linq

// Euler0549_simple.linq
// 2016-06-01

#define FORCE_GC

#if !LINQPAD
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
#endif

//  # benching 'Euler0549_simple' ...
//
//  S(       10):                 39        0,0 ms
//  S(      100):               2012        0,0 ms
//  S(     1000):             136817        0,1 ms
//  S(    10000):           10125843        1,3 ms
//  S(   100000):          793183093       16,1 ms
//  S(  1000000):        64938007616      216,5 ms
//  S( 10000000):      5494373412573    3.538,4 ms
//  S(100000000):    #45a57e1327a02a   66.756,3 ms
//
//  # benching 'Euler0549_lpf_decomp' ...
//
//  S(       10):                 39        0,0 ms
//  S(      100):               2012        0,0 ms
//  S(     1000):             136817        0,0 ms
//  S(    10000):           10125843        0,3 ms
//  S(   100000):          793183093        1,9 ms
//  S(  1000000):        64938007616       19,5 ms
//  S( 10000000):      5494373412573      196,1 ms
//  S(100000000):    #45a57e1327a02a    1.990,3 ms

class Euler0549_simple_linq: Euler0549
{
    static void Main ()
    {
        bench_and_print(Euler0549_lpf_decomp.s_sum, "Euler0549_lpf_decomp");
        bench_and_print(Euler0549_simple.s_sum, "Euler0549_simple");
    }
}

class Euler0549
{
    public const int EULER_N = 100000000;

    // computes s(p^e) for p prime
    public static int s_for_prime_power (int p, int e)
    {
        int k = 0;

        while (e > p)
        {
            k += p;
            e -= p + 1;

            for (int t = k; (t /= p) % p == 0; )
                --e;
        }

        return (k + Math.Max(0 , e)) * p;
    }

    public delegate long test_func (int n);

    public static void bench_and_print (test_func f, string name, int lo = 10, int powers_of_ten = 8)
    {
        Console.WriteLine("# benching '{0}' ...\n", name);

        f(2*2*3*3);  // get it jitted

        for (int n = lo, e = 0; e < powers_of_ten; n *= 10, ++e)
        {
            Console.Write("S({0,9}): ", n);

#if FORCE_GC
            force_gc();
#endif
            var t = Stopwatch.StartNew();
            var x = f(n);
            t.Stop();

            Console.WriteLine("{0,18} {1,10:N1} ms",
                n == EULER_N ? spoiler(x) : x.ToString(),
                t.Elapsed.TotalMilliseconds  );
        }

        Console.WriteLine();
    }

    // last <decimal digit count> hex digits of the MD5 of the text representation
    public static string spoiler (long n)
    {
        var text = ASCIIEncoding.ASCII.GetBytes(n.ToString());
        var hash = new System.Security.Cryptography.MD5Cng()
            .ComputeHash(text)
            .Select(x => x.ToString("x2")).Aggregate((s,x) => s + x);
        return "#" + hash.Substring(Math.Max(0, hash.Length - text.Length + 1));
    }

    static void force_gc ()
    {
        System.Runtime.GCSettings.LargeObjectHeapCompactionMode = System.Runtime.GCLargeObjectHeapCompactionMode.CompactOnce;
        GC.Collect(9, GCCollectionMode.Forced);
    }
}

class Euler0549_lpf_decomp: Euler0549
{
    public static int[] s;

    public static long s_sum (int limit)
    {
        int small_prime_limit = (int)Math.Sqrt(limit);
        var small_primes = new List<int>();
        int m = 2;

        s = new int[limit + 1];

        for (int half = limit / 2; m <= half; ++m)
        {
            if (s[m] == 0)  // prime
            {
                s[m] = m;

                if (m <= small_prime_limit)
                    small_primes.Add(m);
            }

            int s_m = s[m], threshold = limit / m;

            foreach (int p in small_primes)
            {
                if (p > threshold)  // i.e. if (p * m > limit)
                    break;

                if (m % p != 0)  // not a factor of m yet
                {
                    // p is a new least factor -> s(p) < s(m) -> s(p * m) = s(m)
                    s[p * m] = s_m;
                }
                else // p is already a (least) factor of m
                {
                    int e = 2, q = m;

                    while ((q /= p) % p == 0)
                        ++e;

                    s[p * m] = Math.Max(s_for_prime_power(p, e), s[q]);
                    break;
                }
            }
        }

        for ( ; m <= limit; ++m)
            if (s[m] == 0)
                s[m] = m;

        long sum = 0;
        for (int i = 2; i <= limit; ++i)
            sum += s[i];

        return sum;
    }
}

class Euler0549_simple: Euler0549
{
    public static long s_sum (int n)
    {
        long result = 0;

        for (int i = 2; i <= n; ++i)
            result += s(factorisation(i));

        return result;
    }

    public static int s (List<int> factors, int lo = 0)
    {
        int p = factors[lo], hi = lo;

        while (++hi < factors.Count && factors[hi] == p)
            ;

        int result = s_for_prime_power(p, hi - lo);

        if (hi < factors.Count)
            result = Math.Max(result, s(factors, hi));

        return result;
    }

    static ushort[] g_small_primes = U16.SmallPrimesUpTo(1 << 16, 6542).ToArray();

    static List<int> factorisation (int n)
    {
        var result = new List<int>();
        int sqrt_n = (int)Math.Sqrt(n);

        foreach (var prime in g_small_primes)
        {
            if (prime > sqrt_n)
                break;

            if (n % prime == 0)
            {
                do result.Add(prime); while ((n /= prime) % prime == 0);

                sqrt_n = (int)Math.Sqrt(n);
            }
        }

        if (n != 1 || result.Count == 0)
            result.Add(n);

        return result;
    }
}


class U16
{
    public static List<ushort> SmallPrimesUpTo (int n, int initial_allocation_size = 0)
    {
        var result = new List<ushort>(initial_allocation_size);

        if (n < 2)
            return result;

        result.Add(2);  // needs to be pulled out of thin air because of the mod 2 wheel

        if (n < 3)
            return result;

        result.Add(3);  // needs to be pulled out of thin air because of the mod 3 decoder

        int sqrt_n_halved = (int)Math.Sqrt(n) >> 1;
        int max_bit = (n - 1) >> 1;
        var odd_composite = new bool[max_bit + 1];

        for (int i = 5 >> 1; i <= sqrt_n_halved; ++i)
            if (!odd_composite[i])
                for (int p = (i << 1) + 1, j = p * p >> 1; j <= max_bit; j += p)
                    odd_composite[j] = true;

        for (int i = 5 >> 1, d = 1; i <= max_bit; i += d, d ^= 3)
            if (!odd_composite[i])
                result.Add((ushort)((i << 1) + 1));

        return result;
    }
}