text:
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t α β γ δ
@variables x(t) y(t)
D = Differential(t)

eqs = [D(x) ~ α*x - β*x*y,
       D(y) ~ -δ*y + γ*x*y]

sys = ODESystem(eqs)

u0 = [x => 1.0,
      y => 1.0]

p  = [α => 1.5,
      β => 1.0,
      δ => 3.0,
      γ => 1.0]

tspan = (0.0,10.0)
prob = ODEProblem(sys,u0,tspan,p)
sol = solve(prob,Tsit5())

sys2 = liouville_transform(sys)
@variables trJ

u0 = [x => 1.0,
      y => 1.0,
      trJ => 1.0]

prob = ODEProblem(sys2,u0,tspan,p,jac=true)
sol = solve(prob,Tsit5())
@test sol[end,end] ≈ 1.0742818931017244

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
module ImageTools
#
# graythresh - Based on the original Paper: N. Otsu, "A Threshold Selection
# Method from Gray-Level Histograms" 1979. Calculates a threshold of a grayscale
# image, which can be used downstream to convert a grayscale image to binary
# image. One input argument is required - image (expected image format obtained
# by Images.imread, see Images package  for more information). E.g. to convert
# an array to image use grayim(array) or colorim(array).
#
# im2bw - converts image into a binary image. Input: image (expected image
# format obtained by Images.imread,  see Images package  for more information)
# and threshold (floating point number). E.g. to convert an array to image use
# grayim(array) or colorim(array).
#
#
# Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
#
#  MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

export Morph
export graythresh,
im2bw,
imshow,
imread,
bwlabel,
jet,
hsv,
label2rgb

include("imagetools/morph.jl")

import Tk
import Images
import Images: properties, data, raw, label_components
import FixedPointNumbers
import ImageView
import Colors
import BinDeps

function graythresh(img)

  #Check whether the input image is empty
  if isempty(img)
    error("Image is empty");
  end
  #Check whether the input image is gray

  if Images.colorspace(img) != "Gray"
    error("Input Image should be grayscale")
  end

  if isa(raw(img),Array{UInt16,2})
    #Convert image to 8bit and return it's raw values to compute the histogram
    image_array = map(Images.BitShift(Uint8,8),raw(img));
  elseif isa(raw(img),Array{Uint8,2})
    image_array = raw(img);
  else
    warn("Input Image is neither Uint8 or Uint16");
    image_array = raw(img);
  end
  #image_array = raw(img)

  expected_number_of_bins = 2^8;

  (range, counts) = hist(reshape(image_array,length(image_array)),expected_number_of_bins);

  #get real number of bins to be used below
  number_of_bins = size(counts)[1];
  #Next the algorithm implementation according to the original paper
  p = counts / sum(counts);
  ω = cumsum(p);
  μ = cumsum(p .* (1:number_of_bins));
  μ_t = μ[end];

  sigma_b_squared = (μ_t * ω - μ).^2 ./ (ω .* (1 - ω));

  # Find the location of the maximum value of sigma_b_squared.
  # The maximum may extend over several bins, so average together the
  # locations.  If maxval is NaN, meaning that sigma_b_squared is all NaN,
  # then return 0.
  maxvalue = maximum(sigma_b_squared);
  if isfinite(maxvalue)
    index = findfirst(sigma_b_squared, maxvalue);

    # Normalize the threshold to the range [0, 1].
    threshold = (index - 1) / (expected_number_of_bins - 1);
  else
    threshold = 0.0;
  end

  return threshold
end

function im2bw(img, threshold)

if Images.colorspace(img) != "Gray"
    error("Input Image should be grayscale")
  end

  if typeof(threshold) != Float64 && typeof(threshold) != Float32 && typeof(threshold) != Float16
    error("Threshold must be of types Float64, Float32 or Float16");
  end
  if threshold > 1 || threshold < 0
    error("Threshold must be between 0 and 1");
  end

  m,n = size(img);
  blacknWhite = falses(n,m)
  imageData = data(img);

  for i = 1:n
    for j = 1:m
      if imageData[j,i] > threshold
        blacknWhite[i,j] = true;
      end
    end
  end

  return Images.Image(blacknWhite, colorspace = "Binary")
end

function imshow(image)
  # wrapper for ImageView.view, need to add convertion of image to array and back
  ImageView.view(image)
end

function imread(path)
  # wrapper for Images.imread with added image retrieval from url functionality
  if (ismatch(r"http://.*", path) || ismatch(r"https://.*", path) || ismatch(r"ftp://.*", path)|| ismatch(r"smb://.*", path))
    imageContainer = nothing
    if VERSION < v"0.4-"
      imageContainer = download(path);
    else
      # By default, do silent download.
      url = path
      path = tempname()
      cmd = BinDeps.download_cmd(url, path)
      run(pipe(cmd, DevNull, DevNull));
      imageContainer = path
    end
  else
    imageContainer = path;
  end
  image = Images.imread(imageContainer);
  return image
end

function bwlabel(inputImage,connectivity)
  # wrapper for Images.label_components

  #define pixel connectivity matrices for appropriate inputs
  if connectivity == 4
    connectivityMatrix = bool([0 1 0;1 1 1; 0 1 0])
  elseif connectivity == 8
    connectivityMatrix = bool([1 1 1;1 1 1; 1 1 1]);
  else
    error("Invalid pixel connectivity. Pixel connectivity can be either 4 or 8");
  end

  labeledImage = label_components(inputImage,connectivityMatrix);

  return labeledImage;
end

function jet(numberOfColors::Int64)
  #Generate matlab-like jet colormap with specified amount of colors
  return jetColorMap = Colors.RGB{Float64}[
    Colors.RGB(
      clamp(min(4*ellement - 1.5, -4*ellement + 4.5) ,0.0,1.0),
      clamp(min(4*ellement - 0.5, -4*ellement + 3.5) ,0.0,1.0),
      clamp(min(4*ellement + 0.5, -4*ellement + 2.5) ,0.0,1.0))
    for ellement in linspace(0.0,1.0,numberOfColors)];

end

function hsv(numberOfColors::Int64)
  #Generate matlab-like hsv colormap with specified amount of colors
  return hsvColorMap = linspace(Colors.HSV(0,1,1),Colors.HSV(330,1,1),numberOfColors);

end

function label2rgb(labeledMatrix, inputColorMap = "jet",backgroundColor = [1 1 1], isShuffled = "noshuffle")

  #check if the labeledMatrix contains only integers
  if !isinteger(labeledMatrix)
    error("labeledMatrix is not an Integer matrix");
  end

  #change indexing of the labeled Matrix to start from 1
  labeledMatrix = labeledMatrix +1;
  numberOfEllements =  maximum(labeledMatrix); # number of connected ellements in the Matrix

  if inputColorMap == "jet"
    currentColormap =  jet(numberOfEllements);
  elseif inputColorMap == "hsv"
    currentColormap =  hsv(numberOfEllements);
    #check if custom colormap is loaded
  elseif is(typeof(inputColorMap),Array{Colors.RGB,1})

    if length(inputColorMap) == maximum(labeledMatrix)
      currentColormap = inputColorMap;
    else
      error("Input Colormap size is inconsistent with the Labeled Image");
    end
  else
    error("Invalid Input Colormap");
  end

  #convert input backgroundColor to RGB type
  backgroundColor = Colors.RGB(backgroundColor[1],backgroundColor[2],backgroundColor[3]);

  if  isShuffled == "shuffle"
    #shuffle the colormap
    currentColormap  = shuffle(currentColormap);
  elseif isShuffled == "noshuffle"
    #do nothing
  else
    error("isShuffled can be either shuffle or noshuffle ");
  end

  if (currentColormap[1] != backgroundColor)
    currentColormap = [backgroundColor,currentColormap[1:end-1]];
  end

  indexedImage = Images.ImageCmap(labeledMatrix, currentColormap);
  outputRGB = convert(Images.Image, indexedImage);

  return outputRGB;

end

end #End of ImageTools

ext: jl


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
steps(x) = x < -3/2 ? -1 : (x < 3/2 ? 0 : 1)
x = Float64[-4, -1, 2, -3, 0, 3, -2, 1, 4]
Xraw = (x = x, )
yraw = steps.(x);
idxsort = sortperm(x)
xsort = x[idxsort]
ysort = yraw[idxsort]
step(xsort, ysort, label="truth", where="mid")
plot(x, yraw, ls="none", marker="o", label="data")
xlim(-4.5, 4.5)
legend()

savefig(joinpath(@OUTPUT, "s1.svg")) # hide


text:
Base.@deprecate contiguous_view ContiguousView
Base.@deprecate strided_view StridedView
Base.@deprecate unsafe_view unsafe_aview


text:
module hpipm_julia
    using PyCall
    @pyimport numpy as np
    @pyimport hpipm_python as hp

    export hpipm_data, hpipm_dims, hpipm_solver
    include("py2ju_obj.jl")
    include("ju2py_obj.jl")
    include("hpipm_data.jl")
    include("hpipm_dims.jl")
    include("hpipm_solver.jl")
end


text:
export  Parser,
        next_token, peek_token, backup_token,
        peek_token_type,
        expect

"Only evaluate and return `default` expression if `nullable` is null"
macro get(nullable, default)
    quote
        if $(esc(nullable)) == nothing
            $(esc(default))
        else
            $(esc(nullable))
        end
    end
end

"Keeps track of current state of parsing."
mutable struct Parser
    lexer::Lexer
    ahead_token::Union{Token, Nothing}
    backup_token::Union{Token, Nothing}
    function Parser(lexer::Lexer)
        null = nothing
        new(lexer, null, null)
    end
end

function show(io::IO, p::Parser)
    print(io, "Parser($(p.lexer)")

    if p.ahead_token != nothing
        print(", $(get(p.ahead_token))")
    end

    if p.backup_token != nothing
        print(", $(get(p.backup_token))")
    end
    print(")")
end

"Move to next token from lexer"
function next_token(p::Parser)
    t = @get(p.ahead_token, next_token(p.lexer))
    p.backup_token = t
    p.ahead_token  = nothing
    t
end

function peek_token(p::Parser)
    t = @get(p.ahead_token, next_token(p.lexer))
    p.ahead_token = t
    t
end

peek_token_type(p::Parser) = token_type(peek_token(p))

function backup_token(p::Parser)
    p.ahead_token = p.backup_token
end

"""
    expect(parser, token_type)

Verify that next token is of type `token_type`, and if so, get next token.
"""
function expect(p::Parser, kind::TokenType)
    if  token_type(peek_token(p)) == kind
        next_token(p)
    else
        error("Expected $(Token(kind)) but got $(peek_token_type(p))")
    end
end


text:
import Base.Operators.(*)

import Base: scale!, scale, norm, vecdot

import Base: A_mul_B!, At_mul_B,  Ac_mul_B,  A_mul_Bc,  At_mul_Bt,  Ac_mul_Bc,  At_mul_Bt,
                       At_mul_B!, Ac_mul_B!, A_mul_Bc!, At_mul_Bt!, Ac_mul_Bc!, At_mul_Bt!

cublas_size(t::Char, M::CudaVecOrMat) = (size(M, t=='N' ? 1:2), size(M, t=='N' ? 2:1))

###########
#
# BLAS 1
#
###########

#######
# SCAL
#######
scale!{T<:CublasFloat}(x::CudaArray{T}, k::Number) = CUBLAS.scal!(length(x), k, x, 1)
scale{T<:CublasFloat}(x::CudaArray{T}, k::Number) = CUBLAS.scal!(length(x), k, copy(x), 1)

#######
# DOT
#######
function dot{T <: CublasFloat, TI<:Integer}(x::CudaVector{T}, rx::Union{UnitRange{TI},Range{TI}}, y::CudaVector{T}, ry::Union{UnitRange{TI},Range{TI}})
    if length(rx) != length(ry)
        throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))"))
    end
    if minimum(rx) < 1 || maximum(rx) > length(x)
        throw(BoundsError(x, rx))
    end
    if minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError(y, ry))
    end
    dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end

At_mul_B{T<:CublasReal}(x::CudaVector{T}, y::CudaVector{T}) = [CUBLAS.dot(x, y)]
At_mul_B{T<:CublasComplex}(x::CudaVector{T}, y::CudaVector{T}) = [CUBLAS.dotu(x, y)]
Ac_mul_B{T<:CublasComplex}(x::CudaVector{T}, y::CudaVector{T}) = [CUBLAS.dotc(x, y)]

vecdot{T<:CublasReal}(x::CudaVector{T}, y::CudaVector{T}) = dot(x, y)
vecdot{T<:CublasComplex}(x::CudaVector{T}, y::CudaVector{T}) = dotc(x, y)

#######
# NRM2
#######
norm(x::CudaArray) = nrm2(x)

