[GSoC] Introduction - Factorization of linear ordinary differential operators

[Abhinav Baid]

Hi all,

I'm Abhinav, a computer science undergraduate at Birla Institute of Technology and Science, Pilani. I'm interested in doing the Factorization of LODOs project as I'm quite fascinated by the related math. I've already downloaded the van Hoeij thesis and did the SPAD tutorial on the Axiom wiki [1].
I've taken a course on Differential Equations in which differential operators were introduced as a method for solving differential equations. I already know C and Java. Last year, I did a GSoC project under lmonade organisation itself for FLINT by implementing the LLL algorithm [2].
I've tried to implement the function suggested by Waldek as a coding task last year (function which given a linear ordinary differential operator returns true if the operator has constant coefficients and false otherwise). The source code for the same is available at github [3]. Could you please review it? Also, I've implemented the integrate_sols procedure from ISSAC'97 Abramov/Hoeij paper [4] in SPAD and think that it can be added as an operation to the LinearOrdinaryDifferentialOperator1 domain [5]. I'd be glad to hear any comments on this too.
I'd like to contribute to FriCAS and hence, would like ideas for any other code samples I could write.

Thanks and regards,
Abhinav.

[1] http://axiom-wiki.newsynthesis.org/ProgrammingSPAD
[2] http://flintlib.org/news.html
[3] https://github.com/fandango-/spad/blob/master/task.spad
[4] http://www.math.fsu.edu/~hoeij/papers/comments/issac1997.html
[5] https://github.com/fandango-/spad/blob/master/is.spad

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I'm Abhinav, a computer science undergraduate at Birla Institute of
> Technology and Science, Pilani. I'm interested in doing the Factorization
> of LODOs project as I'm quite fascinated by the related math.

Welcome.

> I've already
> downloaded the van Hoeij thesis and did the SPAD tutorial on the Axiom wiki
> [1].
> I've taken a course on Differential Equations in which differential
> operators were introduced as a method for solving differential equations. I
> already know C and Java. Last year, I did a GSoC project under lmonade
> organisation itself for FLINT by implementing the LLL algorithm [2].
> I've tried to implement the function suggested by Waldek as a coding task
> last year (function which given a linear ordinary differential operator
> returns true if the operator has constant coefficients and false
> otherwise). The source code for the same is available at github [3]. Could
> you please review it?

Your code looks fine. A small remark: actual computation can
be done on a single line, without extra variable.

> Also, I've implemented the integrate_sols procedure
> from ISSAC'97 Abramov/Hoeij paper [4] in SPAD and think that it can be
> added as an operation to the LinearOrdinaryDifferentialOperator1 domain
> [5]. I'd be glad to hear any comments on this too.

I took a quck look at it: the code looks good. We need to think
a bit where to put it.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

On Thursday, January 29, 2015 at 8:42:54 AM UTC+5:30, Waldek Hebisch wrote:

Abhinav Baid wrote:
>
> I'm Abhinav, a computer science undergraduate at Birla Institute of
> Technology and Science, Pilani. I'm interested in doing the Factorization
> of LODOs project as I'm quite fascinated by the related math.

Welcome.
> I've already
> downloaded the van Hoeij thesis and did the SPAD tutorial on the Axiom wiki
> [1].
> I've taken a course on Differential Equations in which differential
> operators were introduced as a method for solving differential equations. I
> already know C and Java. Last year, I did a GSoC project under lmonade
> organisation itself for FLINT by implementing the LLL algorithm [2].
> I've tried to implement the function suggested by Waldek as a coding task
> last year (function which given a linear ordinary differential operator
> returns true if the operator has constant coefficients and false
> otherwise). The source code for the same is available at github [3]. Could
> you please review it?

Your code looks fine.  A small remark: actual computation can
be done on a single line, without extra variable.

Thanks for your comment. I've made the required change [1].

> Also, I've implemented the integrate_sols procedure
> from ISSAC'97 Abramov/Hoeij paper [4] in SPAD and think that it can be
> added as an operation to the LinearOrdinaryDifferentialOperator1 domain
> [5]. I'd be glad to hear any comments on this too.

I took a quck look at it: the code looks good.  We need to think
a bit where to put it.

Another place (which, on second thought, I believe would be more apt) where the function could be added
is the RationalLODE package, as it deals with LODO/LODE having rational function coefficients.

[1] https://github.com/fandango-/spad/blob/master/task.spad
 

--
                              Waldek Hebisch
heb...@math.uni.wroc.pl

 
Thanks,
Abhinav.

[Waldek Hebisch]

This one point of view. However Abramov and van Hoeij claim that
the method works for Ore algebras. FriCAS has quite general
implementation of Ore algebras, so it make sense to make
routine more general. ATM important ingredient, that is
"rational" solver is implemented only for differential operators
with rational coefficients, but there are methods to handle
operators with more general coefficients. So in the future
we should be able to use your routine in more general context.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Oh, the situation seems to be similar to the one in Maple. However, in
Maple, the function is included in the diffop package. So, wouldn't it
be better to include this in RationalLODE instead of waiting for a
rational solver for Ore algebras?

Thanks,
Abhinav.

[Waldek Hebisch]

I do not want to wait. One trick is to have a general function
which takes needed operations as parameters. So one can
have something like:

IntegrateSols(F, L) : Exports == Implementation where
F : Join(Field, CharacteristicZero, RetractableTo Integer,
RetractableTo Fraction Integer)
L : UnivariateSkewPolynomialCategory(F)
U ==> Union(F, "failed")
SF ==> (L, F) -> Record(particular : U, basis : List F)
Exports ==> with
integrate_sols : (L, L, L -> L, SF) -> Union(L, List(L))
++ integrate_sols(op, D, adjoint, rat_solve) ...

or put extra requirements on L like

L : UnivariateSkewPolynomialCategory(F) with
adjoint : % -> %
D : () -> %


Then in some other place (like RationalLODE) we can place a wrapper
which just passes extra parameters.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Ralf Hemmecke]

> I do not want to wait. One trick is to have a general function
> which takes needed operations as parameters. So one can
> have something like:
>
> IntegrateSols(F, L) : Exports == Implementation where
> F : Join(Field, CharacteristicZero, RetractableTo Integer,
> RetractableTo Fraction Integer)
> L : UnivariateSkewPolynomialCategory(F)
> U ==> Union(F, "failed")
> SF ==> (L, F) -> Record(particular : U, basis : List F)
> Exports ==> with
> integrate_sols : (L, L, L -> L, SF) -> Union(L, List(L))
> ++ integrate_sols(op, D, adjoint, rat_solve) ...
>
> or put extra requirements on L like
>
> L : UnivariateSkewPolynomialCategory(F) with

> adjoint : % -> % (*)
> D : () -> %

Actually, either way is fine, but I would somehow prefer the second
variant. The bad thing is anyway that if (as in (*)) only the name of
the function is listed, i.e., without its axioms, then it's pretty
unclear what the exact requirements of L are. So instead of listing just
the function signatures, it should be similar to the introduction of
DifferentialRing.
https://github.com/fricas/fricas/blob/master/src/algebra//catdef.spad#L269

Another point, is the return type of integrate_sols. Isn't it a bit
unhandy to go for Union(L, List L)?

Wouldn't List(L) or Union(L, Record(ltilde: L, r: L)) be better? Or
maybe Record(ltilde: L, r: Union(L, failed))?

Ralf

[Abhinav Baid]

Oh, it was my mistake. Sorry. How is the current implementation? [1]
In RationalLODE, we could just use integrate_sols(l, ratDsolve) with it.

[1] https://github.com/fandango-/spad/blob/master/is.spad

Thanks,
Abhinav.

[Abhinav Baid]

Hi Ralf,

On 02/02/2015 03:39 PM, Ralf Hemmecke wrote:
>> I do not want to wait. One trick is to have a general function
>> which takes needed operations as parameters. So one can
>> have something like:
>>
>> IntegrateSols(F, L) : Exports == Implementation where
>> F : Join(Field, CharacteristicZero, RetractableTo Integer,
>> RetractableTo Fraction Integer)
>> L : UnivariateSkewPolynomialCategory(F)
>> U ==> Union(F, "failed")
>> SF ==> (L, F) -> Record(particular : U, basis : List F)
>> Exports ==> with
>> integrate_sols : (L, L, L -> L, SF) -> Union(L, List(L))
>> ++ integrate_sols(op, D, adjoint, rat_solve) ...
>>
>> or put extra requirements on L like
>>
>> L : UnivariateSkewPolynomialCategory(F) with
>> adjoint : % -> % (*)
>> D : () -> %
> Actually, either way is fine, but I would somehow prefer the second
> variant. The bad thing is anyway that if (as in (*)) only the name of
> the function is listed, i.e., without its axioms, then it's pretty
> unclear what the exact requirements of L are. So instead of listing just
> the function signatures, it should be similar to the introduction of
> DifferentialRing.
> https://github.com/fricas/fricas/blob/master/src/algebra//catdef.spad#L269

Thanks for your comments. I've included the function descriptions along
with the signatures.

> Another point, is the return type of integrate_sols. Isn't it a bit
> unhandy to go for Union(L, List L)?
>
> Wouldn't List(L) or Union(L, Record(ltilde: L, r: L)) be better? Or
> maybe Record(ltilde: L, r: Union(L, failed))?

I don't know, actually. Union(L, List(L)) seemed to best fit the
description of the function as mentioned in the paper.
However, Union(L, Record(ltilde: L, r: L)) also seems good as it offers
extra information about what the 2 elements actually mean in the latter
case. So, I've changed the return type to it.
> Ralf
>
Thanks,
Abhinav.

[Ralf Hemmecke]

>> Wouldn't List(L) or Union(L, Record(ltilde: L, r: L)) be better? Or
>> maybe Record(ltilde: L, r: Union(L, failed))?
> I don't know, actually. Union(L, List(L)) seemed to best fit the
> description of the function as mentioned in the paper.
> However, Union(L, Record(ltilde: L, r: L)) also seems good as it offers
> extra information about what the 2 elements actually mean in the latter
> case.

Now suppose you use this function and you are only interested in the
operator, i.e. you don't care whether or not there is a corresponding r.

You would then call

z: Union(L, Record(ltilde: L, r: L)) := integrate_sols(...)
l: L := if z case L then z::L else (z::Record(ltilde: L, r: L)).ltilde

That is why I suggested

Record(ltilde: L, r: Union(L, failed)).

It also better says what actually fails if it does fail.

Actually just having

List L

would work similarly, but I somehow don't like to return a list if it
can only have on or two elements. That's somehow misuse.

Ralf

[Waldek Hebisch]

