From 2ba11da41e25ff279c6acdc95898b9ad15ace3b7 Mon Sep 17 00:00:00 2001 From: alegresor Date: Fri, 4 Aug 2023 12:34:34 -0600 Subject: [PATCH] update lattice orderings --- Project.toml | 2 +- docs/src/tutorial.md | 170 +++++++++++++++++++++---------------------- src/digitalseqb2g.jl | 23 +++--- src/latticeseqb2.jl | 59 ++++++++------- 4 files changed, 128 insertions(+), 126 deletions(-) diff --git a/Project.toml b/Project.toml index 0af352e..7f322fe 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QMCGenerators" uuid = "fcae8fc2-2969-4f4b-b6c9-a95746d81969" authors = ["alegresor agsorokin3@gmail.com"] -version = "1.3.0" +version = "1.3.1" [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index a6d6e8f..8a11985 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -159,9 +159,9 @@ Next(ls,4) Next(ls,4) # output 4×5 Matrix{Float64}: + 0.375 0.875 0.375 0.875 0.875 0.875 0.375 0.875 0.375 0.375 0.625 0.125 0.625 0.125 0.125 - 0.375 0.875 0.375 0.875 0.875 0.125 0.625 0.125 0.625 0.625 ``` ```jldoctest tut_ls; output = false @@ -423,13 +423,12 @@ Next(ls,4) ### Linear Ordering -By default, digital sequences are generated in Gray code order. One may generate the first $2^m$ points in linear order via +By default, digital sequences are generated in Gray code order. One may generate the first 8 points in linear order via ```jldoctest tut_ds_order -m = 3 -n = 2^m +n = 8 ds = DigitalSeqB2G(4) -FirstLinear(ds,m) +FirstLinear(ds,n) # output 8×4 Matrix{Float64}: 0.0 0.0 0.0 0.0 @@ -461,19 +460,18 @@ Next(ds,n) Similarly, Lattices are by default generated in extensible ordering. A linear ordering is also available ```jldoctest tut_ls_order -m = 3 -n = 2^m +n = 8 ls = LatticeSeqB2(4) -FirstLinear(ls,m) +FirstLinear(ls,n) # output 8×4 Matrix{Float64}: 0.0 0.0 0.0 0.0 - 0.125 0.625 0.125 0.625 - 0.25 0.25 0.25 0.25 - 0.375 0.875 0.375 0.875 0.5 0.5 0.5 0.5 - 0.625 0.125 0.625 0.125 + 0.25 0.25 0.25 0.25 0.75 0.75 0.75 0.75 + 0.125 0.625 0.125 0.625 + 0.625 0.125 0.625 0.125 + 0.375 0.875 0.375 0.875 0.875 0.375 0.875 0.375 ``` @@ -487,9 +485,9 @@ Next(ls,n) 0.5 0.5 0.5 0.5 0.75 0.75 0.75 0.75 0.25 0.25 0.25 0.25 + 0.375 0.875 0.375 0.875 0.875 0.375 0.875 0.375 0.625 0.125 0.625 0.125 - 0.375 0.875 0.375 0.875 0.125 0.625 0.125 0.625 ``` @@ -498,7 +496,7 @@ Linear order for randomized sequences has expected syntax ```jldoctest tut_ds_order_rand ds = DigitalSeqB2G(LinearMatrixScramble(4,11)) rds = RandomDigitalShift(ds,1,17) -FirstLinear(rds,2) +FirstLinear(rds,4) # output 4×4 Matrix{Float64}: 0.844967 0.66901 0.686087 0.362606 @@ -509,7 +507,7 @@ FirstLinear(rds,2) ```jldoctest tut_ds_order_rand rds = RandomDigitalShift(ds,2,17) -xs = FirstRLinear(rds,2) +xs = FirstRLinear(rds,4) xs[1] # output 4×4 Matrix{Float64}: @@ -531,7 +529,7 @@ xs[2] ```jldoctest tut_ds_order_rand rds = RandomOwenScramble(ds,2,17) -xs = FirstRLinear(rds,2) +xs = FirstRLinear(rds,4) xs[1] # output 4×4 Matrix{Float64}: @@ -555,13 +553,13 @@ For lattices a similar API holds ```jldoctest tut_ds_order_rand rls = RandomShift(LatticeSeqB2(4),2,17) -xs = FirstRLinear(rls,2) +xs = FirstRLinear(rls,4) xs[1] # output 4×4 Matrix{Float64}: 0.447971 0.593933 0.190141 0.331784 - 0.697971 0.843933 0.440141 0.581784 0.947971 0.0939328 0.690141 0.831784 + 0.697971 0.843933 0.440141 0.581784 0.197971 0.343933 0.940141 0.081784 ``` @@ -570,8 +568,8 @@ xs[2] # output 4×4 Matrix{Float64}: 0.868892 0.697591 0.828991 0.475479 - 0.118892 0.947591 0.0789912 0.725479 0.368892 0.197591 0.328991 0.975479 + 0.118892 0.947591 0.0789912 0.725479 0.618892 0.447591 0.578991 0.225479 ``` @@ -580,23 +578,23 @@ xs[2] For digital sequences, we sometimes want the binary representation of points. We can get the binary representations as integers and then convert them to their floating point values as follows ```jldoctest tut_ds_binary -ds = DigitalSeqB2G(4) +ds = DigitalSeqB2G(3) xb = NextBinary(ds,4) # output -4×4 Matrix{UInt64}: - 0x0000000000000000 0x0000000000000000 … 0x0000000000000000 - 0x0000000080000000 0x0000000080000000 0x0000000080000000 - 0x00000000c0000000 0x0000000040000000 0x0000000040000000 - 0x0000000040000000 0x00000000c0000000 0x00000000c0000000 +4×3 Matrix{UInt64}: + 0x0000000000000000 0x0000000000000000 0x0000000000000000 + 0x0000000080000000 0x0000000080000000 0x0000000080000000 + 0x00000000c0000000 0x0000000040000000 0x0000000040000000 + 0x0000000040000000 0x00000000c0000000 0x00000000c0000000 ``` ```jldoctest tut_ds_binary BinaryToFloat64(xb,ds) # output -4×4 Matrix{Float64}: - 0.0 0.0 0.0 0.0 - 0.5 0.5 0.5 0.5 - 0.75 0.25 0.25 0.25 - 0.25 0.75 0.75 0.75 +4×3 Matrix{Float64}: + 0.0 0.0 0.0 + 0.5 0.5 0.5 + 0.75 0.25 0.25 + 0.25 0.75 0.75 ``` This is also compatible with randomized digital sequences @@ -606,20 +604,20 @@ Reset!(ds) rds_single = RandomDigitalShift(ds,1,11) xb = NextBinary(rds_single,4) # output -4×4 Matrix{UInt64}: - 0x001293891708f524 0x0016da39893c2ad5 … 0x00167b863bac2dff - 0x000293891708f524 0x0006da39893c2ad5 0x00067b863bac2dff - 0x000a93891708f524 0x001eda39893c2ad5 0x001e7b863bac2dff - 0x001a93891708f524 0x000eda39893c2ad5 0x000e7b863bac2dff +4×3 Matrix{UInt64}: + 0x001293891708f524 0x0016da39893c2ad5 0x001b9c72840fa196 + 0x000293891708f524 0x0006da39893c2ad5 0x000b9c72840fa196 + 0x000a93891708f524 0x001eda39893c2ad5 0x00139c72840fa196 + 0x001a93891708f524 0x000eda39893c2ad5 0x00039c72840fa196 ``` ```jldoctest tut_ds_binary BinaryToFloat64(xb,rds_single) # output -4×4 Matrix{Float64}: - 0.58051 0.714139 0.862848 0.702579 - 0.0805097 0.214139 0.362848 0.202579 - 0.33051 0.964139 0.612848 0.952579 - 0.83051 0.464139 0.112848 0.452579 +4×3 Matrix{Float64}: + 0.58051 0.714139 0.862848 + 0.0805097 0.214139 0.362848 + 0.33051 0.964139 0.612848 + 0.83051 0.464139 0.112848 ``` ```jldoctest tut_ds_binary Reset!(ds) @@ -627,88 +625,88 @@ rds_multiple = RandomDigitalShift(ds,2,11) xbs = NextRBinary(rds_multiple,4) xbs[1] # output -4×4 Matrix{UInt64}: - 0x001293891708f524 0x001b9c72840fa196 … 0x001d6f2b8e3bf612 - 0x000293891708f524 0x000b9c72840fa196 0x000d6f2b8e3bf612 - 0x000a93891708f524 0x00139c72840fa196 0x00156f2b8e3bf612 - 0x001a93891708f524 0x00039c72840fa196 0x00056f2b8e3bf612 +4×3 Matrix{UInt64}: + 0x001293891708f524 0x001b9c72840fa196 0x001b92b4f3a3bb3f + 0x000293891708f524 0x000b9c72840fa196 0x000b92b4f3a3bb3f + 0x000a93891708f524 0x00139c72840fa196 0x001392b4f3a3bb3f + 0x001a93891708f524 0x00039c72840fa196 0x000392b4f3a3bb3f ``` ```jldoctest tut_ds_binary xbs[2] # output -4×4 Matrix{UInt64}: - 0x0016da39893c2ad5 0x00167b863bac2dff … 0x001cd7a908e348e1 - 0x0006da39893c2ad5 0x00067b863bac2dff 0x000cd7a908e348e1 - 0x000eda39893c2ad5 0x001e7b863bac2dff 0x0014d7a908e348e1 - 0x001eda39893c2ad5 0x000e7b863bac2dff 0x0004d7a908e348e1 +4×3 Matrix{UInt64}: + 0x0016da39893c2ad5 0x00167b863bac2dff 0x00006a9dbda72567 + 0x0006da39893c2ad5 0x00067b863bac2dff 0x00106a9dbda72567 + 0x000eda39893c2ad5 0x001e7b863bac2dff 0x00086a9dbda72567 + 0x001eda39893c2ad5 0x000e7b863bac2dff 0x00186a9dbda72567 ``` ```jldoctest tut_ds_binary; output = false BinaryToFloat64(xbs,rds_multiple) # output 2-element Vector{Matrix{Float64}}: - [0.5805097055341872 0.8628475741689268 0.8616585501275508 0.9198205736171892; 0.08050970553418724 0.36284757416892677 0.36165855012755077 0.41982057361718916; 0.33050970553418724 0.6128475741689268 0.6116585501275508 0.6698205736171892; 0.8305097055341872 0.11284757416892677 0.11165855012755077 0.16982057361718916] - [0.7141387634631778 0.702578655765535 0.013014669814919055 0.9013257192221112; 0.21413876346317784 0.20257865576553502 0.513014669814919 0.40132571922211124; 0.46413876346317784 0.952578655765535 0.26301466981491906 0.6513257192221112; 0.9641387634631778 0.452578655765535 0.763014669814919 0.15132571922211124] + [0.5805097055341872 0.8628475741689268 0.8616585501275508; 0.08050970553418724 0.36284757416892677 0.36165855012755077; 0.33050970553418724 0.6128475741689268 0.6116585501275508; 0.8305097055341872 0.11284757416892677 0.11165855012755077] + [0.7141387634631778 0.702578655765535 0.013014669814919055; 0.21413876346317784 0.20257865576553502 0.513014669814919; 0.46413876346317784 0.952578655765535 0.26301466981491906; 0.9641387634631778 0.452578655765535 0.763014669814919] ``` ```jldoctest tut_ds_binary Reset!(ds) NextBinary(RandomOwenScramble(ds,1,11),4) # output -4×4 Matrix{UInt64}: - 0x001a5eef4cd370c6 0x00028467570b42be … 0x0001a6384b638476 - 0x000ff20fe5ead1ba 0x001cf89a1a9a2567 0x00151c3bba0d1d72 - 0x000691102ee5e771 0x0009fddc7fb98782 0x000b360939d54a00 - 0x001689e7a7bc5106 0x0014908fa8a77dff 0x001a7c08956e0175 +4×3 Matrix{UInt64}: + 0x001a5eef4cd370c6 0x00028467570b42be 0x0003828fa0349158 + 0x000ff20fe5ead1ba 0x001cf89a1a9a2567 0x001e02e9de330269 + 0x000691102ee5e771 0x0009fddc7fb98782 0x000bb5aca8044288 + 0x001689e7a7bc5106 0x0014908fa8a77dff 0x00158be9ceaeef56 ``` Getting binary points with linear ordering is also supported. ```jldoctest tut_ds_binary Reset!(ds) # resets rds_single and rds_multiple as well -FirstLinearBinary(ds,2) +FirstLinearBinary(ds,4) # output -4×4 Matrix{UInt64}: - 0x0000000000000000 0x0000000000000000 … 0x0000000000000000 - 0x0000000080000000 0x0000000080000000 0x0000000080000000 - 0x0000000040000000 0x00000000c0000000 0x00000000c0000000 - 0x00000000c0000000 0x0000000040000000 0x0000000040000000 +4×3 Matrix{UInt64}: + 0x0000000000000000 0x0000000000000000 0x0000000000000000 + 0x0000000080000000 0x0000000080000000 0x0000000080000000 + 0x0000000040000000 0x00000000c0000000 0x00000000c0000000 + 0x00000000c0000000 0x0000000040000000 0x0000000040000000 ``` ```jldoctest tut_ds_binary -FirstLinearBinary(rds_single,2) +FirstLinearBinary(rds_single,4) # output -4×4 Matrix{UInt64}: - 0x001293891708f524 0x0016da39893c2ad5 … 0x00167b863bac2dff - 0x000293891708f524 0x0006da39893c2ad5 0x00067b863bac2dff - 0x001a93891708f524 0x000eda39893c2ad5 0x000e7b863bac2dff - 0x000a93891708f524 0x001eda39893c2ad5 0x001e7b863bac2dff +4×3 Matrix{UInt64}: + 0x001293891708f524 0x0016da39893c2ad5 0x001b9c72840fa196 + 0x000293891708f524 0x0006da39893c2ad5 0x000b9c72840fa196 + 0x001a93891708f524 0x000eda39893c2ad5 0x00039c72840fa196 + 0x000a93891708f524 0x001eda39893c2ad5 0x00139c72840fa196 ``` ```jldoctest tut_ds_binary -xbs = FirstRLinearBinary(rds_multiple,2) +xbs = FirstRLinearBinary(rds_multiple,4) xbs[1] # output -4×4 Matrix{UInt64}: - 0x001293891708f524 0x001b9c72840fa196 … 0x001d6f2b8e3bf612 - 0x000293891708f524 0x000b9c72840fa196 0x000d6f2b8e3bf612 - 0x001a93891708f524 0x00039c72840fa196 0x00056f2b8e3bf612 - 0x000a93891708f524 0x00139c72840fa196 0x00156f2b8e3bf612 +4×3 Matrix{UInt64}: + 0x001293891708f524 0x001b9c72840fa196 0x001b92b4f3a3bb3f + 0x000293891708f524 0x000b9c72840fa196 0x000b92b4f3a3bb3f + 0x001a93891708f524 0x00039c72840fa196 0x000392b4f3a3bb3f + 0x000a93891708f524 0x00139c72840fa196 0x001392b4f3a3bb3f ``` ```jldoctest tut_ds_binary xbs[2] # output -4×4 Matrix{UInt64}: - 0x0016da39893c2ad5 0x00167b863bac2dff … 0x001cd7a908e348e1 - 0x0006da39893c2ad5 0x00067b863bac2dff 0x000cd7a908e348e1 - 0x001eda39893c2ad5 0x000e7b863bac2dff 0x0004d7a908e348e1 - 0x000eda39893c2ad5 0x001e7b863bac2dff 0x0014d7a908e348e1 +4×3 Matrix{UInt64}: + 0x0016da39893c2ad5 0x00167b863bac2dff 0x00006a9dbda72567 + 0x0006da39893c2ad5 0x00067b863bac2dff 0x00106a9dbda72567 + 0x001eda39893c2ad5 0x000e7b863bac2dff 0x00186a9dbda72567 + 0x000eda39893c2ad5 0x001e7b863bac2dff 0x00086a9dbda72567 ``` ```jldoctest tut_ds_binary Reset!(ds) -FirstLinearBinary(RandomOwenScramble(ds,1,11),2) +FirstLinearBinary(RandomOwenScramble(ds,1,11),4) # output -4×4 Matrix{UInt64}: - 0x001a5eef4cd370c6 0x00028467570b42be … 0x0001a6384b638476 - 0x000ff20fe5ead1ba 0x001cf89a1a9a2567 0x00151c3bba0d1d72 - 0x001691102ee5e771 0x0011fddc7fb98782 0x001b360939d54a00 - 0x000689e7a7bc5106 0x000c908fa8a77dff 0x000a7c08956e0175 +4×3 Matrix{UInt64}: + 0x001a5eef4cd370c6 0x00028467570b42be 0x0003828fa0349158 + 0x000ff20fe5ead1ba 0x001cf89a1a9a2567 0x001e02e9de330269 + 0x001691102ee5e771 0x0011fddc7fb98782 0x0013b5aca8044288 + 0x000689e7a7bc5106 0x000c908fa8a77dff 0x000d8be9ceaeef56 ``` These may be converted to floats as before. diff --git a/src/digitalseqb2g.jl b/src/digitalseqb2g.jl index 2933af1..359e062 100644 --- a/src/digitalseqb2g.jl +++ b/src/digitalseqb2g.jl @@ -53,7 +53,7 @@ mutable struct DigitalSeqB2G const n::Int64 # maximum number of supported points const recipd::Union{BigFloat,Float64} # multiplication factor k::Int64 # index in the sequence - cur::Vector{UInt64} + cur::Vector{UInt64} # current binary point end function DigitalSeqB2G(s::Int64,Cs::Matrix{BigInt}) @@ -87,9 +87,9 @@ function Reset!(seq::DigitalSeqB2G) return end -function NextBinaryLow(i1::Int64,s::Int64,xb::Matrix{UInt64},n::Int64,k::Int64,Csr::Matrix{UInt64},cur::Vector{UInt64}) +function NextBinaryLow(i1::UInt64,s::Int64,xb::Matrix{UInt64},n::Int64,k::Int64,Csr::Matrix{UInt64},cur::Vector{UInt64}) b::Int64 = 0 - for i::Int64=i1:n + for i::UInt64=i1:n k += 1 b = 0; while ~Bool((k>>b)&1) b+= 1 end for j=1:s xb[i,j] = cur[j] ⊻= Csr[j,b+1] end @@ -99,7 +99,7 @@ end function NextBinary(seq::DigitalSeqB2G,n::Int64) (seq.k+n)>=seq.n && throw(DomainError(n,"Generating $n more points will exceed the maximum supported points $(seq.n)")) xb::Matrix{UInt64} = zeros(UInt64,n,seq.s) - NextBinaryLow(seq.k==-1 ? 2 : 1,seq.s,xb,n,seq.k==-1 ? 0 : seq.k,seq.Csr,seq.cur) + NextBinaryLow(seq.k==-1 ? UInt64(2) : UInt64(1),seq.s,xb,n,seq.k==-1 ? 0 : seq.k,seq.Csr,seq.cur) seq.k += n seq.cur .= xb[n,:] xb @@ -109,7 +109,7 @@ Next(seq::DigitalSeqB2G,n::Int64) = BinaryToFloat64(NextBinary(seq,n),seq) function FirstLinearBinaryLow(s::Int64,xb::Matrix{UInt64},n::Int64,k::Int64,Csr::Matrix{UInt64},cur::Vector{UInt64}) b::Int64 = 0 - for i::Int64=1:n-1 + for i::UInt64=1:n-1 igc = (i⊻(i>>1))+1 k += 1 b = 0; while ~Bool((k>>b)&1) b+= 1 end @@ -117,8 +117,8 @@ function FirstLinearBinaryLow(s::Int64,xb::Matrix{UInt64},n::Int64,k::Int64,Csr: end end -function FirstLinearBinary(seq::DigitalSeqB2G,m::Int64) - n = 2^m +function FirstLinearBinary(seq::DigitalSeqB2G,n::Int64) + @assert ispow2(n) n>=seq.n && throw(DomainError(n,"Generating $n more points will exceed the maximum supported points $(seq.n)")) @assert seq.k==-1 xb::Matrix{UInt64} = zeros(UInt64,n,seq.s) @@ -127,7 +127,7 @@ function FirstLinearBinary(seq::DigitalSeqB2G,m::Int64) xb end -FirstLinear(seq::DigitalSeqB2G,m::Int64) = BinaryToFloat64(FirstLinearBinary(seq,m),seq) +FirstLinear(seq::DigitalSeqB2G,n::Int64) = BinaryToFloat64(FirstLinearBinary(seq,n),seq) struct RandomDigitalShift name::String @@ -171,9 +171,8 @@ function NextRBinary(rds::RandomDigitalShift,n::Int64) xr end -function FirstRLinearBinary(rds::RandomDigitalShift,m::Int64) - n = 2^m - xu = FirstLinearBinary(rds.seq,m) +function FirstRLinearBinary(rds::RandomDigitalShift,n::Int64) + xu = FirstLinearBinary(rds.seq,n) xr = [zeros(UInt64,n,rds.seq.s) for k=1:rds.r] NextRLowRandomDigitalShift(rds.tdiff,rds.seq.s,rds.r,n,xu,xr,rds.rshifts) xr @@ -278,4 +277,4 @@ end NextRBinary(rds::RandomOwenScramble,n::Int64) = OwenScramble(rds,NextBinary(rds.seq,n),n) -FirstRLinearBinary(rds::RandomOwenScramble,m::Int64) = OwenScramble(rds,FirstLinearBinary(rds.seq,m),2^m) +FirstRLinearBinary(rds::RandomOwenScramble,n::Int64) = OwenScramble(rds,FirstLinearBinary(rds.seq,n),n) diff --git a/src/latticeseqb2.jl b/src/latticeseqb2.jl index ace3d34..2a543e0 100644 --- a/src/latticeseqb2.jl +++ b/src/latticeseqb2.jl @@ -4,18 +4,20 @@ mutable struct LatticeSeqB2 const z::Vector{UInt64} const m::Int64 # 2^m points supported const n::Int64 # 2^m - const scale::Union{BigFloat,Float64} # 2^(-m) + const scale::Float64 # 2^(-m) k::Int64 # index in the sequence + v::UInt64 # current bitreversed binary end function LatticeSeqB2(s::Int64,z::Vector{BigInt},m::Int64) @assert m > 0 && s <= length(z) + @assert m<=53 n = 2^m - scale = m>53 ? BigFloat(2)^(-m) : Float64(2)^(-m) + scale = Float64(2)^(-m) t = maximum(ndigits.(z[1:s],base=2)) @assert t<=64 z = convert.(UInt64,z[1:s]) - LatticeSeqB2("Lattice Seq B2",s,z,m,n,scale,-1) + LatticeSeqB2("Lattice Seq B2",s,z,m,n,scale,-1,UInt64(0)) end LatticeSeqB2(s::Int64,path::String,m::Int64) = LatticeSeqB2(s,readdlm(download(joinpath("https://bitbucket.org/dnuyens/qmc-generators/raw/cb0f2fb10fa9c9f2665e41419097781b611daa1e/LATSEQ/",path)),BigInt)[:,1],m) @@ -24,45 +26,49 @@ LatticeSeqB2(s::Int64) = LatticeSeqB2(s,DEFAULT_LATTICESEQB2_GVECTOR,20) function Reset!(seq::LatticeSeqB2) seq.k = -1 + seq.v = 0 return end -function NextLow(s::Int64,xb::Matrix{Float64},k::Int64,n::Int64,z::Vector{UInt64}) - idx::Int64,p2::Int64,frac::Float64 = 0,0,0 - for i=1:n - idx = k+i - p2 = 2^ceil(Int64,log2(idx+1)) - frac = (2*(p2-idx)-1)/p2 - for j=1:s - xb[i,j] = convert(Float64,(frac*z[j])%1) - end - end +function NextLow(i1::UInt64,v::UInt64,xb::Matrix{Float64},s::Int64,n::Int64,k::Int64,z::Vector{UInt64},scale::Float64,m::Int64) + b::Int64=0; frac::Float64=0 + for i::UInt64=i1:n + k += 1 + b = 0; while ~Bool((k>>b)&1) b+= 1 end + v ⊻= UInt64(1)<<(m-1-b) + frac = scale*v + for j=1:s xb[i,j] = convert(Float64,(frac*z[j])%1) end + end + v end function Next(seq::LatticeSeqB2,n::Int64) (seq.k+n)>=seq.n && throw(DomainError(n,"Generating $n more points will exceed the maximum supported points $(seq.n)")) xb::Matrix{Float64} = zeros(Float64,n,seq.s) - NextLow(seq.s,xb,seq.k,n,seq.z) + seq.v = NextLow(seq.k==-1 ? UInt64(2) : UInt64(1),seq.v,xb,seq.s,n,seq.k==-1 ? 0 : seq.k,seq.z,seq.scale,seq.m) seq.k += n xb end -function FirstLinearLow(xb::Matrix{Float64},s::Int64,n::Int64,z::Vector{UInt64}) - frac::Float64 = 0 - for i=1:n - frac = (i-1)/n - for j=1:s - xb[i,j] = convert(Float64,(frac*z[j])%1) - end +function FirstLinearLow(xb::Matrix{Float64},s::Int64,n::Int64,k::Int64,z::Vector{UInt64},scale::Float64,m::Int64) + b::Int64=0; v::UInt64=0; frac::Float64=0 + for i::UInt64=1:n-1 + igc = (i⊻(i>>1))+1 + k += 1 + b = 0; while ~Bool((k>>b)&1) b+= 1 end + v ⊻= UInt64(1)<<(m-1-b) + frac = scale*v + for j=1:s xb[igc,j] = convert(Float64,(frac*z[j])%1) end end end -function FirstLinear(seq::LatticeSeqB2,m::Int64) +function FirstLinear(seq::LatticeSeqB2,n::Int64) @assert seq.k == -1 - n = 2^m + @assert ispow2(n) (seq.k+n)>=seq.n && throw(DomainError(n,"Generating $n more points will exceed the maximum supported points $(seq.n)")) xb::Matrix{Float64} = zeros(Float64,n,seq.s) - FirstLinearLow(xb,seq.s,n,seq.z) + FirstLinearLow(xb,seq.s,n,seq.k==-1 ? 0 : seq.k,seq.z,seq.scale,seq.m) + Reset!(seq) xb end @@ -98,9 +104,8 @@ function NextR(rls::RandomShift,n::Int64) xr end -function FirstRLinear(rls::RandomShift,m::Int64) - n = 2^m - xu = FirstLinear(rls.seq,m) +function FirstRLinear(rls::RandomShift,n::Int64) + xu = FirstLinear(rls.seq,n) xr = [zeros(Float64,n,rls.seq.s) for k=1:rls.r] NextRLowRandomShift(rls.seq.s,rls.r,n,xu,xr,rls.rshifts) xr