############
#
# BLAS 2
#
############

#########
# GEMV
##########
function gemv_wrapper!{T<:CublasFloat}(y::CudaVector{T}, tA::Char, A::CudaMatrix{T}, x::CudaVector{T},
                                       alpha = one(T), beta = zero(T))
    mA, nA = cublas_size(tA, A)
    if nA != length(x)
        throw(DimensionMismatch("second dimension of A, $nA, does not match length of x, $(length(x))"))
    end
    if mA != length(y)
        throw(DimensionMismatch("first dimension of A, $mA, does not match length of y, $(length(y))"))
    end
    if mA == 0
        return y
    end
    if nA == 0
        return scale!(y, 0)
    end
    gemv!(tA, alpha, A, x, beta, y)
end

A_mul_B!{T<:CublasFloat}(y::CudaVector{T}, A::CudaMatrix{T}, x::CudaVector{T}) = gemv_wrapper!(y, 'N', A,  x)
At_mul_B!{T<:CublasFloat}(y::CudaVector{T}, A::CudaMatrix{T}, x::CudaVector{T}) = gemv_wrapper!(y, 'T', A, x)
Ac_mul_B!{T<:CublasFloat}(y::CudaVector{T}, A::CudaMatrix{T}, x::CudaVector{T}) = gemv_wrapper!(y, 'T', A, x)
Ac_mul_B!{T<:CublasComplex}(y::CudaVector{T}, A::CudaMatrix{T}, x::CudaVector{T}) = gemv_wrapper!(y, 'C', A, x)

function (*){T<:CublasFloat}(A::CudaMatrix{T}, x::CudaVector{T})
    A_mul_B!(similar(x, T, size(A,1)), A, x)
end

function At_mul_B{T<:CublasFloat}(A::CudaMatrix{T}, x::CudaVector{T})
    At_mul_B!(similar(x, T, size(A,2)), A, x)
end

function Ac_mul_B{T<:CublasFloat}(A::CudaMatrix{T}, x::CudaVector{T})
    Ac_mul_B!(similar(x, T, size(A,2)), A, x)
end

############
#
# BLAS 3
#
############

########
# GEMM
########
function gemm_wrapper!{T <: CublasFloat}(C::CudaVecOrMat{T}, tA::Char, tB::Char,
                                     A::CudaVecOrMat{T},
                                     B::CudaVecOrMat{T},
                                     alpha = one(T),
                                     beta = zero(T))
    mA, nA = cublas_size(tA, A)
    mB, nB = cublas_size(tB, B)

    if nA != mB
        throw(DimensionMismatch("A has dimensions ($mA,$nA) but B has dimensions ($mB,$nB)"))
    end

    if C === A || B === C
        throw(ArgumentError("output matrix must not be aliased with input matrix"))
    end

    if mA == 0 || nA == 0 || nB == 0
        if size(C) != (mA, nB)
            throw(DimensionMismatch("C has dimensions $(size(C)), should have ($mA,$nB)"))
        end
        return scale!(C, 0)
    end

    gemm!(tA, tB, alpha, A, B, beta, C)
end

# Mutating
A_mul_B!{T <: CublasFloat}(C::CudaMatrix{T}, A::CudaMatrix{T}, B::CudaMatrix{T}) = gemm_wrapper!(C, 'N', 'N', A, B)
At_mul_B!(C::CudaMatrix, A::CudaMatrix, B::CudaMatrix) = gemm_wrapper!(C, 'T', 'N', A, B)
At_mul_Bt!(C::CudaMatrix, A::CudaMatrix, B::CudaMatrix) = gemm_wrapper!(C, 'T', 'T', A, B)
Ac_mul_B!{T<:CublasReal}(C::CudaMatrix{T}, A::CudaMatrix{T}, B::CudaMatrix{T}) = At_mul_B!(C, A, B)
Ac_mul_B!(C::CudaMatrix, A::CudaMatrix, B::CudaMatrix) = gemm_wrapper!(C, 'C', 'N', A, B)

function A_mul_B!{T}(C::CudaMatrix{T}, A::CudaVecOrMat{T}, B::CudaVecOrMat{T})
    gemm_wrapper!(C, 'N', 'N', A, B)
end

# Non mutating

# A_mul_Bx
function (*){T <: CublasFloat}(A::CudaMatrix{T}, B::CudaMatrix{T})
    A_mul_B!(similar(B, T,(size(A,1), size(B,2))), A, B)
end

function A_mul_Bt{T}(A::CudaMatrix{T}, B::CudaMatrix{T})
    A_mul_Bt!(similar(B, T, (size(A,1), size(B,1))), A, B)
end

function A_mul_Bc{T}(A::CudaMatrix{T}, B::CudaMatrix{T})
    A_mul_Bc!(similar(B, T,(size(A,1),size(B,1))),A, B)
end

# At_mul_Bx
function At_mul_B{T}(A::CudaMatrix{T}, B::CudaMatrix{T})
    At_mul_B!(similar(B, T, (size(A,2), size(B,2))), A, B)
end

function At_mul_Bt{T}(A::CudaMatrix{T}, B::CudaMatrix{T})
    At_mul_Bt!(similar(B, T, (size(A,2), size(B,1))), A, B)
end

# Ac_mul_Bx
function Ac_mul_B{T}(A::CudaMatrix{T}, B::CudaMatrix{T})
    Ac_mul_B!(similar(B, T, (size(A,2), size(B,2))), A, B)
end

function Ac_mul_Bt{T,S}(A::CudaMatrix{T}, B::CudaMatrix{S})
    Ac_mul_Bt(similar(B, T, (size(A,2), size(B,1))), A, B)
end

function Ac_mul_Bc{T,S}(A::CudaMatrix{T}, B::CudaMatrix{S})
    Ac_mul_Bc!(similar(B, T, (size(A,2), size(B,1))), A, B)
end


text:
import InteractiveUtils
InteractiveUtils.versioninfo(stdout; verbose=true)

include("testsetup.jl")

