diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..2251642
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+Manifest.toml
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index dc2faec..cef11b8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "Sobol"
 uuid = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
-version = "1.5.0"
+version = "2.0.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
diff --git a/README.md b/README.md
index fe9c3da..b488bb3 100644
--- a/README.md
+++ b/README.md
@@ -96,10 +96,9 @@ where `lb` and `ub` are arrays (or other iterables) of length `N`, giving
 the lower and upper bounds of the hypercube, respectively.   For example,
 `SobolSeq([-1,0,0],[1,3,2])` generates points in the box [-1,1]×[0,3]×[0,2].  (Although the generated points are `Float64` by default, in general the precision is determined by that of `lb` and `ub`.)
 
-If you know in advance the number `n` of points that you plan to
-generate, some authors suggest that better uniformity can be attained
-by first skipping the initial portion of the LDS. In particular,
-we skip 2ᵐ−1 for the largest `m` where `2ᵐ-1 ≤ n` (see [Joe and Kuo, 2003][joe03] for a similar suggestion).  This
+If you know in advance the number `n` of points that you plan to generate,
+some authors suggest that skipping the initial portion of the LDS.
+In particular, we skip 2ᵐ points for the largest `m` where `2ᵐ ≤ n` (see [Joe and Kuo, 2003][joe03] for a similar suggestion).  This
 facility is provided by:
 ```julia
 skip(s, n)
@@ -111,8 +110,15 @@ Skipping exactly `n` elements is also possible:
 skip(s, n, exact=true)
 ```
 
