Problem #231 says:
The binomial coefficient \( ^{10}C_{3} = 120\).
\(120 = 2^{3} \ 3 \ 5 = 2\ 2\ 2\ 3\ 5\), and \(2 + 2 + 2 + 3 + 5 = 14\).
So the sum of the terms in the prime factorisation of \(^{10}C_{3}\) is 14.
Find the sum of the terms in the prime factorisation of \(^{20000000}C_{15000000}\).
There’s a simple formula to calculate this. I don’t remember where I found it I wrote this solution a while ago and never posted it. Uses a bunch of prime numbers and the sum of factors to get the answer without brute force.The other way I suppose one could do this is there is a definition for the binomial coefficient as a product of a bunch of stuff. If we could eliminate all rounding error that approach could lead to a working solution.
Solution provided in F# and uses a Sieve of Atkins that seems to be slightly faster than other prime number generators I’ve been using.
#light
open System
open System.Collections
open System.Collections.Generic
open System.Diagnostics
open System.Threading
open Microsoft.FSharp.Collections
open Microsoft.FSharp.Math
let Atkins limit =
let is_prime = new Dictionary<uint64,bool>()
let keys = [5UL..limit]
for i in keys do
is_prime.Add(i,false)
let sqrt_limit = uint64 (Math.Sqrt (float limit))
let xys = [1UL..sqrt_limit]
for x in xys do
for y in xys do
let mutable n = 4UL*(x*x) + (y*y)
if (n <= limit && (n % 12UL = 1UL || n % 12UL = 5UL)) then
is_prime.[n] <- not is_prime.[n]
n <- 3UL*(x*x) + (y*y)
if (n <= limit && (n % 12UL = 7UL)) then
is_prime.[n] <- not is_prime.[n]
n <- 3UL*(x*x) - (y*y)
if ((x>y) && n <= limit && (n % 12UL = 11UL)) then
is_prime.[n] <- not is_prime.[n]
for i in keys do
if is_prime.[i] then
let max = limit / (i*i)
for k in [1UL..max] do
is_prime.[k*(i*i)] <- false
[
yield 2UL;yield 3UL;
for i in keys do
if is_prime.[i] then
yield i
]
let primes = Atkins 20000000UL |> Seq.toArray
let SumOfFactors n p =
let nt = ref n
[
while (!nt > 0UL) do
nt := !nt / p
yield !nt
] |> Seq.sum
let Problem231 n k =
primes
|> PSeq.map (fun i -> ((SumOfFactors n i) - (SumOfFactors k i) - (SumOfFactors (n-k) i) )* i)
|> Seq.sum
Console.WriteLine "Primes Finished"
let sw = new Stopwatch()
sw.Start()
Console.WriteLine (Problem231 20000000UL 15000000UL)
sw.Stop()
Console.WriteLine(sw.Elapsed)
Console.ReadLine() |> ignore
Tags: .net4, f#, fsharp, project euler
Posted in Project Euler | Comments (0)
Problem #124 wants us to:
The radical of n, rad(n), is the product of distinct prime factors of n. For example, \(504 = 2^3 × 3^2 × 7\), so rad(504) = 2 × 3 × 7 = 42.
Let E(k) be the kth element in the sorted n column; for example, E(4) = 8 and E(6) = 9.
If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).
This problem is exceedingly easy once we calculate the radicals for the numbers 1-100000. Just create a simple class to remember the n, and rad values and populate those fields and the problem is 9/10ths of the way solved. From here we can use the power of Linq, which I’m liking more and more as time goes on, and do an OrderBy on the rad values. Now this seems like it would be insufficient to solve the problem because if the radical values are the same we need to sort by the n value. Well Linq has us covered; simply add a ThenBy after the OrderBy and you can sort by a second value. This means that the overwhelming majority of the code in the solution is prime number generation, finding distinct factors, and calculation of the radicals.
A second solution without linq would be to sort by the radical value and get the radical value of the 10000th item. Then find the index of the first number with that radical and store it. Then get a list of only the numbers with that radical and sort this new list by the n value. Using the previously stored index you can then calculate what item in the new list should be the answer. But this is alot more complicated than a simple Linq expression.
A third method would be to write a custom IComparer for the class that holds the radical + n values. But this is a lot more work than is needed and would make the solution much less readable in my opinion. Besides an up to date .Net implementation (Microsoft.Net or Mono) should have Linq available.
Solution below written in c#/mono and requires Linq and .Net 4. Runs in ~5 seconds.
using System;
using System.Collections.Generic;
using System.Linq;
namespace Problem124
{
class MainClass
{
class N
{
public ulong n = 0;
public ulong rad = 0;
public N (ulong num)
{
this.n = num;
this.rad = rad(num);
}
}
public static ulong max = 100000;
public static List<ulong> primes = GeneratePrimes(max);
public static void Main (string[] args)
{
List<N> NCases = new List<N>();
for (ulong n = 1; n<=max;n++)
{
NCases.Add(new N(n));
}
N answer = NCases.OrderBy(x => x.rad)
.ThenBy(x => x.n).ToList()[10000-1];
Console.WriteLine(answer.n);
}
public static ulong rad(ulong n)
{
List<ulong> factors = DistinctFactors(n);
ulong rad = 1;
foreach(ulong f in factors)
rad *= f;
return rad;
}
public static List<ulong> DistinctFactors (ulong n)
{
List<ulong> factors = new List<ulong>();
if(primes.Contains(n))
{
factors.Add(n);
return factors;
}
ulong sqrt = (ulong)Math.Sqrt((double)n);
foreach (ulong p in primes)
{
if (p > sqrt)
break;
if (n%p == 0)
{
factors.Add(p);
while (n%p == 0)
n = n/p;
}
}
if (n!= 1ul)
factors.Add(n);
return factors;
}
static List<ulong> GeneratePrimes(ulong num)
{
List<ulong> RetValue = new List<ulong>();
RetValue.Add((ulong)2);
RetValue.Add((ulong)3);
ulong Stepper = 5;
ulong Check = 1;
while (Stepper <= num)
{
foreach (ulong i in RetValue)
{
if (Stepper % i == 0)
{
Check = 0;
break;
}
if (Math.Sqrt(Stepper) <= i)
{
break;
}
}
if (Check == 1)
{
RetValue.Add(Stepper);
}
Check = 1;
Stepper++;
Stepper++;
}
return RetValue;
}
}
}
Tags: .net4, c#, csharp, mono, project euler
Posted in Project Euler | Comments (0)
Problem #119 is looking for:
The number 512 is interesting because it is equal to the sum of its digits raised to some power: 5 + 1 + 2 = 8, and 83 = 512. Another example of a number with this property is 614656 = 284.
We shall define an to be the nth term of this sequence and insist that a number must contain at least two digits to have a sum.
You are given that a2 = 512 and a10 = 614656.
Find a30.
This is another problem on Project Euler that is much easier if you turn it inside out. The problem wants you to calculate the digit sum \(dsum\) of a number \(n\) and do something like \(log_{dsum}(n) = exp \) then check to see that \(dsum^{exp} == n\) which is doable but really slow. The faster method is to loop through two sets of numbers \(dbase\) and \(dexp\) a base and an exponent and find the digit sum of the resulting number and compare it to the base or something like \(DigitSum(dbase^{dexp}) == dbase\). This second method is much faster.
Solution written in c#/mono and uses the big integer library System.Numerics so requires .Net 4 . Calculates a bunch of extra numbers because of the loose boundaries but still runs in less than 4 seconds.
using System;
using System.Numerics;
using System.Collections.Generic;
namespace Problem119
{
class MainClass
{
public static void Main (string[] args)
{
List<BigInteger> An = new List<BigInteger>();
int Max = 30;
for (BigInteger dbase = 2; dbase <= 500;dbase++)
{
for (int dexp = 2; dexp <= 50;dexp++)
{
BigInteger n = BigInteger.Pow(dbase,dexp);
if (n.ToString().Length == 1)
continue;
if (DigitalSum(n) == dbase)
{
if (!An.Contains(n))
{
An.Add(n);
}
}
}
}
An.Sort();
Console.WriteLine(An[Max-1]);
}
public static BigInteger DigitalSum (BigInteger i)
{
BigInteger digits = (BigInteger) i;
BigInteger digitsum = 0;
while (digits != BigInteger.Zero)
{
BigInteger temp = (BigInteger)(digits%10);
digitsum += temp;
digits = (digits -temp) / 10;
}
return digitsum;
}
}
}
Tags: c#, csharp, mono, project euler
Posted in Project Euler | Comments (0)
Problem #133 says:
A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k; for example, R(6) = 111111.
Let us consider repunits of the form \(R(10^{n})\).
Although R(10), R(100), or R(1000) are not divisible by 17, R(10000) is divisible by 17. Yet there is no value of n for which \(R(10^{n})\) will divide by 19. In fact, it is remarkable that 11, 17, 41, and 73 are the only four primes below one-hundred that can be a factor of R(10n).
Find the sum of all the primes below one-hundred thousand that will never be a factor of R(10n).
This solution works based on the Multiplicative Order of 10 mod p; where p is a prime number (for this problem a prime under 100,000). For the number to never be a divisor of a RepUnit in the form \(R(10^{n})\) then the above multiplicative order has to be something that can not be written in the form \(2^{i} * 5^{i}\). This method expands to any repunit no matter how big the value of n in \(R(10^{n})\) is.
The other method is to calculate a sufficiently large power of 10 such as \(10^{20}\) and use that in pythons pow(10,\(10^{20}\),p) which is faster but if you don’t choose a large enough power of ten you’ll get bad results after a certain point.
The second method is faster but I use the first in the below solution.
#!/usr/bin/python
from math import sqrt
#prime generator code aquired from
#http://code.activestate.com/recipes/366178-a-fast-prime-number-list-generator/
def Primes(n):
if n==2: return [2]
elif n<2: return []
s=range(3,n+1,2)
mroot = n ** 0.5
half=(n+1)/2-1
i=0
m=3
while m <= mroot:
if s[i]:
j=(m*m-3)/2
s[j]=0
while j<half:
s[j]=0
j+=m
i=i+1
m=2*i+3
return [2]+[x for x in s if x]
primes = Primes(100000)
def Factor(n):
global primes
pmax = int(sqrt(n))+1
factors = []
for prime in primes:
if prime > pmax:
break
if (n % prime ==0):
while (n% prime ==0):
factors = factors + [prime]
n = n/prime
if (n != 1):
factors += [n]
return factors
#n is the number to check
#form is the list of factors
#only true if n contains only factors found in form
def OfForm(n,form):
factors = Factor(n)
if (factors == []):
return False
for factor in factors:
if (factor not in form):
return False
return True
#a^k = 1 (mod n).
def MultiplicativeOrder (a,n):
k = 1
if (a%n ==0):
return 0
r = pow(a,k,n)
while (r != 1):
k= k+1
r = pow(a,k,n)
return k
primelist = Primes(100000)
summation = 0
for p in primelist:
if not OfForm(MultiplicativeOrder (10,p),[2,5]):
print p
summation += p
print summation
Tags: c#, csharp, mono, project euler
Posted in Project Euler | Comments (0)
Problem #122 says:
We shall define m(k) to be the minimum number of multiplications to compute \(n^{k}\); for example m(15) = 5.
For 1 ≤ k ≤ 200, find ∑ m(k).
This is a problem about finding Addition Chains that lead to the desired number in the fewest steps. Using the library code I posted yesterday the solution can be written as 5 lines of code (not counting main() and other fluff). If you compile the Addition Chain library so that it doesn’t allow duplicates then the program will give the wrong answer; This is because duplicate paths are needed to ensure that an optimal chain can be found.
Other than that the solution runs in about 45 seconds an decent computer. Solution provided in c#/mono.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using AdditionChain;
namespace Problem122
{
class MainClass
{
public static void Main (string[] args)
{
AdditionChain ac = new AdditionChain(1);
int answer = 0;
for (int i = 1; i<= 200; i++)
answer += ac.NodeLevel(ac.GetNodeByValue(i));
Console.WriteLine(answer);
}
}
}
Tags: c#, programming, project euler
Posted in Project Euler | Comments (0)