Abhinav Baid wrote:
>
> How is the current implementation? [1]
> In RationalLODE, we could just use integrate_sols(l, ratDsolve) with it.
>
> [1] https://github.com/fandango-/spad/blob/master/is.spad

Two remarks:

- actually we do not need D. Namely D is just a shorter way to get
'monomial(1, 1)$L'.
- please keep line length below 80

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Waldek Hebisch]

Abhinav Baid wrote:
> I'd like to contribute to FriCAS and hence, would like ideas for any other
> code samples I could write.

Typically differential operators are written using ordinary derivative D.
But sametimes it is easier to write then in terms of different
derivative, in particular delta = xD (van Hoej uses this a lot).
Some time ago I wrote a domain implementing differantial operators
written in terms of delta (unfortunately in the output it is shown
as D, but it behaves as delta should behave). You can find this
domain at

http://math.uni.wroc.pl/~hebisch/fricas/lodo3.spad

It can be used like this:

F := Fraction(Integer)
uP := UnivariatePolynomial('x, F)
rF := Fraction(uP)
L := LinearOrdinaryDifferentialOperator3(F, uP, rF)

You can think how to convert operator from
LinearOrdinaryDifferentialOperator1 to LinearOrdinaryDifferentialOperator3
and back.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

On 02/03/2015 01:07 AM, Ralf Hemmecke wrote:
>>> Wouldn't List(L) or Union(L, Record(ltilde: L, r: L)) be better? Or
>>> maybe Record(ltilde: L, r: Union(L, failed))?
>> I don't know, actually. Union(L, List(L)) seemed to best fit the
>> description of the function as mentioned in the paper.
>> However, Union(L, Record(ltilde: L, r: L)) also seems good as it offers
>> extra information about what the 2 elements actually mean in the latter
>> case.
> Now suppose you use this function and you are only interested in the
> operator, i.e. you don't care whether or not there is a corresponding r.
>
> You would then call
>
> z: Union(L, Record(ltilde: L, r: L)) := integrate_sols(...)
> l: L := if z case L then z::L else (z::Record(ltilde: L, r: L)).ltilde
>
> That is why I suggested
>
> Record(ltilde: L, r: Union(L, failed)).
>
> It also better says what actually fails if it does fail.

I've changed the return type to the above now.

> Actually just having
>
> List L
>
> would work similarly, but I somehow don't like to return a list if it
> can only have on or two elements. That's somehow misuse.

Yes, even I feel that way. That's why I was initially using Union(L,
List L) instead of List L.
> Ralf
>
Thanks,
Abhinav.

[Abhinav Baid]

On 02/03/2015 09:29 AM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>> How is the current implementation? [1]
>> In RationalLODE, we could just use integrate_sols(l, ratDsolve) with it.
>>
>> [1] https://github.com/fandango-/spad/blob/master/is.spad
>
> Two remarks:
>
> - actually we do not need D. Namely D is just a shorter way to get
> 'monomial(1, 1)$L'.
> - please keep line length below 80
>

Thanks for the remarks. I've made changes to the code according to them.

[Abhinav Baid]

I've tried to write functions for these conversions here [1].
Are they fine?

[1] https://github.com/fandango-/spad/blob/master/task2.spad

Thanks,
Abhinav.

[Waldek Hebisch]

You put GPL heare in the file. Note that FriCAS uses modified BSD
licence, see:

http://fricas.sf.net/copyright.txt

or src/etc/copyright.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Waldek Hebisch]

Close. But compare:

apply(D()$L1, 0, x)

and

apply(D()$L3, 0, x)

We want to have the same operator, so applying result of convertion
should give the same result as original operator.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

On 02/05/2015 03:00 AM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>>
>> On 02/03/2015 09:29 AM, Waldek Hebisch wrote:
>>> Abhinav Baid wrote:
>>>> How is the current implementation? [1]
>>>> In RationalLODE, we could just use integrate_sols(l, ratDsolve) with it.
>>>>
>>>> [1] https://github.com/fandango-/spad/blob/master/is.spad
>>>
>>> Two remarks:
>>>
>>> - actually we do not need D. Namely D is just a shorter way to get
>>> 'monomial(1, 1)$L'.
>>> - please keep line length below 80
>>>
>> Thanks for the remarks. I've made changes to the code according to them.
> You put GPL heare in the file. Note that FriCAS uses modified BSD
> licence, see:
>
> http://fricas.sf.net/copyright.txt
>
> or src/etc/copyright.
>

I've removed the GPL license.

[Abhinav Baid]

Oh, I should've paid closer attention to the warning you gave before
(about LODO3 having the output as plain D).
I've changed it now. [1] Is this okay?

[Ralf Hemmecke]

On 02/06/2015 09:59 AM, Abhinav Baid wrote:
> [1] https://github.com/fandango-/spad/blob/master/task2.spad

Hmmmm.... it says.

convertOp : L3 -> L1
++ convertOp(l3) converts an operator from LODO3 to LODO1.

Wouldn't be the following implementation meet that specification?

convertOp(l3 : L3) == 0$L1

Right. It does not, because it's not exactly clear what LODO3 and LODO1
actually are. I prefer long names instead of abbreviations. And
furthermore parameters are missing in the specification.
You specification also misses the exact input/output conditions.

Ralf

[Waldek Hebisch]

