[fricas-devel] Giving a name to a part of an expression

[Slawomir Kolodynski]

Suppose I have an expression like *(y-m)*sqrt(x)/s* . What I would like to
do is to give a name *e* to the *sqrt(x)/s* part and do some kind of
transformation on this expression so that I get *(y-m)*e *or equivalent as
the result*. *To do that I define a rule

substE := rule (('y-'m)*sqrt('x)/'s == ('y-'m)*'e)

However, when I try to apply this rule to the expression

substE (y-m)*sqrt(x)/s

I get the *(y-m)*sqrt(x)/s *expression back instead of *(y-m)*e. *It looks
like the left hand side of the equality in the rule substE does not pattern
match itself.

Curiously, without sqrt a similar rule works as expected:

substF := rule (('y-'m)*'x/'s == ('y-'m)*'e)

substF ((y-m)*x/s)

gives e*y - e*m

Can you explain what is the difference here and how to approach the goal of
substituting a subexpression with a symbol?

Thanks,

Slawomir

[Waldek Hebisch]

Post by Slawomir Kolodynski
Suppose I have an expression like *(y-m)*sqrt(x)/s* . What I would like to
do is to give a name *e* to the *sqrt(x)/s* part and do some kind of
transformation on this expression so that I get *(y-m)*e *or equivalent as
the result*. *To do that I define a rule
substE := rule (('y-'m)*sqrt('x)/'s == ('y-'m)*'e)
However, when I try to apply this rule to the expression
substE (y-m)*sqrt(x)/s
I get the *(y-m)*sqrt(x)/s *expression back instead of *(y-m)*e. *It looks
like the left hand side of the equality in the rule substE does not pattern
match itself.

Yes. In expression numerator is a sum, that is what you really
have is:

y*sqrt(x) - m*sqrt(x)

In pattern (to which left hand side of the rule is converted) we
have product. Product does not match sum... I need to check
is this is just a bug or an unavoidable misfeature.

BTW: In your other example both pattern and expression is a sum
so the pattern matches.

Concerning your general question: in FriCAS expression is a
quotient of polynomials in kernels. Operations that identify
a useful subexpression usually can also identify complementary
part. For example, if your subexpression is denominator, then
numerator gives complement. If your subepression is a single
term of polynomial, than subtracting this term gives complement.
If your subexpression is a single factor of a term, than
complenent is product of other terms.

Internally FriCAS make only little use of pattern matcher.
Instead, various routines analyze structure of expressions
using lower level operations.

One consequence of internal representation used by FriCAS
is that frequently pattern needed to match given expression
is different than user input would suggest. This make
pattern matching somewhat unreliable (and that is why
FriCAS normally uses other approches). But of course
expression should match to itself...

[Waldek Hebisch]

Post by Waldek Hebisch

Post by Slawomir Kolodynski
Suppose I have an expression like *(y-m)*sqrt(x)/s* . What I would like to
do is to give a name *e* to the *sqrt(x)/s* part and do some kind of
transformation on this expression so that I get *(y-m)*e *or equivalent as
the result*. *To do that I define a rule
substE := rule (('y-'m)*sqrt('x)/'s == ('y-'m)*'e)
However, when I try to apply this rule to the expression
substE (y-m)*sqrt(x)/s
I get the *(y-m)*sqrt(x)/s *expression back instead of *(y-m)*e. *It looks
like the left hand side of the equality in the rule substE does not pattern
match itself.

Yes. In expression numerator is a sum, that is what you really
y*sqrt(x) - m*sqrt(x)
In pattern (to which left hand side of the rule is converted) we
have product. Product does not match sum... I need to check
is this is just a bug or an unavoidable misfeature.

AFAICS the problem is in convertion from polynomials to patterns.
Looks like a bug, but attempting to fix it showed other
problems with rules, so I need to dig deeper.

[Slawomir Kolodynski]

Thank you for looking into this.

Instead, various routines analyze structure of expressions using lower

level operations.

Could you give some examples of those lower level operations? Pattern
matching and rewrite rules seem very useful to me , but I am ok with doing
lower level operations if those fail. I know about numerator, denominator
and monomials so for example I can do

e1 := (y-m)*x/s
(monomials numerator e1)(1)

to get x*y

However when I try to do something similar with

e2 := (y-m)*sqrt(x)/s
numerator e2

to get (y-m)*sqrt(x), I don't know how then to extract the components of
this product. Also, I could not find a way to extract an argument of a
function, for example if I have an expression

sin(a*x+b)

how to get the "a*x+b" part?
Thanks,

Slawomir

[Bill Page]

Maybe this would help:

http://axiom-wiki.newsynthesis.org/ManipulatingExpressions

Post by Slawomir Kolodynski
Thank you for looking into this.

Instead, various routines analyze structure of expressions using lower level operations.

Could you give some examples of those lower level operations? Pattern matching and rewrite rules seem very useful to me , but I am ok with doing lower level operations if those fail. I know about numerator, denominator and monomials so for example I can do
e1 := (y-m)*x/s
(monomials numerator e1)(1)
to get x*y
However when I try to do something similar with
e2 := (y-m)*sqrt(x)/s
numerator e2
to get (y-m)*sqrt(x), I don't know how then to extract the components of this product. Also, I could not find a way to extract an argument of a function, for example if I have an expression
sin(a*x+b)
how to get the "a*x+b" part?
Thanks,
Slawomir

[Waldek Hebisch]

Post by Waldek Hebisch

Instead, various routines analyze structure of expressions using lower

level operations.
Could you give some examples of those lower level operations? Pattern
matching and rewrite rules seem very useful to me , but I am ok with doing
lower level operations if those fail. I know about numerator, denominator
and monomials so for example I can do
e1 := (y-m)*x/s
(monomials numerator e1)(1)
to get x*y
However when I try to do something similar with
e2 := (y-m)*sqrt(x)/s
numerator e2
to get (y-m)*sqrt(x), I don't know how then to extract the components of
this product.

Actually, it is a sum.

Post by Waldek Hebisch
Also, I could not find a way to extract an argument of a
function, for example if I have an expression
sin(a*x+b)
how to get the "a*x+b" part?

Crucial notion is "kernel". Expression above is a single kernel
so you can retract it:

(11) -> ks := retract(sin(a*x + b))@Kernel(Expression(Integer))

(11) sin(a x + b)
Type: Kernel(Expression(Integer))

'argument' gives you list of arguments of the kernel, 'operator'
gives the operator:
(12) -> argument(ks)

(12) [a x + b]
Type: List(Expression(Integer))
(13) -> operator(ks)

(13) sin
Type: BasicOperator

You can use 'kernels' to get list of top level kernels contained in an
expression:
(14) -> kl := kernels(e2)

+-+
(14) [\|x , y, s, m]
Type: List(Kernel(Expression(Integer)))

Using a kernel you can convert numerator to univariate polynomial
with respect to kernel:
(15) -> nu := univariate(numer(e2), kl(1))

(15) (y - m)?
Type: SparseUnivariatePolynomial(SparseMultivariatePolynomial(Integer,Kernel(Expression(Integer))))

And then extract coefficients:
(16) -> coefficient(nu, 1)

(16) y - m
Type: SparseMultivariatePolynomial(Integer,Kernel(Expression(Integer)))

There is also 'tower', which gives all kernels (top level and nested):
(17) -> tower(e2)

+-+
(17) [m, s, x, y, \|x ]
Type: List(Kernel(Expression(Integer)))

(Compared to 'kernels' we also get 'x' above).

Pattern matcher is implemented on top of 'numer', 'denom', 'isPlus' and
'isTimes', 'isExpt' and kernels handling. For example:
(18) -> ne2 := numer(e2)::Expression(Integer)

+-+
(18) (y - m)\|x
Type: Expression(Integer)
(19) -> isTimes(ne2)

(19) "failed"
Type: Union("failed",...)
(20) -> isPlus(ne2)

+-+ +-+
(20) [y\|x , - m\|x ]
Type: Union(List(Expression(Integer)),...)
(21) -> isTimes(isPlus(ne2).1)

+-+
(21) [\|x , y]
Type: Union(List(Expression(Integer)),...)

As you can see structure seen by pattern matcher is somewhat different
from printed form. Note that 'isPlus' essentially gives you list
of monomials, 'isTimes' gives factors of single monomial.

[Themos Tsikas]

Hello

(1) -> f:=(y-m)*sqrt(x)/s

+-+
(y - m)\|x
(1) -----------
s
Type:
Expression(Integer)
(2) -> k_sx:= first kernels f

+-+
(2) \|x
Type:
Kernel(Expression(Integer))
(3) -> subst (f,[k_sx],[e*s])

(3) e y - e m
Type:
Expression(Integer)
(4) -> %::POLY INT

(4) e y - e m
Type:
Polynomial(Integer)
(5) -> factor %

(5) e(y - m)
Type:
Factored(Polynomial(Integer))

Regards
Themos Tsikas

Post by Slawomir Kolodynski
Suppose I have an expression like *(y-m)*sqrt(x)/s* . What I would like
to do is to give a name *e* to the *sqrt(x)/s* part and do some kind of
transformation on this expression so that I get *(y-m)*e *or equivalent
as the result*. *To do that I define a rule
substE := rule (('y-'m)*sqrt('x)/'s == ('y-'m)*'e)
However, when I try to apply this rule to the expression
substE (y-m)*sqrt(x)/s
I get the *(y-m)*sqrt(x)/s *expression back instead of *(y-m)*e. *It
looks like the left hand side of the equality in the rule substE does not
pattern match itself.
substF := rule (('y-'m)*'x/'s == ('y-'m)*'e)
substF ((y-m)*x/s)
gives e*y - e*m
Can you explain what is the difference here and how to approach the goal
of substituting a subexpression with a symbol?
Thanks,
Slawomir

[Slawomir Kolodynski]

I had the code based on Bill Page and Waldek Hebisch suggestions that was
doing what I wanted but this last hint by Themos simplified it a lot. I
will put it here as Marduk seems to be doing a similar exercise in the "Simple
Gaussian integral fails thread".

The goal is to see if we can make FriCAS compute the integral that gives
the Black-Scholes formula for the price of a call option.
So, we define

C:=exp(-r*T)/(s*sqrt(2*%pi*T))
f(x)==(S*exp(x)-K)
g(x) == exp(-(x-T*m)^2/(2*T*s^2))
m := r-s^2/2

The integral integrate(f(x)*g(x),x=log(K/S)..%plusInfinity) fails, but we
can do

I := C*integrate(f(x)*g(x),x=log(K/S)..y)

Now we need to take the limit as y goes to infinity. For that to work we
first have to tell FriCAS that the expression sqrt(1/(2*T))/s is positive.
We do it by replacing this expression by c^2.

withPos:=subst(I,[sqrt(1/(2*T))::Kernel(Expression(Integer))],[c^2*s])

Nowe we can take the limit:

lim:=limit(withPos,y=%plusInfinity)

and replace c^2 back by sqrt(1/(2*T))/s

BS := (rule ('c^2==sqrt(1/(2*T))/s)) (lim::Expression(Integer))