+Note, however, that generally it is not recommended to skip (arbitrary)
+number of points since this can ruin the digital net property and lead
+to worse convergence (as illustrated in [Owen, 2021](https://arxiv.org/abs/2008.08051) for the
+common proactice of dropping the initial point).
+
 `skip` returns `s`, so you can simply do `s = skip(SobolSeq(N))` or similar.
 
+Note, however, that generally it is not recommended to skip arbitrary
+
 ## Example
 
 Here is a simple example, generating 1024 points in two dimensions and
@@ -123,9 +129,12 @@ the [0,1]×[0,1] unit square!
 using Sobol
 using PyPlot
 s = SobolSeq(2)
-p = reduce(hcat, next!(s) for i = 1:1024)'
+p = Matrix{Float64}(undef, 2, 1024)
+for x in eachcol(p)
+   next!(s, x)
+end
 subplot(111, aspect="equal")
-plot(p[:,1], p[:,2], "r.")
+plot(p[1, :], p[2, :], "r.")
 ```
 ![plot of 1024 points of a 2d Sobol sequence](sobol1024.png "1024 points of a 2d Sobol sequence")
 
diff --git a/sobol1024.png b/sobol1024.png
index ba253a1..1a58c37 100644
Binary files a/sobol1024.png and b/sobol1024.png differ
diff --git a/src/Sobol.jl b/src/Sobol.jl
index 2bbb87b..0317683 100644
--- a/src/Sobol.jl
+++ b/src/Sobol.jl
@@ -52,38 +52,43 @@ function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat})
     if s.n == typemax(s.n)
         return rand!(x)
     end
-
-    s.n += one(s.n)
-    c = UInt32(trailing_zeros(s.n))
-    sb = s.b
-    sx = s.x
-    sm = s.m
-    for i=1:ndims(s)
-        @inbounds b = sb[i]
-        # note: ldexp on Float64(sx[i]) is exact, independent of precision of x[i]
-        if b >= c
-            @inbounds sx[i] = sx[i] ⊻ (sm[i,c+1] << (b-c))
-            @inbounds x[i] = ldexp(Float64(sx[i]), ((~b) % Int32))
-        else
-            @inbounds sx[i] = (sx[i] << (c-b)) ⊻ sm[i,c+1]
-            @inbounds sb[i] = c
-            @inbounds x[i] = ldexp(Float64(sx[i]), ((~c) % Int32))
+    
+    if iszero(s.n)
+        # Initial point
+        fill!(x, 0.0)
+    else
+        # Subsequent points
+        c = UInt32(trailing_zeros(s.n))
+        sb = s.b
+        sx = s.x
+        sm = s.m
+        for i=1:ndims(s)
+            @inbounds b = sb[i]
+            # note: ldexp on Float64(sx[i]) is exact, independent of precision of x[i]
+            if b >= c
+                @inbounds sx[i] = sx[i] ⊻ (sm[i,c+1] << (b-c))
+                @inbounds x[i] = ldexp(Float64(sx[i]), ((~b) % Int32))
+            else
+                @inbounds sx[i] = (sx[i] << (c-b)) ⊻ sm[i,c+1]
+                @inbounds sb[i] = c
+                @inbounds x[i] = ldexp(Float64(sx[i]), ((~c) % Int32))
+            end
         end
     end
+
+    s.n += one(s.n)
+
     return x
 end
 next!(s::SobolSeq) = next!(s, Array{Float64,1}(undef, ndims(s)))
 
 # if we know in advance how many points (n) we want to compute, then
-# adopt a suggestion similar to the Joe and Kuo paper, which in turn
-# is taken from Acworth et al (1998), of skipping a number of
-# points one less than the largest power of 2 smaller than n+1.
+# we can adopt a suggestion similar to Joe and Kuo (2003), which in turn
+# is taken from Acworth et al (1998), of skipping a number of initial
+# points. The number of points is the largest power of 2 smaller than n+1.
 # if exactly n points are to be skipped, use the keyword exact=true.
-# (Ackworth and Joe and Kuo seem to suggest skipping exactly
-#  a power of 2, but skipping 1 less seems to produce much better
-#  results: issue #21.)
 #
-# skip!(s, n) skips 2^m - 1 such that 2^m < n ≤ 2^(m+1)
+# skip!(s, n) skips 2^m such that 2^m ≤ n < 2^(m+1)
 # skip!(s, n, exact=true) skips m = n
 
 function skip!(s::SobolSeq, n::Integer, x; exact=false)
@@ -91,7 +96,7 @@ function skip!(s::SobolSeq, n::Integer, x; exact=false)
         n == 0 && return s
         throw(ArgumentError("$n is not non-negative"))
     end
-    nskip = exact ? n : (1 << floor(Int,log2(n+1)) - 1)
+    nskip = exact ? n : prevpow(2, n)
     for unused=1:nskip; next!(s,x); end
     return s
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1555d98..8deff24 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,8 +21,8 @@ using Sobol, Test
             for line in eachline(exp_results_file)
                 values = [parse(Float64, item) for item in split(line)]
                 if length(values) > 0
-                    @test x == values
                     next!(s, x)
+                    @test x == values
                 end
             end
         end
@@ -32,7 +32,7 @@ end
 using Base.Iterators: take
 @testset "iterators" begin
     # issue #8
-    @test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.5,0.75,0.25,0.375,0.875]
+    @test [x[1] for x in collect(take(Sobol.SobolSeq(1),5))] == [0.0,0.5,0.75,0.25,0.375]
 end
 
 @testset "scaled" begin
@@ -44,8 +44,13 @@ end
     @test s isa ScaledSobolSeq{3,Int}
     @test eltype(s) == Vector{Float64}
     @test eltype(SobolSeq(Float32.(lb),Float32.(ub))) == Vector{Float32}
-    @test first(s) == [0,1.5,1]
-    @test first(SobolSeq((x for x in lb), (x for x in ub))) == [0,1.5,1]
+    @test next!(s) == [-1,0,0]
+    @test next!(s) == [0,1.5,1]
+
+    s = SobolSeq((x for x in lb), (x for x in ub))
+    @test next!(s) == [-1,0,0]
+    @test next!(s) == [0,1.5,1]
+
     @test SobolSeq(N,lb,ub) isa ScaledSobolSeq{3,Int}
     @test_throws BoundsError SobolSeq(2,lb,ub)
 end
@@ -71,10 +76,10 @@ end
     p = map(_->next!(s), 1:16)
     sort!(map(x -> sum((x .≥ 0.5) .* (2 .^ (0:3))), p)) == 0:15
 
-    for n in (7,8)
+    for n in (7,8,9)
         s = skip(SobolSeq(2), n)
         p = [next!(s) for _ = 1:n]
-        s2 = skip(SobolSeq(2), 7, exact=true)
+        s2 = skip(SobolSeq(2), n == 7 ? 4 : 8, exact=true)
         p2 = [next!(s2) for _ = 1:n]
         @test p == p2
     end