Pasted

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