@time @testset "VectorizationBase.jl" begin
    # Write your own tests here.
    # Aqua.test_all(VectorizationBase, ambiguities = VERSION < v"1.6-DEV")
    println("Aqua.test_all")
    @time Aqua.test_all(VectorizationBase)
    # @test isempty(detect_unbound_args(VectorizationBase))
    # @test isempty(detect_ambiguities(VectorizationBase))

    W = Int(@inferred(VectorizationBase.pick_vector_width(Float64)))
    @test @inferred(VectorizationBase.pick_integer(Val(W))) == (VectorizationBase.register_size() == VectorizationBase.simd_integer_register_size() ? Int64 : Int32)

    @test first(A) === A[1]
    @test W64S == W64
    println("Struct-Wrapped Vec")
    @time @testset "Struct-Wrapped Vec" begin
        @test data(zero(Vec{4,Float64})) === (VE(0.0),VE(0.0),VE(0.0),VE(0.0)) === data(Vec{4,Float64}(0.0))
        @test data(one(Vec{4,Float64})) === (VE(1.0),VE(1.0),VE(1.0),VE(1.0)) === data(Vec{4,Float64}(1.0)) === data(data(Vec{4,Float64}(1.0)))
        v = Vec((VE(1.0),VE(2.0),VE(3.0),VE(4.0)))
        @test v === Vec{4,Float64}(1, 2, 3, 4) === conj(v) === v' === Vec{4,Float64}(v)
        @test length(v) == 4 == first(size(v))
        @test eltype(v) == Float64
        for i in 1:4
            @test i == VectorizationBase.extractelement(v, i-1)
            # @test i === Vec{4,Int}(v)[i] # should use fptosi (ie, vconvert defined in SIMDPirates).
        end
        @test zero(v) === zero(typeof(v))
        @test one(v) === one(typeof(v))
        # @test Vec{W32,Float32}(one(Vec{W32,Float64})) === Vec(one(Vec{W32,Float32})) === one(Vec{W32,Float32}) # conversions should be tested in SIMDPirates
        @test Vec{1,Int}(1) === 1

        vu = Vec(collect(1.0:16.0)...) + 2
        @test_throws ErrorException vu.data

        @test vu(1,1) === VectorizationBase.data(vu)[1](1)
        @test vu(2,1) === VectorizationBase.data(vu)[1](2)
        @test vu(1,2) === VectorizationBase.data(vu)[2](1)
        @test vu(2,2) === VectorizationBase.data(vu)[2](2)
        if W64 == 8
            @test VectorizationBase.data(vu)[1] === Vec(3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0)
            @test VectorizationBase.data(vu)[2] === Vec(11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0)
        elseif W64 == 4
            @test VectorizationBase.data(vu)[1] === Vec(3.0, 4.0, 5.0, 6.0)
            @test VectorizationBase.data(vu)[2] === Vec(7.0, 8.0, 9.0, 10.0)
            @test VectorizationBase.data(vu)[3] === Vec(11.0, 12.0, 13.0, 14.0)
            @test VectorizationBase.data(vu)[4] === Vec(15.0, 16.0, 17.0, 18.0)
            @test Vec(1.2, 3.4, 3.4) === Vec(1.2, 3.4, 3.4, 0.0)
        elseif W64 == 2
            @test VectorizationBase.data(vu)[1] === Vec(3.0, 4.0)
            @test VectorizationBase.data(vu)[2] === Vec(5.0, 6.0)
            @test VectorizationBase.data(vu)[3] === Vec(7.0, 8.0)
            @test VectorizationBase.data(vu)[4] === Vec(9.0, 10.0)
            @test VectorizationBase.data(vu)[5] === Vec(11.0, 12.0)
            @test VectorizationBase.data(vu)[6] === Vec(13.0, 14.0)
            @test VectorizationBase.data(vu)[7] === Vec(15.0, 16.0)
            @test VectorizationBase.data(vu)[8] === Vec(17.0, 18.0)
        end

    end

    println("alignment.jl")
    @time @testset "alignment.jl" begin

        @test all(i -> VectorizationBase.align(i) == VectorizationBase.register_size(), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i) == 2VectorizationBase.register_size(), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i) == 10VectorizationBase.register_size(), (1:VectorizationBase.register_size()) .+ 9VectorizationBase.register_size())

        @test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid},   Int(VectorizationBase.register_size())), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid},  2Int(VectorizationBase.register_size())), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid}, 20Int(VectorizationBase.register_size())), (1:VectorizationBase.register_size()) .+ 19VectorizationBase.register_size())

        @test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), (1:VectorizationBase.register_size()) .+ 29VectorizationBase.register_size())

        @test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), (1:VectorizationBase.register_size()) .+ 29VectorizationBase.register_size())

        @test reinterpret(Int, VectorizationBase.align(pointer(A))) % VectorizationBase.register_size() === 0

        @test all(i -> VectorizationBase.aligntrunc(i) == 0, 0:VectorizationBase.register_size()-1)
        @test all(i -> VectorizationBase.aligntrunc(i) == VectorizationBase.register_size(), VectorizationBase.register_size():2VectorizationBase.register_size()-1)
        @test all(i -> VectorizationBase.aligntrunc(i) == 9VectorizationBase.register_size(), (0:VectorizationBase.register_size()-1) .+ 9VectorizationBase.register_size())

        @test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), (1:VectorizationBase.register_size()) .+ 29VectorizationBase.register_size())

        @test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), 1:VectorizationBase.register_size())
        @test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), 1+VectorizationBase.register_size():2VectorizationBase.register_size())
        @test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), (1:VectorizationBase.register_size()) .+ 29VectorizationBase.register_size())

        a = Vector{Float64}(undef, 0)
        ptr = pointer(a)
        @test UInt(VectorizationBase.align(ptr, 1 << 12)) % (1 << 12) == 0
    end

    println("masks.jl")
    @time @testset "masks.jl" begin
        # @test Mask{8,UInt8}(0x0f) === @inferred Mask(0x0f)
        # @test Mask{16,UInt16}(0x0f0f) === @inferred Mask(0x0f0f)
        @test EVLMask{8,UInt8}(0xff,8) === mask(Val(8), 0)
        @test EVLMask{8,UInt8}(0xff,8) === mask(Val(8), 8)
        @test EVLMask{8,UInt8}(0xff,8) === mask(Val(8), 16)
        @test EVLMask{8,UInt8}(0xff,8) === mask(Val(8), VectorizationBase.StaticInt(0))
        @test EVLMask{16,UInt16}(0xffff,16) === mask(Val(16), 0)
        @test EVLMask{16,UInt16}(0xffff,16) === mask(Val(16), 16)
        @test EVLMask{16,UInt16}(0xffff,16) === mask(Val(16), 32)
        @test all(w -> VectorizationBase.mask_type(w) === UInt8, 1:8)
        @test all(w -> VectorizationBase.mask_type(w) === UInt16, 9:16)
        @test all(w -> VectorizationBase.mask_type(w) === UInt32, 17:32)
        @test all(w -> VectorizationBase.mask_type(w) === UInt64, 33:64)
        @test all(w -> VectorizationBase.mask_type(w) === UInt128, 65:128)
        if VectorizationBase.register_size() == 64 # avx512
            # @test VectorizationBase.mask_type(Float16) === UInt32
            @test VectorizationBase.mask_type(Float32) === UInt16
            @test VectorizationBase.mask_type(Float64) === UInt8
            # @test VectorizationBase.max_mask(Float16) === 0xffffffff # 32
            @test data(VectorizationBase.max_mask(Float32)) === 0xffff     # 16
            @test data(VectorizationBase.max_mask(Float64)) === 0xff       # 8
        elseif VectorizationBase.register_size() == 32 # avx or avx2
            # @test VectorizationBase.mask_type(Float16) === UInt16
            @test VectorizationBase.mask_type(Float32) === UInt8
            @test VectorizationBase.mask_type(Float64) === UInt8
            # @test VectorizationBase.max_mask(Float16) === 0xffff     # 16
            @test data(VectorizationBase.max_mask(Float32)) === 0xff       # 8
            @test data(VectorizationBase.max_mask(Float64)) === 0x0f       # 4
        elseif VectorizationBase.register_size() == 16 # sse
            # @test VectorizationBase.mask_type(Float16) === UInt8
            @test VectorizationBase.mask_type(Float32) === UInt8
            @test VectorizationBase.mask_type(Float64) === UInt8
            # @test VectorizationBase.max_mask(Float16) === 0xff       # 8
            @test data(VectorizationBase.max_mask(Float32)) === 0x0f       # 4
            @test data(VectorizationBase.max_mask(Float64)) === 0x03       # 2
        end
        @test all(w -> bitstring(VectorizationBase.mask(Val( 8), w)) == reduce(*, ( 8 - i < w ? "1" : "0" for i in 1:8 )), 1:8 )
        @test all(w -> bitstring(VectorizationBase.mask(Val(16), w)) == reduce(*, (16 - i < w ? "1" : "0" for i in 1:16)), 1:16)
        @test all(w -> VectorizationBase.mask(Float64, w) === VectorizationBase.mask(@inferred(VectorizationBase.pick_vector_width(Float64)), w), 1:W64)

        @test VectorizationBase.vbroadcast(Val(8), true) === Vec(true, true, true, true, true, true, true, true)

        @test !VectorizationBase.vall(Mask{8}(0xfc))
        @test !VectorizationBase.vall(Mask{4}(0xfc))
        @test VectorizationBase.vall(EVLMask{8}(0xff,8))
        @test !VectorizationBase.vall(EVLMask{8}(0x1f,5))
        @test VectorizationBase.vall(Mask{4}(0xcf))

        @test VectorizationBase.vany(Mask{8}(0xfc))
        @test VectorizationBase.vany(Mask{4}(0xfc))
        @test !VectorizationBase.vany(Mask{8}(0x00))
        @test !VectorizationBase.vany(Mask{4}(0xf0))

        @test VectorizationBase.vall(Mask{8}(0xfc) + Mask{8}(0xcf) == Vec(0x01,0x01,0x02,0x02,0x01,0x01,0x02,0x02))
        @test VectorizationBase.vall(Mask{4}(0xfc) + Mask{4}(0xcf) == Vec(0x01,0x01,0x02,0x02))
        @test VectorizationBase.vall(Mask{8}(0xcf) + EVLMask{8}(0x1f,5) == Vec(0x02, 0x02, 0x02, 0x02, 0x01, 0x00, 0x01, 0x01))

        @test VectorizationBase.vall(Mask{8}(0xec) != Mask{8}(0x13))
        @test VectorizationBase.vall((!Mask{8}(0xac) & Mask{8}(0xac)) == Mask{8}(0x00))
        @test !VectorizationBase.vany((!Mask{8}(0xac) & Mask{8}(0xac)))
        @test VectorizationBase.vall((!Mask{8}(0xac) | Mask{8}(0xac)) == Mask{8}(0xff))
        @test VectorizationBase.vall((!Mask{8}(0xac) | Mask{8}(0xac)))

        @test VectorizationBase.vall(VectorizationBase.splitint(0xb53a5d6426a9d29d, Int8) == Vec{8,Int8}(-99, -46, -87, 38, 100, 93, 58, -75))
        # other splitint tests for completeness sake
        @test VectorizationBase.splitint(0xb53a5d6426a9d29d, Int64) === 0xb53a5d6426a9d29d
        @test VectorizationBase.splitint(0xff, UInt16) === 0x00ff
        @test !VectorizationBase.vany(
            VectorizationBase.splitint(0x47766b9a9509d175acd77ff497236795, Int8) != Vec{16,Int8}(-107, 103, 35, -105, -12, 127, -41, -84, 117, -47, 9, -107, -102, 107, 118, 71)
        )

        @test (EVLMask{8}(0x1f,5) | EVLMask{8}(0x03,3)) === EVLMask{8}(0x1f,5)
        @test (Mask{8}(0x1f) | EVLMask{8}(0x03,3)) === Mask{8}(0x1f)
        @test (EVLMask{8}(0x1f,5) | Mask{8}(0x03)) === Mask{8}(0x1f)
        @test (Mask{8}(0x1f) | Mask{8}(0x03)) === Mask{8}(0x1f)

        @test (EVLMask{8}(0x1f,5) & EVLMask{8}(0x03,3)) === EVLMask{8}(0x03,3)
        @test (Mask{8}(0x1f) & EVLMask{8}(0x03,3)) === Mask{8}(0x03)
        @test (EVLMask{8}(0x1f,5) & Mask{8}(0x03)) === Mask{8}(0x03)
        @test (Mask{8}(0x1f) & Mask{8}(0x03)) === Mask{8}(0x03)

        @test (Mask{8}(0xac) | false) === Mask{8}(0xac)
        @test (Mask{8}(0xac) | true) === Mask{8}(0xff)
        @test (false | Mask{8}(0xac)) === Mask{8}(0xac)
        @test (true | Mask{8}(0xac)) === Mask{8}(0xff)
        @test (Mask{8}(0xac) & false) === Mask{8}(0x00)
        @test (Mask{8}(0xac) & true) === Mask{8}(0xac)
        @test (false & Mask{8}(0xac)) === Mask{8}(0x00)
        @test (true & Mask{8}(0xac)) === Mask{8}(0xac)
        @test (Mask{8}(0xac) ⊻ false) === Mask{8}(0xac)
        @test (Mask{8}(0xac) ⊻ true) === Mask{8}(0x53)
        @test (false ⊻ Mask{8}(0xac)) === Mask{8}(0xac)
        @test (true ⊻ Mask{8}(0xac)) === Mask{8}(0x53)

        @test (Mask{4}(0x05) | true) === Mask{4}(0x0f)
        @test (Mask{4}(0x05) | false) === Mask{4}(0x05)
        @test (true | Mask{4}(0x05)) === Mask{4}(0x0f)
        @test (false | Mask{4}(0x05)) === Mask{4}(0x05)
        for T ∈ [UInt8, UInt16, UInt32]
            Ws = T === UInt8 ? [2,4,8] : [8sizeof(T)]
            for W ∈ Ws
                u1 = rand(T)
                u2 = rand(T)
                m = ~(typemax(T) << W)
                @test (Mask{W}(u1) & Mask{W}(u2)) === Mask{W}( (u1 & u2) & m )
                @test (Mask{W}(u1) | Mask{W}(u2)) === Mask{W}( (u1 | u2) & m )
                @test (Mask{W}(u1) ⊻ Mask{W}(u2)) === Mask{W}( (u1 ⊻ u2) & m )
            end
        end
        @test convert(Bool, Mask{8}(0xec)) === Vec(false,false,true,true,false,true,true,true) === VectorizationBase.ifelse(convert(Bool, Mask{8}(0xec)), vbroadcast(Val(8),true), vbroadcast(Val(8),false))

        @test (MM{8}(2) ∈ 3:8) === Mask{8}(0x7e)

        fbitvector1 = falses(20);
        fbitvector2 = falses(20);
        mu = VectorizationBase.VecUnroll((Mask{4}(0x0f),Mask{4}(0x0f)))
        GC.@preserve fbitvector begin
            vstore!(stridedpointer(fbitvector1), mu, (VectorizationBase.MM(StaticInt{8}(), 1),))
            vstore!(stridedpointer(fbitvector2), mu, (VectorizationBase.MM(StaticInt{8}(), 1),), Mask{8}(0x7e))
        end
        @test all(fbitvector1[1:8])
        @test !any(fbitvector1[9:end])
        @test !fbitvector2[1]
        @test all(fbitvector2[2:7])
        @test !any(fbitvector2[8:end])
    end

    # @testset "number_vectors.jl" begin
    # # eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :(size(A)), 8)) # doesn't work?
    # @test VectorizationBase.length_loads(A, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() -> 13*17)()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, 13*17, 8)) == divrem(length(A), 8)
    # @test VectorizationBase.size_loads(A,1, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() -> 13   )()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, 13   , 8)) == divrem(size(A,1), 8)
    # @test VectorizationBase.size_loads(A,2, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() ->    17)()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__,    17, 8)) == divrem(size(A,2), 8)
    # end

    println("vector_width.jl")
    @time @testset "vector_width.jl" begin
        for T ∈ (Float32,Float64)
            @test @inferred(VectorizationBase.pick_vector_width(T)) * @inferred(VectorizationBase.static_sizeof(T)) === @inferred(VectorizationBase.register_size(T)) === @inferred(VectorizationBase.register_size())
        end
        for T ∈ (Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64)
            @test @inferred(VectorizationBase.pick_vector_width(T)) * @inferred(VectorizationBase.static_sizeof(T)) === @inferred(VectorizationBase.register_size(T)) === @inferred(VectorizationBase.simd_integer_register_size())
        end

        @test @inferred(VectorizationBase.pick_vector_width(Float64, Int32, Float64, Float32, Float64)) * VectorizationBase.static_sizeof(Float64) === @inferred(VectorizationBase.register_size())
        @test @inferred(VectorizationBase.pick_vector_width(Float64, Int32)) * VectorizationBase.static_sizeof(Float64) === @inferred(VectorizationBase.register_size())

        @test @inferred(VectorizationBase.pick_vector_width(Float32, Float32)) * VectorizationBase.static_sizeof(Float32) === @inferred(VectorizationBase.register_size())
        @test @inferred(VectorizationBase.pick_vector_width(Float32, Int32)) * VectorizationBase.static_sizeof(Float32) === @inferred(VectorizationBase.simd_integer_register_size())

        @test all(VectorizationBase._ispow2, 0:1)
        @test all(i -> !any(VectorizationBase._ispow2, 1+(1 << (i-1)):(1 << i)-1 ) && VectorizationBase._ispow2(1 << i), 2:9)
        @test all(i ->  VectorizationBase.intlog2(1 << i) == i, 0:(Int == Int64 ? 53 : 30))
        FTypes = (Float32, Float64)
        Wv = ntuple(i -> @inferred(VectorizationBase.register_size()) >> (i+1), Val(2))
        for (T, N) in zip(FTypes, Wv)
            W = @inferred(VectorizationBase.pick_vector_width(T))
            # @test Vec{Int(W),T} == VectorizationBase.pick_vector(W, T) == VectorizationBase.pick_vector(T)
            @test W == @inferred(VectorizationBase.pick_vector_width(W, T))
            @test W === @inferred(VectorizationBase.pick_vector_width(W, T)) == @inferred(VectorizationBase.pick_vector_width(T))
            while true
                W >>= VectorizationBase.One()
                W == 0 && break
                W2, Wshift2 = @inferred(VectorizationBase.pick_vector_width_shift(W, T))
                @test W2 == VectorizationBase.One() << Wshift2 == @inferred(VectorizationBase.pick_vector_width(W, T)) == VectorizationBase.pick_vector_width(Val(Int(W)),T)  == W
                @test StaticInt(W) === VectorizationBase.pick_vector_width(Val(Int(W)), T) === VectorizationBase.pick_vector_width(W, T)
                for n in W+1:2W
                    W3, Wshift3 = VectorizationBase.pick_vector_width_shift(StaticInt(n), T)
                    @test W2 << 1 == W3 == 1 << (Wshift2+1) == 1 << Wshift3 == VectorizationBase.pick_vector_width(StaticInt(n), T) == VectorizationBase.pick_vector_width(Val(n),T) == W << 1
                    # @test VectorizationBase.pick_vector(W, T) == VectorizationBase.pick_vector(W, T) == Vec{Int(W),T}
                end
            end
        end

        # @test VectorizationBase.nextpow2(0) == 1
        @test all(i -> VectorizationBase.nextpow2(i) == i, 0:2)
        for j in 1:10
            l, u = (1<<j)+1, 1<<(j+1)
            @test all(i -> VectorizationBase.nextpow2(i) == u, l:u)
        end

    end

    println("Memory")
    @time @testset "Memory" begin
        C = rand(40,20,10) .> 0;
        mtest = vload(stridedpointer(C), ((MM{16})(9), 2, 3));
        v1 = C[9:24,2,3];
        @test tovector(mtest) == v1
        @test [vload(stridedpointer(C), (1+w, 2+w, 3)) for w ∈ 1:W64] == getindex.(Ref(C), 1 .+ (1:W64), 2 .+ (1:W64), 3)
        vstore!(stridedpointer(C), !mtest, ((MM{16})(17), 3, 4))
        @test .!v1 == C[17:32,3,4] == tovector(vload(stridedpointer(C), ((MM{16})(17), 3, 4)))

        dims = (41,42,43) .* 3;
        # dims = (41,42,43);
        A = reshape(collect(Float64(0):Float64(prod(dims)-1)), dims);

        P = PermutedDimsArray(A, (3,1,2));
        O = OffsetArray(P, (-4, -2, -3));
        indices = (
            StaticInt{1}(), StaticInt{2}(), 2, MM{W64}(2), MM{W64,2}(3), MM{W64,-1}(W64+2), Vec(ntuple(i -> 2i + 1, Val(W64))...)#,
            # VectorizationBase.LazyMulAdd{2,-1}(MM{W64}(3))#, VectorizationBase.LazyMulAdd{2,-2}(Vec(ntuple(i -> 2i + 1, Val(W64))...))
        )
        println("LazyMulAdd Loads/Stores")
        @time @testset "LazyMulAdd Loads/Stores" begin
            max_const = 2
            for _i ∈ indices, _j ∈ indices, _k ∈ indices, im ∈ 1:max_const, jm ∈ 1:max_const, km ∈ 1:max_const, B ∈ (A, P, O)
                i = @inferred(VectorizationBase.lazymul(StaticInt(im), _i))
                j = @inferred(VectorizationBase.lazymul(StaticInt(jm), _j))
                k = @inferred(VectorizationBase.lazymul(StaticInt(km), _k))
                iv = tovector(i); jv = tovector(j); kv = tovector(k)
                if B === C
                    off = 9 - iv[1] % 8
                    iv += off
                    i += off
                end
                # @show typeof(B), i, j, k (im, _i), (jm, _j), (km, _k)
                x = getindex.(Ref(B), iv, jv, kv)
                GC.@preserve B begin
                    # @show i,j,k, typeof(B)
                    v = @inferred(vload(stridedpointer(B), (i, j, k)))
                end
                @test x == tovector(v)
                if length(x) > 1
                    m = Mask{W64}(rand(UInt8))
                    mv = tovector(m)
                    x .*= mv
                    GC.@preserve B begin
                        v = @inferred(vload(stridedpointer(B), (i, j, k), m))
                    end
                    @test x == tovector(v)
                end
                for store! ∈ (vstore!, VectorizationBase.vnoaliasstore!)
                    y = isone(length(x)) ? randn() : randnvec(length(x))
                    GC.@preserve B store!(stridedpointer(B), y, (i, j, k))
                    x = getindex.(Ref(B), iv, jv, kv)
                    # @show i, j, k typeof.((i, j, k)), store!, typeof(B) y
                    @test x == tovector(y)
                    if length(x) > 1
                        z = Vec(ntuple(_ -> Core.VecElement(randn()), length(x)))
                        GC.@preserve B store!(stridedpointer(B), z, (i, j, k), m)
                        y = getindex.(Ref(B), iv, jv, kv)
                        @test y == ifelse.(mv, tovector(z), x)
                    end
                end
            end
        end
        println("VecUnroll Loads/Stores")
        @time @testset "VecUnroll Loads/Stores" begin
            for AU ∈ 1:3, B ∈ (A, P, O), i ∈ (StaticInt(1),2,StaticInt(2)), j ∈ (StaticInt(1),3,StaticInt(3)), k ∈ (StaticInt(1),4,StaticInt(4))
                for AV ∈ 1:3
                    v1 = randnvec(); v2 = randnvec(); v3 = randnvec();
                    GC.@preserve B begin
                        if AU == AV
                            vstore!(stridedpointer(B), VectorizationBase.VecUnroll((v1,v2,v3)), VectorizationBase.Unroll{AU,W64,3,AV,W64,zero(UInt)}((i, j, k)))
                            vu = @inferred(vload(stridedpointer(B), VectorizationBase.Unroll{AU,W64,3,AV,W64,zero(UInt)}((i, j, k))))
                        else
                            vstore!(stridedpointer(B), VectorizationBase.VecUnroll((v1,v2,v3)), VectorizationBase.Unroll{AU,1,3,AV,W64,zero(UInt)}((i, j, k)))
                            vu = @inferred(vload(stridedpointer(B), VectorizationBase.Unroll{AU,1,3,AV,W64,zero(UInt)}((i, j, k))))
                        end
                    end
                    @test v1 === VectorizationBase.data(vu)[1]
                    @test v2 === VectorizationBase.data(vu)[2]
                    @test v3 === VectorizationBase.data(vu)[3]

                    ir = 0:(AV == 1 ? W64-1 : 0); jr = 0:(AV == 2 ? W64-1 : 0); kr = 0:(AV == 3 ? W64-1 : 0)
                    x1 = getindex.(Ref(B), i .+ ir, j .+ jr, k .+ kr)
                    if AU == 1
                        ir = ir .+ length(ir)
                    elseif AU == 2
                        jr = jr .+ length(jr)
                    elseif AU == 3
                        kr = kr .+ length(kr)
                    end
                    x2 = getindex.(Ref(B), i .+ ir, j .+ jr, k .+ kr)
                    if AU == 1
                        ir = ir .+ length(ir)
                    elseif AU == 2
                        jr = jr .+ length(jr)
                    elseif AU == 3
                        kr = kr .+ length(kr)
                    end
                    x3 = getindex.(Ref(B), i .+ ir, j .+ jr, k .+ kr)

                    @test x1 == tovector(VectorizationBase.data(vu)[1])
                    @test x2 == tovector(VectorizationBase.data(vu)[2])
                    @test x3 == tovector(VectorizationBase.data(vu)[3])

                end
                v1 = randnvec(); v2 = randnvec(); v3 = randnvec(); v4 = randnvec(); v5 = randnvec()
                GC.@preserve B begin
                    vstore!(VectorizationBase.vsum, stridedpointer(B), VectorizationBase.VecUnroll((v1,v2,v3,v4,v5)), VectorizationBase.Unroll{AU,1,5,0,1,zero(UInt)}((i, j, k)))
                end
                ir = 0:(AU == 1 ? 4 : 0); jr = 0:(AU == 2 ? 4 : 0); kr = 0:(AU == 3 ? 4 : 0)
                xvs = getindex.(Ref(B), i .+ ir, j .+ jr, k .+ kr)
                @test xvs ≈ map(VectorizationBase.vsum, [v1,v2,v3,v4,v5])
            end
        end
        x = Vector{Int}(undef, 100);
        i = MM{1}(0)
        for j ∈ 1:25
            VectorizationBase.__vstore!(pointer(x), j, (i * VectorizationBase.static_sizeof(Int)), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.register_size())
            i += 1
        end
        for j ∈ 26:50
            VectorizationBase.__vstore!(pointer(x), j, (VectorizationBase.static_sizeof(Int) * i), Mask{1}(0xff), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.register_size())
            i += 1
        end
        for j ∈ 51:75
            VectorizationBase.__vstore!(pointer(x), j, VectorizationBase.lazymul(i, VectorizationBase.static_sizeof(Int)), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.register_size())
            i += 1
        end
        for j ∈ 76:100
            VectorizationBase.__vstore!(pointer(x), j, VectorizationBase.lazymul(VectorizationBase.static_sizeof(Int), i), Mask{1}(0xff), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.register_size())
            i += 1
        end
        @test x == 1:100

        xf64 = rand(100);
        GC.@preserve xf64 begin
            vxtu = @inferred(vload(stridedpointer(xf64), (MM{4W64}(1),)));
            @test vxtu isa VectorizationBase.VecUnroll
            vxtutv = tovector(vxtu);
            vxtutvmult = 3.5 .* vxtutv;
            @inferred(vstore!(stridedpointer(xf64), 3.5 * vxtu, (MM{4W64}(1),)));
            @test tovector(@inferred(vload(stridedpointer(xf64), (MM{4W64}(1),)))) == vxtutvmult
            mbig = Mask{4W64}(rand(UInt32)); # TODO: update if any arches support >512 bit vectors
            mbigtv = tovector(mbig);
            @test tovector(@inferred(vload(stridedpointer(xf64), (MM{4W64}(1),), mbig))) == ifelse.(mbigtv, vxtutvmult, 0.0)
            @inferred(vstore!(stridedpointer(xf64), -11 * vxtu, (MM{4W64}(1),), mbig));
            @test tovector(@inferred(vload(stridedpointer(xf64), (MM{4W64}(1),)))) == ifelse.(mbigtv, -11 .* vxtutv, vxtutvmult)
        end
        colors = [(R = rand(), G = rand(), B = rand()) for i ∈ 1:100];
        colormat = reinterpret(reshape, Float64, colors)
        sp = stridedpointer(colormat)
        GC.@preserve colors begin
            @test tovector(@inferred(vload(sp, VectorizationBase.Unroll{1,1,3,2,8,zero(UInt)}((1,9))))) == vec(colormat[:,9:16]')
            vu = @inferred(vload(sp, VectorizationBase.Unroll{1,1,3,2,8,zero(UInt)}((1,41))))
            @inferred(vstore!(sp, vu, VectorizationBase.Unroll{1,1,3,2,8,zero(UInt)}((1,1))))
        end
        @test vec(colormat[:,41:48]) == vec(colormat[:,1:8])
    end

    println("Grouped Strided Pointers")
    @time @testset "Grouped Strided Pointers" begin
        M, K, N = 4, 5, 6
        A = rand(M, K); B = rand(K, N); C = rand(M, N);
        GC.@preserve A B C begin
            fs = [identity, adjoint]
            for fA ∈ fs, fB ∈ fs, fC ∈ fs
                At = fA === identity ? A : copy(A')'
                Bt = fB === identity ? B : copy(B')'
                Ct = fC === identity ? C : copy(C')'
                gsp, pres = @inferred(VectorizationBase.grouped_strided_pointer((At,Bt,Ct), Val{(((1,1),(3,1)),((1,2),(2,1)),((2,2),(3,2)))}()))
                if fA === fC
                    @test sizeof(gsp.strides) == 2sizeof(Int)
                end
                ai = fA === identity
                bi = fB === identity
                ci = fC === identity
                # Test to confirm that redundant strides are not stored in the grouped strided pointer
                @test sizeof(gsp) == sizeof(Int) * (6 - (ai & ci) - ((!ai) & bi) - ((!bi) & (!ci)))
                @test sizeof(gsp.offsets) == 0
                pA, pB, pC = @inferred(VectorizationBase.stridedpointers(gsp))
                @test pA === stridedpointer(At)
                @test pB === stridedpointer(Bt)
                @test pC === stridedpointer(Ct)
            end
        end
    end

    println("Adjoint VecUnroll")
    @time @testset "Adjoint VecUnroll" begin
        W = W64
        while W > 1
            A = rand(W,W); B = similar(A);
            GC.@preserve A B begin
                vut = @inferred(vload(stridedpointer(A), VectorizationBase.Unroll{2,1,W,1,W}((1,1))))
                vu = @inferred(VectorizationBase.transpose_vecunroll(vut))
                @test vu === @inferred(vload(stridedpointer(A'), VectorizationBase.Unroll{2,1,W,1,W}((1,1))))
                @test vu === @inferred(vload(stridedpointer(A), VectorizationBase.Unroll{1,1,W,2,W}((1,1))))
                vstore!(stridedpointer(B), vu, VectorizationBase.Unroll{2,1,W,1,W}((1,1)))
            end
            @test A == B'
            W >>= 1
        end
        W = 2W64
        while W > 1
            A = rand(Float32,W,W); B = similar(A);
            GC.@preserve A B begin
                vut = @inferred(vload(stridedpointer(A), VectorizationBase.Unroll{2,1,W,1,W}((1,1))))
                vu = @inferred(VectorizationBase.transpose_vecunroll(vut))
                @test vu === @inferred(vload(stridedpointer(A'), VectorizationBase.Unroll{2,1,W,1,W}((1,1))))
                @test vu === @inferred(vload(stridedpointer(A), VectorizationBase.Unroll{1,1,W,2,W}((1,1))))
                vstore!(stridedpointer(B), vu, VectorizationBase.Unroll{2,1,W,1,W}((1,1)))
            end
            @test A == B'
            W >>= 1
        end
    end

    println("Unary Functions")
    @time @testset "Unary Functions" begin
        for T ∈ (Float32,Float64)
            for f ∈ [floatmin,floatmax,typemin,typemax]
                @test f(Vec{Int(pick_vector_width(T)),T}) === Vec(ntuple(_ -> f(T), pick_vector_width(T))...)
            end
            v = let W = VectorizationBase.pick_vector_width(T)
                VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> (randn(T)), W)...),
                    Vec(ntuple(_ -> (randn(T)), W)...),
                    Vec(ntuple(_ -> (randn(T)), W)...)
                ))
            end
            x = tovector(v)
            for f ∈ [
                -, abs, inv, floor, ceil, trunc, round, VectorizationBase.relu, abs2,
                Base.FastMath.abs2_fast, Base.FastMath.sub_fast, sign
            ]
                # @show T, f
                @test tovector(@inferred(f(v))) == map(f, x)
            end
            # Don't require exact, but `eps(T)` seems like a reasonable `rtol`, at least on AVX512 systems:
            # function relapprox(x::AbstractVector{T},y) where {T}
            #     t = max(norm(x),norm(y)) * eps(T)
            #     n = norm(x .- y)
            #     n / t
            # end
            # function randapprox(::Type{T}) where {T}
            #     x = Vec(ntuple(_ -> 10randn(T), VectorizationBase.pick_vector_width(T))...)
            #     via = @fastmath inv(x)
            #     vir = inv(x)
            #     relapprox(tovector(via), tovector(vir))
            # end
            # f32t = map(_ -> randapprox(Float32), 1:1_000_000);
            # f64t = map(_ -> randapprox(Float64), 1:1_000_000);
            # summarystats(f64t)
            # summarystats(f32t)
            # for now, I'll use `4eps(T)` if the systems don't have AVX512, but should check to set a stricter bound.
            # also put `sqrt ∘ abs` in here
            let rtol = eps(T) * (Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f))) ? 1 : 4) # more accuracte
                @test isapprox(tovector(@inferred(Base.FastMath.inv_fast(v))), map(Base.FastMath.inv_fast, x), rtol = rtol)
                let f = sqrt ∘ abs
                    if T === Float32
                        @test isapprox(tovector(@inferred(f(v))), map(f, x), rtol = rtol)
                    elseif T === Float64 # exact with `Float64`
                        @test tovector(@inferred(f(v))) == map(f, x)
                    end
                end
            end
            for f ∈ [floor, ceil, trunc, round]
                @test tovector(@inferred(f(Int32, v))) == map(y -> f(Int32,y), x)
                @test tovector(@inferred(f(Int64, v))) == map(y -> f(Int64,y), x)
            end
            invtol = Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f))) ? 2^-14 : 1.5*2^-12 # moreaccurate with AVX512
            @test isapprox(tovector(@inferred(VectorizationBase.inv_approx(v))), map(VectorizationBase.inv_approx, x), rtol = invtol)
        end

        int = Bool(VectorizationBase.has_feature(Val(:x86_64_avx512dq))) ? Int : Int32
        int2 = Bool(VectorizationBase.has_feature(Val(:x86_64_avx2))) ? Int : Int32
        vi = VectorizationBase.VecUnroll((
            Vec(ntuple(_ -> rand(int), Val(W64))...),
            Vec(ntuple(_ -> rand(int), Val(W64))...),
            Vec(ntuple(_ -> rand(int), Val(W64))...)
        )) % int2
        xi = tovector(vi)
        for f ∈ [-, abs, inv, floor, ceil, trunc, round, sqrt ∘ abs, sign]
            @test tovector(@inferred(f(vi))) == map(f, xi)
        end
        let rtol = eps(Float64) * (Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f))) ? 1 : 4) # more accuracte
            @test isapprox(tovector(@inferred(Base.FastMath.inv_fast(vi))), map(Base.FastMath.inv_fast, xi), rtol = rtol)
        end
        # vpos = VectorizationBase.VecUnroll((
        #     Vec(ntuple(_ -> Core.VecElement(rand()), Val(W64))),
        #     Vec(ntuple(_ -> Core.VecElement(rand()), Val(W64))),
        #     Vec(ntuple(_ -> Core.VecElement(rand()), Val(W64)))
        # ))
        # for f ∈ [sqrt]
        #     @test tovector(f(vpos)) == map(f, tovector(vpos))
        # end
    end
    println("Binary Functions")
    @time @testset "Binary Functions" begin
        # TODO: finish getting these tests to pass
        # for I1 ∈ (Int32,Int64,UInt32,UInt64), I2 ∈ (Int32,Int64,UInt32,UInt64)
        let WI = Int(VectorizationBase.pick_vector_width(Int64))
            for I1 ∈ (Int32,Int64), I2 ∈ (Int32,Int64,UInt32)
                # TODO: No longer skip these either.
                sizeof(I1) > sizeof(I2) && continue
                vi1 = VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> Core.VecElement(rand(I1)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(I1)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(I1)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(I1)), Val(WI)))
                ))
                srange = one(I2):(Bool(VectorizationBase.has_feature(Val(:x86_64_avx512dq))) ? I2(8sizeof(I1)-1) : I2(31))
                vi2 = VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> Core.VecElement(rand(srange)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(srange)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(srange)), Val(WI))),
                    Vec(ntuple(_ -> Core.VecElement(rand(srange)), Val(WI)))
                ))
                i = rand(srange); j = rand(I1);
                m1 = VectorizationBase.VecUnroll((MM{WI}(I1(7)), MM{WI}(I1(1)), MM{WI}(I1(13)), MM{WI}(I1(32%last(srange)))));
                m2 = VectorizationBase.VecUnroll((MM{WI,2}(I2(3)), MM{WI,2}(I2(8)), MM{WI,2}(I2(39%last(srange))), MM{WI,2}(I2(17))));
                @test typeof(m1 + I1(11)) === typeof(m1)
                @test typeof(m2 + I2(11)) === typeof(m2)
                xi1 = tovector(vi1); xi2 = tovector(vi2);
                xi3 =  mapreduce(tovector, vcat, VectorizationBase.data(m1));
                xi4 =  mapreduce(tovector, vcat, VectorizationBase.data(m2));
                I3 = promote_type(I1,I2);
                # I4 = sizeof(I1) < sizeof(I2) ? I1 : (sizeof(I1) > sizeof(I2) ? I2 : I3)
                for f ∈ [
                    +, -, *, ÷, /, %, <<, >>, >>>, ⊻, &, |, fld, mod,
                    VectorizationBase.rotate_left, VectorizationBase.rotate_right, copysign, maxi, mini, maxi_fast, mini_fast
                ]
                    # for f ∈ [+, -, *, div, ÷, /, rem, %, <<, >>, >>>, ⊻, &, |, fld, mod, VectorizationBase.rotate_left, VectorizationBase.rotate_right, copysign, max, min]
                    # @show f, I1, I2
                    # if (!Bool(VectorizationBase.has_feature(Val(:x86_64_avx512dq)))) && (f === /) && sizeof(I1) === sizeof(I2) === 8
                    #     continue
                    # end
                    check_within_limits(tovector(@inferred(f(vi1, vi2))),  trunc_int.(f.(size_trunc_int.(xi1, I3), size_trunc_int.(xi2, I3)), I3));
                    check_within_limits(tovector(@inferred(f(j, vi2))), trunc_int.(f.(size_trunc_int.(j, I3), size_trunc_int.(xi2, I3)), I3));
                    check_within_limits(tovector(@inferred(f(vi1, i))), trunc_int.(f.(size_trunc_int.(xi1, I3), size_trunc_int.(i, I3)), I3));
                    check_within_limits(tovector(@inferred(f(m1, i))), trunc_int.(f.(size_trunc_int.(xi3, I3), size_trunc_int.(i, I3)), I3));
                    check_within_limits(tovector(@inferred(f(m1, vi2))), trunc_int.(f.(size_trunc_int.(xi3, I3), size_trunc_int.(xi2, I3)), I3));
                    check_within_limits(tovector(@inferred(f(m1, m2))), trunc_int.(f.(size_trunc_int.(xi3, I3), size_trunc_int.(xi4, I3)), I3));
                    check_within_limits(tovector(@inferred(f(m1, m1))), trunc_int.(f.(size_trunc_int.(xi3, I1), size_trunc_int.(xi3, I1)), I1));
                    check_within_limits(tovector(@inferred(f(m2, i))), trunc_int.(f.(size_trunc_int.(xi4, I3), size_trunc_int.(i, I3)), I2));
                    check_within_limits(tovector(@inferred(f(m2, vi2))), trunc_int.(f.(size_trunc_int.(xi4, I3), size_trunc_int.(xi2, I3)), I2));
                    check_within_limits(tovector(@inferred(f(m2, m2))), trunc_int.(f.(size_trunc_int.(xi4, I3), size_trunc_int.(xi4, I3)), I2));
                    check_within_limits(tovector(@inferred(f(m2, m1))), trunc_int.(f.(size_trunc_int.(xi4, I3), size_trunc_int.(xi3, I3)), I3));
                    if !((f === VectorizationBase.rotate_left) || (f === VectorizationBase.rotate_right))
                        check_within_limits(tovector(@inferred(f(j, m1))), trunc_int.(f.(j, xi3), I1));
                        # @show 12
                        # check_within_limits(tovector(@inferred(f(j, m2))), trunc_int.(f.(size_trunc_int.(j, I1), size_trunc_int.(xi4, I1)), I1));
                    end
                end
                @test tovector(@inferred(vi1 ^ i)) ≈ xi1 .^ i
                @test @inferred(VectorizationBase.vall(@inferred(1 - MM{WI}(1)) == (1 - Vec(ntuple(identity, Val(WI))...)) ))
            end
            vf1 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> Core.VecElement(randn()), Val(WI))),
                Vec(ntuple(_ -> Core.VecElement(randn()), Val(WI)))
            ))
            vf2 = Vec(ntuple(_ -> Core.VecElement(randn()), Val(WI)))
            xf1 = tovector(vf1); xf2 = tovector(vf2); xf22 = vcat(xf2,xf2)
            a = randn();
            for f ∈ [+, -, *, /, %, max, min, copysign, rem, Base.FastMath.max_fast, Base.FastMath.min_fast, Base.FastMath.div_fast, Base.FastMath.rem_fast]
                # @show f
                @test tovector(@inferred(f(vf1, vf2))) ≈ f.(xf1, xf22)
                @test tovector(@inferred(f(a, vf1))) ≈ f.(a, xf1)
                @test tovector(@inferred(f(a, vf2))) ≈ f.(a, xf2)
                @test tovector(@inferred(f(vf1, a))) ≈ f.(xf1, a)
                @test tovector(@inferred(f(vf2, a))) ≈ f.(xf2, a)
            end

            vi2 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(WI))),
                Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(WI))),
                Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(WI))),
                Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(WI)))
            ))
            vones, vi2f, vtwos = promote(1.0, vi2, 2f0); # promotes a binary function, right? Even when used with three args?
            @test vones === VectorizationBase.VecUnroll((vbroadcast(Val(WI), 1.0),vbroadcast(Val(WI), 1.0),vbroadcast(Val(WI), 1.0),vbroadcast(Val(WI), 1.0)));
            @test vtwos === VectorizationBase.VecUnroll((vbroadcast(Val(WI), 2.0),vbroadcast(Val(WI), 2.0),vbroadcast(Val(WI), 2.0),vbroadcast(Val(WI), 2.0)));
            @test VectorizationBase.vall(VectorizationBase.collapse_and(vi2f == vi2))
            W32 = StaticInt(WI)*StaticInt(2)
            vf2 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> Core.VecElement(randn(Float32)), W32)),
                Vec(ntuple(_ -> Core.VecElement(randn(Float32)), W32))
            ))
            vones32, v2f32, vtwos32 = promote(1.0, vf2, 2f0); # promotes a binary function, right? Even when used with three args?
            if Bool(VectorizationBase.has_feature(Val(:x86_64_avx2)))
                @test vones32 === VectorizationBase.VecUnroll((vbroadcast(W32, 1f0),vbroadcast(W32, 1f0)))
                @test vtwos32 === VectorizationBase.VecUnroll((vbroadcast(W32, 2f0),vbroadcast(W32, 2f0)))
                @test vf2 === v2f32
            else
                @test vones32 === VectorizationBase.VecUnroll((vbroadcast(W32, 1.0),vbroadcast(W32, 1.0)))
                @test vtwos32 === VectorizationBase.VecUnroll((vbroadcast(W32, 2.0),vbroadcast(W32, 2.0)))
                @test convert(Float64, vf2) === v2f32
            end
            i = rand(1:31)
            m1 = VectorizationBase.VecUnroll((MM{WI}(7), MM{WI}(1), MM{WI}(13), MM{WI}(18)))
            @test tovector(clamp(m1, 2:i)) == clamp.(tovector(m1), 2, i)
            @test tovector(mod(m1, 1:i)) == mod1.(tovector(m1), i)

            @test VectorizationBase.vdivrem.(1:30, 1:30') == divrem.(1:30, 1:30')
            @test VectorizationBase.vcld.(1:30, 1:30') == cld.(1:30, 1:30')
            @test VectorizationBase.vrem.(1:30, 1:30') == rem.(1:30, 1:30')

            @test gcd(Vec(42,64,0,-37), Vec(18,96,-38,0)) === Vec(6,32,38,37)
            @test lcm(Vec(24,16,42,0),Vec(18,12,18,17)) === Vec(72, 48, 126, 0)
        end
    end
    println("Ternary Functions")
    @time @testset "Ternary Functions" begin
        for T ∈ (Float32, Float64)
            v1, v2, v3, m = let W = @inferred(VectorizationBase.pick_vector_width(T))
                v1 = VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> randn(T), W)...),
                    Vec(ntuple(_ -> randn(T), W)...)
                ))
                v2 = VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> randn(T), W)...),
                    Vec(ntuple(_ -> randn(T), W)...)
                ))
                v3 = VectorizationBase.VecUnroll((
                    Vec(ntuple(_ -> randn(T), W)...),
                    Vec(ntuple(_ -> randn(T), W)...)
                ))
                _W = Int(@inferred(VectorizationBase.pick_vector_width(T)))
                m = VectorizationBase.VecUnroll((Mask{_W}(rand(UInt16)),Mask{_W}(rand(UInt16))))
                v1, v2, v3, m
            end
            x1 = tovector(v1); x2 = tovector(v2); x3 = tovector(v3);
            a = randn(T); b = randn(T)
            a64 = Float64(a); b64 = Float64(b); # test promotion
            mv = tovector(m)
            for f ∈ [
                muladd, fma, clamp, VectorizationBase.vmuladd_fast, VectorizationBase.vfma_fast,
                VectorizationBase.vfmadd, VectorizationBase.vfnmadd, VectorizationBase.vfmsub, VectorizationBase.vfnmsub,
                VectorizationBase.vfmadd_fast, VectorizationBase.vfnmadd_fast, VectorizationBase.vfmsub_fast, VectorizationBase.vfnmsub_fast,
                VectorizationBase.vfmadd231, VectorizationBase.vfnmadd231, VectorizationBase.vfmsub231, VectorizationBase.vfnmsub231
            ]
                @test tovector(@inferred(f(v1, v2, v3))) ≈ map(f, x1, x2, x3)
                @test tovector(@inferred(f(v1, v2, a64))) ≈ f.(x1, x2, a)
                @test tovector(@inferred(f(v1, a64, v3))) ≈ f.(x1, a, x3)
                @test tovector(@inferred(f(a64, v2, v3))) ≈ f.(a, x2, x3)
                @test tovector(@inferred(f(v1, a64, b64))) ≈ f.(x1, a, b)
                @test tovector(@inferred(f(a64, v2, b64))) ≈ f.(a, x2, b)
                @test tovector(@inferred(f(a64, b64, v3))) ≈ f.(a, b, x3)

                @test tovector(@inferred(VectorizationBase.ifelse(f, m, v1, v2, v3))) ≈ ifelse.(mv, f.(x1, x2, x3), x3)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, v1, v2, a64))) ≈ ifelse.(mv, f.(x1, x2, a), a)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, v1, a64, v3))) ≈ ifelse.(mv, f.(x1, a, x3), x3)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, a64, v2, v3))) ≈ ifelse.(mv, f.(a, x2, x3), x3)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, v1, a64, b64))) ≈ ifelse.(mv, f.(x1, a, b), b)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, a64, v2, b64))) ≈ ifelse.(mv, f.(a, x2, b), b)
                @test tovector(@inferred(VectorizationBase.ifelse(f, m, a64, b64, v3))) ≈ ifelse.(mv, f.(a, b, x3), x3)
            end
        end
        let WI = Int(VectorizationBase.pick_vector_width(Int64))
            vi64 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> rand(Int64), Val(WI))...),
                Vec(ntuple(_ -> rand(Int64), Val(WI))...),
                Vec(ntuple(_ -> rand(Int64), Val(WI))...)
            ))
            vi32 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> rand(Int32), Val(WI))...),
                Vec(ntuple(_ -> rand(Int32), Val(WI))...),
                Vec(ntuple(_ -> rand(Int32), Val(WI))...)
            ))
            xi64 = tovector(vi64); xi32 = tovector(vi32);
            @test tovector(@inferred(VectorizationBase.ifelse(vi64 > vi32, vi64, vi32))) == ifelse.(xi64 .> xi32, xi64, xi32)
            @test tovector(@inferred(VectorizationBase.ifelse(vi64 < vi32, vi64, vi32))) == ifelse.(xi64 .< xi32, xi64, xi32)
            @test tovector(@inferred(VectorizationBase.ifelse(true, vi64, vi32))) == ifelse.(true, xi64, xi32)
            @test tovector(@inferred(VectorizationBase.ifelse(false, vi64, vi32))) == ifelse.(false, xi64, xi32)
            vu64_1 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...)
            ))
            vu64_2 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...)
            ))
            vu64_3 = VectorizationBase.VecUnroll((
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...),
                Vec(ntuple(_ -> rand(UInt64), Val(WI))...)
            ))
            xu1 = tovector(vu64_1); xu2 = tovector(vu64_2); xu3 = tovector(vu64_3);
            for f ∈ [clamp, muladd, VectorizationBase.ifmalo, VectorizationBase.ifmahi, VectorizationBase.vfmadd, VectorizationBase.vfnmadd, VectorizationBase.vfmsub, VectorizationBase.vfnmsub]
                @test tovector(@inferred(f(vu64_1,vu64_2,vu64_3))) == f.(xu1, xu2, xu3)
            end
        end
    end
    println("Special functions")
    @time @testset "Special functions" begin
        if VERSION ≥ v"1.6.0-DEV.674" && Bool(VectorizationBase.has_feature(Val(Symbol("x86_64_sse4.1"))))
            erfs = [0.1124629160182849, 0.22270258921047847, 0.3286267594591274, 0.42839235504666845, 0.5204998778130465, 0.6038560908479259, 0.6778011938374184, 0.7421009647076605, 0.7969082124228322, 0.8427007929497149, 0.8802050695740817, 0.9103139782296353, 0.9340079449406524, 0.9522851197626487, 0.9661051464753108, 0.976348383344644, 0.9837904585907745, 0.9890905016357308, 0.9927904292352575, 0.9953222650189527, 0.997020533343667, 0.9981371537020182, 0.9988568234026434, 0.999311486103355, 0.999593047982555, 0.9997639655834707, 0.9998656672600594, 0.9999249868053346, 0.9999589021219005, 0.9999779095030014, 0.9999883513426328, 0.9999939742388483]
            if Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f)))
                v = VectorizationBase.verf(Vec{8, Float64}(0.1:0.1:0.8...,))
                @test [v(i) for i in 1:8] ≈ erfs[1:8]
                v = VectorizationBase.verf(Vec{16, Float32}(0.1:0.1:1.6...,))
                @test [v(i) for i in 1:16] ≈ erfs[1:16]
            end
            if Bool(VectorizationBase.has_feature(Val(:x86_64_avx)))
                v = VectorizationBase.verf(Vec{4, Float64}(0.1:0.1:0.4...,))
                @test [v(i) for i in 1:4] ≈ erfs[1:4]
                v = VectorizationBase.verf(Vec{8, Float32}(0.1:0.1:0.8...,))
                @test [v(i) for i in 1:8] ≈ erfs[1:8]
            end
            if Bool(VectorizationBase.has_feature(Val(Symbol("x86_64_sse4.1"))))
                v = VectorizationBase.verf(Vec{2, Float64}(0.1:0.1:0.2...,))
                @test [v(i) for i in 1:2] ≈ erfs[1:2]
            end
        end
    end
    println("Non-broadcasting operations")
    @time @testset "Non-broadcasting operations" begin
        v1 = Vec(ntuple(_ -> Core.VecElement(randn()), Val(W64))); vu1 = VectorizationBase.VecUnroll((v1, Vec(ntuple(_ -> Core.VecElement(randn()), Val(W64)))));
        v2 = Vec(ntuple(_ -> Core.VecElement(rand(-100:100)), Val(W64))); vu2 = VectorizationBase.VecUnroll((v2, Vec(ntuple(_ -> Core.VecElement(rand(-100:100)), Val(W64)))));
        @test @inferred(VectorizationBase.vsum(2.3, v1)) ≈ @inferred(VectorizationBase.vsum(v1)) + 2.3 ≈ @inferred(VectorizationBase.vsum(VectorizationBase.addscalar(v1, 2.3))) ≈ @inferred(VectorizationBase.vsum(VectorizationBase.addscalar(2.3, v1)))
        @test @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(vu1))) + 2.3 ≈ @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(VectorizationBase.addscalar(vu1, 2.3)))) ≈ @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(VectorizationBase.addscalar(2.3, vu1))))
        @test @inferred(VectorizationBase.vsum(v2)) + 3 == @inferred(VectorizationBase.vsum(VectorizationBase.addscalar(v2, 3))) == @inferred(VectorizationBase.vsum(VectorizationBase.addscalar(3, v2)))
        @test @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(vu2))) + 3 == @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(VectorizationBase.addscalar(vu2, 3)))) == @inferred(VectorizationBase.vsum(VectorizationBase.collapse_add(VectorizationBase.addscalar(3, vu2))))
        @test @inferred(VectorizationBase.vprod(v1)) * 2.3 ≈ @inferred(VectorizationBase.vprod(VectorizationBase.mulscalar(v1, 2.3))) ≈ @inferred(VectorizationBase.vprod(VectorizationBase.mulscalar(2.3, v1)))
        @test @inferred(VectorizationBase.vprod(v2)) * 3 == @inferred(VectorizationBase.vprod(VectorizationBase.mulscalar(3, v2)))
        @test @inferred(VectorizationBase.vall(v1 + v2 == VectorizationBase.addscalar(v1, v2)))
        @test 4.0 == @inferred(VectorizationBase.addscalar(2.0, 2.0))

        v3 = Vec(0, 1, 2, 3); vu3 = VectorizationBase.VecUnroll((v3, v3 - 1))
        v4 = Vec(0.0, 1.0, 2.0, 3.0)
        v5 = Vec(0f0, 1f0, 2f0, 3f0, 4f0, 5f0, 6f0, 7f0)
        @test @inferred(VectorizationBase.vmaximum(v3)) === @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v3, 2)))
        @test @inferred(VectorizationBase.vmaximum(v3 % UInt)) === @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v3 % UInt, 2 % UInt)))
        @test @inferred(VectorizationBase.vmaximum(v4)) === @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v4, prevfloat(3.0))))
        @test @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v4, nextfloat(3.0)))) == nextfloat(3.0)
        @test @inferred(VectorizationBase.vmaximum(v5)) === @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v5, prevfloat(7f0)))) === VectorizationBase.vmaximum(VectorizationBase.maxscalar(prevfloat(7f0), v5))
        @test @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(v5, nextfloat(7f0)))) == @inferred(VectorizationBase.vmaximum(VectorizationBase.maxscalar(nextfloat(7f0), v5))) == nextfloat(7f0)

        @test VectorizationBase.maxscalar(v3, 2) === Vec(2, 1, 2, 3)
        @test VectorizationBase.maxscalar(v3, -1) === v3
        @test VectorizationBase.vmaximum(VectorizationBase.maxscalar(v3 % UInt, -1 % UInt)) === -1 % UInt
        @test VectorizationBase.maxscalar(v4, 1e-16) === Vec(1e-16, 1.0, 2.0, 3.0)
        @test VectorizationBase.maxscalar(v4, -1e-16) === v4
        @test VectorizationBase.vmaximum(VectorizationBase.collapse_max(vu3)) == 3
        @test VectorizationBase.vmaximum(VectorizationBase.collapse_max(VectorizationBase.maxscalar(vu3,2))) == 3
        @test VectorizationBase.vmaximum(VectorizationBase.collapse_max(VectorizationBase.maxscalar(vu3,4))) == 4
        @test VectorizationBase.vminimum(VectorizationBase.collapse_min(vu3)) == -1
        @test VectorizationBase.vminimum(VectorizationBase.collapse_min(VectorizationBase.minscalar(vu3,0))) == -1
        @test VectorizationBase.vminimum(VectorizationBase.collapse_min(VectorizationBase.minscalar(vu3,-2))) == VectorizationBase.vminimum(VectorizationBase.collapse_min(VectorizationBase.minscalar(-2,vu3))) == -2
    end
    println("broadcasting")
    @time @testset "broadcasting" begin
        @test VectorizationBase.vzero(Val(1), UInt32) === VectorizationBase.vzero(StaticInt(1), UInt32) === VectorizationBase.vzero(UInt32) === zero(UInt32)
        @test VectorizationBase.vzero(Val(1), Int) === VectorizationBase.vzero(StaticInt(1), Int) === VectorizationBase.vzero(Int) === 0
        @test VectorizationBase.vzero(Val(1), Float32) === VectorizationBase.vzero(StaticInt(1), Float32) === VectorizationBase.vzero(Float32) === 0f0
        @test VectorizationBase.vzero(Val(1), Float64) === VectorizationBase.vzero(StaticInt(1), Float64) === VectorizationBase.vzero(Float64) === 0.0
        @test VectorizationBase.vzero() === VectorizationBase.vzero(W64S, Float64)
        @test VectorizationBase.vbroadcast(StaticInt(2)*W64S, one(Int64)) === VectorizationBase.vbroadcast(StaticInt(2)*W64S, one(Int32))
        @test VectorizationBase.vbroadcast(StaticInt(2)*W64S, one(UInt64)) === VectorizationBase.vbroadcast(StaticInt(2)*W64S, one(UInt32)) === one(VectorizationBase.Vec{2W64,UInt64}) === oneunit(VectorizationBase.Vec{2W64,UInt64})

        @test VectorizationBase.vall(VectorizationBase.vbroadcast(W64S, pointer(A)) == vbroadcast(W64S, first(A)))
        @test VectorizationBase.vbroadcast(W64S, pointer(A,2)) === Vec{W64}(A[2]) === Vec(A[2])

        @test zero(VectorizationBase.VecUnroll((VectorizationBase.vbroadcast(W64S, pointer(A)), VectorizationBase.vbroadcast(W64S, pointer(A,2))))) === VectorizationBase.VecUnroll((VectorizationBase.vzero(W64S, Float64), VectorizationBase.vzero()))

        @test VectorizationBase.VecUnroll{2,W64,Float64}(first(A)) === VectorizationBase.VecUnroll{2,W64,Float64}(VectorizationBase.vbroadcast(W64S, pointer(A))) === VectorizationBase.VecUnroll((VectorizationBase.vbroadcast(W64S, pointer(A)),VectorizationBase.vbroadcast(W64S, pointer(A)),VectorizationBase.vbroadcast(W64S, pointer(A)))) === VectorizationBase.VecUnroll{2}(VectorizationBase.vbroadcast(W64S, pointer(A)))
    end
    println("CartesianVIndex")
    @time @testset "CartesianVIndex" begin
        @test VectorizationBase.maybestaticfirst(CartesianIndices(A)) === VectorizationBase.CartesianVIndex(ntuple(_ -> VectorizationBase.One(), ndims(A)))
        @test VectorizationBase.maybestaticlast(CartesianIndices(A)) === VectorizationBase.CartesianVIndex(size(A))
        @test VectorizationBase.CartesianVIndex((StaticInt(1),2,VectorizationBase.CartesianVIndex((StaticInt(4), StaticInt(7))), CartesianIndex(12,14), StaticInt(2), 1)) === VectorizationBase.CartesianVIndex((StaticInt(1),2,StaticInt(4),StaticInt(7),12,14,StaticInt(2),1))
        @test Tuple(VectorizationBase.CartesianVIndex((StaticInt(1),2,VectorizationBase.CartesianVIndex((StaticInt(4), StaticInt(7))), CartesianIndex(12,14), StaticInt(2), 1))) === (StaticInt(1),2,StaticInt(4),StaticInt(7),12,14,StaticInt(2),1)
        @test length(VectorizationBase.CartesianVIndex((StaticInt(1),2,VectorizationBase.CartesianVIndex((StaticInt(4), StaticInt(7))), CartesianIndex(12,14), StaticInt(2), 1))) === 8
        @test VectorizationBase.static_length(VectorizationBase.CartesianVIndex((StaticInt(1),2,VectorizationBase.CartesianVIndex((StaticInt(4), StaticInt(7))), CartesianIndex(12,14), StaticInt(2), 1))) === StaticInt{8}()
        @test VectorizationBase.CartesianVIndex((StaticInt(-4), StaticInt(7))):VectorizationBase.CartesianVIndex((StaticInt(14), StaticInt(73))) === CartesianIndices((StaticInt(-4):StaticInt(14),StaticInt(7):StaticInt(73)))
        @test VectorizationBase.maybestaticfirst(CartesianIndices(A)):VectorizationBase.maybestaticlast(CartesianIndices(A)) == CartesianIndices(A)
        @test VectorizationBase.maybestaticfirst(CartesianIndices(A)):VectorizationBase.maybestaticlast(CartesianIndices(A)) === CartesianIndices(map(i -> VectorizationBase.One():i, size(A)))
    end
    println("Promotion")
    @time @testset "Promotion" begin
        vi2 = VectorizationBase.VecUnroll((
            Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(W64))),
            Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(W64))),
            Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(W64))),
            Vec(ntuple(_ -> Core.VecElement(rand(1:M-1)), Val(W64)))
        ))
        vones, vi2f, vtwos = @inferred(promote(1.0, vi2, 2f0)); # promotes a binary function, right? Even when used with three args?
        @test vones === VectorizationBase.VecUnroll((vbroadcast(Val(W64), 1.0),vbroadcast(Val(W64), 1.0),vbroadcast(Val(W64), 1.0),vbroadcast(Val(W64), 1.0)));
        @test vtwos === VectorizationBase.VecUnroll((vbroadcast(Val(W64), 2.0),vbroadcast(Val(W64), 2.0),vbroadcast(Val(W64), 2.0),vbroadcast(Val(W64), 2.0)));
        @test @inferred(VectorizationBase.vall(VectorizationBase.collapse_and(vi2f == vi2)))
        vf2 = VectorizationBase.VecUnroll((
            Vec(ntuple(_ -> Core.VecElement(randn(Float32)), StaticInt(W32))),
            Vec(ntuple(_ -> Core.VecElement(randn(Float32)), StaticInt(W32)))
        ))
        vones32, v2f32, vtwos32 = @inferred(promote(1.0, vf2, 2f0)); # promotes a binary function, right? Even when used with three args?
        @test vones32 === VectorizationBase.VecUnroll((vbroadcast(StaticInt(W32), 1f0),vbroadcast(StaticInt(W32), 1f0)))
        @test vtwos32 === VectorizationBase.VecUnroll((vbroadcast(StaticInt(W32), 2f0),vbroadcast(StaticInt(W32), 2f0)))
        @test vf2 === v2f32

        vm = if Bool(VectorizationBase.has_feature(Val(:x86_64_avx512dq)))
            VectorizationBase.VecUnroll((
                MM{W64}(rand(Int)),MM{W64}(rand(Int)),MM{W64}(rand(Int)),MM{W64}(rand(Int))
            ))
        else
            VectorizationBase.VecUnroll((
                MM{W64}(rand(Int32)),MM{W64}(rand(Int32)),MM{W64}(rand(Int32)),MM{W64}(rand(Int32))
            ))
        end
        @test tovector(@inferred(vm > vi2)) == (tovector(vm) .> tovector(vi2))

        m = Mask{2W64}(rand(UInt));
        v64f = Vec(ntuple(_ -> randn(), Val{2W64}())...)
        v32i = Vec(ntuple(_ -> rand(Int32), Val{2W64}())...)
        mtv = tovector(m); v64ftv = tovector(v64f); v32itv = tovector(v32i);
        vum = @inferred(muladd(v64f, v32i, m))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(v64ftv, v32itv, mtv)

        vum = @inferred(muladd(v64f, m, v32i))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(v64ftv, mtv, v32itv)

        vum = @inferred(muladd(v32i, v64f, m))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(v32itv, v64ftv, mtv)

        vum = @inferred(muladd(v32i, m, v64f))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(v32itv, mtv, v64ftv)

        vum = @inferred(muladd(m, v64f, v32i))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(mtv, v64ftv, v32itv)

        vum = @inferred(muladd(m, v32i, v64f))
        @test vum isa VectorizationBase.VecUnroll
        @test tovector(vum) ≈ muladd.(mtv, v32itv, v64ftv)

        vx = convert(Vec{16,Int64}, 1)
        @test typeof(vx) === typeof(zero(vx)) === Vec{16,Int64}

        vxf32 = Vec(ntuple(_ -> randn(Float32), VectorizationBase.pick_vector_width(Float32))...)
        xf32, yf32 = promote(vxf32, 1.0)
        @test xf32 === vxf32
        @test yf32 === vbroadcast(VectorizationBase.pick_vector_width(Float32), 1f0)
        vxi32 = Vec(ntuple(_ -> rand(Int32), VectorizationBase.pick_vector_width(Int32))...)
        xi32, yi32 = promote(vxi32, one(Int64))
        @test xi32 === vxi32
        @test yi32 === vbroadcast(VectorizationBase.pick_vector_width(Int32), one(Int32))

    end
    println("Lazymul")
    @time @testset "Lazymul" begin
        # partially covered in memory
        for i ∈ (-5, -1, 0, 1, 4, 8), j ∈ (-5, -1, 0, 1, 4, 8)
            @test @inferred(VectorizationBase.lazymul(StaticInt(i), StaticInt(j))) === StaticInt(i*j)
        end
        fi = VectorizationBase.LazyMulAdd{8,0}(MM{8}(StaticInt(16)))
        si = VectorizationBase.LazyMulAdd{2}(240)
        @test @inferred(VectorizationBase.vadd_fast(fi, si)) === VectorizationBase.LazyMulAdd{2,128}(MM{8,4}(240))
    end
    # TODO: Put something here.
    # @time @testset "Arch Functions" begin
    #     @test VectorizationBase.dynamic_register_size() == @inferred(VectorizationBase.register_size()) == @inferred(VectorizationBase.register_size())
    #     @test VectorizationBase.dynamic_integer_register_size() == @inferred(VectorizationBase.simd_integer_register_size()) == @inferred(VectorizationBase.ssimd_integer_register_size())
    #     @test VectorizationBase.dynamic_register_count() == @inferred(VectorizationBase.register_count()) == @inferred(VectorizationBase.sregister_count())
    #     @test VectorizationBase.dynamic_fma_fast() == VectorizationBase.fma_fast()
    #     @test VectorizationBase.dynamic_has_opmask_registers() == VectorizationBase.has_opmask_registers()
    # end
    println("Static Zero and One")
    @time @testset "Static Zero and One" begin
        vx = randnvec(W64)
        vu = VectorizationBase.VecUnroll((vx,randnvec(W64)))
        vm = MM{16}(24);
        for f ∈ [+,Base.FastMath.add_fast]
            @test f(vx, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vx) === vx
            @test f(vu, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vu) === vu
            @test f(vm, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vm) === vm
        end
        for f ∈ [-,Base.FastMath.sub_fast]
            @test f(vx, VectorizationBase.Zero()) ===  vx
            @test f(VectorizationBase.Zero(), vx) === -vx
            @test f(vu, VectorizationBase.Zero()) ===  vu
            @test f(VectorizationBase.Zero(), vu) === -vu
            @test f(vm, VectorizationBase.Zero()) ===  vm
            @test f(VectorizationBase.Zero(), vm) === -vm
        end
        for f ∈ [*,Base.FastMath.mul_fast]
            @test f(vx, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vx) === VectorizationBase.Zero()
            @test f(vu, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vu) === VectorizationBase.Zero()
            @test f(vm, VectorizationBase.Zero()) === f(VectorizationBase.Zero(), vm) === VectorizationBase.Zero()
            @test f(vx, VectorizationBase.One()) === f(VectorizationBase.One(), vx) === vx
            @test f(vu, VectorizationBase.One()) === f(VectorizationBase.One(), vu) === vu
            @test f(vm, VectorizationBase.One()) === f(VectorizationBase.One(), vm) === vm
        end
        vnan = NaN * vx
        for f ∈ [fma, muladd, VectorizationBase.vfma_fast, VectorizationBase.vmuladd_fast]
            @test f(vnan, VectorizationBase.Zero(), vx) === vx
            @test f(VectorizationBase.Zero(), vnan, vx) === vx
        end
    end

    @inline function vlog(x1::VectorizationBase.AbstractSIMD{W,Float64}) where {W} # Testing if an assorted mix of operations
        x2 = reinterpret(UInt64, x1)
        x3 = x2 >>> 0x0000000000000020
        greater_than_zero = x1 > zero(x1)
        alternative = VectorizationBase.ifelse(x1 == zero(x1), -Inf, NaN)
        isinf = x1 == Inf
        x5 = x3 + 0x0000000000095f62
        x6 = x5 >>> 0x0000000000000014
        x7 = x6 - 0x00000000000003ff
        x8 = convert(Float64, x7 % Int)
        x9 = x5 << 0x0000000000000020
        x10 = x9 & 0x000fffff00000000
        x11 = x10 + 0x3fe6a09e00000000
        x12 = x2 & 0x00000000ffffffff
        x13 = x11 | x12
        x14 = reinterpret(Float64, x13)

        x15 = x14 - 1.0
        x16 = x15 * x15
        x17 = 0.5 * x16
        x18 = x14 + 1.0
        x19 = x15 / x18
        x20 = x19 * x19
        x21 = x20 * x20
        x22 = vfmadd(x21, 0.15313837699209373, 0.22222198432149784)
        x23 = vfmadd(x21, x22, 0.3999999999940942)
        x24 = x23 * x21
        x25 = vfmadd(x21, 0.14798198605116586, 0.1818357216161805)
        x26 = vfmadd(x21, x25, 0.2857142874366239)
        x27 = vfmadd(x21, x26, 0.6666666666666735)
        x28 = x27 * x20
        x29 = x24 + x17
        x30 = x29 + x28
        x31 = x8 * 1.9082149292705877e-10
        x32 = vfmadd(x19, x30, x31)
        x33 = x15 - x17
        x34 = x33 + x32
        x35 = vfmadd(x8, 0.6931471803691238, x34)
        x36 = VectorizationBase.ifelse(greater_than_zero, x35, alternative)
        VectorizationBase.ifelse(isinf, Inf, x36)
    end

    println("Defining log")
    @time @testset "Defining log." begin
        vx = Vec(ntuple(_ -> rand(), VectorizationBase.StaticInt(3) * VectorizationBase.pick_vector_width(Float64))...)
        check_within_limits(tovector(@inferred(vlog(vx))), log.(tovector(vx)))
    end

    println("Saturated add")
    @time @testset "Saturated add" begin
        @test VectorizationBase.saturated_add(0xf0, 0xf0) === 0xff
        @test VectorizationBase.saturated_add(2_000_000_000 % Int32, 1_000_000_000 % Int32) === typemax(Int32)
        v = Vec(ntuple(_ -> rand(typemax(UInt)>>1+one(UInt):typemax(UInt)), VectorizationBase.pick_vector_width(UInt))...)
        @test VectorizationBase.saturated_add(v, v) === vbroadcast(VectorizationBase.pick_vector_width(UInt), typemax(UInt))
    end

    println("Special Functions")
    @time @testset "Special Functions" begin
        let T = Float64
            min_non_denormal = nextfloat(abs(reinterpret(T, typemax(Base.uinttype(T)) & (~Base.exponent_mask(T)))))
            l2mnd = log2(min_non_denormal)
            xx = collect(range(0.8l2mnd, 0.8abs(l2mnd), length = 2^20));
            test_acc(exp2, exp2, T, xx, 3)

            lemnd = log(min_non_denormal)
            xx .= range(0.8lemnd, 0.8abs(lemnd), length = 2^20);
            test_acc(exp, exp, T, xx, 3)

            l10mnd = log10(min_non_denormal)
            xx .= range(0.8l10mnd, 0.8abs(l10mnd), length = 2^20);
            test_acc(exp10, exp10, T, xx, 3)

            # xx .= exp2.(range(0.8l2mnd, 0.8abs(l2mnd), length = 2^20));
            # test_acc(VectorizationBase.vlog2, log2, T, xx, 7)
        end
    end
    # end
