chakravala / Reduce.jl
1
__precompile__()
2 14
module Reduce
3
using ForceImport, SyntaxTree, LinearAlgebra
4
!isdefined(Reduce,:linefilter!) && (linefilter! = SyntaxTree.linefilter)
5

6
#   This file is part of Reduce.jl. It is licensed under the MIT license
7
#   Copyright (C) 2017 Michael Reed
8

9
include(joinpath(@__DIR__,"../deps/svn.jl"))
10

11
struct PSL <: Base.AbstractPipe
12
    input::Pipe
13
    output::Pipe
14
    process::Base.Process
15 4
    function PSL()
16
        # Setup pipes and reduce process
17 14
        input = Pipe()
18 14
        output = Pipe()
19 14
        rsl = red()
20 12
        dirf = joinpath(@__DIR__,"..","deps")
21 0
        if !Sys.iswindows()
22 14
            try
23 14
                process = _spawn(red(), input, output)
24
            catch
25 14
                process = _spawn(redsys(dirf), input, output)
26
            end
27
        else
28 0
            process = _spawn(redsys(dirf), input, output)
29
        end
30
        # Close the unneeded ends of Pipes
31 14
        close(input.out)
32 14
        close(output.in)
33 14
        return new(input, output, process)
34
    end
35
end
36

37 14
Base.kill(rs::PSL) = kill(rs.process)
38 14
Base.process_exited(rs::PSL) = process_exited(rs.process)
39

40
export error, ReduceError
41
import Base: error
42

43
struct ReduceError <: Exception
44 14
    errstr::String
45
end
46

47 14
Base.showerror(io::IO, err::ReduceError) = print(io,"Reduce: "*chomp(err.errstr))
48

49 4
function ReduceCheck(output) # check for REDUCE errors
50 14
    occursin(r"(([*]{5})|([+]{3}) )|( ?  \(Y or N\))",output) && throw(ReduceError(output))
51
end
52

53 4
function ReduceWarn(output) # check for REDUCE warnings
54 14
    if occursin(r"[*]{3}",output)
55 12
        @info "REDUCE: $(chomp(chomp(output)))"
56 14
        join(split(output,r"[*]{3}.*\n"))
57
    else
58 14
        output
59
    end
60
end
61

62 4
function PipeClogged(tf::Bool,c::Int,info::String)
63 14
    @warn "Reduce pipe clogged by $info, $(tf ? "success" : "failure") after $c tries"
64
end
65

66 14
clear(rs::PSL) = (write(rs.input,";\n"); readavailable(rs.output))
67 14
clears = (()->(c=true; return (tf=c)->(c≠tf && (c=tf); return c)))()
68

69
const EOT = Char(4) # end of transmission character
70
EOTstr = "symbolic write(int2id $(Int(EOT)))"
71

72 4
function Base.write(rs::PSL, input::String)
73 14
    clears() && clear(rs)
74 14
    write(rs.input,"$input; $EOTstr;\n")
75
end
76

77
const SOS = "[0-9]+: " # REDUCE terminal prompt
78
const RES = Regex("\n($EOT\n$SOS)|(\n$SOS\n$EOT)|(\n$SOS$EOT\n)|($EOT\n)")
79

80 4
function Base.read(rs::PSL) # get result and strip prompts/EOT char
81 14
    out = String(readuntil(rs.output,EOT))*String(readavailable(rs.output))
82 0
    Sys.iswindows() && (out = replace(out,r"\r" => ""))
83 14
    out = replace(replace(out,r"\$\n\n" => "\n\n"),RES=>"")
84 14
    out = replace(out,Regex(SOS) => "") |> chomp |> chomp |> join
85 12
    ReduceCheck(out)
86 14
    return ReduceWarn(out)
87
end
88

89 14
readsp(rs::PSL) = split(read(rs),"\n\n\n")
90

91
include("rexpr.jl") # load RExpr features
92
include("parser.jl") # load parser generator
93
include("repl.jl") # load repl features
94
include("switch.jl") # load switch operators
95

96 14
module Algebra
97
using Reduce, LinearAlgebra
98
import AbstractTensors
99
import AbstractTensors: conj, inv, dot, PROD, SUM, -, /, 
100
import AbstractTensors: sqrt, abs, exp, expm1, log, log1p, sin, cos, sinh, cosh, ^
101
const *,+ = PROD,SUM
102

103
export expm1, log1p, 
104
for T  (RExpr,Expr,Symbol,String)
105
    @eval begin
106 0
        Base.expm1(z::$T) = exp(z)-1
107 0
        Base.log1p(z::$T) = 1+log(z)
108 0
        dot(a::$T,b) = conj(a)*b
109 0
        dot(a,b::$T) = conj(a)*b
110 0
        dot(a::$T,b::$T) = conj(a)*b
111
    end
112
end
113 0
dot(a::Expr,b::Symbol) = conj(a)*b
114 0
dot(a::Symbol,b::Expr) = conj(a)*b
115