>
> On 02/06/2015 09:59 AM, Abhinav Baid wrote:
> > [1] https://github.com/fandango-/spad/blob/master/task2.spad
>
> Hmmmm.... it says.
>
> convertOp : L3 -> L1
> ++ convertOp(l3) converts an operator from LODO3 to LODO1.
>
> Wouldn't be the following implementation meet that specification?
>
> convertOp(l3 : L3) == 0$L1

No. The word 'convert' means that value is preserved.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Ralf Hemmecke]

On 02/06/2015 12:35 PM, Waldek Hebisch wrote:
> No. The word 'convert' means that value is preserved.

;-) If people convert from catholicism to buddhism or any other religion
or direction, what is the preserved value?

Anyway. It's not clear. What's the "preserved value". Abhinav would
otherwise have come up with the right conversion in the first place.
Because here you probably mean that for all r from R:

L(r) = convertOp(L)(r)

(for both directions L1->L3 and L3->L1).

But as long as that is not specified that conversion is useless, because
everyone would have to guess what it actually does in order to actually
use it in some situation.

Just my 2 cents.
Ralf

[abhinav baid]

Um, is it better now?

> Ralf
>
> --
> You received this message because you are subscribed to a topic in the
> Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/fricas-devel/2XuIRtc981E/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> fricas-devel...@googlegroups.com.
> To post to this group, send email to fricas...@googlegroups.com.
> Visit this group at http://groups.google.com/group/fricas-devel.
> For more options, visit https://groups.google.com/d/optout.
>


--
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> Oh, I should've paid closer attention to the warning you gave before
> (about LODO3 having the output as plain D).
> I've changed it now. [1] Is this okay?
>
> [1] https://github.com/fandango-/spad/blob/master/task2.spad
>

Good. Just one more remark: to convert from 'xD' to 'D' we need to
divide by 'x'. Probably the simplest way to ensure that we
can do division is to assume that R is a field. Your
current code will fail if 'x' is not invertible in R. If
you want such behaviour, then IntegralDomain is not needed.
In Ring there is 'recip' which will return "failed" if 'x'
is not invertible and give you inverse otherwise.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Thanks for your remark. I've made the change.

[Abhinav Baid]

Hi all,

Is there something else which I could do now - maybe related to the
project itself? I've been looking for some small functions which could
be useful for the factoring routine but aren't yet implemented without
much success. So, I've limited myself to reading the thesis. Any other
code sample ideas that I could work on?

Thanks,
Abhinav.

[Waldek Hebisch]

Maybe routine for point 5.3.2 (page 93) in the thesis? That
is given nontrivial solution to

fr + lf = 0

find a right factor of f.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Um, I have a doubt here. I feel very silly right now, but I'm not able to finish this routine.
The problem I'm having is this:
I'm able to find the basis of formal solutions at a regular point of f, but after that the thesis mentions computing the matrix of map r in this basis. I think I'm not quite sure about what this means.
For concreteness, take the example mentioned in the thesis
f=∂4+6x∂3+2(x2-1)x4∂2-2(3x2-1)x5∂+1x8f = \partial^4 + \frac{6}{x}\partial^3 + \frac{2(x^2-1)}{x^4}\partial^2 - \frac{2(3x^2-1)}{x^5}\partial + \frac{1}{x^8}
r=−x5∂3-x4∂2+2x3∂+x∂r = \minusx^5\partial^3-x^4\partial^2+2x^3\partial+x\partial
Using these, I'm able to find the basis of formal solutions at x=1 as
       1  4   1  3   1  2   1     23  1  4   2  3    2   1     5  1  2       1
   [- -- x  + - x  - - x  + - x + --, - x  - - x  + x  + - x - -, - x  - x + -,
      24      6      4      6     24  6      3           3     6  2          2
      1  4   7  3     2   3      5
    - - x  + - x  - 2x  + - x - --]
      4      6            2     12
but am not able to figure out what should I do to make operator r act on this basis.
Just applying it to every element in the list and taking the coefficients doesn't work as the matrix mentioned in the thesis is 4x4. Could you please help me? What should I be do to the basis using r so as to get the 4x4 matrix?

Thanks,
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 02/18/2015 10:57 AM, Waldek Hebisch wrote:
> >
> > Maybe routine for point 5.3.2 (page 93) in the thesis? That
> > is given nontrivial solution to
> >

> > fr + lf =3D 0

> >
> > find a right factor of f.
> >
> Um, I have a doubt here. I feel very silly right now, but I'm not able
> to finish this routine.

Doubts are OK -- just ask here.

> The problem I'm having is this:

> I'm able to find the basis of formal solutions at a regular point of f,=20
> but after that the thesis mentions computing the matrix of map r in this=20

> basis. I think I'm not quite sure about what this means.
> For concreteness, take the example mentioned in the thesis

> f = \partial^4 +
> \frac{6}{x}\partial^3 + \frac{2(x^2-1)}{x^4}\partial^2 -
> \frac{2(3x^2-1)}{x^5}\partial + \frac{1}{x^8}

> r= -x^5\partial^3-x^4\partial^2+2x^3\partial+x\partial

> Using these, I'm able to find the basis of formal solutions at x=1 as
> 1 4 1 3 1 2 1 23 1 4 2 3 2 1 5

