Distribution of Random hyperharmonic series

(*The random hyperharmonic series is the infinite series S = Sum(1,inf,d(i)/i^pow),
where integer pow > 1, and d(i) are independent, identically distributed 
random variables with property P(d(i)=0) = P(d(i)=1) = 0.5. The fact that 
for even pows the hyperharmonic series converge to some value at [0, ζ(pow)] with prob 1.
To build cumulative function F(x) = P(S < x) the tail of a series replaced by an
appropriate normal random variable. Complexity ~ O(2^k), although faster calculations 
are possible, the program shows connections between number theory and probability theory.
Expressions employed here are quite simple, therefore many improvements can be made easily*)

#r "System.Windows.Forms.DataVisualization.dll"
open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting

// some primitive actions
let inline square x = x * x
let inline double x = x + x
let inline half x = x*0.5
let inline inverse x = 1.0/x
let sign x = if x=0.0 then x else x/(abs x)

//Binomial coeffs
type binom =
     static member fact n = 
            let rec loop acc = function
                | n when n = 0 -> acc
                | n -> loop (n*acc) (n-1)
            loop 1 n
   
     static member comb n m = 
            binom.fact n / (binom.fact m * binom.fact (n-m))

//Bernoulli numbers
let rec bernoulli = function
    | n when n = 0 -> 1.0
    | n when n = 1 -> -0.5
    | n when n % 2 = 1 -> 0.0
    | n -> -([0..(n-1)] 
                  |> List.map (fun k -> 
                          (binom.comb (n+1) k |> float, bernoulli k))
                  |> List.sumBy (fun x -> fst x * snd x)   
            ) / (float (n+1))

//optional upper limit for a sum                
type UpperLimit =
   | Infinity
   | Integer of int

//Hyperharmonic series of even power   
type HypHarmonic(pow) =
   let m = pow/2
   let repl = 
       let rec loop f n =
           if n = 1 then f
           else loop (f>>f) (n/2)
       loop square (pow/2) 
 
   member h.pow = pow
   //terms
   member h.a i = i |> repl |> inverse
 
   static member (*) (a:HypHarmonic,b:HypHarmonic) = HypHarmonic(a.pow * b.pow)
   static member pi2 = double System.Math.PI 
   static member sign n = let k = n % 2 in 1 - double k
   //sum
   member h.sum = function
          | Infinity -> let s = HypHarmonic.sign (m+1) |> float in
                        let p = pown HypHarmonic.pi2 pow in
                        let f = binom.fact pow |> float in
                        s*p*(bernoulli pow) / f |> half
          | Integer i -> [float i..(-1.0)..1.0] |> List.sumBy h.a
   //remainder of a sum
   member h.rest n = h.sum Infinity - h.sum n

