WIP: more accurate KSDist function

This is a work-in-progress implementation of Simard and L'Ecuyer (2011) meta-algorithm <http://www.jstatsoft.org/v39/i11/paper>.
The current commit contains debug code and duplicated functions in order to check the differences between the current and the proposed implementation.

I've used KolmogorovSmirnovDist.c as a reference, which is an implementation of the meta-algorithm. <http://simul.iro.umontreal.ca/ksdir/>
At some inputs, there seems to be serious dispersion between the current one and the reference.
For example when n = 5 and x = 2.256281e-01 the CDFs are:

reference | current | proposed
-------------|-----------|--------------
8.852118e-02 | 9.711784e-02 | 9.711784e-02

In this case, the `cdf_durbin` function is in charge, which is reused in the proposed implementation.

For other example when n = 1000000 and x = -3.366834e-02 the CDFs are completely different:

reference | current | proposed
-------------|----------|-------------
1.000000e+00 | 0.000000e+00 | 1.000000e+00

What I'm wondering are:

* in the first example, is the current implementation wrong?
* in the second example, returning zero looks to be sensible, am I right?

Here is the complete log of comparison: [out.txt](https://dl.dropboxusercontent.com/u/25281592/out.txt)
(**CAUTION - it's 696KB**)

And here is the way for quick checking and reproducing it:

```
# check out my branch
curl -O -O http://simul.iro.umontreal.ca/ksdir/KolmogorovSmirnovDist.{c,h}
clang -fPIC -shared -dynamiclib KolmogorovSmirnovDist.c -o libksd.dylib
julia check.jl > out.txt
```

check.jl:
```julia
using Distributions

function _cdf(d, x)
    ccall((:KScdf, "libksd"), Float64, (Int, Float64), d.n, x)
end

begin
    ns = vcat([1:50], [100, 200, 500, 1000, 2000, 10_000, 50_000, 100_000, 500_000, 1_000_000])
    xs = linspace(-0.1, 1.1, 200)
    n_all = 0
    n_eq1 = 0
    n_eq2 = 0
    println("n\tx\treference\tcurrent\tproposed\tcurrent_ok\tproposed_ok")
    for n in ns
        ksd = Distributions.KSDist(n)
        for x in xs
            v = _cdf(ksd, x)
            v1 = Distributions.cdf(ksd, x)
            v2 = Distributions.cdf2(ksd, x)
            eq1 = isapprox(v, v1)
            eq2 = isapprox(v, v2)
            @printf "%d\t%12.6e\t%12.6e\t%12.6e\t%12.6e\t%s\t%s\n" n x v v1 v2 (eq1 ? "o" : "x") (eq2 ? "o" : "x")

            n_all += 1
            if eq1
                n_eq1 += 1
            end
            if eq2
                n_eq2 += 1
            end
        end
    end

    println(STDERR, "1: $n_eq1/$n_all")
    println(STDERR, "2: $n_eq2/$n_all")
end
```

Thank you.