> [- -- x + - x - - x + - x + --, - x - - x + x + - x - -,

> 24 6 4 6 24 6 3 3 6

<snip>

> but am not able to figure out what should I do to make operator r act on
> this basis.
> Just applying it to every element in the list and taking the
> coefficients doesn't work as the matrix mentioned in the thesis is 4x4.
> Could you please help me? What should I be do to the basis using r so as
> to get the 4x4 matrix?

The basis got mangled in the mail but it looks wrong. First, you
should have four solutions and it looked as what you posted had
only three. Second, convenient form of basis is such initial
condition is of form (0, ..., 1, 0, ..., 0), that is i-th
element of basis (counting from 0) has i-th derivative at x=1
equal to 1 and other derivatives of order < d (d = 4 in this example)
equal to 0.

I do not know how you produced the basis above, but for this
task convenient routine is 'ode' from UTSODE. You may look
at 'candidates' in ODERAT to see how it is used.

For experimantation 'seriesSolve' may be more convenient:

op := D(y(x), x, 4) + (6/x)*D(y(x), x, 3) + 2*(x^2 - 1)/x^4*D(y(x), x, 2) - 2*(3*x^2 - 1)/x^5*D(y(x), x) + 1/x^8*y(x)

(85)
8 (iv) 7 ,,, 6 4 ,, 5 3 ,
x y (x) + 6x y (x) + (2x - 2x )y (x) + (- 6x + 2x )y (x) + y(x)

-----------------------------------------------------------------------
8
x
Type: Expression(Integer)

b1 := seriesSolve(op, y, x = 1, [1, 0, 0, 0])
Compiling function %CS with type List(UnivariateTaylorSeries(
Expression(Integer),x,1)) -> UnivariateTaylorSeries(Expression(
Integer),x,1)

1 4 7 5 7 6 7
(86) 1 - -- (x - 1) + -- (x - 1) - -- (x - 1) + O((x - 1) )
24 60 30
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

xx := taylor(x, x = 1)

(20) 1 + (x - 1)
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

apr(t) == -xx^5*D(t, x, 3) - xx^4*D(t, x, 2) + (2*xx^3 + x)*D(t, x)

apr(b1)

(88)
3 2 13 3 37 4 171 5
(x - 1) - - (x - 1) + -- (x - 1) - -- (x - 1) + --- (x - 1)
2 6 12 40
+
1393 6 7
- ---- (x - 1) + O((x - 1) )
240
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

As you see first four coefficients (that is 0, 1, -3/2, 13/6)
give you the first column of the matrix. Similarly

b2 := seriesSolve(op, y, x = 1, [0, 1, 0, 0])

1 4 11 5 137 6 7
(89) (x - 1) + - (x - 1) - -- (x - 1) + --- (x - 1) + O((x - 1) )
6 40 360
Type: UnivariateTaylorSeries(Expression(Integer),x,1)
(90) -> apr(b2)

(90)
1 2 5 3 23 4 159 5
3 + 3(x - 1) + - (x - 1) - - (x - 1) + -- (x - 1) - --- (x - 1)
2 3 8 40
+
1243 6 7
---- (x - 1) + O((x - 1) )
240
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

Now coefficient give you second column. Next:

b3 := seriesSolve(op, y, x = 1, [0, 0, 2, 0]);
apr(b3)

(92)
2 1 3 3 4 61 5
- 2 - 2(x - 1) + 2(x - 1) + - (x - 1) - - (x - 1) + -- (x - 1)
3 4 60
+
53 6 7
- -- (x - 1) + O((x - 1) )
45
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

b4 := seriesSolve(op, y, x = 1, [0, 0, 0, 6]);

(94)
2 3 41 4 5 741 6
- 6 - 3(x - 1) + 7(x - 1) - -- (x - 1) + 14(x - 1) - --- (x - 1)
4 40
+
7
O((x - 1) )
Type: UnivariateTaylorSeries(Expression(Integer),x,1)

Note that 'seriesSolve' (and 'ode') takes values of derivatives,
while above I (and van Hoej) used coefficients, so I had to specifiy
factorial(i) as values of i-th derivative. Alternatively, one
can choose basis so that derivatives are 1, then one would have
to divide i-th coefficient by factorial(i). This would lead to
different matrix, but the same eigenvalues.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Hi Waldek,

Thanks for the detailed explanation. Actually, there were 4 solutions in
the basis, but they indeed got mangled in the post. I did use the 'ode'
function to produce the basis above. But, instead of a Taylor series
centered at 1, I converted the first 4 terms to UnivariatePolynomial
using UTS2UP. Hence, I was confused as to what coefficients I should be
taking. I thought we had to take the coefficients of 1,x,x^2,x^3 after
applying r, but now understand that we should instead take the
coefficients of 1,(x-1),(x-1)^2,(x-1)^3. Thanks for clearing this up.
I've coded up the routine here [1].
Could you please review it?

Thanks,
Abhinav.

[1] https://github.com/fandango-/spad/blob/master/task3.spad

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I've coded up the routine here [1].
> Could you please review it?
>
> Thanks,
> Abhinav.
>
> [1] https://github.com/fandango-/spad/blob/master/task3.spad

Hmm, it look that you are still converting Taylor series to
polynomials and computing coefficients of polynomials. Did
you push correct version? For applying r to elements of
basis you will need a little routine which expands each
coefficient of r into series and computes combination
of coefficients of r times derivaties of basis element.

Another remark: 'zerosOf' may be quite expensive, because
it produces _all_ zeros. Big advantage of van Hoej methods
is that usually it is enough to work with single zero.
So use 'zeroOf' instead of computing all zeros and discarding
the other ones. Actually, before computing zero you should
factor the characteristic polynomial and compute zero of
a single factor. For various reasons factorizer should be
passed as an extra parameter to 'find_right_factor'.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

On 02/22/2015 06:47 PM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>> I've coded up the routine here [1].
>> Could you please review it?
>>
>> Thanks,
>> Abhinav.
>>
>> [1] https://github.com/fandango-/spad/blob/master/task3.spad
> Hmm, it look that you are still converting Taylor series to
> polynomials and computing coefficients of polynomials. Did
> you push correct version? For applying r to elements of
> basis you will need a little routine which expands each
> coefficient of r into series and computes combination
> of coefficients of r times derivaties of basis element.

Oh, sorry. I convert the coefficients of r to Taylor Series now. [1]

> Another remark: 'zerosOf' may be quite expensive, because
> it produces _all_ zeros. Big advantage of van Hoej methods
> is that usually it is enough to work with single zero.
> So use 'zeroOf' instead of computing all zeros and discarding
> the other ones.

Yes, I was aware of the zeroOf function. But for the example mentioned
in the thesis, the zeroOf function wasn't giving a radical solution even
though 2 existed (repeated twice), but zerosOf was. So, I decided to go
with zerosOf.

> Actually, before computing zero you should
> factor the characteristic polynomial and compute zero of
> a single factor. For various reasons factorizer should be
> passed as an extra parameter to 'find_right_factor'.
>

Okay, I now take the factorizer as an extra parameter. Surprisingly,
now, zeroOf returns a radical result, so I use that as well.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 02/22/2015 06:47 PM, Waldek Hebisch wrote:
>
> > Hmm, it look that you are still converting Taylor series to
> > polynomials and computing coefficients of polynomials. Did
> > you push correct version? For applying r to elements of
> > basis you will need a little routine which expands each
> > coefficient of r into series and computes combination
> > of coefficients of r times derivaties of basis element.
> Oh, sorry. I convert the coefficients of r to Taylor Series now. [1]
> > Another remark: 'zerosOf' may be quite expensive, because
> > it produces _all_ zeros. Big advantage of van Hoej methods
> > is that usually it is enough to work with single zero.
> > So use 'zeroOf' instead of computing all zeros and discarding
> > the other ones.
> Yes, I was aware of the zeroOf function. But for the example mentioned
> in the thesis, the zeroOf function wasn't giving a radical solution even
> though 2 existed (repeated twice), but zerosOf was. So, I decided to go
> with zerosOf.

I see. However, note that in algebraic case you may easily get
equal but differently looking answers.

> > Actually, before computing zero you should
> > factor the characteristic polynomial and compute zero of
> > a single factor. For various reasons factorizer should be
> > passed as an extra parameter to 'find_right_factor'.
> >
> Okay, I now take the factorizer as an extra parameter. Surprisingly,
> now, zeroOf returns a radical result, so I use that as well.

Well, using 'zeroOf' and 'zerosOf' without factoring may lead
to subtly wrong answers. In this case 'zerosOf' compensated
for lack of factoring, but that was just an accident due to
simple form of data. 'zerosOf' will produce answer in terms
of radicals if polynomial is simple enough, otherwise it will
produce "general" root of polynomial which displays as a new
symbol and behaves in proper way in subseqent calculations.


> [1] https://github.com/fandango-/spad/blob/master/task3.spad

Now it looks OK. One extra remark: in other parts it is natural
to use factorizer which takes SparseUnivariatePolynomial(F).

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[abhinav baid]

Sorry, I'm not sure I understand. What do you mean by other parts?
Does it refer to other spad code which takes a factorizer as a
parameter? I could use the InnerEigenPackage code to get a
characteristic polynomial of type SparseUnivariatePolynomial(F) and
thus, change the factorizer domain now. Is that what you want to
suggest?

> --
> Waldek Hebisch
> heb...@math.uni.wroc.pl

>
> --
> You received this message because you are subscribed to a topic in the
> Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/fricas-devel/2XuIRtc981E/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> fricas-devel...@googlegroups.com.
> To post to this group, send email to fricas...@googlegroups.com.
> Visit this group at http://groups.google.com/group/fricas-devel.
> For more options, visit https://groups.google.com/d/optout.
>

Thanks,
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 2/26/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> > Abhinav Baid wrote:
> >>
> >
> >> [1] https://github.com/fandango-/spad/blob/master/task3.spad
> >
> > Now it looks OK. One extra remark: in other parts it is natural
> > to use factorizer which takes SparseUnivariatePolynomial(F).
> >
> Sorry, I'm not sure I understand. What do you mean by other parts?

First step of computing exponential parts is to compute Newton
polynomials associated with slopes of Newton polygon and factoring
them.

> Does it refer to other spad code which takes a factorizer as a
> parameter? I could use the InnerEigenPackage code to get a
> characteristic polynomial of type SparseUnivariatePolynomial(F) and
> thus, change the factorizer domain now. Is that what you want to
> suggest?

Eventually yes. Just there is no rush to do this now. For such
change you would need to fetch new (developement) FriCAS version.
If you want to do this now, then fine. If you prefer to work
with releases that is fine too.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[abhinav baid]

On 3/1/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> On 2/26/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
>> > Abhinav Baid wrote:
>> >>
>> >
>> >> [1] https://github.com/fandango-/spad/blob/master/task3.spad
>> >
>> > Now it looks OK. One extra remark: in other parts it is natural
>> > to use factorizer which takes SparseUnivariatePolynomial(F).
>> >
>> Sorry, I'm not sure I understand. What do you mean by other parts?
>
> First step of computing exponential parts is to compute Newton
> polynomials associated with slopes of Newton polygon and factoring
> them.
>

Thanks for the clarification.

>> Does it refer to other spad code which takes a factorizer as a
>> parameter? I could use the InnerEigenPackage code to get a
>> characteristic polynomial of type SparseUnivariatePolynomial(F) and
>> thus, change the factorizer domain now. Is that what you want to
>> suggest?
>
> Eventually yes. Just there is no rush to do this now. For such
> change you would need to fetch new (developement) FriCAS version.
> If you want to do this now, then fine. If you prefer to work
> with releases that is fine too.
>

I got the latest version and changed the code. I have the development
repository cloned anyway, so it wasn't a problem.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I got the latest version and changed the code. I have the development
> repository cloned anyway, so it wasn't a problem.

Good. How far are you with reading van Hoej thesis? You probaly
should start thinking about steps and subalgorithms needed for
factorization.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[abhinav baid]

On 3/4/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> I got the latest version and changed the code. I have the development
>> repository cloned anyway, so it wasn't a problem.
>
> Good. How far are you with reading van Hoej thesis?

I've done the first 2 chapters and skimmed through chapter 3 which
seems to contain the main algorithm (along with chapter 5). It seems
to be that chapter 4 is not required for the algorithm implementation.
Can it be skipped or am I missing something?

> You probaly
> should start thinking about steps and subalgorithms needed for
> factorization.

Yes, I was contemplating about making a list of all the functions that
would be needed.

[Waldek Hebisch]

> Abhinav Baid wrote:
>
> On 3/4/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> >
> > Good. How far are you with reading van Hoej thesis?
> I've done the first 2 chapters and skimmed through chapter 3 which
> seems to contain the main algorithm (along with chapter 5). It seems
> to be that chapter 4 is not required for the algorithm implementation.
> Can it be skipped or am I missing something?

Directly Chapter 4 is not needed for factorization. Its main
use is for finding elementary solutions and computing differential
Galois groups. Ideas from Chapter 4 may be useful for hardest
cases of factorization. Namely, main van Hoej method will
factor operator which have first order factor. It may find
higher order factor, but there is no warranty. So for higher
order factors van Hoej method should be combined with some
other method. However, for GSOC project van Hoej method
already is enough work. So yes, Chapter 4 can be skipped.

Actually, main part is in Chapters 2 and 3. This already
gives quite powerful method. Method from Chapter 5
is interesting because help factor operators which are
hard when using method from Chaper 3.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Hi Waldek,

Sorry for being inactive for so many days. My laptop screen broke and I
had a hard time finding a replacement. But now, it's all okay. As the
application period is quite near, I was thinking about starting to work
on the proposal. I'll also include the details of the sub-algorithms
necessary there itself. Is that fine?

Thanks,
Abhinav.

[Waldek Hebisch]

OK. Thanks for information (I was wondering what is happening).

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[abhinav baid]

I've submitted a proposal to Melange though it doesn't include a
timeline yet. Could you please review it and suggest any modifications
necessary?

Thanks,
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I've submitted a proposal to Melange though it doesn't include a
> timeline yet. Could you please review it and suggest any modifications
> necessary?

Few remarks:

- you write about reconstructing factors via Pade aproximation.
This is in fact called Hermite-Pade problem (Pade would be in case
of first order factor). FriCAS contains routine to solve
Hermite-Pade problem. In fact 'guessHolo' routine (maybe after
a small modification) should do.

- concerning 'nm_mult' routine I do not undertand what is the
purpose. If you need computations on operatiors written in
terms of xD, then the domain I provided and convertion routines
you wrote should do. In particular multiplication of operators
written in term of 'xD' should be just '*'.

- what 'ramification_of' should do? Is it the substitution
from end of page 19 in van Hoej thesis?

- could you say more precisely what 'coefs_operator' should do?

- in chapter 2 van Hoej writes about exponential parts, then
in chapter 3 about very similar notion of generalized exponents.
Later developements showed that generalized exponents are more
important. I think you should provide routine to compute
generalized exponents (it should probably just take part of
result returned by 'factor_riccati').

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[abhinav baid]

Hi Waldek,

Thanks for your remarks.

On 3/18/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> I've submitted a proposal to Melange though it doesn't include a
>> timeline yet. Could you please review it and suggest any modifications
>> necessary?
>
> Few remarks:
>
> - you write about reconstructing factors via Pade aproximation.
> This is in fact called Hermite-Pade problem (Pade would be in case
> of first order factor). FriCAS contains routine to solve
> Hermite-Pade problem. In fact 'guessHolo' routine (maybe after
> a small modification) should do.
>

Thanks for mentioning this. I added it to the proposal.

> - concerning 'nm_mult' routine I do not undertand what is the
> purpose. If you need computations on operatiors written in
> terms of xD, then the domain I provided and convertion routines
> you wrote should do. In particular multiplication of operators
> written in term of 'xD' should be just '*'.
>

Oh, yes. I've removed that function from the proposal now.

> - what 'ramification_of' should do? Is it the substitution
> from end of page 19 in van Hoej thesis?
>

Yes, it is. I've tried to make it more clear in the proposal.

> - could you say more precisely what 'coefs_operator' should do?
>

I've modified the proposal to reflect its purpose. Is it fine now?

> - in chapter 2 van Hoej writes about exponential parts, then
> in chapter 3 about very similar notion of generalized exponents.
> Later developements showed that generalized exponents are more
> important. I think you should provide routine to compute
> generalized exponents (it should probably just take part of
> result returned by 'factor_riccati').
>

Okay. I was confused about whether or not I should've added a separate
function for generalized exponents. But, as it seems to be important,
I've included it now.

Thanks,
Abhinav.

[abhinav baid]

Are there any other changes you'd like me to make - any corrections,
modifications or additions? Is the proposal okay right now? If yes,
then I can start with the timeline part as well.

Thanks,
Abhinav.

[Waldek Hebisch]

I am in a hurry now, so just quick remarks:

1) After my remarks you changed function list, but it you should
also make small changes to text before so that it matches
2) Some functions still needs clearer specification (I will write
later which)
3) You should think a bit about types of arguments and return
values. In particular I did not see any mention of streams
which are likely to be useful during lifting.
4) You should work on timeline without delay. Changes to earlier
parts are not going to be so large to invalidate timeline.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