//calculate exact distribution for first "apr" terms
//and approximate distribution for the remainder
let distribApprox pow apr x =
  let k = (float (pown 2.0 apr)) in
  let n = Integer apr in                    //number of terms to sum
  let h = HypHarmonic pow in                 //series of even power pow S(n) = 1/1^pow + 1/2^pow + .. +1/n^pow
  let M (h:HypHarmonic) = h.rest n |> half   //remainder of the h's series. It's expectation  Eh
  let m = M h
  let s = M (square h) |> half |> sqrt      //remainder of the h's series. It's stdev  sqrt(Vh)

  //Normal distribution
  let rec erfc x =
      if x < 0.0 then 2.0 - erfc (-x)
      elif x<0.5 then 1.0 - erf x
      elif x>=10.0 then 0.0
      else 
         let P =x*(x*(x*(x*(x*(x*(x*(0.5641877825507397413087057563)
                 + 9.675807882987265400604202961)
                  + 77.08161730368428609781633646)
                   + 368.5196154710010637133875746)
                    + 1143.262070703886173606073338)
                     + 2320.439590251635247384768711)
                      + 2898.0293292167655611275846)
                       + 1826.3348842295112592168999                      
         let Q = x*(x*(x*(x*(x*(x*(x*(1.0 + 17.14980943627607849376131193)
                 + 137.1255960500622202878443578)
                  + 661.7361207107653469211984771)
                   + 2094.384367789539593790281779)
                    + 4429.612803883682726711528526)
                      + 6089.5424232724435504633068)
                       + 4958.82756472114071495438422)
                        + 1826.3348842295112595576438       
         exp (-(square x))*P/Q
  and erf x =
      let z = abs x
      let S = x |> sign |> float
      if z >= 10.0 then S
      else
         if z<0.5 then
          let Xsq = square x 
          let P =Xsq*(Xsq*(Xsq*(Xsq*(Xsq*(Xsq*0.007547728033418631287834
                    + 0.288805137207594084924010)
                     + 14.3383842191748205576712)
                      + 38.0140318123903008244444)
                       + 3017.82788536507577809226)
                        + 7404.07142710151470082064)
                         + 80437.3630960840172832162
          let Q = Xsq*(Xsq*(Xsq*(Xsq*(Xsq + 38.0190713951939403753468)
                     + 658.070155459240506326937)
                      + 6379.60017324428279487120)
                       + 34216.5257924628539769006)
                        + 80437.3630960840172826266
          S*1.1283791670955125738961589031*z*P/Q
         else
            S*(1.0-erfc z) 
            
  let pnormStd t = 
      let sqrt2 = 1.41421356237309504880
      half (erf (t / sqrt2)+1.0) 
  
  let standartize m s x = (x - m) / s in
  //cumulative normal distribution N(m,s)
  let pnorm x m s = x |> standartize m s |> pnormStd 
  
  //builds mediate distribution table for sum S(i-1) + a(i)
  let combine L i =
    let stat = L 
                |> List.fold (fun acc an -> 
                                  let b = h.a (float i) + an in                                
                                  match b with 
                                    | b when b + h.rest (Integer i) < x -> acc        //never reach x -> ignore
                                    | b when b < x -> b::fst acc, snd acc             //lt x -> append to table
                                    | _ -> fst acc, 1.0/float(pown 2.0 i) + snd acc)  //gt x -> add to prob estim
                             (L,0.0)
    //exact table, crude prob estimation
    fst stat, snd stat
   
  //build comlete distribution table for sum S(apr)
  let stat = [1..apr] 
                |> List.fold 
                        (fun acc i -> let c = combine (fst acc) i
                                      (fst c), snd c + snd acc
                        ) ([0.0],0.0)
  
  //compute sum of exact and normal distributions
  fst stat
    |> List.sumBy (fun v -> (1. - pnorm x (m+v) s))  //refine
    |> (*) (inverse k)
    |> (+) (snd stat) 

//comulative distribution for Sum(1,inf, d(i)/i^pow), 
//where P(d(i)=0) = P(d(i)=1) = 0.5 
//tol - converge criteria 10^(-tol) 
let distrib pow tol x = 
   let rec loop cur prev n =
      if abs (prev-cur)>pown 0.1 tol then loop (distribApprox pow n x) cur (n+1)
      else cur
   1. - loop 0.0 1.0 1         

let h = 0.001 //grid step
let mesh = [|0.0..h..1.7|] //grid
//integral distribution
let intF = mesh |> Array.Parallel.map (distrib 2 8) 

//differential distribution
let diffF = [| for i in 0..(intF.Length-1) do 
                  if i=0 then 
                     yield 0.0
                  else
                     yield (intF.[i]-intF.[i-1])/h 
            |]


// Create a chart 
let chart = new Chart(Dock = DockStyle.Fill)
let form = new Form(Visible = true, Width = 700, Height = 500)
chart.ChartAreas.Add(new ChartArea("MainArea"))
form.Controls.Add(chart)

// Create series and add it to the chart
let seriesInt = new Series(ChartType = SeriesChartType.Line)
chart.Series.Add(seriesInt)
let seriesDiff = new Series(ChartType = SeriesChartType.Line)
chart.Series.Add(seriesDiff)
// Specify data for the series using data-binding
seriesInt.Points.DataBindXY(mesh,intF)
seriesDiff.Points.DataBindXY(mesh,diffF)