login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

# Edit detail for SandBoxFSpace revision 5 of 5

 1 2 3 4 5 Editor: Bill Page Time: 2015/03/06 21:48:22 GMT+0 Note:

changed:
-)lib FSPECF
)lib BOP,BOP1,FSPECF

fricas
)lib BOP,BOP1,FSPECF
)library cannot find the file BOP,BOP1,FSPECF.

)abbrev category ES ExpressionSpace
++ Category for domains on which operators can be applied
++ Author: Manuel Bronstein
++ Date Created: 22 March 1988
++ Date Last Updated: 27 May 1994
++ Description:
++ An expression space is a set which is closed under certain operators;
++ Keywords: operator, kernel, expression, space.
ExpressionSpace() : Category == Defn where
N   ==> NonNegativeInteger
K   ==> Kernel %
OP  ==> BasicOperator
SY  ==> Symbol
PAREN  ==> '%paren
BOX    ==> '%box
import from Set(K) import from List(K) import from List(Set(K)) import from List(NonNegativeInteger) import from List(Symbol)
DUMMYVAR := '%dummyVar -- the 7 functions not provided are: -- kernels minPoly definingPolynomial -- coerce: K -> % eval: (%, List K, List %) -> % -- subst: (%, List K, List %) -> % -- eval: (%, List Symbol, List(List % -> %)) -> %
allKernels : % -> Set K listk : % -> List K allk : List % -> Set K unwrap : (List K, %) -> % okkernel : (OP, List %) -> % mkKerLists : List Equation % -> Record(lstk : List K, lstv : List %)
oppren := operator(PAREN)\$CommonOperators() opbox := operator(BOX)\$CommonOperators()
box(x : %) == box [x] paren(x : %) == paren [x] belong? op == has?(op, 'any) and (is?(op, PAREN) or is?(op, BOX)) listk f == parts allKernels f tower(f : %) : List(K) == sort! listk f allk l == reduce("union", [allKernels f for f in l], set()) tower(lf : List(%)) : List(K) == sort! parts allk(lf) kernels(lf : List(%)) : List(K) == reduce(setUnion\$List(K), [kernels(f) for f in lf], [])\$List(List(K)) operators f == [operator k for k in listk f] height f == reduce("max", [height k for k in kernels f], 0) freeOf?(x : %, s : SY) == not member?(s, [name k for k in listk x]) distribute x == unwrap([k for k in listk x | is?(k, oppren)], x) box(l : List %) == opbox l paren(l : List %) == oppren l freeOf?(x : %, k : %) == not member?(retract k, listk x) kernel(op : OP, arg : %) == kernel(op, [arg]) elt(op : OP, x : %) == op [x] elt(op : OP, x : %, y : %) == op [x, y] elt(op : OP, x : %, y : %, z : %) == op [x, y, z] elt(op : OP, x : %, y : %, z : %, t : %) == op [x, y, z, t] elt(op : OP, x : %, y : %, z : %, t : %, s : %) == op [x, y, z, t, s] elt(op : OP, x : %, y : %, z : %, t : %, s : %, r : %) == op [x, y, z, t, s, r] elt(op : OP, x : %, y : %, z : %, t : %, s : %, r : %, q : %) == op [x, y, z, t, s, r, q] elt(op : OP, x : %, y : %, z : %, t : %, s : %, r : %, q : %, p : %) == op [x, y, z, t, s, r, q, p] elt(op : OP, x : %, y : %, z : %, t : %, s : %, r : %, q : %, p : %, o : %) == op [x, y, z, t, s, r, q, p, o] eval(x : %, s : SY, f : List % -> %) == eval(x, [s], [f]) eval(x : %, s : OP, f : List % -> %) == eval(x, [name s], [f]) eval(x : %, s : SY, f : % -> %) == eval(x, [s], [(y : List %) : % +-> f(first y)]) eval(x : %, s : OP, f : % -> %) == eval(x, [s], [(y : List %) : % +-> f(first y)]) subst(x : %, e : Equation %) == subst(x, [e])
eval(x : %, ls : List OP, lf : List(% -> %)) == eval(x, ls, [y +-> f(first y) for f in lf]\$List(List % -> %))
eval(x : %, ls : List SY, lf : List(% -> %)) == eval(x, ls, [y +-> f(first y) for f in lf]\$List(List % -> %))
eval(x : %, ls : List OP, lf : List(List % -> %)) == eval(x, [name s for s in ls]\$List(SY), lf)
map(fn, k) == (l := [fn x for x in argument k]\$List(%)) = argument k => k::% (operator k) l
operator op == is?(op, PAREN) => oppren is?(op, BOX) => opbox error concat("Unknown operator 1: ",string(name(op)))\$String
mainKernel x == empty?(l := kernels x) => "failed" n := height(k := first l) for kk in rest l repeat if height(kk) > n then n := height kk k := kk k
-- takes all the kernels except for the dummy variables, which are second -- arguments of rootOf's, integrals, sums and products which appear only in -- their first arguments allKernels f == s := brace(l := kernels f) for k in l repeat t := (u := property(operator k, DUMMYVAR)) case None => arg := argument k s0 := remove!(retract(second arg)@K, allKernels first arg) arg := rest rest arg n := (u::None) pretend N if n > 1 then arg := rest arg union(s0, allk arg) allk argument k s := union(s, t) s
kernel(op : OP, args : List %) == not belong? op => error concat("Unknown operator 2: ",string(name(op)))\$String okkernel(op, args)
okkernel(op, l) == kernel(op, l, 1 + reduce("max", [height f for f in l], 0))\$K :: %
elt(op : OP, args : List %) == not belong? op => error concat("Unknown operator 3: ",string(name(op)))\$String ((u := arity op) case N) and (#args ~= u::N) => error "Wrong number of arguments" (v := evaluate(op, args)\$BasicOperatorFunctions1(%)) case % => v::% okkernel(op, args)
retract f == (k := mainKernel f) case "failed" => error "not a kernel" k::K::% ~= f => error "not a kernel" k::K
retractIfCan f == (k := mainKernel f) case "failed" => "failed" k::K::% ~= f => "failed" k
is?(f : %, s : SY) == (k := retractIfCan f) case "failed" => false is?(k::K, s)
is?(f : %, op : OP) == (k := retractIfCan f) case "failed" => false is?(k::K, op)
unwrap(l, x) == for k in reverse! l repeat x := eval(x, k, first argument k) x
distribute(x, y) == ky := retract y unwrap([k for k in listk x | is?(k, '%paren) and member?(ky, listk(k::%))], x)
-- in case of conflicting substitutions e.g. [x = a, x = b], -- the first one prevails. -- this is not part of the semantics of the function, but just -- a feature of this implementation. eval(f : %, leq : List Equation %) == rec := mkKerLists leq eval(f, rec.lstk, rec.lstv)
subst(f : %, leq : List Equation %) == rec := mkKerLists leq subst(f, rec.lstk, rec.lstv)
mkKerLists leq == lk := empty()\$List(K) lv := empty()\$List(%) for eq in leq repeat (k := retractIfCan(lhs eq)@Union(K, "failed")) case "failed" => error "left hand side must be a single kernel" if not member?(k::K, lk) then lk := concat(k::K, lk) lv := concat(rhs eq, lv) [lk, lv]
if % has RetractableTo Integer then intpred? : (%, Integer -> Boolean) -> Boolean
even? x == intpred?(x, even?) odd? x == intpred?(x, odd?)
intpred?(x, pred?) == (u := retractIfCan(x)@Union(Integer, "failed")) case Integer and pred?(u::Integer)
)abbrev package ES1 ExpressionSpaceFunctions1 ++ Lifting of maps from expression spaces to kernels over them ++ Author: Manuel Bronstein ++ Date Created: 23 March 1988 ++ Date Last Updated: 19 April 1991 ++ Description: ++ This package allows a map from any expression space into any object ++ to be lifted to a kernel over the expression set, using a given ++ property of the operator of the kernel. -- should not be exposed ExpressionSpaceFunctions1(F : ExpressionSpace, S : Type) : with map : (F -> S, Symbol, Kernel F) -> S ++ map(f, p, k) uses the property p of the operator ++ of k, in order to lift f and apply it to k.
== add -- prop contains an evaluation function List S -> S map(F2S, prop, k) == args := [F2S x for x in argument k]\$List(S) (p := property(operator k, prop)) case None => ((p::None) pretend (List S -> S)) args error "Operator does not have required property"
)abbrev package ES2 ExpressionSpaceFunctions2 ++ Lifting of maps from expression spaces to kernels over them ++ Author: Manuel Bronstein ++ Date Created: 23 March 1988 ++ Date Last Updated: 19 April 1991 ++ Description: ++ This package allows a mapping E -> F to be lifted to a kernel over E; ++ This lifting can fail if the operator of the kernel cannot be applied ++ in F; Do not use this package with E = F, since this may ++ drop some properties of the operators. ExpressionSpaceFunctions2(E : ExpressionSpace, F : ExpressionSpace) : with map : (E -> F, Kernel E) -> F ++ map(f, k) returns \spad{g = op(f(a1), ..., f(an))} where ++ \spad{k = op(a1, ..., an)}. == add map(f, k) == (operator(operator k)\$F) [f x for x in argument k]\$List(F)
)abbrev category FS FunctionSpace ++ Category for formal functions ++ Author: Manuel Bronstein ++ Date Created: 22 March 1988 ++ Date Last Updated: 14 February 1994 ++ Description: ++ A space of formal functions with arguments in an arbitrary ++ ordered set. ++ Keywords: operator, kernel, function. FunctionSpace(R : Comparable) : Category == Definition where OP ==> BasicOperator O ==> OutputForm SY ==> Symbol N ==> NonNegativeInteger Z ==> Integer K ==> Kernel % Q ==> Fraction R PR ==> Polynomial R MP ==> SparseMultivariatePolynomial(R, K) QF==> PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, MP, %) FSPECF ==> FunctionalSpecialFunction
SPECIALDISP ==> '%specialDisp SPECIALEQUAL ==> '%specialEqual SPECIALINPUT ==> '%specialInput
import from BasicOperatorFunctions1(%) import from K import from Integer import from NonNegativeInteger import from String import from List(OutputForm) import from List(%)
ODD := 'odd EVEN := 'even SPECIALDIFF := '%specialDiff
-- these are needed in Ring only, but need to be declared here -- because of compiler bug: if they are declared inside the Ring -- case, then they are not visible inside the IntegralDomain case. smpIsMult : MP -> Union(Record(coef:Z, var:K),"failed") smpret : MP -> Union(PR, "failed") smpeval : (MP, List K, List %) -> % smpsubst : (MP, List K, List %) -> % smpderiv : (MP, SY) -> % smpunq : (MP, List SY, Boolean) -> % kerderiv : (K, SY) -> % kderiv : K -> List % opderiv : (OP, N) -> List(List % -> %) smp2O : MP -> O bestKernel : List K -> K worse? : (K, K) -> Boolean diffArg : (List %, OP, N) -> List % dispdiff : List % -> Record(name : O, sub : O, arg : List O, orig_op : OP, level : N) ddiff : List % -> O diffEval : List % -> % dfeval : (List %, K) -> % smprep : (List SY, List N, List(List % -> %), MP) -> % diffdiff : (List %, SY) -> % subs : (% -> %, K) -> % symsub : (SY, Z) -> SY kunq : (K, List SY, Boolean) -> % pushunq : (List SY, List %) -> List % notfound : (K -> %, List K) -> (K -> %)
equaldiff : (K, K)->Boolean debugA : (List % , List %, Boolean) -> Boolean opdiff := operator('%diff)\$CommonOperators opquote := operator('%quote)\$CommonOperators
ground? x == retractIfCan(x)@Union(R,"failed") case R ground x == retract x coerce(x : SY) : % == kernel(x)@K :: % retract(x : %) : SY == symbolIfCan(retract(x)@K)::SY applyQuote(s : SY, x : %) == applyQuote(s, [x]) applyQuote(s, x, y) == applyQuote(s, [x, y]) applyQuote(s, x, y, z) == applyQuote(s, [x, y, z]) applyQuote(s, x, y, z, t) == applyQuote(s, [x, y, z, t]) applyQuote(s : SY, l : List %) == opquote concat(s::%, l) belong? op == has?(op, 'any) and (is?(op, '%diff) or is?(op, '%quote)) subs(fn, k) == kernel(operator k, [fn x for x in argument k]\$List(%))
operator op == is?(op, '%diff) => opdiff is?(op, '%quote) => opquote error concat("Unknown operator 4: ",string(name(op)))\$String
if R has ConvertibleTo InputForm then INP==>InputForm import from MakeUnaryCompiledFunction(%, %, %) indiff : List % -> INP pint : List INP-> INP differentiand : List % -> %
differentiand l == eval(first l, retract(second l)@K, third l) pint l == convert concat(convert('D)@INP, l) indiff l == r2 := convert([convert("::"::SY)@INP,convert(third l)@INP,convert('Symbol)@INP]@List INP)@INP pint [convert(differentiand l)@INP, r2] eval(f : %, s : SY) == eval(f, [s]) eval(f : %, s : OP, g : %, x : SY) == eval(f, [s], [g], x)
eval(f : %, ls : List OP, lg : List %, x : SY) == eval(f, ls, [compiledFunction(g, x) for g in lg])
setProperty(opdiff, SPECIALINPUT, indiff@(List % -> InputForm) pretend None)
variables(lx : List(%)) == l := empty()\$List(SY) for k in tower lx repeat if ((s := symbolIfCan k) case SY) then l := concat(s::SY, l) reverse! l
variables(x : %) == variables([x])
retractIfCan(x:%):Union(SY, "failed") == (k := retractIfCan(x)@Union(K,"failed")) case "failed" => "failed" symbolIfCan(k::K)
if R has Ring then
import from UserDefinedPartialOrdering(SY) import from MP import from SparseUnivariatePolynomial(%)
-- cannot use new()\$Symbol because of possible re-instantiation gendiff := '%%0
characteristic() == characteristic()\$R coerce(k : K) : % == k::MP::% symsub(sy, i) == new()\$Symbol numerator x == numer(x)::% eval(x : %, s : SY, n : N, f : % -> %) == eval(x, [s], [n], [(y : List %) : % +-> f(first(y))]) eval(x : %, s : SY, n : N, f : List % -> %) == eval(x, [s], [n], [f]) eval(x : %, l : List SY, f : List(List % -> %)) == eval(x, l, new(#l, 1), f)
elt(op : OP, args : List %) == unary? op and ((od? := has?(op, ODD)) or has?(op, EVEN)) and smaller?(leadingCoefficient(numer first args), 0) => x := op(- first args) od? => -x x elt(op, args)\$ExpressionSpace_&(%)
eval(x : %, s : List SY, n : List N, l : List(% -> %)) == eval(x, s, n, [y +-> f(first(y)) for f in l]\$List(List % -> %))
-- op(arg)^m ==> func(arg)^(m quo n) * op(arg)^(m rem n) smprep(lop, lexp, lfunc, p) == (v := mainVariable p) case "failed" => p::% -- symbolIfCan(k := v::K) case SY => p::% k := v::K g := (op := operator k) (arg := [eval(a, lop, lexp, lfunc) for a in argument k]\$List(%)) q := map(y +-> eval(y::%, lop, lexp, lfunc), univariate(p, k))\$SparseUnivariatePolynomialFunctions2(MP, %) (n := position(name op, lop)) < minIndex lop => q g a : % := 0 f := eval((lfunc.n) arg, lop, lexp, lfunc) e := lexp.n while q ~= 0 repeat m := degree q qr := divide(m, e) t1 := f ^ (qr.quotient)::N t2 := g ^ (qr.remainder)::N a := a + leadingCoefficient(q) * t1 * t2 q := reductum q a
dispdiff l == s := second(l)::O t := third(l)::O a := argument(k := retract(first l)@K) is?(k, opdiff) => rec := dispdiff a i := position(s, rec.arg) rec.arg.i := t [rec.name, hconcat(rec.sub, hconcat(","::SY::O, (i+1-minIndex a)::O)), rec.arg, rec.orig_op, (zero?(rec.level) => 0; rec.level + 1)] i := position(second l, a) m := [x::O for x in a]\$List(O) m.i := t [name(operator k)::O, hconcat(","::SY::O, (i+1-minIndex a)::O), m, operator(k), (empty? rest a => 1; 0)]
ddiff l == rec := dispdiff l opname := zero?(rec.level) => sub(rec.name, rec.sub) differentiate(rec.name, rec.level) (f := property(rec.orig_op, '%diffDisp)) case None => ((f::None) pretend (List OutputForm -> OutputForm)) (cons(opname, rec.arg)) prefix(opname, rec.arg)
Dec_Rec ==> Record(orig_k : K, sub : List(Integer), oarg : List %, arg : List %, dummies : List(%))
decode_diff : (List %) -> Dec_Rec
encode_diff(ker : K, nsub : List(Integer), args : List(%), nd : List(%)) : % == li : List(Integer) := [] lk : List(Integer) := [] i := first(nsub) k : Integer := 1 nsub := rest(nsub) while not(empty?(nsub)) repeat ii := first(nsub) nsub := rest(nsub) ii = i => k := k + 1 lk := cons(k, lk) li := cons(i, li) i := ii k := 1 lk := cons(k, lk) li := cons(i, li) nargs := copy(args) pts : List % := [] for i in reverse(li) for dm in nd repeat pts := cons(nargs(i), pts) nargs(i) := dm res := kernel(operator(ker), nargs) for i in li for k in lk for pt in pts for dm in reverse(nd) repeat for kk in 2..k repeat res := kernel(opdiff, [res, dm, dm]) res := kernel(opdiff, [res, dm, pt]) res
insert_sub(lsub : List(Integer), j : Integer) : List(Integer) == nsub : List(Integer) := [] to_insert : Boolean := true for i in lsub repeat if to_insert and not(i > j) then nsub := cons(j, nsub) to_insert := false nsub := cons(i, nsub) if to_insert then nsub := cons(j, nsub) reverse!(nsub)
pos_diff(l : List %, r1 : Dec_Rec, i : Integer) : % == nsub := insert_sub(r1.sub, i) nd := r1.dummies if not(member?(i, r1.sub)\$List(Integer)) then k : NonNegativeInteger := 0 ii := first(nsub) for j in nsub repeat i = j => break if not(ii = j) then k := k + 1 nd1 := first(nd, k)\$List(%) nd2 := rest(nd, k) dm := (new()\$Symbol)::% nd := concat(nd1, cons(dm, nd2)) encode_diff(r1.orig_k, nsub, r1.arg, nd)
diffdiff(l, x) == r1 := decode_diff(l) args := r1.arg res : % := 0 for i in minIndex args .. maxIndex args for b in args repeat (db := differentiate(b, x)) = 0 => "iterate" res := res + db*pos_diff(l, r1, i) res
dfeval(l, g) == eval(differentiate(first l, symbolIfCan(g)::SY), g, third l)
diffEval l == k : K g := retract(second l)@K ((u := retractIfCan(first l)@Union(K, "failed")) case "failed") or (u case K and symbolIfCan(k := u::K) case SY) => dfeval(l, g) op := operator k (ud := derivative op) case "failed" => -- possible trouble -- make sure it is a dummy var dumm : % := symsub(gendiff, 1)::% ss := subst(l.1, l.2 = dumm) -- output(nl::OutputForm)\$OutputPackage -- output("fixed"::OutputForm)\$OutputPackage nl := [ss, dumm, l.3] kernel(opdiff, nl) (n := position(second l, argument k)) < minIndex l => dfeval(l, g) d := ud::List(List % -> %) eval((d.n)(argument k), g, third l)
diffArg(l, op, i) == n := i - 1 + minIndex l z := copy l z.n := g := symsub(gendiff, n)::% [kernel(op, z), g, l.n]
opderiv(op, n) == -- one? n => (n = 1) => g := symsub(gendiff, n)::% [x +-> kernel(opdiff, [kernel(op, g), g, first x])] [x +-> kernel(opdiff, diffArg(x, op, i)) for i in 1..n]
kderiv k == --output("kderiv k: ",k::OutputForm)\$OutputPackage zero?(n := #(args := argument k)) => empty() op := operator k grad := (u := derivative op) case "failed" => opderiv(op, n) u::List(List % -> %) --output("kderiv grad: ",grad::OutputForm)\$OutputPackage if #grad ~= n then grad := opderiv(op, n) gs:=[g args for g in grad] --output("kderiv gs: ",gs::OutputForm)\$OutputPackage return gs
smpderiv(p, x) == map((s : R) : R +-> retract differentiate(s::PR, x), p)::% + +/[differentiate(p, k)::% * kerderiv(k, x) for k in variables p]
coerce(p : PR) : % == map(s +-> s::%, s+-> s::%, p)\$PolynomialCategoryLifting( IndexedExponents SY, SY, R, PR, %)
worse?(k1, k2) == (u := less?(name operator k1,name operator k2)) case "failed" => k1 < k2 u::Boolean
bestKernel l == empty? rest l => first l a := bestKernel rest l worse?(first l, a) => a first l
smp2O p == (r := retractIfCan(p)@Union(R,"failed")) case R =>r::R::OutputForm a := userOrdered?() => bestKernel variables p mainVariable(p)::K outputForm(map((x : MP) : % +-> x::%, univariate(p, a))\$SparseUnivariatePolynomialFunctions2(MP, %), a::OutputForm)
smpsubst(p, lk, lv) == map(x +-> match(lk, lv, x, notfound((z : K) : % +-> subs(s +-> subst(s, lk, lv), z), lk))\$ListToMap(K, %), y +-> y::%, p)\$PolynomialCategoryLifting(IndexedExponents K, K, R, MP, %)
smpeval(p, lk, lv) == map(x +-> match(lk, lv, x, notfound((z : K) : % +-> map(s +-> eval(s, lk, lv), z), lk))\$ListToMap(K, %), y +-> y::%, p)\$PolynomialCategoryLifting(IndexedExponents K, K, R, MP, %)
-- this is called on k when k is not a member of lk notfound(fn, lk) == k +-> empty? setIntersection(tower(f := k::%), lk) => f fn k
if R has ConvertibleTo InputForm then pushunq(l, arg) == empty? l => [eval a for a in arg] [eval(a, l) for a in arg]
kunq(k, l, givenlist?) == givenlist? and empty? l => k::% is?(k, opquote) and (member?(s := retract(first argument k)@SY, l) or empty? l) => interpret(convert(concat(convert(s)@InputForm, [convert a for a in pushunq(l, rest argument k) ]@List(InputForm)))@InputForm)\$InputFormFunctions1(%) (operator k) pushunq(l, argument k)
smpunq(p, l, givenlist?) == givenlist? and empty? l => p::% map(x +-> kunq(x, l, givenlist?), y +-> y::%, p)\$PolynomialCategoryLifting(IndexedExponents K, K, R, MP, %)
smpret p == "or"/[symbolIfCan(k) case "failed" for k in variables p] => "failed" map(x +-> symbolIfCan(x)::SY::PR, y +-> y::PR, p)\$PolynomialCategoryLifting(IndexedExponents K, K, R, MP, PR)
isExpt(x : %, op : OP) == (u := isExpt x) case "failed" => "failed" v := (u::Record(var : K, exponent : Z)).var is?(v, op) and #argument(v) = 1 => u "failed"
isExpt(x : %, sy : SY) == (u := isExpt x) case "failed" => "failed" v := (u::Record(var : K, exponent : Z)).var is?(v, sy) and #argument(v) = 1 => u "failed"
if R has RetractableTo Z then smpIsMult p == -- (u := mainVariable p) case K and one? degree(q := univariate(p, u::K)) (u := mainVariable p) case K and (degree(q := univariate(p, u::K))=1) and zero?(leadingCoefficient reductum q) and ((r := retractIfCan(leadingCoefficient q)@Union(R,"failed")) case R) and (n := retractIfCan(r::R)@Union(Z, "failed")) case Z => [n::Z, u::K] "failed"
evaluate(opdiff, diffEval)
debugA(a1, a2, t) == -- uncomment for debugging -- output(hconcat [a1::OutputForm, a2::OutputForm, t::OutputForm])\$OutputPackage t
decode_diff(l : List %) : Dec_Rec == da := second(l) pt := third(l) a := argument(k := retract(first l)@K) is?(k, opdiff) => rec := decode_diff(a) i := position(da, rec.arg) rec.arg.i := pt nd := rec.dummies )if false -- This would work with derivatives in canonical order i0 := first(rec.sub) i < i0 => print(l::OutputForm) error "Malformed formal derivative" if i > i0 then nd := cons(da, nd) )else -- For now no order if not(member?(i, rec.sub)\$List(Integer)) then nd := cons(da, nd) )endif [rec.orig_k, cons(i, rec.sub), rec.oarg, rec.arg, nd] i := position(da, a) a1 := copy(a) a1.i := pt [k, [i], a, a1, [da]]
equaldiff(k1, k2) == a1 := argument k1 a2 := argument k2 -- check the operator res := operator k1 = operator k2 not res => debugA(a1, a2, res) -- check the evaluation point res := (a1.3 = a2.3) not res => debugA(a1, a2, res) a1.2 = a2.2 => a1.1 = a2.1 r1 := decode_diff(a1) r2 := decode_diff(a2) not(operator(r1.orig_k) = operator(r2.orig_k)) => false not(r1.sub =\$List(Integer) r2.sub) => false not(r1.arg = r2.arg) => false od : List K := [retract(dk)@K for dk in r1.dummies] ok1 : % := (r1.orig_k)::% sk1 : % := subst(ok1, od, r2.dummies) sk1 = (r2.orig_k)::%
setProperty(opdiff, SPECIALEQUAL, equaldiff@((K, K) -> Boolean) pretend None) setProperty(opdiff, SPECIALDIFF, diffdiff@((List %, SY) -> %) pretend None) setProperty(opdiff, SPECIALDISP, ddiff@(List % -> OutputForm) pretend None)
if not(R has IntegralDomain) then -- SPECIALDIFF contains a map (List %, Symbol) -> % -- it is used when the usual chain rule does not apply, -- for instance with implicit algebraics. kerderiv(k, x) == (v := symbolIfCan(k)) case SY => v::SY = x => 1 0 (fn := property(operator k, SPECIALDIFF)) case None => ((fn::None) pretend ((List %, SY) -> %)) (argument k, x) +/[g * differentiate(y, x) for g in kderiv k for y in argument k]
mainKernel x == mainVariable numer x kernels(x : %) : List(K) == variables numer x retract(x : %) : R == retract numer x retract(x : %) : PR == smpret(numer x)::PR retractIfCan(x:%):Union(R, "failed") == retract numer x retractIfCan(x:%):Union(PR, "failed") == smpret numer x eval(x : %, lk : List K, lv : List %) == smpeval(numer x, lk, lv) subst(x : %, lk : List K, lv : List %) == smpsubst(numer x, lk, lv) differentiate(x : %, s : SY) == smpderiv(numer x, s) coerce(x : %) : OutputForm == smp2O numer x
if R has ConvertibleTo InputForm then eval(f : %, l : List SY) == smpunq(numer f, l, true) eval f == smpunq(numer f, empty(), false)
eval(x : %, s : List SY, n : List N, f : List(List % -> %)) == smprep(s, n, f, numer x)
isPlus x == (u := isPlus numer x) case "failed" => "failed" [p::% for p in u::List(MP)]
isTimes x == (u := isTimes numer x) case "failed" => "failed" [p::% for p in u::List(MP)]
isExpt x == (u := isExpt numer x) case "failed" => "failed" r := u::Record(var : K, exponent : NonNegativeInteger) [r.var, r.exponent::Z]
isPower x == (u := isExpt numer x) case "failed" => "failed" r := u::Record(var : K, exponent : NonNegativeInteger) [r.var::%, r.exponent::Z]
if R has ConvertibleTo Pattern Z then convert(x : %) : Pattern(Z) == convert numer x
if R has ConvertibleTo Pattern Float then convert(x : %) : Pattern(Float) == convert numer x
if R has RetractableTo Z then isMult x == smpIsMult numer x
if R has CommutativeRing then r : R * x : % == r::MP::% * x
if R has IntegralDomain then
import from Fraction(PR) import from MP
par : % -> %
mainKernel x == mainVariable(x)\$QF kernels(x : %) : List(K) == variables(x)\$QF univariate(x : %, k : K) == univariate(x, k)\$QF isPlus x == isPlus(x)\$QF isTimes x == isTimes(x)\$QF isExpt x == isExpt(x)\$QF isPower x == isPower(x)\$QF denominator x == denom(x)::% coerce(q : Q) : % == (numer q)::MP / (denom q)::MP coerce(q : Fraction PR) : % == (numer q)::% / (denom q)::% coerce(q : Fraction Polynomial Q) == (numer q)::% / (denom q)::% retract(x : %) : PR == retract(retract(x)@Fraction(PR)) retract(x : %) : Fraction(PR) == smpret(numer x)::PR / smpret(denom x)::PR retract(x : %) : R == (retract(numer x)@R exquo retract(denom x)@R)::R
alg_ker_set(lx : List %) : List(K) == resl : List(K) := [] ak1 : List(K) := [] for x in lx repeat for k in kernels x repeat not(is?(k, 'nthRoot) or is?(k, 'rootOf)) => "iterate" ak1 := cons(k, ak1) while not(empty?(ak1)) repeat ak := ak1 ak1 := [] for k in ak repeat needed := true for k1 in resl while needed repeat if EQ(k1, k)\$Lisp then needed := false for k1 in resl while needed repeat if k1 = k then needed := false not(needed) => "iterate" resl := cons(k, resl) ak1 := cons(k, ak1) arg := argument(k) for k1 in kernels(arg.1) repeat if (is?(k1, 'nthRoot) or is?(k1, 'rootOf)) then ak1 := cons(k1, ak1) resl
algtower(lx : List %) : List K == reverse!(sort! alg_ker_set(lx))
algtower(x : %) : List K == algtower([x])
coerce(x : %) : OutputForm == -- one?(denom x) => smp2O numer x ((denom x) = 1) => smp2O numer x smp2O(numer x) / smp2O(denom x)
retractIfCan(x:%):Union(R, "failed") == (n := retractIfCan(numer x)@Union(R, "failed")) case "failed" or (d := retractIfCan(denom x)@Union(R, "failed")) case "failed" or (r := n::R exquo d::R) case "failed" => "failed" r::R
eval(f : %, l : List SY) == smpunq(numer f, l, true) / smpunq(denom f, l, true)
if R has ConvertibleTo InputForm then eval f == smpunq(numer f, empty(), false) / smpunq(denom f, empty(), false)
eval(x : %, s : List SY, n : List N, f : List(List % -> %)) == smprep(s, n, f, numer x) / smprep(s, n, f, denom x)
-- cannot use new()\$Symbol because of possible re-instantiation gendiff := '%%1
opderiv' : (OP, N) -> List(List % -> %) opderiv'(op, n) == -- one? n => (n = 1) => g := symsub(gendiff, n)::% [x +-> conjugate(kernel(opdiff, [kernel(conjugate op, g), g, first x]))\$FSPECF(R,%)] [x +-> conjugate(kernel(opdiff, diffArg(x, conjugate op, i)))\$FSPECF(R,%) for i in 1..n]
kderiv' : K -> List % kderiv' k == --output("kderiv' k: ",k::OutputForm)\$OutputPackage zero?(n := #(args := argument k)) => empty() op := operator k grad := (u := derivative conjugate op) case "failed" => opderiv'(op, n) u::List(List % -> %) --output("kderiv' grad: ",grad::OutputForm)\$OutputPackage if #grad ~= n then grad := opderiv'(op, n) gs:=[g args for g in grad] --output("kderiv' gs: ",gs::OutputForm)\$OutputPackage return gs
-- Wirtinger's chain rule kerderiv(k, x) == --output("kerderiv(k,x): ",paren([k::OutputForm,x::OutputForm]))\$OutputPackage if (v := symbolIfCan(k)) case SY then if v::SY = x then k1:=1 else k1:=0 --output("kerderiv k1+k2: ",k1::OutputForm)\$OutputPackage return k1 if is?(k,'conjugate) then if (v := retractIfCan(argument(k)(1))@Union(SY,"failed")) case SY then k1:=0 --output("kerderiv k1+k2: ",k1::OutputForm)\$OutputPackage return k1
(fn := property(operator k, SPECIALDIFF)) case None => ((fn::None) pretend ((List %, SY) -> %)) (argument k, x) k1:= +/[(g=0=>0;g * differentiate(y, x)) _ for g in kderiv k for y in argument k] --output("kerderiv k1: ",k1::OutputForm)\$OutputPackage
gs:=kderiv' k ys:=argument k --output("kerderiv ys: ",ys::OutputForm)\$OutputPackage k2:= +/[(g'=0=>0; g' * differentiate(conjugate(y)\$FSPECF(R,%), x)) _ for g' in gs for y in ys] -- eval(g',[coerce('%conjugate)@%],[conjugate(coerce(x)@%)\$FSPECF(R,%)]) --output("kerderiv k2: ",k2::OutputForm)\$OutputPackage --output("kerderiv k1+k2: ",(k1+k2)::OutputForm)\$OutputPackage return k1+k2
differentiate(f : %, x : SY) == (smpderiv(numer f, x) * denom(f)::% - numer(f)::% * smpderiv(denom f, x)) / (denom(f)::% ^ 2)
eval(x : %, lk : List K, lv : List %) == smpeval(numer x, lk, lv) / smpeval(denom x, lk, lv)
subst(x : %, lk : List K, lv : List %) == smpsubst(numer x, lk, lv) / smpsubst(denom x, lk, lv)
par x == (r := retractIfCan(x)@Union(R, "failed")) case R => x paren x
convert(x : Factored %) : % == par(unit x) * */[par(f.factor) ^ f.exponent for f in factors x]
retractIfCan(x:%):Union(PR, "failed") == (u := retractIfCan(x)@Union(Fraction PR,"failed")) case "failed" => "failed" retractIfCan(u::Fraction(PR))
retractIfCan(x:%):Union(Fraction PR, "failed") == (n := smpret numer x) case "failed" => "failed" (d := smpret denom x) case "failed" => "failed" n::PR / d::PR
coerce(p : Polynomial Q) : % == map(x +-> x::%, y +-> y::%, p)\$PolynomialCategoryLifting(IndexedExponents SY, SY, Q, Polynomial Q, %)
if R has RetractableTo Z then coerce(x : Fraction Z) : % == numer(x)::MP / denom(x)::MP
isMult x == (u := smpIsMult numer x) case "failed" or (v := retractIfCan(denom x)@Union(R, "failed")) case "failed" or (w := retractIfCan(v::R)@Union(Z, "failed")) case "failed" => "failed" r := u::Record(coef : Z, var : K) (q := r.coef exquo w::Z) case "failed" => "failed" [q::Z, r.var]
if R has ConvertibleTo Pattern Z then convert(x : %) : Pattern(Z) == convert(numer x) / convert(denom x)
if R has ConvertibleTo Pattern Float then convert(x : %) : Pattern(Float) == convert(numer x) / convert(denom x)
)abbrev package FS2 FunctionSpaceFunctions2 ++ Lifting of maps to function spaces ++ Author: Manuel Bronstein ++ Date Created: 22 March 1988 ++ Date Last Updated: 3 May 1994 ++ Description: ++ This package allows a mapping R -> S to be lifted to a mapping ++ from a function space over R to a function space over S; FunctionSpaceFunctions2(R, A, S, B) : Exports == Implementation where R, S : Join(Ring, Comparable) A : FunctionSpace R B : FunctionSpace S
K ==> Kernel A P ==> SparseMultivariatePolynomial(R, K)
Exports ==> with map : (R -> S, A) -> B ++ map(f, a) applies f to all the constants in R appearing in \spad{a}.
Implementation ==> add smpmap : (R -> S, P) -> B
smpmap(fn, p) == map(x +-> map(z +-> map(fn, z), x)\$ExpressionSpaceFunctions2(A, B), y +-> fn(y)::B, p)\$PolynomialCategoryLifting(IndexedExponents K, K, R, P, B)
if R has IntegralDomain then if S has IntegralDomain then map(f, x) == smpmap(f, numer x) / smpmap(f, denom x) else map(f, x) == smpmap(f, numer x) * (recip(smpmap(f, denom x))::B) else map(f, x) == smpmap(f, numer x)
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-- SPAD files for the functional world should be compiled in the -- following order: -- -- op kl FSPACE expr funcpkgs
Compiling FriCAS source code from file
using old system compiler.
ES abbreviates category ExpressionSpace
------------------------------------------------------------------------
initializing NRLIB ES for ExpressionSpace
compiling into NRLIB ES
;;; *** |ExpressionSpace| REDEFINED Time: 0.01 SEC.
ES- abbreviates domain ExpressionSpace& ------------------------------------------------------------------------ initializing NRLIB ES- for ExpressionSpace& compiling into NRLIB ES- importing Set Kernel S importing List Kernel S importing List Set Kernel S importing List NonNegativeInteger importing List Symbol compiling exported box : S -> S Time: 0.02 SEC.
compiling exported paren : S -> S Time: 0 SEC.
compiling exported belong? : BasicOperator -> Boolean Time: 0 SEC.
compiling local listk : S -> List Kernel S Time: 0.01 SEC.
compiling exported tower : S -> List Kernel S Time: 0.01 SEC.
compiling local allk : List S -> Set Kernel S Time: 0 SEC.
compiling exported tower : List S -> List Kernel S Time: 0 SEC.
compiling exported kernels : List S -> List Kernel S Time: 0.01 SEC.
compiling exported operators : S -> List BasicOperator Time: 0.01 SEC.
compiling exported height : S -> NonNegativeInteger Time: 0 SEC.
compiling exported freeOf? : (S,Symbol) -> Boolean Time: 0 SEC.
compiling exported distribute : S -> S Time: 0 SEC.
compiling exported box : List S -> S Time: 0 SEC.
compiling exported paren : List S -> S Time: 0 SEC.
compiling exported freeOf? : (S,S) -> Boolean Time: 0.01 SEC.
compiling exported kernel : (BasicOperator,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S,S,S,S) -> S Time: 0.01 SEC.
compiling exported elt : (BasicOperator,S,S,S,S,S,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S,S,S,S,S,S) -> S Time: 0 SEC.
compiling exported elt : (BasicOperator,S,S,S,S,S,S,S,S,S) -> S Time: 0.01 SEC.
compiling exported eval : (S,Symbol,List S -> S) -> S Time: 0.01 SEC.
compiling exported eval : (S,BasicOperator,List S -> S) -> S Time: 0.01 SEC.
compiling exported eval : (S,Symbol,S -> S) -> S Time: 0.01 SEC.
compiling exported eval : (S,BasicOperator,S -> S) -> S Time: 0.01 SEC.
compiling exported subst : (S,Equation S) -> S Time: 0 SEC.
compiling exported eval : (S,List BasicOperator,List S -> S) -> S Time: 0.01 SEC.
compiling exported eval : (S,List Symbol,List S -> S) -> S Time: 0.01 SEC.
compiling exported eval : (S,List BasicOperator,List List S -> S) -> S Time: 0 SEC.
compiling exported map : (S -> S,Kernel S) -> S Time: 0.01 SEC.
compiling exported operator : BasicOperator -> BasicOperator Time: 0.01 SEC.
compiling exported mainKernel : S -> Union(Kernel S,failed) Time: 0.01 SEC.
compiling local allKernels : S -> Set Kernel S ****** comp fails at level 3 with expression: ****** error in function allKernels
(SEQ (|:=| |s| | << | (|brace| (|:=| |l| (|kernels| |f|))) | >> |) (REPEAT (IN |k| |l|) (SEQ (|:=| |t| (SEQ (|:=| |u| (|property| (|operator| |k|) DUMMYVAR)) (|exit| 1 (IF (|case| |u| (|None|)) (SEQ (|:=| |arg| (|argument| |k|)) (|:=| |s0| (|remove!| (@ (|retract| (|second| |arg|)) (|Kernel| S)) (|allKernels| (|first| |arg|)))) (|:=| |arg| (|rest| (|rest| |arg|))) (|:=| |n| (|pretend| (|::| |u| (|None|)) (|NonNegativeInteger|))) (IF (> |n| 1) (|:=| |arg| (|rest| |arg|)) |noBranch|) (|exit| 1 (|union| |s0| (|allk| |arg|)))) (|allk| (|argument| |k|)))))) (|exit| 1 (|:=| |s| (|union| |s| |t|))))) (|exit| 1 |s|)) ****** level 3 ****** \$x:= (brace (:= l (kernels f))) \$m:= \$EmptyMode \$f:= ((((|f| # #) (|opbox| #) (|oppren| #) (< #) ...)))
>> Apparent user error: NoValueMode is an unknown mode