end # @testset VectorizationBase

            # ptr_A = pointer(A)
            # vA = VectorizationBase.stridedpointer(A)
            # Att = copy(A')'
            # vAtt = VectorizationBase.stridedpointer(Att)
            # @test eltype(vA) == Float64
            # @test Base.unsafe_convert(Ptr{Float64}, vA) === ptr_A === pointer(vA)
            # @test vA == VectorizationBase.stridedpointer(vA)
            # @test all(i -> A[i+1] === VectorizationBase.vload(ptr_A + 8i) === VectorizationBase.vload(vA, (i,)) === Float64(i), 0:15)
            # VectorizationBase.vstore!(vA, 99.9, (3,))
            # @test 99.9 === VectorizationBase.vload(ptr_A + 8*3) === VectorizationBase.vload(vA, (VectorizationBase.StaticInt(3),)) === VectorizationBase.vload(vA, (3,0)) === A[4,1]
            # VectorizationBase.vstore!(vAtt, 99.9, (3,1))
            # @test 99.9 === VectorizationBase.vload(vAtt, (3,1)) === VectorizationBase.vload(vAtt, (VectorizationBase.StaticInt(3),1)) === Att[4,2]
            # VectorizationBase.vnoaliasstore!(ptr_A+8*4, 999.9)
            # @test 999.9 === VectorizationBase.vload(ptr_A + 8*4) === VectorizationBase.vload(pointer(vA), 4*sizeof(eltype(A))) === VectorizationBase.vload(vA, (4,))
            # @test vload(vA, (7,2)) == vload(vAtt, (7,2)) == A[8,3]
            # @test vload(VectorizationBase.subsetview(vA, Val(1), 7), (2,)) == vload(VectorizationBase.subsetview(vAtt, Val(1), 7), (2,)) == A[8,3]
            # @test vload(VectorizationBase.subsetview(vA, Val(2), 2), (7,)) == vload(VectorizationBase.subsetview(vAtt, Val(2), 2), (7,)) == A[8,3]
            #     @test vload(VectorizationBase.double_index(vA, Val(0), Val(1)), (2,)) == vload(VectorizationBase.double_index(vA, Val(0), Val(1)), (VectorizationBase.StaticInt(2),)) == A[3,3]
            #     @test vload(VectorizationBase.double_index(vAtt, Val(0), Val(1)), (1,)) == vload(VectorizationBase.double_index(vAtt, Val(0), Val(1)), (VectorizationBase.StaticInt(1),)) == A[2,2]
            #     B = rand(5, 5)
            # vB = VectorizationBase.stridedpointer(B)
            # @test vB[1, 2] == B[2, 3] == vload(VectorizationBase.stridedpointer(B, 2, 3))
            # @test vB[3] == B[4] == vload(VectorizationBase.stridedpointer(B, 4))
            # @test vload(Vec{4,Float64}, vB) == Vec{4,Float64}(ntuple(i->B[i], Val(4)))


