BitBasis with number of qubits greater than 64 (#49)

* update

* fix tests

* update

* update shift

* update

Author: GiggleLiu

src/BitBasis.jl

@@ -1,12 +1,14 @@
 module BitBasis
 
 export DitStr, BitStr, @bit_str, @lbit_str, @dit_str, @ldit_str, BitStr64, LongBitStr, DitStr64, LongDitStr, bit_literal
-export onehot, buffer
+export onehot, buffer, indicator
 export bitarray, basis, packbits, bfloat, bfloat_r, bint, bint_r, flip
 export anyone, allone, bmask, baddrs, readbit, setbit, controller
 export swapbits, ismatch, neg, breflect, btruncate
+export LongLongUInt
 
 include("utils.jl")
+include("longlonguint.jl")
 include("DitStr.jl")
 include("BitStr.jl")
 include("bit_operations.jl")

src/BitStr.jl

@@ -36,7 +36,7 @@ julia> bit"1101"[2]
 """
 const BitStr{N,T} = DitStr{2,N,T}
 const BitStr64{N} = BitStr{N,Int64}
-const LongBitStr{N} = BitStr{N,BigInt}
+const LongBitStr{N} = BitStr{N,LongLongUInt{C}} where C
 const BitStr{N}(x::T) where {N,T<:IntStorage} = DitStr{2,N}(x)
 const BitStr(x::Union{Vector,Tuple}) = DitStr{2}(x)
 
@@ -99,7 +99,7 @@ end
 Long bit string version of `@bit_str` macro.
 """
 macro lbit_str(str)
-    return _parse_dit(Val(2), BigInt, str)
+    return _parse_dit(Val(2), LongLongUInt, str)
 end
 
 """

src/DitStr.jl

@@ -1,4 +1,4 @@
-const UIntStorage = Union{UInt8,UInt16,UInt32,UInt64,UInt128}
+const UIntStorage = Union{UInt8,UInt16,UInt32,UInt64,UInt128,LongLongUInt}
 const IntStorage = Union{Int8,Int16,Int32,Int64,Int128,BigInt,UIntStorage}
 
 ########## DitStr #########
@@ -37,23 +37,27 @@ function DitStr{D,T}(vector::Union{AbstractVector,Tuple}) where {D,T}
     val = zero(T)
     D_power_k = one(T)
     for k in 1:length(vector)
-        0 <= vector[k] <= D-1 || error("expect 0 or 1, got $(vector[k])")
-        val += vector[k] * D_power_k
-        D_power_k *= D
+        0 <= vector[k] <= D-1 || error("expect 0-$(D-1), got $(vector[k])")
+        val = accum(Val{D}(), val, vector[k], D_power_k)
+        D_power_k = _lshift(Val{D}(), D_power_k, 1)
     end
     return DitStr{D,length(vector),T}(val)
 end
+# val += x * y
+accum(::Val{D}, val, x, y) where D = val + x * y
+accum(::Val{2}, val, x, y) = iszero(x) ? val : val ⊻ y
 DitStr{D}(vector::Tuple{T,Vararg{T,N}}) where {N,T,D} = DitStr{D,T}(vector)
 DitStr{D}(vector::AbstractVector{T}) where {D,T} = DitStr{D,T}(vector)
 DitStr{D,N,T}(val::DitStr) where {D,N,T<:Integer} = convert(DitStr{D,N,T}, val)
 DitStr{D,N,T}(val::DitStr{D,N,T}) where {D,N,T<:Integer} = val
 
 const DitStr64{D,N} = DitStr{D,N,Int64}
-const LongDitStr{D,N} = DitStr{D,N,BigInt}
+const LongDitStr{D,N} = DitStr{D,N,LongLongUInt{C}} where C
+LongDitStr{D}(vector::AbstractVector{T}) where {D,T} = DitStr{D,longinttype(length(vector), D)}(vector)
 
 Base.show(io::IO, ditstr::DitStr{D,N,<:Integer}) where {D,N} =
     print(io, string(buffer(ditstr), base = D, pad = N), " ₍$('₀'+D)₎")
-Base.show(io::IO, ditstr::DitStr{D,N,<:BigInt}) where {D,N} =
+Base.show(io::IO, ditstr::DitStr{D,N,<:LongLongUInt}) where {D,N} =
     print(io, join(map(string, [ditstr[end:-1:1]...])), " ₍$('₀'+D)₎")
 
 Base.zero(::Type{DitStr{D,N,T}}) where {D,N,T} = DitStr{D,N,T}(zero(T))
@@ -295,20 +299,29 @@ function parse_dit(::Type{T}, str::String) where {T<:Integer}
 end
 
 function _parse_dit(::Val{D}, ::Type{T}, str::AbstractString) where {D, T<:Integer}
+    TT = T <: LongLongUInt ? longinttype(count(isdigit, str), D) : T
+    _parse_dit_safe(Val(D), TT, str)
+end
+
+function _parse_dit_safe(::Val{D}, ::Type{T}, str::AbstractString) where {D, T<:Integer}
     val = zero(T)
-    k = 1
-    maxk = T <: BigInt ? Inf : log(typemax(T))/log(D)
+    k = 0
+    maxk = max_num_elements(T, D)
     for each in reverse(str)
-        k >= maxk && error("string length is larger than $(maxk), use @ldit_str instead")
+        k >= maxk-1 && error("string length is larger than $(maxk), use @ldit_str instead")
         v = each - '0'
         if 0 <= v < D
-            val += _lshift(Val(D), v, k-1)
+            val += _lshift(Val(D), T(v), k)
             k += 1
         elseif each == '_'
             continue
         else
-            error("expect char in range 0-$(D-1), got $each at $k-th dit")
+            error("expect char in range 0-$(D-1), got $each at $(k+1)-th dit")
         end
     end
-    return DitStr{D,k-1,T}(val)
+    return DitStr{D,k,T}(val)
 end
+
+max_num_elements(::Type{T}, D::Int) where T<:Integer = floor(Int, log(typemax(T))/log(D))
+max_num_elements(::Type{BigInt}, D::Int) = typemax(Int)
+max_num_elements(::Type{LongLongUInt{C}}, D::Int) where {C} = max_num_elements(UInt, D) * C

src/bit_operations.jl

@@ -37,12 +37,19 @@ packbits(arr::AbstractArray) = _packbits(arr)
 _packbits(arr) =
     selectdim(sum(mapslices(x -> x .* (1 .<< (0:size(arr, 1)-1)), arr, dims = 1), dims = 1), 1, 1)
 
+"""
+    indicator(::Type{T}, k) -> T
+
+Return an integer with `k`-th bit set to 1.
+"""
+indicator(::Type{T}, k::Int) where T<:Integer = one(T) << (k-1)
+
 """
     btruncate(b, n)
 
 Truncate bits `b` to given length `n`.
 """
-btruncate(b::Integer, n) = b & (1 << n - 1)
+btruncate(b::T, n) where T<:Integer = b & (one(T) << n - 1)
 
 """
     bfloat(b::Integer; nbits::Int=bit_length(b)) -> Float64
@@ -93,7 +100,7 @@ bmask(args...) = bmask(Int, args...)
 bmask(::Type{T}) where {T<:Integer} = zero(T)
 bmask(::Type{T}, positions::Int...) where {T<:Integer} = bmask(T, positions)
 bmask(::Type{T}, itr) where {T<:Integer} =
-    isempty(itr) ? 0 : reduce(+, one(T) << (b - 1) for b in itr)
+    isempty(itr) ? 0 : reduce(+, indicator(T, b) for b in itr)
 
 @inline function bmask(::Type{T}, range::UnitRange{Int})::T where {T<:Integer}
     ((one(T) << (range.stop - range.start + 1)) - one(T)) << (range.start - 1)
@@ -318,6 +325,6 @@ function controller(::Type{T}, cbits::IntIterator, cvals::IntIterator) where T
     do_mask = bmask(T, cbits...)
     target =
         length(cvals) == 0 ? zero(T) :
-        mapreduce(xy -> (xy[2] == 1 ? one(T) << T(xy[1] - 1) : zero(T)), |, zip(cbits, cvals))
+        mapreduce(xy -> (xy[2] == 1 ? indicator(T, xy[1]) : zero(T)), |, zip(cbits, cvals))
     return b -> ismatch(b, do_mask, target)
 end

src/longlonguint.jl

@@ -0,0 +1,85 @@
+"""
+    LongLongUInt{C} <: Integer
+
+A `LongLongUInt{C}` is an integer with `C` `UInt` numbers to store the value.
+"""
+struct LongLongUInt{C} <: Integer
+    content::NTuple{C, UInt}
+    function LongLongUInt{C}(content::NTuple{C, UInt}) where {C}
+        new{C}(content)
+    end
+    function LongLongUInt{C}(content::NTuple{C}) where {C}
+        new{C}(UInt.(content))
+    end
+    function LongLongUInt(content::NTuple{C}) where {C}
+        LongLongUInt{C}(content)
+    end
+end
+bsizeof(::LongLongUInt{C}) where C = bsizeof(UInt64) * C
+nint(::LongLongUInt{C}) where {C} = C
+Base.Int(x::LongLongUInt{1}) = Int(x.content[1])
+Base.UInt(x::LongLongUInt{1}) = x.content[1]
+Base.zero(::Type{LongLongUInt{C}}) where {C} = LongLongUInt{C}(ntuple(_->UInt(0), C))
+Base.zero(::LongLongUInt{C}) where {C} = zero(LongLongUInt{C})
+# convert from integers
+LongLongUInt{C}(x::T) where {C, T<:Integer} = LongLongUInt{C}(ntuple(i->i==1 ? UInt(x) : zero(UInt), C))
+Base.promote_type(::Type{LongLongUInt{C}}, ::Type{Int}) where {C} = LongLongUInt{C}
+Base.promote_type(::Type{LongLongUInt{C}}, ::Type{UInt}) where {C} = LongLongUInt{C}
+function Base.mod(x::LongLongUInt{C}, D::Int) where {C}
+    D == 2 ? mod(x.content[end], 2) : error("mod only supports 2")
+end
+
+function Base.:(>>)(x::LongLongUInt{C}, y::Int) where {C}
+    y < 0 && return x << -y
+    nshift = y ÷ bsizeof(UInt)
+    nshift_inner = y % bsizeof(UInt)
+    LongLongUInt(ntuple(C) do k
+            right = k<=nshift ? zero(UInt) : x.content[k-nshift]
+            left = k<=(nshift+1) ? zero(UInt) : x.content[k-nshift-1]
+            (right >> nshift_inner) | (left << (bsizeof(UInt) - nshift_inner))
+        end
+    )
+end
+function Base.:(<<)(x::LongLongUInt{C}, y::Int) where C
+    y < 0 && return x >> -y
+    nshift = y ÷ bsizeof(UInt)
+    nshift_inner = y % bsizeof(UInt)
+    LongLongUInt(ntuple(C) do k
+            left = k>C-nshift ? zero(UInt) : x.content[k+nshift]
+            right = k>C-nshift-1 ? zero(UInt) : x.content[k+nshift+1]
+            (left << nshift_inner) | (right >> (bsizeof(UInt) - nshift_inner))
+        end
+    )
+end
+function indicator(::Type{LongLongUInt{C}}, i::Int) where C
+    k = (i-1) ÷ bsizeof(UInt)
+    LongLongUInt{C}(ntuple(j->j==C-k ? indicator(UInt, i-k*bsizeof(UInt)) : zero(UInt), C))
+end
+for OP in [:&, :|, :xor]
+    @eval Base.$OP(x::LongLongUInt{C}, y::LongLongUInt{C}) where {C} = LongLongUInt{C}($OP.(x.content, y.content))
+end
+Base.:(~)(x::LongLongUInt{C}) where {C} = LongLongUInt{C}((~).(x.content))
+function Base.:(+)(x::LongLongUInt{C}, y::LongLongUInt{C}) where {C}
+    return LongLongUInt(_sadd(x.content, y.content, false))
+end
+function _sadd(x::NTuple{1,UInt}, y::NTuple{1,UInt}, c::Bool)
+    return (x[1] + y[1] + c,)
+end
+function _sadd(x::NTuple{C,UInt}, y::NTuple{C,UInt}, c::Bool) where {C}
+    v1, c1 = Base.add_with_overflow(x[C], y[C])
+    if c
+        v2, c2 = Base.add_with_overflow(v1, c)
+        c = c1 || c2
+        return (_sadd(x[1:C-1], y[1:C-1], c)..., v2)
+    else
+        return (_sadd(x[1:C-1], y[1:C-1], c1)..., v1)
+    end
+end
+Base.count_ones(x::LongLongUInt) = sum(count_ones, x.content)
+
+function longinttype(n::Int, D::Int)
+    N = ceil(Int, n * log2(D))
+    C = (N-1) ÷ bsizeof(UInt) + 1
+    return LongLongUInt{C}
+end
+Base.hash(x::LongLongUInt) = hash(x.content)

test/BitStr.jl

@@ -6,7 +6,8 @@ using Test, BitBasis
     @test bit"10_100_11" == bit"1010011"
     @test bit"10_100_11" != bit"01010011"
     @test bit"1011" === bit_literal(Int64.((1, 1, 0, 1))...)
-    @test BitBasis._parse_dit(Val(2), BigInt, "10101010101010101010101010101010010101010101"^10) isa LongBitStr
+    @test BitBasis._parse_dit(Val(2), BigInt, "10101010101010101010101010101010010101010101"^10) isa BitStr
+    @test BitBasis._parse_dit(Val(2), LongLongUInt, "10101010101010101010101010101010010101010101"^10) isa LongBitStr
     @test_throws ErrorException BitBasis._parse_dit(Val(2),
         Int64,
         "10101010101010101010101010101010010101010101"^10,
@@ -44,7 +45,7 @@ end
         @test res == 2
         @test typeof(res) == T
 
-        @test convert(BitStr{2,BigInt}, convert(T, bit"10")) == lbit"10"
+        @test convert(BitStr{2,LongLongUInt{1}}, convert(T, bit"10")) == lbit"10"
         @test convert(BitStr{2,Int64}, convert(T, bit"10")) === bit"10"
     end
 end
@@ -175,3 +176,10 @@ end
     @test bfloat_r(x) === 14 / 2^6
     @test btruncate(x, 3) === bit"000110"
 end
+
+@testset "LongDitStr" begin
+    x = LongBitStr(rand(0:1, 64 * 3))
+    @test bsizeof(x.buf) == 64 * 3
+    x = LongBitStr(rand(0:1, 64 * 3+1))
+    @test bsizeof(x.buf) == 64 * 4
+end

test/DitStr.jl

@@ -24,4 +24,4 @@ using BitBasis, Test
     @test_throws ErrorException BitBasis.parse_dit(Int64, "12341111111111111111111111111111111111111111111111111111111;5")
 
     @test hash(x) isa UInt64
-end
+end

test/longlonguint.jl

@@ -0,0 +1,65 @@
+using Test, BitBasis
+
+@testset "longlonguint" begin
+    x = LongLongUInt((3,))
+    @test Int(x) === 3
+    @test bsizeof(x) == 64
+    @test BitBasis.nint(x) == 1
+
+    x = LongLongUInt((3, 6))
+    @test zero(x) == LongLongUInt((0, 0))
+    @test x << 1 == LongLongUInt((6, 12))
+    @test x >> 1 == LongLongUInt((UInt(1), UInt(3) + UInt(1)<<63))
+    @test count_ones(x) == 4
+    @test ~x == LongLongUInt((~UInt(3), ~UInt(6)))
+
+    y = LongLongUInt((5, 7))
+    @test x & y == LongLongUInt((1, 6))
+    @test x | y == LongLongUInt((7, 7))
+    @test x ⊻ y == LongLongUInt((6, 1))
+    @test x + y == LongLongUInt((8, 13))
+
+    # add with overflow
+    z = LongLongUInt((UInt(17), typemax(UInt)-1))
+    @test z + x == LongLongUInt((21, 4))
+    @test (@allocated z + x) == 0
+
+    # maximum number of elements
+    BitBasis.max_num_elements(LongLongUInt{2}, 2) == 80
+    BitBasis.max_num_elements(LongLongUInt{2}, 4) == 40
+    BitBasis.max_num_elements(UInt, 2) == 64
+    BitBasis.max_num_elements(Int, 2) == 63
+end
+
+@testset "shift" begin
+    x = LongLongUInt((3, 6))
+    @test x << 0 == LongLongUInt((3, 6))
+    @test x << 64 == LongLongUInt((6, 0))
+    @test x << 65 == LongLongUInt((6<<1, 0))
+    @test x >> 64 == LongLongUInt((0, 3))
+    @test x >> 65 == LongLongUInt((0, 1))
+    @test x << -65 == LongLongUInt((0, 1))
+    @test x >> -65 == LongLongUInt((6<<1, 0))
+
+    x = LongLongUInt((3, 6, 7))
+    @test x >> 66 == LongLongUInt((UInt(0), UInt(0), UInt(3) << 62 + UInt(1)))
+    @test x << 63 == LongLongUInt((UInt(1) << 63 + UInt(6) >> 1, UInt(7) >> 1, UInt(1)<<63))
+end
+
+@testset "longinttype" begin
+    @test BitBasis.longinttype(1, 2) == LongLongUInt{1}
+    @test BitBasis.longinttype(64, 2) == LongLongUInt{1}
+    @test BitBasis.longinttype(65, 2) == LongLongUInt{2}
+end
+
+@testset "hash" begin
+    @test hash(LongLongUInt((3, 6))) == hash(LongLongUInt((3, 6)))
+    @test hash(LongLongUInt((3, 6))) != hash(LongLongUInt((6, 3)))
+end
+
+@testset "indicator" begin
+    @test indicator(LongLongUInt{1}, 1) == LongLongUInt((1,))
+    @test indicator(LongLongUInt{1}, 64) == LongLongUInt((UInt(1)<<63,))
+    @test indicator(LongLongUInt{2}, 64) == LongLongUInt((zero(UInt), UInt(1)<<63))
+    @test indicator(LongLongUInt{2}, 65) == LongLongUInt((1, 0))
+end

test/runtests.jl

@@ -30,5 +30,9 @@ end
     include("iterate_control.jl")
 end
 
+@testset "longlonguint" begin
+    include("longlonguint.jl")
+end
+
 DocMeta.setdocmeta!(BitBasis, :DocTestSetup, :(using BitBasis); recursive=true)
 Documenter.doctest(BitBasis; manual=false)