116
include("unary.jl") # load unary operators
117
include("args.jl") # load calculus operators
118

119 4
function init_subtype(name)
120 14
    Expr(:block,[:(Base.$i(r::$name...)=$i(RExpr.(r)...)|>$name) for i  [alg;iops]]...) |> eval
121 14
    Expr(:block,[:($i(r::$name...)=$i(RExpr.(r)...)|>$name) for i  [calculus;cnan;cmat]]...) |> eval
122 14
    Expr(:block,[:(Base.$i(r::$name)=$i(RExpr(r))|>$name) for i  [sbas;[:length]]]...) |> eval
123 14
    Expr(:block,[:($i(r::$name)=$i(RExpr(r))|>$name) for i  [sfun;snan;snum;scom;sint;sran;smat]]...) |> eval
124 14
    :(Base.promote(r::$name,x) = (r,$name(x))) |> eval
125 14
    :(Base.promote(x,r::$name) = (r,$name(x))) |> eval
126 14
    :(Base.signbit(r::$name) = false) |> eval
127
end
128

129
const variables = [
130
    :root_multiplicities,
131
    :requirements,
132
    :assumptions,
133
    :low_pow,
134
    :high_pow,
135
    :euler_gamma,
136
    :khinchin,
137
]
138

139 8
for var  [variables;[:ws]]
140 14
    :($var() = rcall(RExpr($(string(var)))) |> parse) |> eval
141
end
142

143 2
function bypass(op::Symbol,T::ExprSymbol,args)
144 10
    1  args && (@eval $op(x::X) where X<:$T = Base.$op(x))
145 10
    if 2  args
146 10
        @eval begin
147 0
            $op(a::A,b::B) where {A<:ExprSymbol,B<:$T} = Base.$op(a,b)
148 0
            $op(a::A,b::B) where {A<:$T,B<:ExprSymbol} = Base.$op(a,b)
149
        end
150
    end
151
end
152

153
import AbstractTensors: TensorAlgebra, value
154 8
for op  (:exp,:reverse)
155 10
    bypass(op,:TensorAlgebra,(1,))
156
end
157 8
for op  (:+,:*,:/)
158 10
    bypass(op,:TensorAlgebra,(2,))
159
end
160
@eval bypass(:-,:TensorAlgebra,(1,2))
161 0
rank(x::X) where X<:TensorAlgebra = LinearAlgebra.rank(x)
162

163
@doc """
164
    ws()
165

166
Several mechanisms are available for saving and retrieving previously evaluated expressions. The simplest of these refers to the last algebraic expression simplified. When an assignment of an algebraic expression is made, or an expression is evaluated at the top level, (i.e., not inside a compound statement or procedure) the results of the evaluation are automatically saved in a variable `ws` that we shall refer to as the workspace. (More precisely, the expression is assigned to the variable `ws` that is then available for further manipulation.)
167

168
*Example:* If we evaluate the expression `(x+y)^2` at the top level and next wish to differentiate it with respect to `y`, we can simply say
169
``
170
julia> Algebra.df(Algebra.ws(),:y)
171
```
172
to get the desired answer.
173
""" Reduce.Algebra.ws
174

175
@doc """
176
    root_multiplicities()
177

178
Solution multiplicities are stored in the global variable `root_multiplicities` rather than the solution list. The value of this variable is a list of the multiplicities of the solutions for the last call of `solve`. For example,
179
```Julia
180
julia> Algebra.solve(:(x^2==2x-1),:x); Algebra.root_multiplicities()
181
```
182
gives the results
183
```Julia
184
(:(x = 1),)
185
 
186
(2,)
187
```
188
""" Reduce.Algebra.root_multiplicities
189

190
@doc """
191
    requirements()
192

193
The proper design of a variable sequence supplied as a second argument to `solve` is important for the structure of the solution of an equation system. Any unknown in the system not in this list is considered totally free. E.g. the call
194
```Julia
195
Algebra.solve((:(x==2z),:(z==2y)),(:z,))
196
```
197
produces an empty list as a result because there is no function `z = z(x,y)` which fulfills both equations for arbitrary `x` and `y` values. In such a case the share variable `requirements` displays a set of restrictions for the parameters of the system:
198

199
```Julia
200
julia> Algebra.requirements()
201
(:(x - 4y),)
202
```
203

204
The non-existence of a formal solution is caused by a contradiction which disappears only if the parameters of the initial system are set such that all members of the requirements list take the value zero. For a linear system the set is complete: a solution of the requirements list makes the initial system solvable. E.g. in the above case a substitution ``x = 4y`` makes the equation set consistent. For a non-linear system only one inconsistency is detected. If such a system has more than one inconsistency, you must reduce them one after the other. 1 The set shows you also the dependency among the parameters: here one of ``x`` and ``y`` is free and a formal solution of the system can be computed by adding it to the variable list of `solve`. The requirement set is not unique – there may be other such sets.
205
""" Reduce.Algebra.requirements
206

