#! /usr/bin/gawk -f
# Last edited on 2011-03-05 21:10:54 by stolfilocal

BEGIN{

  n = 31;
  for (dev = 3; dev < 1.5*n; dev += 0.5)
    { 
      f0 = folded_bell(0.0,dev,n);
      fm = folded_bell(0.5*n,dev,n);
      printf "%19.15f %19.15f  %19.15f %19.15f  %19.15f\n", dev, z, f0, fm, fm-f0 > "/dev/stderr";
    }
}

function folded_bell(z,dev,   kmax,sum,um,up,uz,k,big)
  {
    # Assumes {z >= 0}:
    # Assumes {n > 0}:
    big = 8.6;
    # Reduce {z} to the range {[0_n)}:
    z = z + n;
    z = z - n*int(z/n);
    while (z >= n) { z = z - n; }

    # Add all fold-over terms that are {10^{-16}} or more:
    kmax = int((big*dev + z)/n);
    sum = 0;

    # Add all terms for {z - kmax*n} to {z + kmax*n}:
    k = kmax;
    while (k > 0) { 
        um = (z - k*n)/dev; 
        up = (z + k*n)/dev; 
        sum += exp(-um*um/2) + exp(-up*up/2);
        k--;
      }
    uz = z/dev;
    sum += exp(-uz*uz/2);
    return sum;
  }