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