207
@doc """
208
    assumptions()
209

210
A system with parameters may have a formal solution, e.g. 
211
```Julia
212
julia> Algebra.solve((:(x==a*z+1),:(0==b*z-y)),(:z,:x))
213
(:(z = y // b), :(x = (a * y + b) // b))
214
```
215
which is not valid for all possible values of the parameters. The variable `assumptions` contains then a list of restrictions: the solutions are valid only as long as none of these expressions vanishes. Any zero of one of them represents a special case that is not covered by the formal solution. In the above case the value is
216
```Julia
217
julia> Algebra.assumptions()
218
(:b,)
219
```
220
which excludes formally the case ``b = 0``; obviously this special parameter value makes the system singular. The set of assumptions is complete for both, linear and non–linear systems.
221
""" Reduce.Algebra.assumptions
222

223
@doc """
224
    euler_gamma()
225

226
Euler's constant, also available as -\$\\psi(1)\$.
227
""" Reduce.Algebra.euler_gamma
228

229
@doc """
230
    khinchin()
231

232
Khinchin's constant, defined as
233

234
\$ \\prod_{n=1}^\\infty \\left( 1 + \\frac{1}{n(n+2)} \\right)^{\\log_2 n}. \$
235
""" Reduce.Algebra.khinchin
236

237
end
238

239
export Algebra, @force, @vars
240

241
macro vars(syms...)
242
    Expr(:block,[Expr(:(=),esc(sym),QuoteNode(sym)) for sym  syms]...)
243
end
244

245 14
Base.write(rs::PSL,r::RExpr) = write(rs,convert(String,r))
246

247
import Base: zero, one
248

249
for T  [:Any,:Expr,:Symbol]
250
    @eval begin
251 0
        zero(::Type{$T}) = 0
252 0
        zero(::$T) = 0
253 0
        one(::Type{$T}) = 1
254 0
        one(::$T) = 1
255
    end
256
end
257

258 0
zero(::Type{RExpr}) = R"0"
259 0
zero(::RExpr) = R"0"
260 0
one(::Type{RExpr}) = R"1"
261 0
one(::RExpr) = R"1"
262

263
import LinearAlgebra: transpose, adjoint
264

265 14
transpose(r::ExprSymbol) = r
266 0
transpose(r::RExpr) = r
267 14
adjoint(r::ExprSymbol) = Algebra.conj(r)
268 0
adjoint(r::RExpr) = Algebra.conj(r)
269

270
## Setup
271

272
const offlist = [:nat,:latex,:exp]
273

274
export load_package, @load_package
275

276
"""
277
    load_package(::Symbol)
278

279
Loads the specified package into REDUCE
280

281
## Examples
282
```julia-repl
283
julia> load_package(:rlfi)
284
```
285
"""
286 4
function load_package(pkg::Union{String,Symbol},pkgs...)
287 14
    "load_package $pkg" |> rcall
288 14
    for extra in pkgs
289 14
        load_package(extra)
290
    end
291 14
    return nothing
292
end
293 4
function load_package(pkgs::Union{Array{String,1},Array{Symbol,1}})
294 14
    for pkg in pkgs
295 14
        load_package(pkg)
296
    end
297
end
298

299
macro load_package(pkg...)
300
    load_package(pkg...)
301
end
302

303
"""
304
    Reduce.Reset()
305

306
Kills the REDUCE process and starts a new instance.
307

308
## Examples
309
```julia-repl
310
julia> Reduce.Reset()
311
Reduce (Free PSL version, revision 4015),  5-May-2017 ...
312
```
313
"""
314 14
Reset() = (kill(rs); Load())
315 14
__init__() = (Load(); repl_init(); atexit(() -> kill(rs)))
316

317
# Server setup
318

319 4
function Load()
320 14
    global rs = PSL()
321

322 0
    offs = ""
323 14
    for o in offlist
324 14
        o != :nat && (offs = offs*"off $o; ")
325
    end
326 14
    write(rs.input,"off nat; $EOTstr;\n")
327 14
    banner = readuntil(rs.output,EOT) |> String
328 14
    readavailable(rs.output)
329 14
    rcsl = occursin(" CSL ",banner)
330 0
    if Sys.iswindows()
331 0
        banner = replace(banner,r"\r" => "")
332
    else
333 14
        ReduceCheck(banner)
334
    end
335 14
    global banner = split(String(banner),'\n')[rcsl ? 1 : end-3]
336 14
    load_package(:rlfi)
337 14
    offs |> RExpr |> rcall
338 14
    rcall(R"on savestructr")
339 14
    return nothing
340
end
341

342
global preload = false
343
try
344
    global preload
345
    (ENV["REDPRE"]  "0") && (preload = true)
346
catch
347
end
348
preload && include("precomp.jl")
349

350
end # module

Read our documentation on viewing source code .

Loading