text:
module WebAPI

using HTTP
using HTTP.URIs
using JSON3
using Sockets

export App, Req
export add_get!, add_post!, add_put!, add_patch!, add_delete!
export serve

include("app.jl")
include("req.jl")
include("routing.jl")
include("serve.jl")

end


text:
### # Running value divided by average i.e. x/mean(x)
### using JuliaDB
### t = table([1,2,3,4,5], [1,1,1,2,2], names=[:t,:x])
###
### @time t2 = pushcol(t, :z,
### 		map(mr->mr.x/mean(select(filter(row->row.t <= mrg.t, t), :x) ),
### 		    t)
### 	);
###


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

sigma = b"expand 32-byte k"

mutable struct ChaCha20
    block::Vector{UInt32}
    pos::UInt32
    xor_table::Vector{UInt8}
    xor_table_pos::Int

    function ChaCha20(key::Vector{UInt8}, nonce::Vector{UInt8}, pos::Int32)
        block_bytes::Vector{UInt8} = vcat(sigma, key, int32_to_bytes(pos), nonce)
        @assert length(block_bytes) == 64

        block::Vector{UInt32} = []
        for i in range(1, length(block_bytes), step=4)
            append!(block, bytes_to_uint32(block_bytes[i:i+3]))
        end
        xor_table = chacha20_round_bytes(block)

        new(block, UInt32(pos), xor_table, 1)
    end
