diff --git a/StraightLinePrograms/.travis.yml b/StraightLinePrograms/.travis.yml new file mode 100644 index 000000000000..06e1b62c2621 --- /dev/null +++ b/StraightLinePrograms/.travis.yml @@ -0,0 +1,16 @@ +## Documentation: http://docs.travis-ci.com/user/languages/julia/ +language: julia +os: + - linux + - osx +julia: +# - 1.0 + - 1.1 + - 1.2 + - 1.3 + - 1.4 + - nightly +notifications: + email: false +git: + depth: 99999999 diff --git a/StraightLinePrograms/Project.toml b/StraightLinePrograms/Project.toml new file mode 100644 index 000000000000..687592ba21bd --- /dev/null +++ b/StraightLinePrograms/Project.toml @@ -0,0 +1,12 @@ +name = "StraightLinePrograms" +uuid = "42eed0b5-a112-4d56-a5b1-e566078f5bf3" +authors = ["Rafael Fourquet "] +version = "0.1.0" + +[deps] +AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +julia = "1" +AbstractAlgebra = "^0.9.1" diff --git a/StraightLinePrograms/README.md b/StraightLinePrograms/README.md new file mode 100644 index 000000000000..6757884cd24f --- /dev/null +++ b/StraightLinePrograms/README.md @@ -0,0 +1,419 @@ +[![Build Status](https://travis-ci.org/rfourquet/StraightLinePrograms.jl.svg?branch=master)](https://travis-ci.org/rfourquet/StraightLinePrograms.jl) + +# StraightLinePrograms + +This is a WIP implementation of straight-line programs (SLP) +This is part of the [Oscar](https://oscar.computeralgebra.de/) project. + + +## Introduction + +The main SLP type is `SLProgram`, to which other types can "compile" (or +"transpile"). The easiest way to create an `SLProgram` is to combine +"generators": + +```julia +julia> using StraightLinePrograms; const SL = StraightLinePrograms + +julia> x, y, z = gens(SLProgram, 3) +3-element Array{SLProgram{Union{}},1}: + x + y + z + +julia> p = (x*y^2 + 1.3*z)^-1 +#1 = ^ y 2 ==> y^2 +#2 = * x #1 ==> (xy^2) +#3 = * 1.3 z ==> (1.3z) +#4 = + #2 #3 ==> ((xy^2) + (1.3z)) +#5 = ^ #4 -1 ==> ((xy^2) + (1.3z))^-1 +return: #5 +``` + +On the right side of the above output is the representation of the computation +so far. It's done via another SLP type (tentatively) called `Free` which +represent "formulas" as trees: + +```julia +julia> X, Y, Z = gens(Free, 3) +3-element Array{Free,1}: + x + y + z + +julia> q = (X*Y^2 + 1.3*Z)^-1 +((xy^2) + (1.3z))^-1 + +julia> f = evaluate(p, [X, Y, Z]) +((xy^2) + (1.3z))^-1 + +julia> evaluate(f, [X, Y, Z]) == f +true + +julia> evaluate(p, Any[x, y, z]) == p +true + +julia> dump(q) # q::Free is a tree +Free + x: StraightLinePrograms.Exp + p: StraightLinePrograms.Plus + xs: Array{StraightLinePrograms.Lazy}((2,)) + 1: StraightLinePrograms.Times + xs: Array{StraightLinePrograms.Lazy}((2,)) + 1: StraightLinePrograms.Input + n: Int64 1 + 2: StraightLinePrograms.Exp + p: StraightLinePrograms.Input + n: Int64 2 + e: Int64 2 + 2: StraightLinePrograms.Times + xs: Array{StraightLinePrograms.Lazy}((2,)) + 1: StraightLinePrograms.Const{Float64} + c: Float64 1.3 + 2: StraightLinePrograms.Input + n: Int64 3 + e: Int64 -1 + gens: Array{Symbol}((3,)) + 1: Symbol x + 2: Symbol y + 3: Symbol z +``` + +Evaluation of SLPs is done via `evaluate`, which can take a vector of +anything which supports the operations used in the SLP (e.g. `*`, `+` and `^` +in this example; `-` is also supported but division not yet). +Note that currently, the `eltype` of the input vector for `SLProgram` +must be a supertype of any intermediate computation (so it's always safe to +pass a `Vector{Any}`). + +```julia +julia> evaluate(p, [2.0, 3.0, 5.0]) +0.04081632653061224 + +julia> evaluate(X*Y^2, ['a', 'b']) +"abb" +``` + +## Returning multiple values + +There is a low-level interface to return multiple values from an `SLProgram`; +for example, to return the second and last intermediate values from `p` +above, we would "assign" these values to positions `#1` and `#2`, +delete all other positions (via the "keep" operation), and return the +resulting array (the one used for intermediate computations): + +```julia +julia> SL.pushop!(p, SL.assign, SL.Arg(2), SL.Arg(1)) + SL.pushop!(p, SL.assign, SL.Arg(5), SL.Arg(2)) + SL.pushop!(p, SL.keep, SL.Arg(2)) + SL.setmultireturn!(p) +#1 = ^ y 2 ==> y^2 +#2 = * x #1 ==> (xy^2) +#3 = * 1.3 z ==> (1.3z) +#4 = + #2 #3 ==> ((xy^2) + (1.3z)) +#5 = ^ #4 -1 ==> ((xy^2) + (1.3z))^-1 +#1 = #2 ==> (xy^2) +#2 = #5 ==> ((xy^2) + (1.3z))^-1 +keep: #1..#2 +return: [#1, #2] + +julia> evaluate(p, [X, Y, Z]) +2-element Array{Free,1}: + (xy^2) + ((xy^2) + (1.3z))^-1 +``` + +## Straight line decisions + +A "decision" is a special operation which allows to stop prematurely the +execution of the program if a condition is `false`, and the program returns +`true` if no condition failed. +Currently, the interface is modeled after GAP's SLPs and defaults to testing +the `AbstractAlgebra.order` of an element. More specifically, +`test(prg, n::Integer)` tests whether the order of the result of `prg` is +equal to `n`, and `dec1 & dec2` chains two programs with a short-circuiting +behavior: + +```julia +julia> p1 = SL.test(x*y^2, 2) +#1 = ^ y 2 ==> y^2 +#2 = * x #1 ==> (xy^2) +test: order(#2) == 2 || return false +return: true + +julia> p2 = SL.test(y, 4) +test: order(y) == 4 || return false +return: true + +julia> p1 & p2 +#1 = ^ y 2 ==> y^2 +#2 = * x #1 ==> (xy^2) +test: order(#2) == 2 || return false +test: order(y) == 4 || return false +return: true + +julia> evaluate(p1 & p2, [X, Y]) +test((xy^2), 2) & test(y, 4) + +julia> using AbstractAlgebra; perm1, perm2 = perm"(1, 4)", perm"(1, 3, 4, 2)"; + +julia> evaluate(p1 & p2, [perm1, perm2]) +true + +julia> evaluate(p1 & p2, [perm2, perm1]) +false +``` + +## GAP's SLPs + +There are two other available SLP types: `GAPSLProgram` and `AtlasSLProgram`, +and related `GAPSLDecision` and `AtlasSLDecision`, which are constructed +similarly as in GAP: + +```julia +julia> prg = GAPSLProgram( [ [1,2,2,3], [3,-1] ], 2 ) +# input: +r = [ g1, g2 ] +# program: +r[3] = r[1]^2*r[2]^3 +r[4] = r[3]^-1 +# return value: +r[4] + +julia> evaluate(prg, [perm1, perm2]) +(1,3,4,2) + +julia> evaluate(prg, [x, y]) +#1 = ^ x 2 ==> x^2 +#2 = ^ y 3 ==> y^3 +#3 = * #1 #2 ==> (x^2y^3) +#4 = ^ #3 -1 ==> (x^2y^3)^-1 +return: #4 + +julia> SLProgram(prg) # direct compilation (with room for optimizations obviously) +#1 = x ==> x +#2 = y ==> y +#3 = ^ #1 2 ==> x^2 +#4 = ^ #2 3 ==> y^3 +#5 = * #3 #4 ==> (x^2y^3) +#3 = #5 ==> (x^2y^3) +keep: #1..#3 +#4 = ^ #3 -1 ==> (x^2y^3)^-1 +keep: #1..#4 +return: #4 + +julia> GAPSLProgram( [ [2,3], [ [3,1,1,4], [1,2,3,1] ] ], 2 ) +# input: +r = [ g1, g2 ] +# program: +r[3] = r[2]^3 +# return values: +[ r[3]*r[1]^4, r[1]^2*r[3] ] + +julia> GAPSLDecision([ [ [ 1, 1, 2, 1 ], 3 ], [ "Order", 1, 2 ], [ "Order", 2, 3 ], [ "Order", 3, 5 ] ] ) +# input: +r = [ g1, g2 ] +# program: +r[3] = r[1]*r[2] +order( r[1] ) == 2 || return false +order( r[2] ) == 3 || return false +order( r[3] ) == 5 || return false +# return value: +true + +julia> SLProgram(ans) +#1 = x ==> x +#2 = y ==> y +#3 = * #1 #2 ==> (xy) +keep: #1..#3 +test: order(#1) == 2 || return false +test: order(#2) == 3 || return false +test: order(#3) == 5 || return false +return: true + +julia> d = AtlasSLDecision("inp 2\nchor 1 2\nchor 2 3\nmu 1 2 3\nchor 3 5") +inp 2 +chor 1 2 +chor 2 3 +mu 1 2 3 +chor 3 5 + +376> evaluate(d, [perm1, perm2]) +false + +julia> GAPSLDecision(d) +# input: +r = [ g1, g2 ] +# program: +order( r[1] ) == 2 || return false +order( r[2] ) == 3 || return false +r[3] = r[1]*r[2] +order( r[3] ) == 5 || return false +# return value: +true + +julia> SLProgram(d) +#1 = x ==> x +#2 = y ==> y +test: order(#1) == 2 || return false +test: order(#2) == 3 || return false +#3 = * #1 #2 ==> (xy) +keep: #1..#3 +test: order(#3) == 5 || return false +return: true +``` + +## AbstractAlgebra's polynomial interface + +
+ +This is the initial API of SLPs which hasn't been updated in a while +and might not work as-is with the current state of the package. + + +Currently, SLPs have a polynomial interface (`SLPoly`). + +## Examples + +```julia +julia> using AbstractAlgebra, StraightLinePrograms, BenchmarkTools; + +julia> S = SLPolyRing(zz, [:x, :y]); x, y = gens(S) +2-element Array{SLPoly{Int64,SLPolyRing{Int64,AbstractAlgebra.Integers{Int64}}},1}: + x + y + +julia> p = 3 + 2x * y^2 # each line of the SLP is shown with current value + #1 = * 2 x ==> (2x) + #2 = ^ y 2 ==> y^2 + #3 = * #1 #2 ==> (2xy^2) + #4 = + 3 #3 ==> (3 + (2xy^2)) + +julia> p.cs # constants used in the program +2-element Array{Int64,1}: + 3 + 2 + +julia> p.lines # each line is a UInt64 encoding an opcode and 2 arguments + Line(0x0500000028000001) + Line(0x8780000020000002) + Line(0x0500000030000004) + Line(0x0300000010000005) + +julia> evaluate(p, [2, 3]) +39 + +julia> p2 = (p*(x*y))^6 +#1 = * 2 x ==> (2x) +#2 = ^ y 2 ==> y^2 +#3 = * #1 #2 ==> (2xy^2) +#4 = + 3 #3 ==> (3 + (2xy^2)) +#5 = * x y ==> (xy) +#6 = * #4 #5 ==> ((3 + (2xy^2))xy) +#7 = ^ #6 6 ==> ((3 + (2xy^2))xy)^6 + +julia> R, (x1, y1) = PolynomialRing(zz, ["x", "y"]); R +Multivariate Polynomial Ring in x, y over Integers + +julia> q = convert(R, p2) +64*x^12*y^18+576*x^11*y^16+2160*x^10*y^14+4320*x^9*y^12+4860*x^8*y^10+2916*x^7*y^8+729*x^6*y^6 + +julia> v = [3, 5]; @btime evaluate($q, $v) + 32.101 μs (634 allocations: 45.45 KiB) +-1458502820125772303 + +julia> @btime evaluate($p2, $v) + 270.341 ns (4 allocations: 352 bytes) +-1458502820125772303 + +julia> res = Int[]; @btime StraightLinePrograms.evaluate!($res, $p2, $v) + 171.013 ns (0 allocations: 0 bytes) +-1458502820125772303 + +julia> res # intermediate computations (first 2 elements are constants) +9-element Array{Int64,1}: + 3 + 2 + 6 + 25 + 150 + 153 + 15 + 2295 + -1458502820125772303 + +julia> f2 = StraightLinePrograms.compile!(p2) # compile to machine code +#3 (generic function with 1 method) + +julia> @btime evaluate($p2, $v) + 31.153 ns (1 allocation: 16 bytes) +-1458502820125772303 + +julia> @btime $f2($v) # use a function barrier for last bit of efficiency + 7.980 ns (0 allocations: 0 bytes) +-1458502820125772303 + +julia> q +64*x^12*y^18+576*x^11*y^16+2160*x^10*y^14+4320*x^9*y^12+4860*x^8*y^10+2916*x^7*y^8+729*x^6*y^6 + +julia> p3 = convert(S, q) # convert back q::Mpoly to an SLPoly + #1 = ^ x 6 ==> x^6 + #2 = ^ x 7 ==> x^7 + #3 = ^ x 8 ==> x^8 + #4 = ^ x 9 ==> x^9 + #5 = ^ x 10 ==> x^10 + #6 = ^ x 11 ==> x^11 + #7 = ^ x 12 ==> x^12 + #8 = ^ y 6 ==> y^6 + #9 = ^ y 8 ==> y^8 +#10 = ^ y 10 ==> y^10 +#11 = ^ y 12 ==> y^12 +#12 = ^ y 14 ==> y^14 +#13 = ^ y 16 ==> y^16 +#14 = ^ y 18 ==> y^18 +#15 = * 64 #7 ==> (64x^12) +#16 = * #15 #14 ==> (64x^12y^18) +#17 = * 576 #6 ==> (576x^11) +#18 = * #17 #13 ==> (576x^11y^16) +#19 = + #16 #18 ==> ((64x^12y^18) + (576x^11y^16)) +#20 = * 2160 #5 ==> (2160x^10) +#21 = * #20 #12 ==> (2160x^10y^14) +#22 = + #19 #21 ==> ((64x^12y^18) + (576x^11y^16) + (2160x^10y^14)) +#23 = * 4320 #4 ==> (4320x^9) +#24 = * #23 #11 ==> (4320x^9y^12) +#25 = + #22 #24 ==> ((64x^12y^18) + (576x^11y^16) + (2160x^10y^14) + (4320x^9y^12)) +#26 = * 4860 #3 ==> (4860x^8) +#27 = * #26 #10 ==> (4860x^8y^10) +#28 = + #25 #27 ==> ((64x^12y^18) + (576x^11y^16) + (2160x^10y^14) + (4320x^9y^12) + (4860x^8y^10)) +#29 = * 2916 #2 ==> (2916x^7) +#30 = * #29 #9 ==> (2916x^7y^8) +#31 = + #28 #30 ==> ((64x^12y^18) + (576x^11y^16) + (2160x^10y^14) + (4320x^9y^12) + (4860x^8y^10) + (2916x^7y^8)) +#32 = * 729 #1 ==> (729x^6) +#33 = * #32 #8 ==> (729x^6y^6) +#34 = + #31 #33 ==> ((64x^12y^18) + (576x^11y^16) + (2160x^10y^14) + (4320x^9y^12) + (4860x^8y^10) + (2916x^7y^8) + (729x^6y^6)) + +julia> @btime evaluate($p3, $v) + 699.465 ns (5 allocations: 1008 bytes) +-1458502820125772303 + +julia> @time f3 = StraightLinePrograms.compile!(p3); + 0.002830 seconds (1.40 k allocations: 90.930 KiB) + +julia> @btime $f3($v) + 80.229 ns (0 allocations: 0 bytes) +-1458502820125772303 + +julia> p4 = convert(S, q, limit_exp=true); # use different encoding + +julia> @btime evaluate($p4, $v) + 766.864 ns (5 allocations: 1008 bytes) +-1458502820125772303 + +julia> @time f4 = StraightLinePrograms.compile!(p4); + 0.002731 seconds (1.74 k allocations: 108.676 KiB) + +julia> @btime $f4($v) + 11.781 ns (0 allocations: 0 bytes) +-1458502820125772303 +``` +
diff --git a/StraightLinePrograms/src/StraightLinePrograms.jl b/StraightLinePrograms/src/StraightLinePrograms.jl new file mode 100644 index 000000000000..cfb12ddb090c --- /dev/null +++ b/StraightLinePrograms/src/StraightLinePrograms.jl @@ -0,0 +1,28 @@ +module StraightLinePrograms + +import Base: +, -, *, ^, parent + + +using AbstractAlgebra: AbstractAlgebra, Ring, RingElement, RingElem, MPolyRing, MPolyElem, elem_type, Generic + +import AbstractAlgebra: base_ring, gen, gens, ngens, nvars, symbols, evaluate, order, addeq! + +export LazyPolyRing, LazyPoly, SLPolyRing, SLPoly, SLProgram +export AbstractSLProgram, GAPSLProgram, GAPSLDecision, AtlasSLProgram, AtlasSLDecision, Free +export slpcst, slpgen, slpgens, compile, gens, evaluate, nsteps, compose, list + +abstract type AbstractSLProgram end + +(::Type{SLP})(p::AbstractSLProgram) where {SLP<:AbstractSLProgram} = + compile(SLP, p) + + +include("straightline.jl") +include("lazy.jl") +include("lazypolys.jl") +include("slpolys.jl") +include("gap.jl") +include("atlas.jl") + + +end # module diff --git a/StraightLinePrograms/src/atlas.jl b/StraightLinePrograms/src/atlas.jl new file mode 100644 index 000000000000..8617c44b5469 --- /dev/null +++ b/StraightLinePrograms/src/atlas.jl @@ -0,0 +1,266 @@ +struct AtlasLine + cmd::Symbol + dst::Int + i::Int + j::Int +end + +const ATLAS_VOID_SLOT = 0 + +AtlasLine(cmd, dst, i) = AtlasLine(cmd, dst, i, ATLAS_VOID_SLOT) + +abstract type AbstractAtlasSL <: AbstractSLProgram end + +struct AtlasSLProgram <: AbstractAtlasSL + code::String # source code + echo::Vector{String} # echo lines + ngens::Int + outputs::Vector{Int} + lines::Vector{AtlasLine} +end + +struct AtlasSLDecision <: AbstractAtlasSL + code::String # source code + echo::Vector{String} # echo lines + ngens::Int + lines::Vector{AtlasLine} +end + +function (::Type{SLP})(code::String) where SLP <: AbstractAtlasSL + codelines = split(code, "\n", keepempty=false) + labels = String[] + outputs = Int[] + lines = AtlasLine[] + echo = String[] + + ngens::Int = 0 + + function getidx!(label, create=false) + # no "inp" line, set defaults + isempty(labels) && push!(labels, "1", "2") + if ngens == 0 + ngens = length(labels) + end + + i = findfirst(==(label), labels) + if i === nothing + if create + push!(labels, label) + lastindex(labels) + else + throw(ArgumentError("label $label does not exist")) + end + else + i + end + end + + for codeline in codelines + codeline0 = split(codeline, '#', limit=2)[1] # remove comments + codeline = split(codeline0, keepempty=false) + isempty(codeline) && continue + if codeline[1] == "echo" + push!(echo, strip(split(codeline0, "echo", limit=2)[2])) + continue + end + + length(codeline) < 2 && error_invalid_line(codeline) + cmd = Symbol(codeline[1]) + + if cmd == :inp + ngens != 0 && throw(ArgumentError("\"inp\" line not at the beginning")) + n = tryparse(Int, codeline[2]) + n === nothing && error_invalid_line(codeline) + if length(codeline) == 2 + isempty(labels) || + throw(ArgumentError("inp must not omit the names")) + append!(labels, string.(1:n)) + else + length(codeline) == 2+n || error_invalid_line(codeline) + append!(labels, codeline[3:end]) + end + continue + end + + if cmd == :oup + SLP <: AtlasSLProgram || + throw(ArgumentError("\"oup\" line only allowed in AtlasSLProgram")) + n = tryparse(Int, codeline[2]) + n === nothing && error_invalid_line(codeline) + if length(codeline) == 2 + isempty(outputs) || + throw(ArgumentError("oup must not omit the names")) + append!(outputs, getidx!.(string.(1:n))) + else + length(codeline) == 2+n || error_invalid_line(codeline) + append!(outputs, getidx!.(codeline[3:end])) + end + continue + end + + !isempty(outputs) && throw(ArgumentError("\"oup\" line not at the end")) + + line = + if cmd == :cjr + check_line_length(codeline, 2) + args = getidx!.(codeline[2:end], (false, false)) + AtlasLine(:cj, args[1], args[1], args[2]) + elseif cmd in [:cj, :com, :mu] + check_line_length(codeline, 3) + args = getidx!.(codeline[2:end], (false, false, true)) + AtlasLine(cmd, args[3], args[1], args[2]) + elseif cmd in [:iv, :cp] + check_line_length(codeline, 2) + args = getidx!.(codeline[2:end], (false, true)) + AtlasLine(cmd, args[2], args[1]) + elseif cmd == :pwr + check_line_length(codeline, 3) + arg = getidx!(codeline[3]) + dst = getidx!(codeline[4], true) + AtlasLine(cmd, dst, parse(Int, codeline[2]), arg) + elseif cmd == :chor + SLP <: AtlasSLDecision || + throw(AtlasSLDecision("\"chor\" line only allowed in AtlasSLDecision")) + check_line_length(codeline, 2) + arg = getidx!(codeline[2]) + n = parse(Int, (codeline[3])) + AtlasLine(cmd, ATLAS_VOID_SLOT, arg, n) + else + error_invalid_line(codeline) + end + push!(lines, line) + end + + if SLP <: AtlasSLProgram + if isempty(outputs) + push!(outputs, getidx!("1"), getidx!("2")) + end + AtlasSLProgram(code, echo, ngens, outputs, lines) + else + AtlasSLDecision(code, echo, ngens, lines) + end +end + +check_line_length(codeline, n) = + length(codeline) == n+1 || throw(ArgumentError( + "wrong number of arguments in $(codeline[1]) line")) + +error_invalid_line(codeline) = + throw(ArgumentError("""invalid line: $(join(codeline, " "))""")) + + +## show + +Base.show(io::IO, p::AbstractAtlasSL) = print(io, p.code) + + +## evaluate + +evaluate(p::AbstractAtlasSL, xs::Vector) = evaluate!(empty(xs), p, xs) + +function evaluate!(res::Vector{S}, p::AbstractAtlasSL, xs::Vector{S}) where S + append!(empty!(res), view(xs, 1:p.ngens)) + + for l in p.lines + cmd, dst, i, j = l.cmd, l.dst, l.i, l.j + a, b = get(res, i, nothing), get(res, j, nothing) + + k = + if cmd == :cj + b^-1 * a * b + elseif cmd == :com + a^-1 * b^-1 * a * b + elseif cmd == :iv + @assert j == ATLAS_VOID_SLOT + a^-1 + elseif cmd == :mu + a * b + elseif cmd == :pwr + b^i + elseif cmd == :cp + @assert j == ATLAS_VOID_SLOT + a + elseif cmd == :chor + @assert dst == ATLAS_VOID_SLOT + if order(a) == j + continue + else + return false + end + else + @assert false "unexpected command" + end + + n = lastindex(res) + if dst == n + 1 + push!(res, k) + else + res[dst] = k + end + end + + if p isa AtlasSLProgram + n = lastindex(res) + # first copy the ouputs at the end of res before moving them + # at the beginning of the array, to avoid overwriting values + # which are still needed + # TODO: this algo can be slightly optimized + for i in p.outputs + push!(res, res[i]) + end + r = length(p.outputs) + for i = 1:r + res[i] = res[i+n] + end + resize!(res, r) + else + true + end +end + + +## compilation to GAPSLProgram + +compile(::Type{SLProgram}, p::AbstractAtlasSL) = + compile(SLProgram, compile(AbstractGAPSL, p)) + +compile(p::AbstractAtlasSL) = compile(SLProgram, p) + +function compile(::Type{SLP}, p::AbstractAtlasSL) where SLP<:AbstractGAPSL + SLP === AbstractGAPSL || + (SLP === GAPSLProgram) == (p isa AtlasSLProgram) || + throw(ArgumentError("cannot compile an $(typeof(p)) into $SLP")) + lines = [] + for l in p.lines + cmd, dst, i, j = l.cmd, l.dst, l.i, l.j + k = + if cmd == :cj + [j, -1, i, 1, j, 1] + elseif cmd == :com + [i, -1, j, -1, i, 1, j, 1] + elseif cmd == :iv + @assert j == ATLAS_VOID_SLOT + [i, -1] + elseif cmd == :mu + [i, 1, j, 1] + elseif cmd == :pwr + [j, i] + elseif cmd == :cp + @assert j == ATLAS_VOID_SLOT + [i, 1] + elseif cmd == :chor + @assert dst == ATLAS_VOID_SLOT + push!(lines, (i, j)) + continue + else + @assert false "unexpected command" + end + push!(lines, [k, dst]) + end + if p isa AtlasSLProgram + push!(lines, [[r, 1] for r in p.outputs]) + GAPSLProgram(lines, p.ngens) + else + GAPSLDecision(lines, p.ngens) + end +end diff --git a/StraightLinePrograms/src/gap.jl b/StraightLinePrograms/src/gap.jl new file mode 100644 index 000000000000..04fc2b587aec --- /dev/null +++ b/StraightLinePrograms/src/gap.jl @@ -0,0 +1,312 @@ +const ReturnList = Vector{Vector{Int}} + +const GAPStraightLine = Union{Vector{Int}, # e.g. [1, 2, 2, -1] + Tuple{Vector{Int},Int}, # e.g. ([1, 2], 2) + Tuple{Int,Int}, # (i, n) for ["Order", i, n] + ReturnList} # return list + +issimpleline(line) = line isa Vector{Int} +isassignline(line) = line isa Tuple{Vector{Int},Int} +isreturnline(line) = line isa ReturnList +isorderline(line) = line isa Tuple{Int,Int} + +abstract type AbstractGAPSL <: AbstractSLProgram end + +struct GAPSLProgram <: AbstractGAPSL + lines::Vector{GAPStraightLine} + ngens::Int + slp::Ref{SLProgram} +end + +struct GAPSLDecision <: AbstractGAPSL + lines::Vector{GAPStraightLine} + ngens::Int + slp::Ref{SLProgram} +end + +function (::Type{SLP})(lines::Vector, ngens::Integer=-1) where {SLP<:AbstractGAPSL} + ls = GAPStraightLine[] + n = length(lines) + ng = 0 + + have = ngens < 0 ? BitSet(0) : BitSet(0:ngens) + + function maxnohave(word) + local ng = 0 + skip = true + for i in word + skip = !skip + skip && continue + if !(i in have) + ng = max(ng, i) + end + end + ng + end + + undef_slot() = throw(ArgumentError("invalid use of undef memory")) + + for (i, line) in enumerate(lines) + pushline!(ls, line) + l = ls[end] + if ngens < 0 && issimpleline(l) && (i < n || SLP === GAPSLDecision) + throw(ArgumentError("ngens must be specified")) + end + if isassignline(l) + ng = max(ng, maxnohave(l[1])) + if ngens < 0 + union!(have, 1:ng) + elseif ng > 0 + undef_slot() + end + push!(have, l[2]) + elseif issimpleline(l) + ng = max(ng, maxnohave(l)) + ngens >= 0 && ng > 0 && undef_slot() + push!(have, maximum(have)+1) + elseif isorderline(l) + SLP <: GAPSLDecision || + throw(ArgumentError("\"Order\" line only allowed in GAPSLDecision")) + # GAP doesn't seem to take an "Order" line into account for determining + # the number of generators + else + @assert isreturnline(l) "unknown line" + SLP <: GAPSLProgram || + throw(ArgumentError("return list only allowed in GAPSLProgram")) + for li in l + ng = max(ng, maxnohave(li)) + ngens >= 0 && ng > 0 && undef_slot() + end + end + end + SLP(ls, ngens < 0 ? ng : ngens, Ref{SLProgram}()) +end + +invalid_list_error(line) = throw(ArgumentError("invalid line or list: $line")) + +check_element(list) = iseven(length(list)) || invalid_list_error(list) + +check_last(lines) = isempty(lines) || !isa(lines[end], ReturnList) || + throw(ArgumentError("return list only allowed in last position")) + +function pushline!(lines::Vector{GAPStraightLine}, line::Vector{Int}) + check_last(lines) + check_element(line) + push!(lines, line) +end + +function pushline!(lines::Vector{GAPStraightLine}, line::Vector{Vector{Int}}) + check_last(lines) + for l in line + check_element(l) + end + push!(lines, line) +end + +function pushline!(lines::Vector{GAPStraightLine}, line::Tuple{Vector{Int},Int}) + check_last(lines) + check_element(line[1]) + push!(lines, line) +end + +function pushline!(lines::Vector{GAPStraightLine}, line::Tuple{Int,Int}) + check_last(lines) # redundant + push!(lines, line) +end + +# catch-all function to parse input in GAP style +function pushline!(lines::Vector{GAPStraightLine}, line::AbstractVector) + newline = + if length(line) == 2 && + line[1] isa AbstractVector && line[2] isa Integer + all(isinteger, line[1]) || invalid_list_error(line) + (Vector{Int}(line[1]), Int(line[2])) # assign line + elseif length(line) == 3 && line[1] == "Order" + isinteger(line[2]) && isinteger(line[3]) || + invalid_list_error(line) + (Int(line[2]), Int(line[3])) + elseif isempty(line) || line[1] isa AbstractVector # return line + all(v -> v isa AbstractVector && all(isinteger, v), line) || + invalid_list_error(line) + Vector{Int}[Vector{Int}(l) for l in line] + elseif line isa AbstractVector # simple line + all(isinteger, line) || invalid_list_error(line) + Vector{Int}(line) + else + invalid_list_error(line) + end + pushline!(lines, newline) +end + + +## show + +expshow(x, i) = i == 1 ? "r[$x]" : "r[$x]^$i" +prodshow(io, x) = join(io, (expshow(x[2*i-1], x[2*i]) for i=1:(length(x)>>1)), '*') + +function Base.show(io::IO, p::AbstractGAPSL) + k = p.ngens + println(io, "# input:") + print(io, "r = [ ") + join(io, ("g$i" for i in 1:k), ", ") + println(io, " ]") + print("# program:") + + for line in p.lines + if issimpleline(line) + k += 1 + print(io, "\nr[$k] = ") + prodshow(io, line) + elseif isassignline(line) + l, k = line + print(io, "\nr[$k] = ") + prodshow(io, l) + elseif isorderline(line) + x, e = line + print(io, "\norder( r[$x] ) == $e || return false") + elseif isreturnline(line) + print("\n# return values:\n[ ") + for (i, l) in enumerate(line) + i == 1 || print(io, ", ") + prodshow(io, l) + end + return print(io, " ]") + else + @assert false "unknown line" + end + end + if p isa GAPSLProgram + print("\n# return value:\nr[$k]") + else + print("\n# return value:\ntrue") + end +end + + +## evaluate + +evaluate(p::AbstractGAPSL, xs::Vector) = evaluate!(empty(xs), p, xs) +# TODO: use compiled version if present? + +expterm(x, i) = i == 1 ? x : x^i +prodlist(res, x) = prod(expterm(res[x[2*i-1]], x[2*i]) for i=1:(length(x)>>1)) + +function evaluate!(res::Vector{S}, p::AbstractGAPSL, xs::Vector{S}) where {S} + append!(empty!(res), view(xs, 1:p.ngens)) + + local k + for line in p.lines + if issimpleline(line) + push!(res, prodlist(res, line)) + k = lastindex(res) + elseif isassignline(line) + k = line[2] + r = prodlist(res, line[1]) + if k == lastindex(res) + 1 + push!(res, r) + else + res[k] = r + end + elseif isorderline(line) + x = res[line[1]] + order(x) == line[2] || return false + else + k = lastindex(res) + 1 + for l in line + push!(res, prodlist(res, l)) + end + if k > 1 + for i in 1:length(line) + res[i] = res[k] + k += 1 + end + end + resize!(res, length(line)) + return res + end + end + if p isa GAPSLProgram + res[k] + else + true + end +end + + +## compile + +# compile to an SLProgram + +compile!(gp::AbstractGAPSL) = gp.slp[] = compile(gp) + +compile(gp::AbstractGAPSL) = compile(SLProgram, gp) + +function compile(::Type{SLProgram}, gp::AbstractGAPSL) + p = SLProgram() + local res::Arg + + for i = 1:gp.ngens + pushop!(p, assign, input(i)) + end + + multi = false + for line in gp.lines + @assert !multi + if issimpleline(line) + res = Arg(p.len + 1) + k = write_list!(p, line) + if res != k + pushop!(p, assign, k, res) + end + pushop!(p, keep, res) + elseif isassignline(line) + list, dst = line + res = Arg(dst) + ptr = Arg(max(p.len, dst)) + k = write_list!(p, list) + if res != k + pushop!(p, assign, k, res) + end + pushop!(p, keep, ptr) + elseif isorderline(line) + x = Arg(line[1]) + ord = pushint!(p, line[2]) + pushop!(p, decision, x, ord) + else + reslist = Arg[] + for l in line + k = write_list!(p, l) + push!(reslist, k) + end + for (i, r) in enumerate(reslist) + @assert i != r.x # can happen? if so, don't assign + pushop!(p, assign, r, Arg(i)) + end + pushop!(p, keep, Arg(length(reslist))) + multi = true + end + end + if gp isa GAPSLDecision + setdecision!(p) + elseif !multi + pushfinalize!(p, res) + end + p +end + +function write_list!(p::SLProgram, list::Vector{Int}) + @assert !isempty(list) # TODO: handle empty lists + n = length(list) >> 1 + local k + for i = 1:n + x = Arg(list[2*i-1]) + e = list[2*i] + if e == 1 + l = x + else + l = pushop!(p, exponentiate, x, intarg(e)) + end + k = i == 1 ? l : pushop!(p, times, k, l) + end + k +end diff --git a/StraightLinePrograms/src/lazy.jl b/StraightLinePrograms/src/lazy.jl new file mode 100644 index 000000000000..f8953c523704 --- /dev/null +++ b/StraightLinePrograms/src/lazy.jl @@ -0,0 +1,582 @@ +## Free + +abstract type Lazy end + +struct Free <: AbstractSLProgram + x::Lazy # must not contain Gen + gens::Vector{Symbol} +end + +Free(x::Lazy) = Free(x, collect(Symbol, slpsyms(maxinput(x)))) + +Free(x) = Free(Const(x), Symbol[]) + +function freegens(n::Integer, syms=slpsyms(n)) + if !isa(syms, Vector{Symbol}) + syms = collect(Symbol, syms) + end + Free[Free(Input(i), syms) for i=1:n] +end + +freegens(syms::AbstractVector{Symbol}) = freegens(length(syms), syms) + +gens(::Type{Free}, n::Integer) = freegens(n) + +ngens(f::Free) = length(f.gens) + +gens(f::Free) = f.gens + +Base.show(io::IO, f::Free) = + show(IOContext(io, :slp_free_gens => gens(f)), f.x) + +evaluate(f::Free, xs) = evaluate(gens(f), f.x, xs) + +# check compatibility and return the biggest of the two gens arrays +function gens(x::Free, xs::Free...) + gx = gens(x) + for y in xs + gy = gens(y) + if length(gy) > length(gx) + gx, gy = gy, gx + end + gy == view(gx, eachindex(gy)) || throw(ArgumentError( + "incompatible symbols")) + end + gx +end + +# TODO: must they also have the same gens? +Base.:(==)(x::Free, y::Free) = x.x == y.x + +compile(::Type{SLProgram}, f::Free) = + compile(SLProgram{constantstype(f.x)}, f) + +function compile(::Type{SLProgram{T}}, f::Free) where T + p = SLProgram{T}() + i = pushlazy!(p, f.x, gens(f)) + pushfinalize!(p, i) +end + +compile(f::Free) = compile(SLProgram, f) + + +### unary/binary ops + ++(x::Free, y::Free) = Free(x.x + y.x, gens(x, y)) +-(x::Free, y::Free) = Free(x.x - y.x, gens(x, y)) +*(x::Free, y::Free) = Free(x.x * y.x, gens(x, y)) + +Base.:(&)(x::Free, y::Free) = Free(x.x & y.x, gens(x, y)) + +-(x::Free) = Free(-x.x, gens(x)) +^(x::Free, e::Integer) = Free(x.x^e, gens(x)) + +Base.literal_pow(::typeof(^), p::Free, ::Val{e}) where {e} = p^e + +test(x::Free, n) = Free(test(x.x, n), gens(x)) + +# TODO: don't splat +list(::Type{Free}, xs) = Free(List(Lazy[x.x for x in xs]), gens(xs...)) + +function compose(p::Free, q::Free; flatten=true) + if flatten + q.x isa List || throw(ArgumentError( + "first argument must return a list")) + gs = gens(q) + evaluate(p, [Free(qi, gs) for qi in q.x.xs])::Free + else + Free(Compose(p.x, q.x), gens(q)) + end +end + +Base.getindex(p::Free, is::Free...) = + Free(Getindex(p.x, Lazy[i.x for i in is]), gens(p, is...)) + + +#### adhoc + ++(x::Free, y) = Free(x.x + y, gens(x)) ++(x, y::Free) = Free(x + y.x, gens(y)) + +-(x::Free, y) = Free(x.x - y, gens(x)) +-(x, y::Free) = Free(x - y.x, gens(y)) + +*(x::Free, y) = Free(x.x * y, gens(x)) +*(x, y::Free) = Free(x * y.x, gens(y)) + +Base.getindex(x::Free, is...) = + getindex(x, map(f -> f isa Free ? f : Free(f), is)...) +Base.getindex(x, y::Free) = Free(x[y.x], gens(y)) + +# special case for integer literals + +function Base.getindex(x::Free, i::Integer) + if x.x isa List + Free(x.x.xs[i], gens(x)) + else + getindex(x, Free(i)) + end +end + +function Base.getindex(x::Free, is::AbstractVector) + T = eltype(is) + if x.x isa List && (T <: Integer || + typeintersect(T, Integer) !== Union{} && + all(x -> isa(x, Integer), is)) + Free(List(x.x.xs[is]), gens(x)) + else + getindex(x, Free(is)) + end +end + + +## Lazy + +gens(l::Lazy) = sort!(unique!(pushgens!(Symbol[], l)::Vector{Symbol})) + +Base.:(==)(k::Lazy, l::Lazy) = false + +Lazy(x::Lazy) = x +Lazy(x::AbstractSLProgram) = compile(Lazy, x) +Lazy(x) = Const(x) + +evaluate(l::Lazy, xs) = evaluate(gens(l), l, xs) +# TODO: remove the 3-arg evaluate methods? + + +### Const + +struct Const{T} <: Lazy + c::T + + # TODO: this is a hack, delete ASAP + Const{T}(x::AbstractSLProgram) where {T} = new{T}(x) + Const(x) = new{typeof(x)}(x) +end + +Const(x::AbstractSLProgram) = Const{typeof(x)}(x) + + +Base.show(io::IO, c::Const) = print(io, c.c) + +pushgens!(gs, c::Const) = gs +constantstype(l::Const{T}) where {T} = T + +Base.:(==)(k::Const, l::Const) = k.c == l.c + +evaluate(gs, c::Const, xs) = c.c + +maxinput(c::Const) = 0 + + +### Input + +struct Input <: Lazy + n::Int +end + +function Base.show(io::IO, i::Input) + syms = io[:slp_free_gens] + print(io, syms[i.n]) +end + +evaluate(gs, i::Input, xs) = xs[i.n] + +maxinput(i::Input) = i.n + +Base.:(==)(i::Input, j::Input) = i.n == j.n + +constantstype(::Input) = Union{} + + +### Gen + +struct Gen <: Lazy + g::Symbol +end + +Base.show(io::IO, g::Gen) = print(io, g.g) + +pushgens!(gs, g::Gen) = + if findfirst(==(g.g), gs) === nothing + push!(gs, g.g) + else + gs + end + +constantstype(l::Gen) = Union{} + +Base.:(==)(k::Gen, l::Gen) = k.g == l.g + +evaluate(gs, g::Gen, xs) = xs[findfirst(==(g.g), gs)] + +maxinput(g::Gen) = throw(ArgumentError("logic error: Gen not allowed in Free")) + + +### Plus + +struct Plus <: Lazy + xs::Vector{Lazy} +end + +Plus(xs::Lazy...) = Plus(collect(Lazy, xs)) + +function Base.show(io::IO, p::Plus) + print(io, '(') + join(io, p.xs, " + ") + print(io, ')') +end + +pushgens!(gs, p::Plus) = foldl(pushgens!, p.xs, init=gs) + +constantstype(p::Plus) = + mapreduce(constantstype, typejoin, p.xs, init=Union{}) + +Base.:(==)(k::Plus, l::Plus) = k.xs == l.xs + +evaluate(gs, p::Plus, xs) = mapreduce(q -> evaluate(gs, q, xs), +, p.xs) + +maxinput(p::Plus) = mapreduce(maxinput, max, p.xs) + + +### Minus + +struct Minus <: Lazy + p::Lazy + q::Lazy +end + +Base.show(io::IO, p::Minus) = print(io, '(', p.p, " - ", p.q, ')') + +pushgens!(gs, l::Minus) = pushgens!(pushgens!(gs, l.p), l.q) + +constantstype(m::Minus) = typejoin(constantstype(m.p), constantstype(m.q)) + +Base.:(==)(k::Minus, l::Minus) = k.p == l.p && k.q == l.q + +evaluate(gs, p::Minus, xs) = evaluate(gs, p.p, xs) - evaluate(gs, p.q, xs) + +maxinput(m::Minus) = max(maxinput(m.p), maxinput(m.q)) + + +### UniMinus + +struct UniMinus <: Lazy + p::Lazy +end + +Base.show(io::IO, p::UniMinus) = print(io, "(-", p.p, ')') + +pushgens!(gs, l::UniMinus) = pushgens!(gs, l.p) + +constantstype(p::UniMinus) = constantstype(p.p) + +Base.:(==)(k::UniMinus, l::UniMinus) = k.p == l.p + +evaluate(gs, p::UniMinus, xs) = -evaluate(gs, p.p, xs) + +maxinput(p::UniMinus) = maxinput(p.p) + + +### Times + +struct Times <: Lazy + xs::Vector{Lazy} +end + +Times(xs::Lazy...) = Times(collect(Lazy, xs)) + +function Base.show(io::IO, p::Times) + print(io, '(') + join(io, p.xs, '*') + print(io, ')') +end + +pushgens!(gs, p::Times) = foldl(pushgens!, p.xs, init=gs) + +constantstype(p::Times) = + mapreduce(constantstype, typejoin, p.xs, init=Union{}) + +Base.:(==)(k::Times, l::Times) = k.xs == l.xs + +evaluate(gs, p::Times, xs) = mapreduce(q -> evaluate(gs, q, xs), *, p.xs) + +maxinput(p::Times) = mapreduce(maxinput, max, p.xs) + + +### Exp + +struct Exp <: Lazy + p::Lazy + e::Int +end + +Base.show(io::IO, p::Exp) = print(io, p.p, '^', p.e) + +pushgens!(gs, l::Exp) = pushgens!(gs, l.p) + +constantstype(p::Exp) = constantstype(p.p) + +Base.:(==)(k::Exp, l::Exp) = k.p == l.p && k.e == l.e + +Base.literal_pow(::typeof(^), p::Lazy, ::Val{e}) where {e} = Exp(p, e) + +evaluate(gs, p::Exp, xs) = evaluate(gs, p.p, xs)^p.e + +maxinput(e::Exp) = maxinput(e.p) + + +### Decision + +struct Decision <: Lazy + ps::Vector{Tuple{Lazy,Int}} +end + +test(p::Lazy, n) = Decision(Tuple{Lazy,Int}[(p, n)]) + +Base.show(io::IO, d::Decision) = + join(io, (sprint(print, "test(", p, ", ", n, ")", context=io) + for (p, n) in d.ps), + " & ") + +constantstype(l::Decision) = + mapreduce(constantstype∘first, typejoin, l.ps, init=Union{}) + +pushgens!(gs, l::Decision) = + foldl((gs, d) -> pushgens!(gs, first(d)), l.ps, init=gs) + +function Base.:(&)(p::Decision, q::Decision) + d = Decision(copy(p.ps)) + append!(d.ps, q.ps) + d +end + +Base.:(==)(p::Decision, q::Decision) = p.ps == q.ps + +function evaluate(gs, p::Decision, xs) + (d, i), rest = Iterators.peel(p.ps) # p.ps should never be empty + res = test(evaluate(gs, d, xs), i) + res === false && return false + for (d, i) in rest + r = test(evaluate(gs, d, xs), i) + r === false && return false + res &= r + end + res +end + +maxinput(p::Decision) = mapreduce(x -> maxinput(x[1]), max, p.ps) + + +### List + +struct List <: Lazy + xs::Vector{Lazy} +end + +function Base.show(io::IO, l::List) + print(io, "list([") + join(io, l.xs, ", ") + print(io, "])") +end + +constantstype(l::List) = + mapreduce(constantstype, typejoin, l.xs, init=Union{}) + +pushgens!(gs, l::List) = foldl(pushgens!, l.xs, init=gs) + +Base.:(==)(p::List, q::List) = p.xs == q.xs + +evaluate(gs, l::List, xs) = + list(eltype(xs)[evaluate(gs, p, xs) for p in l.xs]) + +maxinput(l::List) = mapreduce(maxinput, max, l.xs) + + +### Compose + +struct Compose <: Lazy + p::Lazy + q::Lazy # q must return a list! + + function Compose(p::Lazy, q::Lazy) + q isa List || throw(ArgumentError( + "second argument of Compose must return a list")) + new(p, q) + end +end + +Base.show(io::IO, c::Compose) = print(io, c.p, " ∘ ", c.q) + +constantstype(c::Compose) = typejoin(constantstype(c.p, constantstype(c.q))) + +pushgens!(gs, c::Compose) = pushgens!(gs, c.q) + +Base.:(==)(k::Compose, l::Compose) = k.p == l.p && k.q == l.q + +evaluate(gs, p::Compose, xs) = evaluate(gs, p.p, evaluate(gs, p.q, xs)) + +maxinput(m::Compose) = maxinput(m.q) + + +### Getindex + +struct Getindex <: Lazy + p::Lazy + is::Vector{Lazy} +end + +function Base.show(io::IO, g::Getindex) + print(io, g.p, '[') + join(io, g.is, ", ") + print(io, ']') +end + +pushgens!(gs, l::Getindex) = foldl(pushgens!, l.is, init=pushgens!(gs, l.p)) + +constantstype(m::Getindex) = mapreduce(constantstype, typejoin, m.is, init=constantstype(m.p)) + +Base.:(==)(k::Getindex, l::Getindex) = k.p == l.p && k.is == l.is + +evaluate(gs, p::Getindex, xs) = evaluate(gs, p.p, xs)[(evaluate(gs, i, xs) for i in p.is)...] + +maxinput(m::Getindex) = mapreduce(maxinput, max, m.is, init=maxinput(m.p)) + + +### binary ops + +#### + + ++(x::Lazy, y::Lazy) = Plus(x, y) + +function +(x::Plus, y::Lazy) + p = Plus(copy(x.xs)) + push!(p.xs, y) + p +end + +function +(x::Lazy, y::Plus) + p = Plus(copy(y.xs)) + pushfirst!(p.xs, x) + p +end + +function +(x::Plus, y::Plus) + p = Plus(copy(x.xs)) + append!(p.xs, y.xs) + p +end + + +#### - + +-(p::Lazy, q::Lazy) = Minus(p, q) +-(p::Lazy) = UniMinus(p) + + +#### * + +*(x::Lazy, y::Lazy) = Times(x, y) + +function *(x::Times, y::Lazy) + p = Times(copy(x.xs)) + push!(p.xs, y) + p +end + +function *(x::Lazy, y::Times) + p = Times(copy(y.xs)) + pushfirst!(p.xs, x) + p +end + +function *(x::Times, y::Times) + p = Times(copy(x.xs)) + append!(p.xs, y.xs) + p +end + + +#### ^ + +^(x::Lazy, e::Integer) = Exp(x, Int(e)) + +#### getindex + +Base.getindex(x::Lazy, i::Lazy...) = Getindex(x, collect(i)) + + +### adhoc binary ops + +*(x, y::Lazy) = Const(x) * y +*(x::Lazy, y) = x * Const(y) + ++(x, y::Lazy) = Const(x) + y ++(x::Lazy, y) = x + Const(y) + +-(x, y::Lazy) = Const(x) - y +-(x::Lazy, y) = x - Const(y) + +Base.getindex(x, i::Lazy) = Getindex(Const(x), [i]) +Base.getindex(x::Lazy, is...) = Getindex(x, Lazy[i isa Lazy ? i : Const(i) for i in is]) + + +## compile to SLProgram + +# TODO: this is legacy, only for tests +SLProgram(l::Lazy) = compile(SLProgram, l) +SLProgram{T}(l::Lazy) where {T} = compile(SLProgram{T}, l) + +compile(::Type{SLProgram}, l::Lazy) = compile(SLProgram{constantstype(l)}, l) + +function compile(::Type{SLProgram{T}}, l::Lazy, gs=gens(l)) where T + p = SLProgram{T}() + i = pushlazy!(p, l, gs) + pushfinalize!(p, i) +end + +pushlazy!(p::SLProgram, l::Const, gs) = pushconst!(p, l.c) + +pushlazy!(p::SLProgram, l::Input, gs) = input(l.n) + +pushlazy!(p::SLProgram, l::Gen, gs) = input(findfirst(==(l.g), gs)) + +function pushlazy!(p, l::Union{Plus,Times}, gs) + # TODO: handle isempty(p.xs) ? + op = l isa Plus ? plus : times + x, xs = Iterators.peel(l.xs) + i = pushlazy!(p, x, gs) + for x in xs + j = pushlazy!(p, x, gs) + i = pushop!(p, op, i, j) + end + i +end + +function pushlazy!(p, l::Minus, gs) + i = pushlazy!(p, l.p, gs) + j = pushlazy!(p, l.q, gs) + pushop!(p, minus, i, j) +end + +pushlazy!(p, l::UniMinus, gs) = pushop!(p, uniminus, pushlazy!(p, l.p, gs)) + +pushlazy!(p, l::Exp, gs) = + pushop!(p, exponentiate, pushlazy!(p, l.p, gs), intarg(l.e)) + +function pushlazy!(p, l::Decision, gs) + local k + for (x, i) in l.ps + d = pushlazy!(p, x, gs) + k = pushop!(p, decision, d, pushint!(p, i)) + end + k +end + +function pushlazy!(p, g::Getindex, gs) + length(g.is) == 1 || + throw(ArgumentError("getindex with multiple indices not supported")) + i = pushlazy!(p, g.p, gs) + j = pushlazy!(p, g.is[1], gs) + pushop!(p, getindex_, i, j) +end diff --git a/StraightLinePrograms/src/lazypolys.jl b/StraightLinePrograms/src/lazypolys.jl new file mode 100644 index 000000000000..4c1038a97f91 --- /dev/null +++ b/StraightLinePrograms/src/lazypolys.jl @@ -0,0 +1,42 @@ +## LazyPolyRing + +struct LazyPolyRing{T<:RingElement,R<:Ring} <: MPolyRing{T} + base_ring::R + + LazyPolyRing(r::Ring) = new{elem_type(r),typeof(r)}(r) +end + +base_ring(F::LazyPolyRing) = F.base_ring + + +## LazyPoly + +struct LazyPoly{T<:RingElement,PR<:MPolyRing{T}} <: MPolyElem{T} + parent::PR + p::Lazy +end + +parent(p::LazyPoly) = p.parent + +gen(R::LazyPolyRing, s::Symbol) = LazyPoly(R, Gen(s)) + +(R::LazyPolyRing)(s::Symbol) = gen(R, s) + +(R::LazyPolyRing{T})(c::T) where {T} = LazyPoly(R, Const(c)) + +function check_parent(p::LazyPoly, q::LazyPoly) + par = parent(p) + par === parent(q) || throw(ArgumentError("incompatible parents")) + par +end + +Base.show(io::IO, p::LazyPoly) = show(io, p.p) + + +### ops + ++(p::LazyPoly, q::LazyPoly) = LazyPoly(check_parent(p, q), p.p + q.p) +-(p::LazyPoly, q::LazyPoly) = LazyPoly(check_parent(p, q), p.p - q.p) +-(p::LazyPoly) = LazyPoly(parent(p), -p.p) +*(p::LazyPoly, q::LazyPoly) = LazyPoly(check_parent(p, q), p.p * q.p) +^(p::LazyPoly, e::Integer) = LazyPoly(parent(p), p.p^e) diff --git a/StraightLinePrograms/src/slpolys.jl b/StraightLinePrograms/src/slpolys.jl new file mode 100644 index 000000000000..1dc6219b65d3 --- /dev/null +++ b/StraightLinePrograms/src/slpolys.jl @@ -0,0 +1,371 @@ +## SLPolyRing (SL = straight-line) + +struct SLPolyRing{T<:RingElement,R<:Ring} <: MPolyRing{T} + base_ring::R + S::Vector{Symbol} + + SLPolyRing(r::Ring, s::Vector{Symbol}) = new{elem_type(r),typeof(r)}(r, s) +end + +SLPolyRing(r::Ring, s::Union{AbstractVector{<:AbstractString}, + AbstractVector{<:AbstractChar}}) = + SLPolyRing(r, Symbol.(s)) + +SLPolyRing(r::Ring, n::Integer) = SLPolyRing(r, [Symbol("x$i") for i=1:n]) + +# cf. mpoly.jl in Oscar +SLPolyRing(r::Ring, v::Pair{<:Union{String,Symbol}, + <:AbstractVector{<:Integer}}...) = + SLPolyRing(r, [Symbol(s, n) for (s, ns) in v for n in ns]) + +base_ring(S::SLPolyRing) = S.base_ring + +symbols(S::SLPolyRing) = S.S + +# have to constrain T <: RingElement so that this is more specific than the second +# method taking c::RingElement +(S::SLPolyRing{T})(c::T=zero(base_ring(S))) where {T<:RingElement} = S(Const(c)) + +(S::SLPolyRing{T})(c::RingElement) where {T<:RingElement} = S(Const(base_ring(S)(c))) + +function gen(S::SLPolyRing{T}, i::Integer) where {T} + s = symbols(S)[i] + S(Gen(s)) +end + +gens(S::SLPolyRing) = [S(Gen(s)) for s in symbols(S)] + +ngens(S::SLPolyRing) = length(symbols(S)) +nvars(S::SLPolyRing) = ngens(S) + +function PolynomialRing(R::Ring, s) + S = SLPolyRing(R, s) + S, gens(S) +end + +function PolynomialRing(R::Ring, v::Pair{<:Union{String,Symbol}, + <:AbstractVector{<:Integer}}...) + S = SLPolyRing(R, v...) + + # TODO: enable on Julia 1.5 (required for init keyword) + # rs = Iterators.accumulate(v; init=0:0) do x, a + # last(x)+1:last(x)+length(a[2]) + # end + rs = [] + prev = 0 + for a in v + newprev = prev+length(a[2]) + push!(rs, prev+1:newprev) + prev = newprev + end + + gs = gens(S) + S, (gs[r] for r in rs)... +end + +Base.one(S::SLPolyRing) = S(one(base_ring(S))) +Base.zero(S::SLPolyRing) = S() + +# TODO: merge this with method in AbstractAlgebra +function Base.show(io::IO, p::SLPolyRing) + max_vars = 5 + n = nvars(p) + print(io, "SLP Multivariate Polynomial Ring in ") + if n > max_vars + print(io, n) + print(io, " variables ") + end + for i = 1:min(n - 1, max_vars - 1) + print(io, string(p.S[i]), ", ") + end + if n > max_vars + print(io, "..., ") + end + print(io, string(p.S[n])) + print(io, " over ") + print(IOContext(io, :compact => true), base_ring(p)) +end + + + +## SLPoly + +struct SLPoly{T<:RingElement,SLPR<:SLPolyRing{T}} <: MPolyElem{T} + parent::SLPR + slprogram::SLProgram{T} + + SLPoly(parent, slp::SLProgram) = + new{elem_type(base_ring(parent)),typeof(parent)}(parent, slp) +end + +constants(p::SLPoly) = constants(p.slprogram) +lines(p::SLPoly) = lines(p.slprogram) + + +# create invalid poly +SLPoly(parent::SLPolyRing{T}) where {T} = SLPoly(parent, SLProgram{T}()) + +isvalid(p::SLPoly) = !hasmultireturn(p.slprogram) + +function assert_valid(p::SLPoly) + isvalid(p) || throw(ArgumentError("SLPoly is in an invalid state")) + p +end + +parent(p::SLPoly) = p.parent + +function check_parent(p::SLPoly, q::SLPoly) + p.parent === q.parent || + throw(ArgumentError("incompatible parents")) + p.parent +end + +function (S::SLPolyRing{T})(p::SLPoly{T}) where T <: RingElement + parent(p) != S && throw(ArgumentError("unable to coerce polynomial")) + p +end + +Base.zero(p::SLPoly) = zero(parent(p)) +Base.one(p::SLPoly) = one(parent(p)) + +function Base.copy!(p::SLPoly{T}, q::SLPoly{T}) where {T} + check_parent(p, q) + copy!(p.slprogram, q.slprogram) + p +end + +function Base.copy(q::SLPoly) + p = SLPoly(q.parent) + copy!(p, q) + p +end + +""" + nsteps(p::SLPoly) + +Return the number of steps ("lines") involved in the underlying +straight-line program. +""" +nsteps(p::SLPoly) = nsteps(p.slprogram) + + +## show + +function Base.show(io::IO, p::SLPoly) + io = IOContext(io, :SLPsymbols => symbols(parent(p))) + show(io, p.slprogram) +end + + +## mutating ops + +pushinit!(p::SLPoly) = pushinit!(p.slprogram) + +function pushfinalize!(p::SLPoly, i) + pushfinalize!(p.slprogram, i) + p +end + +pushop!(p::SLPoly, op::Op, i::Arg, j::Arg=Arg(0)) = + pushop!(p.slprogram, op, i, j) + +function combine!(op::Op, p::SLPoly, q::SLPoly) + combine!(op, p.slprogram, q.slprogram) + p +end + +addeq!(p::SLPoly{T}, q::SLPoly{T}) where {T} = combine!(plus, p, q) + +subeq!(p::SLPoly{T}, q::SLPoly{T}) where {T} = combine!(minus, p, q) + +function subeq!(p::SLPoly) + combine!(uniminus, p.slprogram) + p +end + +muleq!(p::SLPoly{T}, q::SLPoly{T}) where {T} = combine!(times, p, q) + +function expeq!(p::SLPoly, e::Integer) + combine!(exponentiate, p.slprogram, e) + p +end + +function permutegens!(p::SLPoly, perm) + permute_inputs!(p.slprogram, perm, + perm isa Union{AbstractArray,AbstractAlgebra.AbstractPerm}) + p +end + + +## unary/binary ops + ++(p::SLPoly{T}, q::SLPoly{T}) where {T} = addeq!(copy(p), q) + +*(p::SLPoly{T}, q::SLPoly{T}) where {T} = muleq!(copy(p), q) + +-(p::SLPoly{T}, q::SLPoly{T}) where {T} = subeq!(copy(p), q) + +-(p::SLPoly) = subeq!(copy(p)) + +^(p::SLPoly, e::Integer) = expeq!(copy(p), e) + +# should be AbstractPerm instead of GroupElem, but we need to support GAP's +# permutations as provided in Oscar +^(p::SLPoly, perm::AbstractAlgebra.GroupElem) = permutegens!(copy(p), perm) + + +## adhoc ops + ++(p::SLPoly{T}, q::T) where {T<:RingElem} = p + parent(p)(q) ++(q::T, p::SLPoly{T}) where {T<:RingElem} = parent(p)(q) + p + +-(p::SLPoly{T}, q::T) where {T<:RingElem} = p - parent(p)(q) +-(q::T, p::SLPoly{T}) where {T<:RingElem} = parent(p)(q) - p + +*(p::SLPoly{T}, q::T) where {T<:RingElem} = p * parent(p)(q) +*(q::T, p::SLPoly{T}) where {T<:RingElem} = parent(p)(q) * p + + +## comparison + +# TODO: two SLPolys migth be considered equal if they "canonical" form +# would be equal (e.g. their conversions to MPoly) +function Base.:(==)(p::SLPoly{T}, q::SLPoly{T}) where {T} + check_parent(p, q) + p.slprogram == q.slprogram +end + + +## evaluate + +evaluate(p::SLPoly{T}, xs::Vector{S}, conv::F=identity + ) where {T<:RingElement,S<:RingElement,F} = + evaluate(p.slprogram, xs, conv) + +function evaluate!(res::Vector{S}, p::SLPoly{T}, xs::Vector{S}, + conv::F=identity + ) where {S,T,F} + evaluate!(res, p.slprogram, xs, conv) +end + + +## conversion Lazy -> SLP + +(R::SLPolyRing{T})(p::LazyPoly{T}) where {T} = R(p.p) + +# TODO: remove this method (this is an ambiguity fix) +(R::SLPolyRing{T})(p::LazyPoly{T}) where {T<:RingElement} = R(p.p) + +function (R::SLPolyRing{T})(p::Lazy) where {T} + pr = compile(SLProgram{T}, p, symbols(R)) + SLPoly(R, pr) +end + + +## conversion MPoly -> SLPoly + +function Base.convert(R::SLPolyRing, p::Generic.MPoly; limit_exp::Bool=false) + # TODO: currently handles only default ordering + symbols(R) == symbols(parent(p)) || + throw(ArgumentError("incompatible symbols")) + q = SLPoly(R) + @assert lastindex(p.coeffs) < cstmark + qcs = constants(q) + @assert isempty(qcs) + # have to use p.length, as p.coeffs and p.exps + # might contain trailing gargabe + resize!(qcs, p.length) + copyto!(qcs, 1, p.coeffs, 1, p.length) + exps = UInt64[] + monoms = [Pair{UInt64,Arg}[] for _ in axes(p.exps, 1)] + for v in reverse(axes(p.exps, 1)) + copy!(exps, view(p.exps, v, 1:p.length)) + unique!(sort!(exps)) + if !isempty(exps) && first(exps) == 0 + popfirst!(exps) + end + isempty(exps) && continue + xref = input(size(p.exps, 1) + 1 - v) + if limit_exp # experimental + # TODO: move to a general SLP optimization pass? + # and check in which case it's an improvement + + # 1) find all exponents which will be computed, i.e. all e÷2 + # for all e in exps, recursively + n = length(exps) - 1 + while length(exps) > n + n = length(exps) + for i in eachindex(exps) + exps[i] == 1 && continue + e0 = exps[i] >> 1 + push!(exps, e0) + end + unique!(sort!(exps)) + end + + # 2) for all e from 1), compute x^2 and store the result in monoms + for e in exps + e1 = e >> 1 + if e == 1 + k = xref + elseif monoms[v][end][1] == 2*e1 + @assert e == 2*e1+1 + k = pushop!(q, times, monoms[v][end][2], xref) + else + m1 = searchsortedfirst(monoms[v], e1, by=first) + k1 = monoms[v][m1][2] + k = pushop!(q, times, k1, k1) + if e1+e1 != e + @assert e1+e1+1 == e + k = pushop!(q, times, k, xref) + end + end + push!(monoms[v], e => k) + end + else + for e in exps + e == 0 && continue + k = pushop!(q, exponentiate, xref, Arg(e)) + push!(monoms[v], e => k) + end + end + end + + k = Arg(0) # TODO: don't use 0 + for t in eachindex(qcs) + i = asconstant(t) + j = 0 + for v in reverse(axes(p.exps, 1)) + e = p.exps[v, t] + if e != 0 + j = monoms[v][searchsortedfirst(monoms[v], e, by=first)][2] + i = pushop!(q, times, i, j) + end + end + if k == Arg(0) + k = i + else + k = pushop!(q, plus, k, i) + end + end + if isempty(qcs) + k = pushconst!(q.slprogram, base_ring(R)()) + end + pushfinalize!(q, k) + q +end + + +## conversion SLPoly -> MPoly + +function Base.convert(R::MPolyRing, p::SLPoly) + symbols(R) == symbols(parent(p)) || + throw(ArgumentError("incompatible symbols")) + assert_valid(p) + evaluate(p, gens(R)) +end + + +## compile! + +compile!(p::SLPoly) = compile!(p.slprogram) diff --git a/StraightLinePrograms/src/straightline.jl b/StraightLinePrograms/src/straightline.jl new file mode 100644 index 000000000000..34dc60389c4e --- /dev/null +++ b/StraightLinePrograms/src/straightline.jl @@ -0,0 +1,839 @@ +## constants + +const opmask = 0xff00000000000000 +const argmask = 0x000000000fffffff +const inputmark = 0x0000000008000000 +const cstmark = 0x0000000004000000 +const intmark = inputmark | cstmark +const typemark = intmark +const payloadmask = 0x0000000003ffffff +const negbit = 0x0000000002000000 +const argshift = 28 # argmask == 2^argshift - 1 + +@assert payloadmask == argmask ⊻ inputmark ⊻ cstmark + +const dectrue = argmask + + +## Op, Line, Arg + +struct Op + x::UInt64 +end + +struct Line + x::UInt64 +end + +struct Arg + x::UInt64 + + function Arg(x::UInt64) + iszero(x & ~argmask ) || + throw(ArgumentError("argument too big")) + new(x) + end +end + +Arg(x::Integer) = Arg(UInt64(x)) + + +## predicates + +const showop = Dict{Op,String}() + +for (i, (op, unary, showchar)) in enumerate([(:assign , false, "->"), + (:uniminus , true, "-"), + (:plus , false, "+"), + (:minus , false, "-"), + (:times , false, "*"), + (:divide , false, "/"), + (:exponentiate , true, "^"), + (:keep , true, "keep"), + (:decision , false, "&"), + (:getindex_ , false, "[]"), + ]) + isop = Symbol(:is, op) + c = UInt64(i) << (2*argshift) + if unary + c |= 0x8000000000000000 + end + @eval begin + const $op = Op($c) + $isop(op::Op) = op === $op + end + showop[Op(c)] = showchar +end + +isquasiunary(op) = (op.x & 0x8000000000000000) != 0 +isunary(op) = isquasiunary(op) & (op != exponentiate) + + +## raw manips + +Line(op::Op, i::Arg, j::Arg) = Line(op.x | i.x << argshift | j.x) + +pack(op::Op, i::Arg, j::Arg) = Line(op, i, j) + +function unpack(line::Line) + line = line.x + op = opmask & line + j = line & argmask + i = (line >> argshift) & argmask + Op(op), Arg(i), Arg(j) +end + + +## show + +function Base.show(io::IO, l::Line) + op, i, j = unpack(l) + print(io, op, " : ", i, " , ", j) +end + +Base.show(io::IO, op::Op) = print(io, showop[op]) + +function Base.show(io::IO, x::Arg) + if isint(x) + print(io, 'i', intidx(x)) + elseif isinput(x) + print(io, '$', inputidx(x)) + elseif isconstant(x) + print(io, '+', constantidx(x)) + else + print(io, ' ', x.x) + end +end + + +## SLProgram + +mutable struct SLProgram{T} <: AbstractSLProgram + cs::Vector{T} # constants + lines::Vector{Line} # instructions + int::Int # number of stored Int at the beginning of lines + len::Int # length of result vector + # (where to store next result for any new line, minus 1) + ret::Arg # result index to return (0 for all) + f::Ref{Function} # compiled execution +end + +SLProgram{T}() where {T} = SLProgram(T[], Line[], 0, 0, Arg(0), Ref{Function}()) +SLProgram() = SLProgram{Union{}}() + +# return an input +function SLProgram{T}(i::Integer) where {T} + p = SLProgram{T}() + pushfinalize!(p, input(i)) +end + +SLProgram(i::Integer) = SLProgram{Union{}}(i) + +function Base.copy!(q::SLProgram, p::SLProgram) + copy!(q.cs, p.cs) + copy!(q.lines, p.lines) + q.int = p.int + q.len = p.len + q.ret = p.ret + q +end + +copy_oftype(p::SLProgram, ::Type{T}) where {T} = copy!(SLProgram{T}(), p) + +Base.copy(p::SLProgram{T}) where {T} = copy_oftype(p, T) + +Base.:(==)(p::SLProgram, q::SLProgram) = + p.cs == q.cs && p.lines == q.lines && p.int == q.int && p.ret == q.ret + +constants(p::SLProgram) = p.cs +_integers(p::SLProgram) = @view p.lines[1:p.int] +integers(p::SLProgram) = @view p.lines[p.int:-1:1] +lines(p::SLProgram) = @view p.lines[1+p.int : end] +linesindices(p::SLProgram) = 1+p.int:lastindex(p.lines) +nsteps(p::SLProgram) = length(p.lines) - p.int + +constantstype(p::SLProgram{T}) where {T} = T + +# return the (max) number of inputs +function ninputs(p::SLProgram) + m0 = isinput(p.ret) ? inputidx(p.ret) % Int : 0 + m1 = mapreduce(max, lines(p); init=0) do line + op, i, j = unpack(line) + max(isinput(i) ? inputidx(i) : 0, + isinput(j) ? inputidx(j) : 0) % Int + end + max(m0, m1) +end + +slpgen(n::Integer) = SLProgram(n) +slpgens(n::Integer) = [SLProgram(i) for i=1:n] + +gens(::Type{SLProgram}, n::Integer) = slpgens(n) +gens(::Type{SLProgram{T}}, n::Integer) where {T} = [SLProgram{T}(i) for i=1:n] + +slpcst(c) = SLProgram(Const(c)) + + +## show + +# old basic version, useful when normal show is broken +function showsimple(io::IO, p::SLProgram) + println("SLProgram:") + if !isempty(constants(p)) + println("with constants:") + for (i, c) in enumerate(constants(p)) + println(io, i, " | ", c) + end + end + if !isempty(integers(p)) + println("with integers:") + for (i, c) in enumerate(integers(p)) + println(io, i, " | ", c.x % Int) + end + end + if !isempty(lines(p)) + println("with lines:") + for (i, c) in enumerate(lines(p)) + println(io, i, " | ", c) + end + end + println("return: ", p.ret == Arg(dectrue) ? "true" : p.ret) +end + +showsimple(p::SLProgram) = showsimple(stdout, p) + +slpsyms(n::Integer) = + n <= 3 ? + (:x, :y, :z)[1:n] : + (Symbol(:x, i) for i = 1:n) + +function Base.show(io::IO, p::SLProgram) + gs = get(io, :SLPsymbols, slpsyms(ninputs(p))) + show(io, evaluate(p, lazygens(gs), Const)) +end + +function Base.show(io::IO, ::MIME"text/plain", p::SLProgram{T}) where T + n = length(lines(p)) + gs = get(io, :SLPsymbols, slpsyms(ninputs(p))) + syms = lazygens(gs) + reslazy = Any[] + if n == 0 && !hasmultireturn(p) + # trivial program, show only result + return show(io, retrieve(integers(p), constants(p), syms, reslazy, p.ret)) + end + + str(x...) = sprint(print, x...; context=io) + strlines = Vector{String}[] + widths = [0, 0, 0, 0, 0] + + ptmp = SLProgram{T}() + copy!(ptmp.cs, p.cs) + copy!(ptmp.lines, _integers(p)) + ptmp.int = p.int + + cs = constants(p) + ints = integers(p) + + showarg(x) = string(retrieve(ints, cs, syms, + ['#'*string(i) for i in eachindex(lines(p))], + x)) + + for line in lines(p) + op, i, j = unpack(line) + k = updatelen!(ptmp, op, i, j) + push!(ptmp.lines, line) + pushfinalize!(ptmp, Arg(k)) + + line = String[] + push!(strlines, line) + if iskeep(op) + push!(line, "keep:") + if i.x == 0 + push!(line, " nothing") + elseif i.x == 1 + push!(line, " #1") + else + push!(line, " #1..#$(i.x)") + end + continue + end + + x = showarg(i) + y = isunary(op) || isassign(op) ? "" : + isexponentiate(op) ? string(getint(j)) : + isquasiunary(op) ? string(j.x) : + showarg(j) + + if isdecision(op) + push!(line, "test:") + push!(line, " order($x) == $y || return false") + continue + end + + push!(line, str('#', string(k), " =")) + # TODO: add a way to evaluate step by step instead of re-evaluating + # from the beginning each time + plazy = evaluate!(reslazy, ptmp, syms, Const) + + strop = isassign(op) ? "" : showop[op] + push!(line, str(strop), x, y, str(plazy)) + widths .= max.(widths, textwidth.(line)) + end + for (k, line) in enumerate(strlines) + k == 1 || println(io) + if line[1] ∈ ["keep:", "test:"] + join(io, line) + continue + end + print(io, ' '^(widths[1]-length(line[1])), line[1]) + for i = 2:4 + print(io, ' '^(1+widths[i]-textwidth(line[i])), line[i]) + end + print(io, " ==> ", line[5]) + end + + print(io, "\nreturn: ") + # TODO: ^--- remove when unnecessary? + if hasmultireturn(p) + print(io, '[') + join(io, + (showarg(Arg(i)) for i in 1:p.len), + ", ") + print(io, ']') + elseif hasdecision(p) + print(io, true) + else + print(io, showarg(p.ret)) + end + +end + + +## building SLProgram + +isregister(i::Arg) = typemark & i.x == 0 +registeridx(i::Arg) = i.x +asregister(i::Integer) = Arg(UInt64(i)) # TODO: sanitize input + +# return #ref for i-th input +function input(i::Integer) + i = Int(i) + @assert i < cstmark + Arg((i % UInt64) | inputmark) +end + +isinput(i::Arg) = typemark & i.x === inputmark +inputidx(i::Arg) = i.x ⊻ inputmark + +isconstant(i::Arg) = typemark & i.x === cstmark +constantidx(i::Arg) = i.x ⊻ cstmark + +asconstant(i::Integer) = Arg(UInt64(i) | cstmark) + +isint(i::Arg) = typemark & i.x === intmark +asint(i::Int) = Arg(i % UInt64 | intmark) +intidx(i::Arg) = i.x ⊻ intmark + +hasmultireturn(p::SLProgram) = p.ret.x == 0 +setmultireturn!(p::SLProgram) = (p.ret = Arg(0); p) + +hasdecision(p::SLProgram) = p.ret.x == argmask +setdecision!(p::SLProgram) = (p.ret = Arg(argmask); p) + + +# make an Arg out of x, considered as "computational" integer (e.g. not an index) +# call to create an integer exponent +function intarg(x::Integer) + y = Int(x) % UInt64 + if x >= 0 + y > (payloadmask ⊻ negbit) && + throw(ArgumentError("positive integer argument too big")) + Arg(y) + else + y | (payloadmask >> 1) == typemax(UInt64) || # >> 1 for signbit + throw(ArgumentError("negative integer argument too small")) + Arg(y & payloadmask) + end +end + +# retrieve a "computational" integer from an Arg +function getint(x::Arg) + y = x.x + iszero(y & ~payloadmask) || + throw(ArgumentError("Arg does not contain an Int")) + if iszero(y & negbit) + y % Int + else + (y | ~payloadmask) % Int + end +end + +# call before mutating, unless p is empty (opposite of pushfinalize!) + +# TODO: pushinit! & pushfinalize! not really nececessary anymore? +pushinit!(p::SLProgram) = p.ret + +function pushconst!(p::SLProgram{T}, c::T) where T + push!(constants(p), c) + l = lastindex(constants(p)) + @assert l < cstmark + asconstant(l) +end + +function pushint!(p::SLProgram, i::Integer) + pushfirst!(p.lines, Line(Int(i) % UInt64)) + p.int += 1 + @assert p.int < cstmark + asint(p.int) +end + +function updatelen!(p, op, i, j) + if isassign(op) && j != Arg(0) + ptr = Int(j.x) + if ptr == p.len + 1 + p.len += 1 + end + 1 <= ptr <= p.len || + throw(ArgumentError("invalid `assign` destination")) + elseif iskeep(op) + ptr = Int(i.x) + ptr <= p.len || throw(ArgumentError("cannot `keep` so many items")) + p.len = ptr + elseif isdecision(op) + ptr = argmask # cf. hasdecision + else + ptr = p.len += 1 + @assert ptr < cstmark + end + ptr +end + +function pushop!(p::SLProgram, op::Op, i::Arg, j::Arg=Arg(0)) + @assert i.x <= argmask && j.x <= argmask + ptr = updatelen!(p, op, i, j) + push!(p.lines, Line(op, i, j)) + Arg(ptr % UInt64) +end + +function pushfinalize!(p::SLProgram, ret::Arg) + p.ret = ret + p +end + +function _combine!(p::SLProgram, q::SLProgram; makeinput=identity) + i1 = pushinit!(p) + i2 = pushinit!(q) + koffset = length(constants(p)) + ioffset = p.int + len = length(lines(p)) + loffset = len - count(lines(p)) do line + op, i, j = unpack(line) + isdecision(op) + # TODO: probably needs also to account for "keep" lines + end + p.int += q.int # must be *after* computing len + prepend!(p.lines, _integers(q)) + append!(p.lines, lines(q)) + append!(constants(p), constants(q)) + + @assert length(constants(p)) < cstmark # TODO: should not be @assert + + for n = len+1:lastindex(lines(p)) + op, i, j = unpack(lines(p)[n]) + if isconstant(i) + i = Arg(i.x + koffset) + elseif isinput(i) + i = makeinput(i) + elseif isint(i) + i = Arg(i.x + ioffset) + else + i = Arg(i.x + loffset) + end + if isconstant(j) + j = Arg(j.x + koffset) + elseif isinput(j) + j = makeinput(j) + elseif isint(j) + j = Arg(j.x + ioffset) + elseif !isquasiunary(op) + # TODO: normalize assign so that j.x is never 0 + if !(isassign(op) && j.x == 0) + j = Arg(j.x + loffset) + end + end + lines(p)[n] = Line(op, i, j) + # TODO: write conditionally only when modifications + end + if isconstant(i2) + i2 = Arg(i2.x + koffset) + elseif isinput(i2) + i2 = makeinput(i2) + elseif hasdecision(q) + elseif hasmultireturn(q) + elseif isint(i2) + i2 = Arg(i2.x + ioffset) + else + i2 = Arg(i2.x + loffset) + end + p.len += q.len + i1, i2 +end + +function combine!(op::Op, p::SLProgram, q::SLProgram) + i = pushop!(p, op, _combine!(p, q)...) + pushfinalize!(p, i) +end + +function combine!(op::Op, p::SLProgram) + i = pushinit!(p) + i = pushop!(p, op, i) + pushfinalize!(p, i) +end + +function combine!(op::Op, p::SLProgram, e::Integer) + i = pushinit!(p) + i = pushop!(p, op, i, intarg(e)) + pushfinalize!(p, i) +end + +applyperm(perm, i, usegetindex) = usegetindex ? perm[i] : perm(i) + +function permute_inputs!(p::SLProgram, perm, usegetindex=true) + for l in linesindices(p) + op, i, j = unpack(p.lines[l]) + changed = false + if isinput(i) + i = input(applyperm(perm, inputidx(i) % Int, usegetindex)) + changed = true + end + if isinput(j) + j = input(applyperm(perm, inputidx(j) % Int, usegetindex)) + changed = true + end + if changed + p.lines[l] = pack(op, i, j) + end + end + p +end + + +### adhoc + +function combine(op::Op, p::SLProgram, q) + r = copy_oftype(p, typejoin(constantstype(p), typeof(q))) + i = pushinit!(r) + j = pushop!(r, op, i, pushconst!(r, q)) + pushfinalize!(r, j) +end + +function combine(op::Op, p, q::SLProgram) + r = copy_oftype(q, typejoin(constantstype(q), typeof(p))) + i = pushinit!(r) + j = pushop!(r, op, pushconst!(r, p), i) + pushfinalize!(r, j) +end + + +## mutating ops + +addeq!(p::SLProgram, q::SLProgram) = combine!(plus, p, q) + +subeq!(p::SLProgram, q::SLProgram) = combine!(minus, p, q) + +function subeq!(p::SLProgram) + combine!(uniminus, p) + p +end + +muleq!(p::SLProgram, q::SLProgram) = combine!(times, p, q) + +function expeq!(p::SLProgram, e::Integer) + combine!(exponentiate, p, e) + p +end + +function testeq!(p::SLProgram, q::SLProgram) + hasdecision(p) && hasdecision(q) || throw(ArgumentError( + "cannot &-combine two programs which are not decisions")) + _combine!(p, q) + setdecision!(p) +end + + +## unary/binary ops + +copy_jointype(p::SLProgram, q::SLProgram) = + copy_oftype(p, typejoin(constantstype(p), constantstype(q))) + ++(p::SLProgram, q::SLProgram) = addeq!(copy_jointype(p, q), q) + +*(p::SLProgram, q::SLProgram) = muleq!(copy_jointype(p, q), q) + +-(p::SLProgram, q::SLProgram) = subeq!(copy_jointype(p, q), q) + +-(p::SLProgram) = subeq!(copy(p)) + +^(p::SLProgram, e::Integer) = expeq!(copy(p), e) + +Base.literal_pow(::typeof(^), p::SLProgram, ::Val{e}) where {e} = p^e + +Base.:&(p::SLProgram, q::SLProgram) = testeq!(copy_jointype(p, q), q) + +function test!(p::SLProgram, x::Integer) + hasdecision(p) && throw(ArgumentError("SLProgram is already a decision")) + i = pushinit!(p) + j = pushint!(p, x) + pushfinalize!(p, pushop!(p, decision, i, j)) +end + +test(p::SLProgram, x::Integer) = test!(copy(p), x) + +function list(::Type{SL}, ps) where {SL<:SLProgram} + idx = Arg[] + q = SL === SLProgram ? SLProgram{Any}() : SL() + # TODO: ^^^--- use typejoin or smthg instead of Any + for (n, p) in enumerate(ps) + _, i = _combine!(q, p) + if !isregister(i) + # we write explicitly in res-vector, so that the indices increase + # strictly and it's always safe in the final loop to copy from + # this `i` into final destination + i = pushop!(q, assign, i) + @assert isregister(i) + end + push!(idx, i) + @assert nsteps(q) >= registeridx(i) + end + for (n, i) in enumerate(idx) + if registeridx(i) > n + pushop!(q, assign, i, Arg(n)) + else + @assert registeridx(i) == n + end + end + if nsteps(q) > length(idx) + pushop!(q, keep, Arg(length(idx))) + end + q +end + +function composewith!(q::SLProgram, p::SLProgram) + hasmultireturn(q) || throw(ArgumentError( + "second argument of `compose` must return a list")) + ninp = q.len # max ninputs that p can use + _, j = _combine!(q, p, makeinput = function (i::Arg) + idx = inputidx(i) + idx <= ninp || throw(ArgumentError( + "inner compose argument does not provide enough input")) + asregister(idx) + end) + if j.x == 0 && p.len != 0 # multireturn + newlen = q.len - ninp + for i = 1:newlen + pushop!(q, assign, asregister(i+ninp), asregister(i)) + end + pushop!(q, keep, asregister(newlen)) + end + pushfinalize!(q, j) +end + + +compose(p::SLProgram, q::SLProgram) = composewith!(copy_jointype(q, p), p) + +getindex!(p::SLProgram, q::SLProgram) = combine!(getindex_, p, q) +Base.getindex(p::SLProgram, q::SLProgram) = getindex!(copy_jointype(p, q), q) + + +### adhoc + ++(p::SLProgram, q) = combine(plus, p, q) ++(p, q::SLProgram) = combine(plus, p, q) +*(p::SLProgram, q) = combine(times, p, q) +*(p, q::SLProgram) = combine(times, p, q) +-(p::SLProgram, q) = combine(minus, p, q) +-(p, q::SLProgram) = combine(minus, p, q) + +function getindex!(p::SLProgram, i::Integer) + if hasmultireturn(p) + pushfinalize!(p, asregister(i)) + else + a = pushinit!(p) + k = pushop!(p, getindex_, a, pushint!(p, i)) + pushfinalize!(p, k) + end +end + +Base.getindex(p::SLProgram, i::Integer) = getindex!(copy(p), i) + +## conversion SLProgram -> Lazy + +lazygens(gs) = Any[Gen(s) for s in gs] +lazygens(n::Integer) = lazygens(slpsyms(n)) + +# strictly convenience function +aslazy(p::SLProgram, gs=lazygens(ninputs(p))) = Lazy(evaluate(p, gs)) + + +## evaluate + +getres(p::SLProgram, xs) = Vector{eltype(xs)}(undef, length(p.lines)) +# length(p.lines) might be an overestimate in some cases, but +# this generally limits the number of allocations + +function evaluate(p::SLProgram{T}, xs::Vector{S}, conv::F=identity + ) where {T,S,F} + if isassigned(p.f) + p.f[](xs)::S + else + evaluate!(getres(p, xs), p, xs, conv) + end +end + +function evaluates(p::SLProgram{T}, xs::Vector{S}, conv::F=identity + ) where {T,S,F} + res = getres(p, xs) + evaluate!(res, p, xs, conv) + res +end + +retrieve(ints, cs, xs, res, i, conv::F=identity) where {F} = + isint(i) ? ints[intidx(i)].x % Int : + isconstant(i) ? conv(cs[constantidx(i)]) : + isinput(i) ? xs[inputidx(i)] : + res[i.x] + +function evaluate!(res::Vector{S}, p::SLProgram{T}, xs::Vector{S}, + conv::F=identity) where {S,T,F} + # TODO: handle isempty(lines(p)) + empty!(res) + + cs = constants(p) + ints = integers(p) + + local decide::Union{S,Bool} + + for line in lines(p) + local r::S + op, i, j = unpack(line) + if iskeep(op) + @assert isregister(i) + resize!(res, i.x) + continue + end + x = retrieve(ints, cs, xs, res, i, conv) + if isexponentiate(op) + r = x^getint(j) # TODO: support bigger j + elseif isassign(op) + dst = j.x % Int + if dst != 0 && dst != lastindex(res) + 1 + res[dst] = x + continue + end + r = x + elseif isuniminus(op) + r = -x + else + y = retrieve(ints, cs, xs, res, j, conv) + if isplus(op) + r = x + y + elseif isminus(op) + r = x - y + elseif istimes(op) + r = x * y + elseif isdivide(op) + r = divexact(x, y) + elseif isgetindex_(op) + r = x[y] + elseif isdecision(op) + t = test(x, y) + t === false && return false + if @isdefined(decide) + decide &= t + else + decide = t + end + continue + else + throw(ArgumentError("unknown operation")) + end + end + push!(res, r) + end + + @assert length(res) == p.len + if hasdecision(p) + decide + elseif hasmultireturn(p) + list(res) + else + retrieve(ints, cs, xs, res, p.ret, conv) + end +end + +test(x, o) = order(x) == o + +list(res) = list(eltype(res), res) + +list(::Type{T}, res) where {T} = res + + +## compile! + +cretrieve(i) = + isinput(i) ? Symbol(:x, inputidx(i)) => inputidx(i) : + isconstant(i) ? Symbol(:c, constantidx(i)) => 0 : + Symbol(:res, i.x) => 0 + +# TODO: handle the "conv" argument like in evaluate! +# (works without, but there can be type-instability) + +# return compiled execution function f, and updates +# p.f[] = f, which is not invalidated when p is mutated +function compile!(p::SLProgram) + res = Expr[] + fn = :(function (xs::Vector) + end) + k = 0 + cs = constants(p) + for k in eachindex(cs) + push!(res, :($(Symbol(:c, k)) = @inbounds $cs[$k])) + end + mininput = 0 + for line in lines(p) + k += 1 + rk = Symbol(:res, k) + op, i, j = unpack(line) + x, idx = cretrieve(i) + mininput = max(mininput, idx) + line = + if isexponentiate(op) + :($rk = $x^$(getint(j))) + elseif isassign(op) + :($rk = $x) + elseif isuniminus(op) + :($rk = -$x) + else + y, idx = cretrieve(j) + mininput = max(mininput, idx) + if isplus(op) + :($rk = $x + $y) + elseif isminus(op) + :($rk = $x - $y) + elseif istimes(op) + :($rk = $x * $y) + elseif isdivide(op) + :($rk = divexact($x, $y)) + end + end + push!(res, line) + end + for k = 1:mininput-1 + pushfirst!(res, :($(Symbol(:x, k)) = @inbounds xs[$k])) + end + if mininput >= 1 + pushfirst!(res, :($(Symbol(:x, mininput)) = xs[$mininput])) + end + append!(fn.args[2].args, res) + p.f[] = eval(fn) +end diff --git a/StraightLinePrograms/test/atlas.jl b/StraightLinePrograms/test/atlas.jl new file mode 100644 index 000000000000..046fb6fe84cf --- /dev/null +++ b/StraightLinePrograms/test/atlas.jl @@ -0,0 +1,182 @@ +@testset "AtlasSLProgram" begin + for s0 in ["", "\n# comment\n \n # comment \n\n # comment", + """ # \n echo "1 2 4" \necho a echo becho""" ] + p0 = AtlasSLProgram(s0) + @test p0.code == s0 + @test p0.ngens == 2 + @test isempty(p0.lines) + @test p0.outputs == 1:2 + @test endswith(p0.code, "echo") == !isempty(p0.echo) + end + + s1 = " inp 3 3 1 2 \n oup 2 2 2" + p1 = AtlasSLProgram(s1) + @test p1.code == s1 + @test p1.ngens == 3 + @test isempty(p1.lines) + @test p1.outputs == [3, 3] + + p2 = AtlasSLProgram("inp 2 \n inp 1 4 \n oup 3 4 1 2") + @test p2.ngens == 3 + @test p2.outputs == [3, 1, 2] + + p3 = AtlasSLProgram("oup 1 2") + @test p3.ngens == 2 + @test p3.outputs == [2] + + @test_throws ArgumentError AtlasSLProgram("inp 1 2") # implicit "oup 1 2" + + p4 = AtlasSLProgram("inp 1 2 \n oup 2 2 2") + @test p4.ngens == 1 + @test p4.outputs == [1, 1] + + for bad in ["inp 1 2", "oup 3", "oup 2 \n inp 2", + "inp 2 \n inp 1", "oup 2 \n oup 1", + "inp 2 1 2 3", "oup 2 1", + "inp 3 \n oup 4", "oup 2 2 3", + "mu 1 2 3 \n inp 3", "oup 2 \n mu 1 2 3"] + @test_throws ArgumentError AtlasSLProgram(bad) + end + + for bad in ["cjr 1", "cjr 1 2 3", + "cj 1 2", "cj 1 2 3 4", + "com 1 2", "com 1 2 3 4", + "mu 1 2", "mu 1 2 3 4", + "pwr 1 2", "pwr 1 2 3 4", + "iv 1", "iv 1 2 3", + "cp 1", "cp 1 2 3"] + @test_throws ArgumentError AtlasSLProgram(bad) + end + q1 = AtlasSLProgram(""" + cjr 1 2 + cj 1 2 3 + com 2 3 4 + mu 3 4 5 + pwr 4 5 6 + iv 6 7 + cp 7 8 + oup 3 6 7 8 + """) + @test q1.ngens == 2 + @test q1.outputs == [6, 7, 8] + @test q1.lines == [AtlasLine(:cj, 1, 1, 2), + AtlasLine(:cj, 3, 1, 2), + AtlasLine(:com, 4, 2, 3), + AtlasLine(:mu, 5, 3, 4), + AtlasLine(:pwr, 6, 4, 5), + AtlasLine(:iv, 7, 6), + AtlasLine(:cp, 8, 7) + ] + @test_throws ArgumentError AtlasSLProgram("unknown 1 2") + @test_throws ArgumentError AtlasSLProgram("chor 1 2") + + d1 = AtlasSLDecision("echo something\ninp 2\nchor 1 2\nchor 2 3\nmu 1 2 3\nchor 3 5") + @test d1.lines == [AtlasLine(:chor, 0, 1, 2), + AtlasLine(:chor, 0, 2, 3), + AtlasLine(:mu, 3, 1, 2), + AtlasLine(:chor, 0, 3, 5)] + @test d1.ngens == 2 + @test d1.echo == ["something"] + @test_throws ArgumentError AtlasSLDecision("chor 1 2\noup 1") +end + +@testset "AtlasSLProgram evaluate" begin + ab = Lazy[Gen(:a), Gen(:b)] + a, b = ab + res = empty(ab) + + p = AtlasSLProgram("cjr 2 1 \noup 1 2") + @test evaluate(p, ab) == [a^-1 * b * a] + @test res === evaluate!(res, p, ab) == [a^-1 * b * a] + g = SL.compile(GAPSLProgram, p) + @test g isa GAPSLProgram + @test evaluate(g, ab) == [a^-1 * b * a] + sl = SL.compile(SLProgram, p) + @test sl isa SLProgram + @test evaluate(sl, ab) == [a^-1 * b * a] + sl = SL.compile(p) + @test sl isa SLProgram + @test evaluate(sl, ab) == [a^-1 * b * a] + + p = AtlasSLProgram("cj 1 2 3\noup 2 1 3") + @test evaluate(p, ab) == [a, b^-1 * a * b] + @test res === evaluate!(res, p, ab) == [a, b^-1 * a * b] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a, b^-1 * a * b] + + p = AtlasSLProgram("com 1 2 3 \noup 1 3") + @test evaluate(p, ab) == [a^-1 * b^-1 * a * b] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a^-1 * b^-1 * a * b] + + p = AtlasSLProgram("iv 1 2") + @test evaluate(p, ab) == [a, a^-1] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a, a^-1] + + p = AtlasSLProgram("mu 1 2 1") + @test evaluate(p, ab) == [a * b, b] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a * b, b] + + p = AtlasSLProgram("inp 1 \n pwr 4 1 2") + @test evaluate(p, ab) == [a, a^4] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a, a^4] + + p = AtlasSLProgram("cp 1 3 \n cp 2 4 \n oup 4") + @test evaluate(p, ab) == [a, b, a, b] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [a, b, a, b] + + @test_throws ArgumentError SL.compile(GAPSLDecision, p) + + # bug with overwriting results + p = AtlasSLProgram("oup 2 2 1") + @test evaluate(p, ab) == [b, a] + g = SL.compile(GAPSLProgram, p) + @test evaluate(g, ab) == [b, a] +end + +@testset "AtlasSLDecision evaluate /compile" begin + p = perm"(1, 2)(3, 4)" + q = Perm([1, 4, 2, 3]) + pq = [p, q] + + # same test as in tests for GAPSLDecision + for (i, j, k, r) in [(2, 3, 5, false), + (2, 3, 3, true), + (2, 2, 3, false), + (1, 3, 3, false)] + d = AtlasSLDecision( + """ + mu 1 2 3 + chor 1 $i + chor 2 $j + chor 3 $k + """) + + @test evaluate(d, pq) == r + @test_throws ArgumentError SL.compile(GAPSLProgram, d) + end + + for (p, q, r) in eachcol(rand(parent(p), 3, 50)) + pqr = [p, q, r] + for (i, j, k) in eachcol(rand(1:6, 3, 10)) + ad = AtlasSLDecision( + """ + inp 3 + chor 1 $i + chor 2 $j + chor 3 $k + """) + gd = GAPSLDecision([ + ["Order", 1, i], + ["Order", 2, j], + ["Order", 3, k]], + 3) + @test evaluate(ad, pqr) == evaluate(gd, pqr) + @test SL.compile(GAPSLDecision, ad).lines == gd.lines + end + end +end diff --git a/StraightLinePrograms/test/gap.jl b/StraightLinePrograms/test/gap.jl new file mode 100644 index 000000000000..0cf46ca858a9 --- /dev/null +++ b/StraightLinePrograms/test/gap.jl @@ -0,0 +1,189 @@ +@testset "GAPSLProgram construction" begin + slines = GAPStraightLine[] + + pushline!(slines, [1, 2]) + @test slines[1] == [1, 2] + @test_throws ArgumentError pushline!(slines, [1, 2, 3]) + + pushline!(slines, ([2, 3], 1)) + @test slines[2] == ([2, 3], 1) + @test_throws MethodError pushline!(slines, ([2, 3], 2, 3)) + + pushline!(slines, [[4, 3], 2]) + @test slines[3] == ([4, 3], 2) + @test_throws ArgumentError pushline!(slines, [[2, 3], 2, 3]) + @test_throws ArgumentError pushline!(slines, [[2, 3, 4], 2, 3]) + + pushline!(slines, [[4, 3], [1, 2, 3, 4]]) + @test slines[4] == [[4, 3], [1, 2, 3, 4]] + @test_throws ArgumentError pushline!(slines, [[2, 3], [2, 3, 4]]) + + # return list not at the end + @test_throws ArgumentError pushline!(slines, [1, 2]) + @test_throws ArgumentError pushline!(slines, ([2, 3], 1)) + @test_throws ArgumentError pushline!(slines, [[4, 3], 2]) + @test_throws ArgumentError pushline!(slines, [[4, 3], [1, 2, 3, 4]]) + + @test_throws ArgumentError GAPSLProgram(slines) + p = GAPSLProgram(slines, 9) + @test p.lines == slines + q = GAPSLProgram([[1, 2], + ([2, 3], 1), + [[4, 3], 2], + [[4, 3], [1, 2, 3, 4]]], 9) + @test q.lines == slines + + r = GAPSLProgram( + [ [[1, 2], 3], + [[3, 2], 2], + [1, 2, 2, 1] ] ) + @test r.ngens == 1 + + @test_throws ArgumentError GAPSLProgram( + [ [[1, 2], 3], + [[3, 2], 2], + [1, 2, 4, 1] ], 2) + + @test_throws ArgumentError GAPSLProgram( + [ [[1, 2], 3], + [[3, 2], 2], + [[1, 2, 4, 1]] ], 2 ) + + r = GAPSLProgram( + [ [[1, 2], 3], + [[1, 2, 4, 2], 2], + [1, 2, 2, 1] ]) + @test r.ngens == 4 + + r = GAPSLProgram( + [ [[2, 2], 3], + [[3, 2], 2], + [1, 2, 2, 1] ]) + @test r.ngens == 2 + + r = GAPSLProgram( + [ [[2, 3], 3], + [[3, 2], 2], + [[1, 2, 2, 1], [4, 1]] ]) + @test r.ngens == 4 +end + +@testset "GAPSLDecision construction" begin + d = GAPSLDecision([ + [ [ 1, 1, 2, 1 ], 3 ], + (1, 2), # "Order" line + (2, 3), + (3, 5) ]) + @test d.ngens == 2 + + d = GAPSLDecision([ + [ [ 1, 1, 2, 1 ], 3 ], + ["Order", 1, 2], + ["Order", 2, 3], + ["Order", 3, 5] ]) + @test d.ngens == 2 + + d = GAPSLDecision([ + [ 1, 1, 2, 1 ], + (1, 2), # "Order" line + (2, 3), + (3, 5) ], 2) + @test d.ngens == 2 + + # ngens must be specified + @test_throws ArgumentError GAPSLDecision([ [1, 1], (1, 1) ]) + @test_throws ArgumentError GAPSLDecision([ [1, 1], ["Order", 1, 1] ]) + # can't have return line + @test_throws ArgumentError GAPSLDecision([ [[1, 1], 2], [[1, 1]] ]) +end + +@testset "GAPSLProgram evaluate / compile!" begin + for xy = (slpgens(3), Lazy[Gen(:x), Gen(:y)]) + x, y = xy + res = empty(xy) + + g = GAPSLProgram([ [2,3], [ 1, 2, 3, 1] ], 2 ) + + @test evaluate(g, xy) == x^2*y^3 + @test evaluate!(res, g, xy) == x^2*y^3 + slp0 = SL.compile(g) + @test slp0 isa SLProgram + @test evaluate(slp0, xy) == x^2*y^3 + slp = SL.compile(SLProgram, g) + @test slp isa SLProgram + @test slp == slp0 + @test evaluate(slp, xy) == x^2*y^3 + slp = SL.compile!(g) + @test slp == slp0 + @test g.slp[] === slp + @test evaluate(slp, xy) == x^2*y^3 + + g = GAPSLProgram([ [2,3], [ 3, 1, 1, 4] ], 2 ) + @test evaluate(g, xy) == y^3*x^4 + @test evaluate!(res, g, xy) == y^3*x^4 + slp = SL.compile!(g) + @test evaluate(slp, xy) == y^3*x^4 + + g = GAPSLProgram([ [2,3], [ [3, 1, 1, 4], [ 1, 2, 3, 1]] ], 2 ) + @test evaluate(g, xy) == [y^3*x^4, x^2*y^3] + @test evaluate!(res, g, xy) == [y^3*x^4, x^2*y^3] + slp = SL.compile!(g) + # convert to Vector{Any} as otherwise we don't get a vector + # but an SLP returning a list + @test evaluate(slp, convert(Vector{Any}, xy)) == [y^3*x^4, x^2*y^3] + + g = GAPSLProgram([ [ [1,1,2,2], 2 ], [2,3,1,1] ] ) + @test evaluate(g, xy) == (x*y^2)^3*x + @test evaluate!(res, g, xy) == (x*y^2)^3*x + slp = SL.compile!(g) + @test evaluate(slp, xy) == (x*y^2)^3*x + + g = GAPSLProgram([ [ [1,-1,2,-2], 2 ], [2,-3,1,1] ] ) + t = (x^-Int(1)*y^-Int(2))^-Int(3)*x + @test evaluate(g, xy) == t + @test evaluate!(res, g, xy) == t + slp = SL.compile!(g) + @test evaluate(slp, xy) == t + + # assignments + g = GAPSLProgram([ [[2, 3], 3] ]) + @test evaluate(g, xy) == y^3 + slp = SL.compile!(g) + @test evaluate(slp, xy) == y^3 + + g = GAPSLProgram([ [[2, 3], 3], [[1, 2, 3, 2], 2] ]) + @test evaluate(g, xy) == x^2*(y^3)^2 + slp = SL.compile!(g) + @test evaluate(slp, xy) == x^2*(y^3)^2 + end +end + +@testset "GAPSLDecision evaluate / compile!" begin + p = perm"(1, 2)(3, 4)" + q = Perm([1, 4, 2, 3]) + pq = [p, q] + + for (i, j, k, r) in [(2, 3, 5, false), + (2, 3, 3, true), + (2, 2, 3, false), + (1, 3, 3, false)] + gd = GAPSLDecision([ + [ [ 1, 1, 2, 1 ], 3 ], + ["Order", 1, i], + ["Order", 2, j], + ["Order", 3, k] ]) + + @test evaluate(gd, pq) == r + + gdc = SL.compile!(gd) + @test evaluate(gdc, pq) == r + + sld = SLProgram() + m = SL.pushop!(sld, SL.times, SL.input(1), SL.input(2)) + SL.pushop!(sld, SL.decision, SL.input(1), SL.pushint!(sld, i)) + SL.pushop!(sld, SL.decision, SL.input(2), SL.pushint!(sld, j)) + dec = SL.pushop!(sld, SL.decision, m, SL.pushint!(sld, k)) + SL.pushfinalize!(sld, dec) + @test evaluate(sld, pq) == r + end +end diff --git a/StraightLinePrograms/test/runtests.jl b/StraightLinePrograms/test/runtests.jl new file mode 100644 index 000000000000..2af8447d536d --- /dev/null +++ b/StraightLinePrograms/test/runtests.jl @@ -0,0 +1,4 @@ +include("setup.jl") +include("straightline.jl") +include("gap.jl") +include("atlas.jl") diff --git a/StraightLinePrograms/test/setup.jl b/StraightLinePrograms/test/setup.jl new file mode 100644 index 000000000000..e445954dfa41 --- /dev/null +++ b/StraightLinePrograms/test/setup.jl @@ -0,0 +1,12 @@ +using Test, StraightLinePrograms, AbstractAlgebra + +using StraightLinePrograms: Const, Exp, Gen, Minus, Plus, Lazy, + Times, UniMinus, pushconst!, pushop!, + Line, Arg, constants, lines, evaluate, evaluate! + +using StraightLinePrograms: pushline!, GAPStraightLine +using StraightLinePrograms: AtlasSLProgram, AtlasLine + +const SL = StraightLinePrograms + +replstr(c) = sprint((io, x) -> show(io, "text/plain", x), c) diff --git a/StraightLinePrograms/test/straightline.jl b/StraightLinePrograms/test/straightline.jl new file mode 100644 index 000000000000..fd02e7cfe858 --- /dev/null +++ b/StraightLinePrograms/test/straightline.jl @@ -0,0 +1,1056 @@ +@testset "LazyPolyRing" begin + F = LazyPolyRing(ZZ) + @test F isa LazyPolyRing{elem_type(ZZ)} + @test F isa MPolyRing{elem_type(ZZ)} + @test base_ring(F) == ZZ +end + +@testset "Lazy" begin + x, y, z = Gen.([:x, :y, :z]) + xyz = Any[x, y, z] + + # Const + c = Const(1) + @test c isa Const{Int} + @test c isa Lazy + @test string(c) == "1" + @test isempty(SL.gens(c)) + @test c == Const(1) == Const(0x1) + @test c != Const(2) + @test evaluate(c, rand(3)) == 1 + @test evaluate(c, xyz) == 1 + + # Gen + g = Gen(:x) + @test g isa Gen + @test g isa Lazy + @test string(g) == "x" + @test SL.gens(g) == [:x] + @test g == x + @test g != y + @test evaluate(g, [2, 3, 4]) == 2 + @test evaluate(g, xyz) == x + @test evaluate(g, [Gen(:a), Gen(:b), Gen(:x)]) == Gen(:a) + + # Plus + p = Plus(c, g) + @test p isa Plus <: Lazy + @test p.xs[1] == c && p.xs[2] == g + @test string(p) == "(1 + x)" + @test SL.gens(p) == [:x] + @test p == 1+x == 0x1+x + @test p != 2+x && p != 1+y + @test evaluate(p, [2]) == 3 + @test evaluate(p, xyz) == 1 + x + + # Minus + m = Minus(p, g) + @test m isa Minus <: Lazy + @test string(m) == "((1 + x) - x)" + @test SL.gens(m) == [:x] + @test m == (1+x)-x + @test m != (1+x)+x && m != x-x && m != (1+x)-y + @test evaluate(m, [2]) == 1 + @test evaluate(m, xyz) == (1+x)-x + + # UniMinus + u = UniMinus(p) + @test u isa UniMinus <: Lazy + @test string(u) == "(-(1 + x))" + @test SL.gens(u) == [:x] + @test u == -(1+x) + @test u != (1+x) && u != -(1+y) + @test evaluate(u, [2]) == -3 + @test evaluate(u, xyz) == -(1 + x) + + # Times + t = Times(g, p) + @test t isa Times <: Lazy + @test string(t) == "(x*(1 + x))" + @test SL.gens(t) == [:x] + @test t == x*(1+x) + @test t != (1+x)*x && t != y*(1+x) && t != x*(1+y) + @test evaluate(t, [2]) == 6 + @test evaluate(t, xyz) == x*(1+x) + + # Exp + e = Exp(p, 3) + @test e isa Exp <: Lazy + @test string(e) == "(1 + x)^3" + @test SL.gens(e) == [:x] + @test e == (1+x)^3 + @test e != (1+x)^4 && e != (1+y)^3 + @test evaluate(e, [2]) == 27 + @test evaluate(e, xyz) == (1+x)^3 + + # + + p1 = e + t + @test p1 isa Plus + @test p1.xs[1] === e + @test p1.xs[2] === t + @test p1 == e+t + @test evaluate(p1, xyz) == (1+x)^3 + x*(1+x) + + p2 = p + e + @test p2 isa Plus + @test p2.xs[1] === p.xs[1] + @test p2.xs[2] === p.xs[2] + @test p2.xs[3] === e + @test p2 == p+e + @test evaluate(p2, xyz) == (1 + x + (1 + x)^3) + + p3 = e + p + @test p3 isa Plus + @test p3.xs[1] === e + @test p3.xs[2] === p.xs[1] + @test p3.xs[3] === p.xs[2] + @test p3 == e+p + + p4 = p + p + @test p4 isa Plus + @test p4.xs[1] === p.xs[1] + @test p4.xs[2] === p.xs[2] + @test p4.xs[3] === p.xs[1] + @test p4.xs[4] === p.xs[2] + @test p4 == p+p + + # - + m1 = e - t + @test m1 isa Minus + @test m1.p === e + @test m1.q === t + @test m1 == e-t + m2 = -e + @test m2 isa UniMinus + @test m2.p === e + @test m2 == -e + + # * + t1 = e * p + @test t1 isa Times + @test t1.xs[1] === e + @test t1.xs[2] === p + @test t1 == e*p + t2 = t * e + @test t2 isa Times + @test t2.xs[1] === t.xs[1] + @test t2.xs[2] === t.xs[2] + @test t2.xs[3] === e + @test t2 == t*e + t3 = e * t + @test t3 isa Times + @test t3.xs[1] === e + @test t3.xs[2] === t.xs[1] + @test t3.xs[3] === t.xs[2] + @test t3 == e*t + t4 = t * t + @test t4 isa Times + @test t4.xs[1] === t.xs[1] + @test t4.xs[2] === t.xs[2] + @test t4.xs[3] === t.xs[1] + @test t4.xs[4] === t.xs[2] + @test t4 == t*t + + # adhoc * + at1 = 3 * p + @test at1 isa Times + at2 = big(3) * p + @test at2 isa Times + @test at1 == at2 + at3 = p * 3 + @test at3 isa Times + at4 = p * big(3) + @test at4 isa Times + @test at3 == at4 + + # adhoc + + ap1 = 3 + p + @test ap1 isa Plus + ap2 = big(3) + p + @test ap2 isa Plus + @test ap1 == ap2 + ap3 = p + 3 + @test ap3 isa Plus + ap4 = p + big(3) + @test ap4 isa Plus + @test ap3 == ap4 + + # adhoc - + am1 = 3 - p + @test am1 isa Minus + am2 = big(3) - p + @test am2 isa Minus + @test am1 == am2 + am3 = p - 3 + @test am3 isa Minus + am4 = p - big(3) + @test am4 isa Minus + @test am3 == am4 + + # ^ + e1 = p^3 + @test e1 isa Exp + @test e1.p === p + @test e1.e == 3 + @test e1 == p^3 + + h = Gen(:y) + q = e1+t4*h + @test gens(q) == [:x, :y] + @test h == y + @test q == e1+t4*h == ((1 + x)^3 + (x*(1 + x)*x*(1 + x)*y)) + @test evaluate(q, [2, 3]) == 135 + @test evaluate(q, xyz) == ((1 + x)^3 + (x*(1 + x)*x*(1 + x)*y)) +end + +@testset "LazyPoly" begin + F = LazyPolyRing(zz) + r = Const(1) + Gen(:x) + p = LazyPoly(F, r) + @test parent(p) === F + @test string(p) == "(1 + x)" + for x in (gen(F, :x), F(:x)) + @test x isa LazyPoly{Int} + @test x.p isa Gen + @test x.p.g == :x + end + c1 = F(2) + @test c1 isa LazyPoly{Int} + @test c1.p isa Const{Int} + @test c1.p.c === 2 + + @test (p+c1).p isa Plus + @test (p-c1).p isa Minus + @test (-p).p isa UniMinus + @test (p*c1).p isa Times + @test (p^3).p isa Exp + @test_throws ArgumentError LazyPolyRing(ZZ)(big(1)) + c1 +end + +@testset "SLPolyRing" begin + S = SLPolyRing(zz, [:x, :y]) + @test S isa SLPolyRing{Int} + @test base_ring(S) == zz + @test symbols(S) == [:x, :y] + + for S2 in (SLPolyRing(zz, ["x", "y"]), + SLPolyRing(zz, ['x', 'y'])) + @test S2 isa SLPolyRing{Int} + @test base_ring(S2) == zz + @test symbols(S2) == [:x, :y] + end + for Sxy in (SL.PolynomialRing(zz, ["x", "y"]), + SL.PolynomialRing(zz, ['x', 'y'])) + S3, (x, y) = Sxy + @test S3 isa SLPolyRing{Int} + @test base_ring(S3) == zz + @test symbols(S3) == [:x, :y] + @test string(x) == "x" + @test string(y) == "y" + @test parent(x) == S3 + @test parent(y) == S3 + end + + S4 = SLPolyRing(zz, 3) + @test S4 isa SLPolyRing{Int} + @test ngens(S4) == 3 + X = gens(S4) + @test length(X) == 3 + @test string.(X) == ["x1", "x2", "x3"] + + S4, X = SL.PolynomialRing(zz, 0x2) + @test S4 isa SLPolyRing{Int} + @test ngens(S4) == 2 + X = gens(S4) + @test length(X) == 2 + @test string.(X) == ["x1", "x2"] + + S5 = SLPolyRing(zz, :x => 1:3, :y => [2, 4]) + @test S5 isa SLPolyRing{Int} + XS = gens(S5) + @test string.(XS) == ["x1", "x2", "x3", "y2", "y4"] + + S5, X, Y = SL.PolynomialRing(zz, :x => 1:3, "y" => [2, 4]) + @test S5 isa SLPolyRing{Int} + XS = gens(S5) + @test string.(XS) == ["x1", "x2", "x3", "y2", "y4"] + @test X == XS[1:3] + @test Y == XS[4:5] + + s1 = one(S) + @test s1 == 1 + @test s1 isa SLPoly{Int} + s0 = zero(S) + @test s0 == 0 + @test s0 isa SLPoly{Int} + + R, (x1, y1) = PolynomialRing(zz, ["x", "y"]) + + x0, y0 = gens(S) + @test ngens(S) == 2 + @test nvars(S) == 2 + @test string(x0) == "x" + @test string(y0) == "y" + @test replstr(x0) == "x" + @test replstr(y0) == "y" + @test string(S(2)) == "2" + @test replstr(S(2)) == "2" + + for x in (gen(S, 1), x0) + @test string(x) == "x" + @test x isa SLPoly{Int} + @test convert(R, x) == x1 + end + for y in (gen(S, 2), y0) + @test string(y) == "y" + @test y isa SLPoly{Int} + @test convert(R, y) == y1 + end + + for t = (2, big(2), 0x2) + @test S(t) isa SLPoly{Int,typeof(S)} + end +end + +@testset "SLPoly" begin + S = SLPolyRing(zz, [:x, :y]) + p = SLPoly(S, SLProgram{Int}()) + @test p isa SLPoly{Int,typeof(S)} <: MPolyElem{Int} + @test parent(p) === S + p = SLPoly(S) + @test p isa SLPoly{Int,typeof(S)} <: MPolyElem{Int} + @test parent(p) === S + + @test S(p) === p + + @test zero(p) == zero(S) + @test zero(p) isa SLPoly{Int} + @test one(p) == one(S) + @test one(p) isa SLPoly{Int} + + # copy + q = SLPoly(S) + # TODO: do smthg more interesting with q + push!(constants(q), 3) + push!(q.slprogram.lines, Line(0)) + copy!(p, q) + p2 = copy(q) + for p1 in (p, p2) + @test constants(p1) == constants(q) && constants(p1) !== constants(q) + @test lines(p1) == lines(q) && lines(p1) !== lines(q) + end + S2 = SLPolyRing(zz, [:z, :t]) + @test_throws ArgumentError copy!(SLPoly(S2, SLProgram{Int}()), p) + @test_throws ArgumentError S2(p) # wrong parent + + # building + p = SLPoly(S) + l1 = pushconst!(p.slprogram, 1) + @test constants(p) == [1] + @test l1 === SL.asconstant(1) + + # currently not supported anymore + # l2 = pushconst!(p, 3) + # @test constants(p) == [1, 3] + # @test l2 === SL.asconstant(2) + + l3 = pushop!(p, SL.plus, l1, SL.input(1)) + @test l3 == Arg(UInt64(1)) + @test lines(p)[1].x == 0x0340000018000001 + l4 = pushop!(p, SL.times, l3, SL.input(2)) + @test l4 == Arg(UInt64(2)) + @test lines(p)[2].x ==0x0500000018000002 + pl = copy(lines(p)) + @test p === SL.pushfinalize!(p, l4) + + @test lines(p) == [Line(0x0340000018000001), Line(0x0500000018000002)] + SL.pushinit!(p) + @test pl == lines(p) + SL.pushfinalize!(p, l4) + @test lines(p) == [Line(0x0340000018000001), Line(0x0500000018000002)] + # p == (1+x)*y + @test SL.evaluate!(Int[], p, [1, 2]) == 4 + @test SL.evaluate!(Int[], p, [0, 3]) == 3 + l5 = SL.pushinit!(p) + l6 = SL.pushop!(p, SL.times, SL.input(1), SL.input(2)) # xy + l7 = SL.pushop!(p, SL.exponentiate, l5, Arg(2)) # ((1+x)y)^2 + l8 = SL.pushop!(p, SL.minus, l6, l7) # xy - ((1+x)y)^2 + SL.pushfinalize!(p, l8) + @test string(p) == "((x*y) - ((1 + x)*y)^2)" + @test SL.evaluate!(Int[], p, [2, 3]) == -75 + @test SL.evaluate!(Int[], p, [-2, -1]) == 1 + + @test SL.evaluate(p, [2, 3]) == -75 + @test SL.evaluate(p, [-2, -1]) == 1 + + # nsteps + @test nsteps(p) == 5 + + # compile! + pf = SL.compile!(p) + @test pf([2, 3]) == -75 + @test pf([-2, -1]) == 1 + res = Int[] + for xy in eachcol(rand(-99:99, 2, 100)) + v = Vector(xy) # TODO: don't require this + @test pf(v) == SL.evaluate(p, v) == SL.evaluate!(res, p, v) + end + + # conversion -> MPoly + R, (x1, y1) = PolynomialRing(zz, ["x", "y"]) + q = convert(R, p) + @test q isa Generic.MPoly + @test parent(q) === R + @test q == -x1^2*y1^2-2*x1*y1^2+x1*y1-y1^2 + R2, (x2, y2) = PolynomialRing(zz, ["y", "x"]) + @test_throws ArgumentError convert(R2, p) + + @test convert(R, S()) == R() + @test convert(R, S(3)) == R(3) + @test convert(R, S(-4)) == R(-4) + + # conversion MPoly -> SLPoly + for _=1:100 + r = rand(R, 1:20, 0:13, -19:19) + @test convert(R, convert(S, r)) == r + @test convert(R, convert(S, r; limit_exp=true)) == r + end + r = R() + @test convert(R, convert(S, r)) == r + + # construction from LazyPoly + L = LazyPolyRing(zz) + x, y = L(:x), L(:y) + q = S(L(1)) + @test string(q) == "1" + @test convert(R, q) == R(1) + @test convert(R, S(x*y^2-x)) == x1*y1^2-x1 + @test convert(R, S(-(x+2*y)^3-4)) == -(x1+2*y1)^3-4 + + # corner cases + @test_throws ArgumentError convert(R, SLPoly(S)) # error can change + @test convert(R, S(L(1))) == R(1) + @test convert(R, S(L(:x))) == x1 + + # mutating ops + X, Y = gens(S) + p = S(x*y-16*y^2) + # SL.addeq!(p, S(x)) # TODO: this bugs + @test p === SL.addeq!(p, S(x*y)) + @test convert(R, p) == 2*x1*y1-16y1^2 + @test p == X*Y-16Y^2+X*Y + + @test p === SL.subeq!(p, S(-16y^2)) + @test convert(R, p) == 2*x1*y1 + @test p == X*Y-16Y^2+X*Y- (-16*Y^2) + + @test p === SL.subeq!(p) + @test convert(R, p) == -2*x1*y1 + # @test p === SL.muleq!(p, p) # TODO: this bugs + @test p === SL.muleq!(p, S(-2*x*y)) + @test convert(R, p) == 4*(x1*y1)^2 + @test p === SL.expeq!(p, 3) + @test convert(R, p) == 64*(x1*y1)^6 + + # permutegens! and ^ + _, (X1, X2, X3) = SL.PolynomialRing(zz, [:x1, :x2, :x3]) + p = X1*X2^2+X3^3 + p0 = copy(p) + perm = [3, 1, 2] + q = SL.permutegens!(p, perm) + @test q === p + @test q == X3*X1^2+X2^3 + @test q == p0^Perm(perm) + @test p0 == X1*X2^2+X3^3 # not mutated + + # binary/unary ops + p = S(x*y - 16y^2) + p = p + S(x*y) + @test p isa SLPoly{Int} + @test convert(R, p) == 2*x1*y1-16y1^2 + p = p - S(-16y^2) + @test p isa SLPoly{Int} + @test convert(R, p) == 2*x1*y1 + p = -p + @test p isa SLPoly{Int} + @test convert(R, p) == -2*x1*y1 + p = p * S(-2*x*y) + @test p isa SLPoly{Int} + @test convert(R, p) == 4*(x1*y1)^2 + p = p^3 + @test p isa SLPoly{Int} + @test convert(R, p) == 64*(x1*y1)^6 + + # adhoc ops + p = S(x*y - 16y^2) + q = convert(R, p) + @test convert(R, 2p) == 2q + @test convert(R, p*3) == q*3 + @test convert(R, 2+p) == 2+q + @test convert(R, p+2) == q+2 + @test convert(R, 2-p) == 2-q + @test convert(R, p-2) == q-2 + + R = ResidueRing(ZZ, 3) + a = R(2) + S = SLPolyRing(R, [:x, :y]) + x, y = gens(S) + @test parent(a*x) == S + @test parent(x*a) == S + @test parent(a+x) == S + @test parent(x+a) == S + @test parent(a-x) == S + @test parent(x-a) == S + + # 3-args evaluate + S = SLPolyRing(zz, [:x, :y]) + x, y = gens(S) + p = 2*x^3+y^2+3 + @test evaluate(p, [2, 3]) == 28 + @test evaluate!(Int[], p, [2, 3]) == 28 + @test evaluate(p, [2, 3], identity) == 28 + @test evaluate!(Int[], p, [2, 3], x -> x) == 28 + @test evaluate(p, [2, 3], x -> -x) == -10 + @test evaluate!(Int[], p, [2, 3], x -> -x) == -10 + + # trivial rings + S = SLPolyRing(zz, Symbol[]) + gs = gens(S) + @test evaluate(S(1), gs) == S(1) + @test evaluate!(empty(gs), S(1), gs) == S(1) + + # evaluate MPoly at SLPolyRing generators + R, (x, y) = PolynomialRing(zz, ["x", "y"]) + S = SLPolyRing(zz, [:x, :y]) + X, Y = gens(S) + p = evaluate(x+y, [X, Y]) + # this is bad to hardcode exactly how evaluation of `x+y` happens, + # we just want to test that this works and looks correct + @test p == 0 + 1*(1*X^1*Y^0) + 1*(1*X^0*Y^1) +end + +@testset "Free" begin + x, y, z = xyz = SL.freegens(3) + + @test xyz == gens(SL.Free, 3) + + xs = Float64[2, 3, 4] + + @test evaluate(x, xs) == 2 + @test evaluate(y, xs) == 3 + @test evaluate(z, xs) == 4 + + @test evaluate(x, xyz) == x + @test evaluate(y, xyz) == y + @test evaluate(z, xyz) == z + + # constant + p = Free(1) + @test evaluate(p, rand(Int, rand(0:20))) === 1 + p = Free([1, 2, 4]) + @test evaluate(p, rand(Int, rand(0:20))) == [1, 2, 4] + + p = 9 + 3*x*y^2 + ((y+z+3-x-3)*2)^-2 * 100 + @test evaluate(p, xs) == 64 + @test evaluate(p, xyz) == p + + a, b = SL.freegens([:a, :bc]) + @test string(a) == "a" && string(b) == "bc" + + q1 = SL.compile(SLProgram, p) + @test evaluate(q1, xs) == 64 + q2 = SL.compile(p) + @test q2 == q1 + @test evaluate(q2, xs) == 64 + q3 = SLProgram(p) + @test q3 == q1 + @test evaluate(q3, xs) == 64 + + x1, x2, x3 = gens(SLProgram, 3) + @test SLProgram(x) == slpgen(1) == x1 + @test SLProgram(y) == slpgen(2) == x2 + @test SLProgram(z) == slpgen(3) == x3 +end + +@testset "SL internals" begin + @test SL.showop == Dict(SL.assign => "->", + SL.plus => "+", + SL.uniminus => "-", + SL.minus => "-", + SL.times => "*", + SL.divide => "/", + SL.exponentiate => "^", + SL.keep => "keep", + SL.decision => "&", + SL.getindex_ => "[]") + @test length(SL.showop) == 10 # tests all keys are distinct + for op in keys(SL.showop) + @test SL.isassign(op) == (op == SL.assign) + @test SL.istimes(op) == (op == SL.times) + # ... + @test (op.x & 0x8000000000000000 != 0) == + SL.isquasiunary(op) == + (op ∈ (SL.uniminus, SL.exponentiate, SL.keep)) + @test SL.isunary(op) == (op ∈ (SL.uniminus, SL.keep)) + end + + # pack & unpack + ops = SL.Op.(rand(UInt64(0):UInt64(0xff), 100) .<< 62) + is = rand(UInt64(0):SL.argmask, 100) + js = rand(UInt64(0):SL.argmask, 100) + @test SL.unpack.(SL.pack.(ops, Arg.(is), Arg.(js))) == tuple.(ops, Arg.(is), Arg.(js)) + + for x = rand(Int64(0):Int(SL.cstmark-1), 100) + if SL.isinput(Arg(x)) + @test SL.input(x) == Arg(x) + else + @test SL.input(x).x ⊻ SL.inputmark == x + end + end +end + +@testset "Arg" begin + @test_throws InexactError Arg(-1) + @test_throws ArgumentError Arg(typemax(Int)) + @test_throws ArgumentError Arg(1 + SL.argmask % Int) + a = Arg(SL.argmask) + @test a.x == SL.argmask + + @test_throws ArgumentError SL.intarg(SL.payloadmask % Int) + x = (SL.payloadmask ⊻ SL.negbit) % Int + @test SL.getint(SL.intarg(x)) == x + @test_throws ArgumentError SL.intarg(x+1) + @test SL.getint(SL.intarg(-x)) == -x + @test SL.getint(SL.intarg(-x-1)) == -x-1 + @test_throws ArgumentError SL.intarg(-x-2) +end + +@testset "SLProgram" begin + x, y, z = Gen.([:x, :y, :z]) + xyz = Any[x, y, z] + + p = SLProgram() + @test p isa SLProgram{Union{}} + @test isempty(p.cs) + @test isempty(p.lines) + @test !isassigned(p.f) + + p = SLProgram{Int}() + @test p isa SLProgram{Int} + @test isempty(p.cs) + @test isempty(p.lines) + @test !isassigned(p.f) + + # construction/evaluate/ninputs/aslazy + p = SLProgram{Int}(1) + @test evaluate(p, [10, 20]) == 10 + @test SL.ninputs(p) == 1 + @test SL.aslazy(p) == Gen(:x) + p = SLProgram(3) + @test evaluate(p, [10, "20", 'c']) == 'c' + @test SL.ninputs(p) == 3 + @test SL.aslazy(p) == Gen(:z) + + p = SLProgram(Const(3)) + @test evaluate(p, [10, 20]) == 3 + @test SL.aslazy(p) == Const(3) + p = SLProgram(Const('c')) + @test evaluate(p, ["10", 20]) == 'c' + @test SL.ninputs(p) == 0 + @test SL.aslazy(p) == Const('c') + + # exponent + p = SLProgram(x*y^-2) + @test SL.aslazy(p) == x*y^Int(-2) + e = (SL.payloadmask ⊻ SL.negbit) % Int + @test x^e == SL.aslazy(SLProgram(x^e)) + @test_throws ArgumentError SLProgram(x^(e+1)) + e = -e-1 + @test x^e == SL.aslazy(SLProgram(x^e)) + @test_throws ArgumentError SLProgram(x^(e-1)) + p = SLProgram(x^2*y^(Int(-3))) + @test evaluate(p, [sqrt(2), 4^(-1/3)]) === 8.0 + @test evaluate(p^-1, [sqrt(2), 4^(-1/3)]) == 0.125 + p = SLProgram(x^0) + @test evaluate(p, [2]) == 1 + @test evaluate(p, xyz) == x^0 + + # assign + p = SLProgram{Int}() + k = SL.pushop!(p, SL.plus, SL.input(1), SL.input(2)) + k = SL.pushop!(p, SL.times, k, SL.input(2)) + @assert length(p.lines) == 2 + k = SL.pushop!(p, SL.assign, k, Arg(1)) + @test k == Arg(1) + SL.pushfinalize!(p, k) + @test evaluate(p, Lazy[x, y]) == (x+y)*y + k = SL.pushop!(p, SL.exponentiate, k, Arg(2)) + @test evaluate(p, Lazy[x, y]) == (x+y)*y + SL.pushfinalize!(p, k) + @test evaluate(p, Lazy[x, y]) == ((x+y)*y)^2 + SL.pushfinalize!(p, Arg(2)) + @test evaluate(p, Lazy[x, y]) == ((x+y)*y) + k = SL.pushop!(p, SL.assign, k, Arg(2)) + @test k == Arg(2) + @test evaluate(p, Lazy[x, y]) == ((x+y)*y)^2 + + # nsteps + @test nsteps(p) == 5 + + # permute_inputs! + x1, x2, x3 = slpgens(3) + p = x1*x2^2+x3^3 + for perm = ([3, 1, 2], Perm([3, 1, 2])) + q = copy(p) + SL.permute_inputs!(q, perm) + @test q == x3*x1^2+x2^3 + end + + # nsteps again + @test nsteps(p) == 4 + @test nsteps(x1) == nsteps(x3) == 0 + + p = SLProgram() + i = SL.pushint!(p, 2) + j = SL.pushint!(p, 4) + k = SL.pushop!(p, SL.plus, i, j) + SL.pushfinalize!(p, k) + @test nsteps(p) == 1 + + # mutating ops + p = SLProgram{Int}(1) + q = SLProgram(Const(6)) + r = SLProgram(2) + + @test p === SL.addeq!(p, q) + @test evaluate(p, [3]) == 9 + @test SL.aslazy(p) == x+6 + + @test p === SL.subeq!(p, r) + @test evaluate(p, [3, 2]) == 7 + @test SL.aslazy(p) == x+6-y + + @test p === SL.subeq!(p) + @test evaluate(p, [3, 2]) == -7 + @test SL.aslazy(p) == -(x+6-y) + + @test p === SL.muleq!(p, r) + @test evaluate(p, [3, 2]) == -14 + @test SL.aslazy(p) == -(x+6-y)*y + + @test p === SL.expeq!(p, 3) + @test evaluate(p, [3, 2]) == -2744 + @test SL.evaluates(p, [3, 2]) == [9, 7, -7, -14, -2744] + @test SL.aslazy(p) == (-(x+6-y)*y)^3 + + @test SL.ninputs(p) == 2 + + p = SLProgram{UInt8}(1) + q = SLProgram(Const(2)) + + SL.addeq!(p, q) + @test p.cs[1] === 0x2 + @test SL.aslazy(p) == x+2 + + SL.muleq!(p, SLProgram(Const(3.0))) + @test p.cs[2] === 0x3 + @test SL.aslazy(p) == (x+2)*3.0 + + SL.subeq!(p, SLProgram(Const(big(4)))) + @test p.cs[3] === 0x4 + @test SL.aslazy(p) == (x+2)*3.0-big(4) + + @test_throws InexactError SL.addeq!(p, SLProgram(Const(1.2))) + @assert length(p.cs) == 4 # p.cs was resized before append! failed + pop!(p.cs) # set back consistent state + @assert length(p.lines) == 3 # p.lines was *not* resized before append! failed + @test SL.aslazy(p) == (x+2)*3.0-big(4) + + p2 = SL.copy_oftype(p, Float64) + @test p2 == p + @test p2.cs == p.cs + @test p2.lines == p.lines + SL.addeq!(p2, SLProgram(Const(1.2))) + @test p2.cs[4] == 1.2 + @test SL.aslazy(p2) == ((((x + 2.0)*3.0) - 4.0) + 1.2) + + p3 = copy(p) + @test p3 == p + @test p3.cs == p.cs + @test p3.lines == p.lines + @test_throws InexactError SL.addeq!(p3, SLProgram(Const(1.2))) + + # unary/binary ops + p = SLProgram{BigInt}(1) + p2 = SLProgram(1) + q = SLProgram(Const(2)) + + r = p+q + @test SL.aslazy(r) == x+2 + @test SL.constantstype(r) === Signed + + r2 = p2+q + @test SL.aslazy(r) == x+2 + @test SL.constantstype(r2) === Int + + r = r*SLProgram(Const(0x3)) + @test SL.aslazy(r) == (x+2)*3 + @test SL.constantstype(r) === Integer + + r2 = r2*SLProgram(Const(0x3)) + @test SL.aslazy(r2) == (x+2)*3 + @test SL.constantstype(r2) === Integer + + r = r-SLProgram(Const(1.2)) + @test SL.aslazy(r) == (x+2)*3-1.2 + @test SL.constantstype(r) === Real + + r = -r + @test SL.aslazy(r) == -((x+2)*3-1.2) + @test SL.constantstype(r) === Real + + r = r^3 + @test SL.aslazy(r) == (-((x+2)*3-1.2))^3 + @test SL.constantstype(r) === Real + + @testset "adhoc" begin + r = p+q + @test typeof(r) == SLProgram{Signed} + @assert evaluate(r, xyz) == x+2 + + @test evaluate(2+r, xyz) == 2+(x+2) + @test evaluate(r+big(4), xyz) == (x+2) + big(4) + + @test evaluate(2*r, xyz) == 2*(x+2) + @test evaluate(r*1.3, xyz) == (x+2) * 1.3 + + @test evaluate(2 - r, xyz) == 2 - (x+2) + @test evaluate(r - 0x12, xyz) == (x+2) - 0x12 + end + + # conversion Lazy -> SLProgram + @test SLProgram(x^2+y) isa SLProgram{Union{}} + p = SL.muleq!(SLProgram(Const(2)), SLProgram{Int}(x^2+y)) + @test p isa SLProgram{Int} + @test evaluate(p, [2, 3]) == 14 + @test SL.aslazy(p) == 2*(x^2 + y) + + l = SL.test(x^2 * y, 2) & SL.test(y^2 * x, 3) + p = SLProgram(l) + P = SymmetricGroup(4) + p1, p2 = P("(1,4,3)"), P("(1, 3)") + @test evaluate(p, [p1, p2]) + @test !evaluate(p, [p2, p1]) + @test SL.aslazy(p) == l + + # multiple return + p = SLProgram{Int}() + inputs = Any[x, y, z] + + @test evaluate(p, inputs) == [] + k1 = SL.pushop!(p, SL.plus, SL.input(1), SL.input(2)) + @test evaluate(p, inputs) == [x+y] + k2 = SL.pushconst!(p, 3) + k3 = SL.pushop!(p, SL.times, k1, k2) + @test evaluate(p, inputs) == [x+y, (x+y)*3] + SL.pushop!(p, SL.assign, k3, k1) + @test evaluate(p, inputs) == [(x+y)*3, (x+y)*3] + SL.pushfinalize!(p, k3) + @test evaluate(p, inputs) == (x+y)*3 + SL.setmultireturn!(p) + @test evaluate(p, inputs) == [(x+y)*3, (x+y)*3] + + # multiple return & list + X, Y = slpgens(2) + pl = SL.list([X*Y, X+1-Y]) + @test evaluate(pl, inputs) == [(x*y), (x+1-y)] + pl = SL.list([X, Y, X+Y]) # first elements don't add a "step"/"line" + @test evaluate(pl, inputs) == [x, y, x+y] + + @test evaluate(X^2*Y+Y^2, [X, Y]) == X^2*Y+Y^2 + + # keep + @test_throws ArgumentError SL.pushop!(p, SL.keep, Arg(3)) + # test we are still in valid state: + @test evaluate(p, inputs) == [(x+y)*3, (x+y)*3] + SL.pushop!(p, SL.keep, Arg(2)) + @test evaluate(p, inputs) == [(x+y)*3, (x+y)*3] + SL.pushop!(p, SL.keep, Arg(1)) + @test evaluate(p, inputs) == [(x+y)*3] + SL.pushop!(p, SL.keep, Arg(0)) + @test evaluate(p, inputs) == [] + + # integers + p = SLProgram() + i = SL.pushint!(p, 123) + j = SL.pushint!(p, -4) + + k = SL.pushop!(p, SL.plus, i, j) + SL.pushfinalize!(p, k) + @test evaluate(p, inputs) == 119 + + k = SL.pushop!(p, SL.times, k, SL.input(1)) + SL.pushfinalize!(p, k) + @test evaluate(p, inputs) == 119 * x + + l = SL.pushint!(p, -2) + SL.pushfinalize!(p, l) + @test evaluate(p, inputs) == -2 + + m = SL.pushop!(p, SL.minus, k, l) + SL.pushfinalize!(p, m) + @test evaluate(p, inputs) == 119 * x - (-2) + + @testset "bug with ints" begin + u = SLProgram() + i = SL.pushint!(u, 1) + j = SL.pushint!(u, 2) + k = SL.pushop!(u, SL.plus, i, j) + SL.pushfinalize!(u, k) + v = SLProgram() + i = SL.pushint!(v, 3) + j = SL.pushint!(v, 4) + k = SL.pushop!(v, SL.plus, i, j) + SL.pushfinalize!(v, k) + w = u + v + @test evaluate(w, []) == 10 + end + + # 3-args evaluate + x, y = slpgens(2) + p = 3*x*y^3 + s = evaluate(p, ["a", "b"], string) + @test s == "3abbb" + s = evaluate!(String[], p, ["a", "b"], string) + @test s == "3abbb" + + p = [1, 3, 2]*x + (1, 2) + @test evaluate(p, [2], length) == 8 +end + +@testset "SL Decision" begin + x, y = SL.lazygens(2) + a, b = ab = gens(SL.Free, 2) + + p = SLProgram() + pushop!(p, SL.decision, SL.input(1), SL.pushint!(p, 3)) + c = pushop!(p, SL.times, SL.input(1), SL.input(2)) + pushop!(p, SL.decision, c, SL.pushint!(p, 2)) + SL.setdecision!(p) + + l = SL.test(x, 3) & SL.test(x*y, 2) + f = SL.test(a, 3) & SL.test(a*b, 2) + @test evaluate(p, Any[x, y]) == l + @test evaluate(l, Any[x, y]) == l + @test evaluate(l, gens(SLProgram, 2)) == p + @test evaluate(p, ab) == f + @test evaluate(f, slpgens(2)) == p + + S = SymmetricGroup(4) + for (x, y) in eachcol(rand(S, 2, 200)) + res = order(x) == 3 && order(x*y) == 2 + @test evaluate(p, [x, y]) == res + @test evaluate(l, [x, y]) == res + end +end + +@testset "SL lists" begin + x, y = SL.freegens(2) + X, Y = slpgens(2) + + q = SL.list([x*y^2, x+1-y]) + @test evaluate(q, [x, y]) == q + @test evaluate(q, Any[x, y]) == [x*y^2, x+1-y] + @test evaluate(q, [2, 3]) == [18, 0] + @test evaluate(evaluate(q, SLProgram[X, Y]), [x, y]) == q + + @test evaluate(SL.list([q, q]), [x, y]) == SL.list([q, q]) + # TODO: list of list of SLProgram hits an assertion error, so + # handle this case more gracefully + + # test 3+ elements + r = SL.list([x, y+1, x+y, y-x]) + @test evaluate(r, [2, 3]) == [2, 4, 5, 1] + @test evaluate(r, [x, y]) == r +end + +@testset "SL compose" begin + x, y = SL.freegens(2) + X, Y = slpgens(2) + + q = SL.compose(x - y, SL.list([y, x])) + p = SL.compose(x - y, SL.list([y, x]), flatten=false) + r = SL.compose(X - Y, SL.list([Y, X])) + + for s = (q, p, r) + @test evaluate(s, [2, 3]) == 1 + @test evaluate(s, [3.0, 1.0]) == -2 + end + @test evaluate(r, [x, y]) == q + + q = SL.compose(SL.list([x+y, x-y]), SL.list([y-x, y+x])) + r = SL.compose(SL.list([X+Y, X-Y]), SL.list([Y-X, Y+X])) + + for s = (q, r) + @test evaluate(s, [2, 3]) == [6, -4] + @test evaluate(s, [3.0, 1.0]) == [2, -6] + end + @test evaluate(r, [x, y]) == q + + p = SL.compose(2.0*X+1.0, SL.list([3*X])) + @test p isa SLProgram{Real} + @test evaluate(p, [x]) == 2*(3*x)+1 +end + +@testset "SL getindex" begin + x, y = SL.freegens(2) + X, Y = slpgens(2) + + p = y[x+1] + @test p.x isa SL.Getindex + c = compile(SLProgram, p) + @test c isa SLProgram{Int} + q = Y[X+1] + @test q == c + @test q isa SLProgram{Int} + + for r = (p, q) + @test evaluate(p, Any[2, [4, 5, 6]]) == 6 + end + + # adhoc + p = y[2] + x[big(1)] + @test evaluate(p, [[1, 2, 3], [10, 11, 12]]) == 12 + p = y[x[3]] + @test evaluate(p, Any[[1, 2, 4], 1:4]) == 4 + + p = SL.slpcst([10, 20, 30])[2] + @test evaluate(p, []) == 20 + p = SL.list([X, Y, SL.slpcst(30)]) + @test p isa SLProgram + @test p[2] isa SLProgram + @test p[3] isa SLProgram + @test evaluate(p[2], [10, 20]) == 20 + @test evaluate(p[0x3], [10, 20]) == 30 + + # multi-indices + p = y[x, 2, 1] + @test evaluate(p, [3, reshape(1:27, 3, 3, 3)]) == 6 + + # integer & integer-array indexing + p = list([x, y, y-x, y+x]) + + @test p[1] == x + @test p[0x2] == y + @test p[big(4)] == y+x + + @test p[[1, 3]] == list([x, y-x]) + @test p[Any[1, 3]] == list([x, y-x]) + @test p[[4]] == list([y+x]) + @test p[Number[4]] == list([y+x]) + + p = Free(Vector{Char})[['a']] + e = evaluate(p, []) + @test e isa Vector{Vector{Char}} + @test e == [['a']] +end diff --git a/examples/GaloisGrp.jl b/examples/GaloisGrp.jl index 87d1185a33e5..92871790e315 100644 --- a/examples/GaloisGrp.jl +++ b/examples/GaloisGrp.jl @@ -7,7 +7,7 @@ module GaloisGrp using Oscar import Base: ^, +, -, * import Oscar: Hecke, AbstractAlgebra, GAP -using StraightLinePrograms +using Oscar.StraightLinePrograms export galois_group, isprimitive, istransitive, transitive_group_identification diff --git a/src/Oscar.jl b/src/Oscar.jl index bf7a890b8aec..85ec9f14fd36 100644 --- a/src/Oscar.jl +++ b/src/Oscar.jl @@ -182,6 +182,8 @@ include("Modules/FreeModules-graded.jl") include("Polymake/Ineq.jl") include("Polymake/NmbThy.jl") +include("../StraightLinePrograms/src/StraightLinePrograms.jl") + if is_dev include("../examples/ModStdNF.jl") include("../examples/PrimDec.jl")