Project Euler 47

http://projecteuler.net/index.php?section=problems&id=47

素数関係のコードをいちいち書いたりコピペするのが面倒なので、別ファイルにします。そのときに、名前空間とモジュールを使います。

// mod1.fs
namespace NS1
module Mod1 =
   let a = 1

// mod2.fs
namespace NS1
module Mod2 =
   let b = 2

名前空間の中でモジュールを定義し、その中で定数や関数を定義します。
これをそれぞれコンパイルします。

>fsc mod1.fs --target:library

>fsc mod2.fs --target:library

そうすると、mod1.dllとmod2.dllが生成されます。
これらを利用したコードを書きます。

// test1.fs
printfn "%d" NS1.Mod1.a     // 1
printfn "%d" NS1.Mod2.b     // 2

これを次のようにコンパイルします。

>fsc test1.fs -r:mod1.dll -r:mod2.dll

test1.exeが生成されます。

open Arithmetic

let count n = Seq.initInfinite ((+) n)

let num_facts = seq {
   for n in count 2 -> List.length (Primes.factorize n)
}

let N = 4
let s = seq {
   let counter = ref 0
   for n, m in Seq.zip (count 2) num_facts do
      if m = N then
         counter := !counter + 1
         if !counter = N then
            yield n - N + 1
      else
         counter := 0
}

printfn "%d" (Seq.head s)

namespace Arithmetic
module Primes =
   let count n = Seq.initInfinite ((+) n)
   
   let rec calc_exp n p =
      if n % p = 0 then
         let t = calc_exp (n / p) p
         (fst t + 1, snd t)
      else
         (0, n)
   
   // prime numbers
   let is_prime n =
      let rec is_prime_core n k =
         let p = if k > 1 then k * 3 - (k &&& 1) - 1 else k + 2
         if p * p > n then true
         else if n % p = 0 then false
         else is_prime_core n (k + 1)
      is_prime_core n 0
   
   let sieve max_n =
      let a = Array.create (max_n + 1) true
      for p in Seq.takeWhile (fun p -> p * p <= max_n)
                  (Seq.filter (fun n -> a.[n]) (count 2)) do
         for k in seq { p*2..p..max_n } do
            a.[k] <- false
      a
   
   let make_primes max_n =
      let a = sieve max_n
      Seq.toList (Seq.filter (fun n -> a.[n]) (seq { 2..max_n }))
   
   // factorization
   let rec fac n k =
      let p = if k > 1 then k * 3 - (k &&& 1) - 1 else k + 2
      if n = 1 then
         []
      else if p * p > n then
         [(n, 1)]
      else
         let t = calc_exp n p
         if fst t > 0 then
            (p, fst t) :: (fac (snd t) (k + 1))
         else
            fac n (k + 1)
   
   let rec pow n e = 
      if e = 0 then 1
      else
         let m = pow n (e / 2)
         if e % 2 = 1 then m * m * n else m * m
   
   let factorize n = fac n 0
   let value f = List.fold (fun x y -> x * (pow (fst y) (snd y))) 1 f
   let rec div_f f g =
      if f = [] then
         []
      else if g = [] then
         f
      else
         let f0 = List.head f
         let g0 = List.head g
         if fst f0 = fst g0 then
            if snd f0 = snd g0 then
               div_f (List.tail f) (List.tail g)
            else
               (fst f0, (snd f0) - (snd g0)) :: (div_f (List.tail f) (List.tail g))
         else
            div_f (List.tail f) g
   
   // divisor
   let rec divisors = function
         | [] -> [[]]
         | (p,e)::fs -> 
               [ for t in (divisors fs) do
                     for e' in [0..e]
                        -> if e' = 0 then t else (p, e')::t ]
   let divs n = divisors (factorize n)