I've tried to improve the proposal based on your remarks and also added
the timeline. Could you please review it now?

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I've tried to improve the proposal based on your remarks and also added
> the timeline. Could you please review it now?

1) You write that factor_riccati takes input that has only 1 slope.
Do you understand that algorithm 'Ricatti solution' may produce
in step 4 operator that has more than 1 slope. So you will
have to handle case of multiple slopes anyway.
2) Algorithm 'Ricatti solution' is enough to produce exponential
parts but not enough for generalized exponents. Namely
'Ricatti solution' looses information about generalized
exponents ei of form ei = ej + k where k is positive integer
and ej is another generalized exponent. That is why van
Hoej also have algorithm 'semi-regular parts' in 2.8.4
When we want only exponential parts some care is needed
to find multiplicities. This affects your milestone 1.

--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Hi Waldek,

On 03/27/2015 10:24 AM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>> I've tried to improve the proposal based on your remarks and also added
>> the timeline. Could you please review it now?
> 1) You write that factor_riccati takes input that has only 1 slope.
> Do you understand that algorithm 'Ricatti solution' may produce
> in step 4 operator that has more than 1 slope. So you will
> have to handle case of multiple slopes anyway.

Yes, but according to step 2, if the Newton polygon has more than 1
slope, then we compute a coprime index 1 factorization and apply
recursion to the right-hand factor. So, steps 3, 4, 5 always assume that
f has only one slope. Am I missing something here? If so, could you
please explain? Should I remove the restriction on f having only 1 slope?