end

function rotate(x::UInt32, n::Int)::UInt32
    (x << n) | (x >> (32 - n))
end

function quaterround(a::UInt32, b::UInt32, c::UInt32, d::UInt32)
    a += b
    d = xor(d, a)
    d = rotate(d, 16)

    c += d
    b = xor(b, c)
    b = rotate(b, 12)

    a += b
    d = xor(d, a)
    d = rotate(d, 8)

    c += d
    b = xor(b, c)
    b = rotate(b, 7)

    a, b, c, d
end

function chacha20_round_bytes(block::Vector{UInt32})::Vector{UInt8}
    x = copy(block)

    for _ in range(1, length=10)
        # column rounds
        x[1], x[5], x[9], x[13] = quaterround(x[1], x[5], x[9], x[13])
        x[2], x[6], x[10], x[14] = quaterround(x[2], x[6], x[10], x[14])
        x[3], x[7], x[11], x[15] = quaterround(x[3], x[7], x[11], x[15])
        x[4], x[8], x[12], x[16] = quaterround(x[4], x[8], x[12], x[16])
        # diagonal rounds
        x[1], x[6], x[11], x[16] = quaterround(x[1], x[6], x[11], x[16])
        x[2], x[7], x[12], x[13] = quaterround(x[2], x[7], x[12], x[13])
        x[3], x[8], x[9], x[14] = quaterround(x[3], x[8], x[9], x[14])
        x[4], x[5], x[10], x[15] = quaterround(x[4], x[5], x[10], x[15])
    end

    for i in range(1, length=16)
        x[i] += block[i]
    end

    r::Vector{UInt8} = []
    for i in x
        r = vcat(r, uint32_to_bytes(i))
    end
    r
end

function translate(chacha20::ChaCha20, plain::Vector{UInt8})
    enc::Vector{UInt8} = []

    for i in range(1, length=length(plain))
        append!(enc, xor(plain[i], chacha20.xor_table[chacha20.xor_table_pos]))
        chacha20.xor_table_pos += 1
        if chacha20.xor_table_pos > length(chacha20.xor_table)
            chacha20.block[13] += 1
            chacha20.xor_table = chacha20_round_bytes(chacha20.block)
            chacha20.xor_table_pos = 1
        end
    end

    enc
end


text:
function gendercount(nds::NameDataset, name::AbstractString, years::Union{T, Vector{T}, AbstractRange{T}}) where T <: Integer
    if any(y-> firstyear(nds) > y || y > lastyear(nds), years)
        @debug "Database only has years from $(nds.firstyear)-$(nds.lastyear)" maxlog = 1
        return (female=missing, male=missing)
    elseif !haskey(nds, name)
        @debug "Database doesn't have name $name"
        return (female=0, male=0)
    else
        return (
            female = sum(nds[name][:female][years]),
            male   = sum(nds[name][:male][years])
            )
    end
end

gendercount(name::AbstractString, years::Union{T, Vector{T}, AbstractRange{T}}) where T <: Integer = gendercount(NAMES, name, years)
gendercount(nds::NameDataset, name::AbstractString) = gendercount(nds, name, firstyear(nds):lastyear(nds))
gendercount(name::AbstractString) = gendercount(NAMES, name)