> 2) Algorithm 'Ricatti solution' is enough to produce exponential
> parts but not enough for generalized exponents. Namely
> 'Ricatti solution' looses information about generalized
> exponents ei of form ei = ej + k where k is positive integer
> and ej is another generalized exponent. That is why van
> Hoej also have algorithm 'semi-regular parts' in 2.8.4
> When we want only exponential parts some care is needed
> to find multiplicities. This affects your milestone 1.
>

Yes, I'd thought about this but forgot to include it in the description.
I was thinking that since most of the steps in Riccati Solution and
semi-regular parts are the same, we could have some argument which
indicates which algorithm should be applied in the same function. Is
such an idea supported in FriCAS? Would it be better to split into 2 -
one for Riccati and the other for semi-regular?

Thanks,
Abhinav.

[Waldek Hebisch]

>
> Hi Waldek,
>
> On 03/27/2015 10:24 AM, Waldek Hebisch wrote:
> > Abhinav Baid wrote:
> >> I've tried to improve the proposal based on your remarks and also added
> >> the timeline. Could you please review it now?
> > 1) You write that factor_riccati takes input that has only 1 slope.
> > Do you understand that algorithm 'Ricatti solution' may produce
> > in step 4 operator that has more than 1 slope. So you will
> > have to handle case of multiple slopes anyway.
> Yes, but according to step 2, if the Newton polygon has more than 1
> slope, then we compute a coprime index 1 factorization and apply
> recursion to the right-hand factor. So, steps 3, 4, 5 always assume that
> f has only one slope. Am I missing something here? If so, could you
> please explain? Should I remove the restriction on f having only 1 slope?

Yes. In steps 3, 4, 5 f indeed has only one slope, but you need to
recurse from them to step 1. So step 1 is natural entry point and
multiple slopes are needed here.

> > 2) Algorithm 'Ricatti solution' is enough to produce exponential
> > parts but not enough for generalized exponents. Namely
> > 'Ricatti solution' looses information about generalized
> > exponents ei of form ei = ej + k where k is positive integer
> > and ej is another generalized exponent. That is why van
> > Hoej also have algorithm 'semi-regular parts' in 2.8.4
> > When we want only exponential parts some care is needed
> > to find multiplicities. This affects your milestone 1.
> >
> Yes, I'd thought about this but forgot to include it in the description.
> I was thinking that since most of the steps in Riccati Solution and
> semi-regular parts are the same, we could have some argument which
> indicates which algorithm should be applied in the same function. Is
> such an idea supported in FriCAS? Would it be better to split into 2 -
> one for Riccati and the other for semi-regular?

Thanks for explanation. Yes, you can merge both algorithms into
a single functiom.


--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Abhinav Baid]

Thanks for the prompt reply. I've removed the conditions on slope and
Newton polynomial and added a note about the fact that the 2 algorithms
will be merged in factor_riccati (which may be renamed).