function proportionfemale(name, years=1880:2018)
    (f, m) = gendercount(name, years)
    total = f+m
    (ismissing(total) || total == 0) && return missing
    return f / total
end

proportionmale(name, years=1880:2018) = 1. - proportionfemale(name, years)

function gender(name::AbstractString, year; threshold::AbstractFloat=0.5)
    0.5 ≤ threshold ≤ 1. || raise(ArgumentError("Threshold must be between 0.5 and 1"))
    pf = proportionfemale(name, year)

    ismissing(pf) && return missing

    if pf > threshold
        return :female
    elseif pf < 1. - threshold
        return :male
    else
        return missing
    end
end

gender(name::AbstractString; threshold::AbstractFloat=0.5) = gender(name, 1880:2017, threshold=threshold)


text:
"""
    Entity{S,D,I}
Immutable data structure to represent entities (vehicle, pedestrian, ...).
Entities are defined by a state, a definition, and an id.
The state of an entity usually models changing values while the definition and the id should not change.

# Constructor

`Entity(state, definition, id)`

Copy constructor that keeps the definition and id but changes the state (a new object is still created):

`Entity(entity::Entity{S,D,I}, s::S)`

# Fields

- `state::S`
- `def::D`
- `id::I`
"""
struct Entity{S,D,I} # state, definition, identification
    state::S
    def::D
    id::I
end
Entity(entity::Entity{S,D,I}, s::S) where {S,D,I} = Entity(s, entity.def, entity.id)