Re: [fricas-devel] Factorization of linear ordinary

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I've added all the functions required to implement factor, but don't use
> guessHolo yet as I couldn't get how to use it for a list of power
> series, instead of a list of coefficients, and so use some custom PadĂŠ
> functions which may be quite inefficient as they make use of Polynomial
> data type [1]. The factor function does work for some input operators,
> but for others, it seems to be stuck inside the loop in
> try_factorization. Could you please see if there's some mistake and
> suggest what changes I'll have to make?

Some remarks:

1) In factor_global you unconditionally add infinity as a singularity.
However, it may happen that infinity is a regular point.
2) In factor_global you use 'zerosOf' to find singularities.
However, given irreducible factor of denominator its zeros
can not be distingushed using algebraic operations, so
normally it is enough to work with a single zero (the result
for other are conjugate via action of Galois group). Also
'zeroOf' tries to express zero in terms of radicals. But
this may lead to troubles here, so it is better to use
'rootOf'.
3) Computing 'lcm' and then factoring is likely to be more
expensive than using already known factors. In particular,
you can use 'gcdBasis' routine to produce list of relatively
prime factors and factor each separately. Also, it makes
sense to use 'squareFree' first to simplify computations.
4) In 'compute_bound' you seem to ignore ramified exponential
parts. This looks wrong.
5) Infinity is used differently in bounds: other singularities
bound denomiantors of a_i. Infinity bounds differences
between degrees of mumerators and denominators of a_i (so
given bound on denominators infinity gives bound for degree
of numarators).

Concerning solving Hermite-Pade problem: if you have a solution
at your disposal (which is easy to obtain from a factor), then
'guessHolo' should be quite efficient. More precisely, given
first order factor delta - r, the solution is exponential of
the integral of r/x which is an easy power series computations.
Using substitution S_e you can do this with Taylor series.
Given Taylor series is solution 'guessHolo' will produce the
factor (of course then you need to do S_{-e} to undo effect of S_e).
The above is a bit different from what van Hoej wrote.
If you want to follow van Hoej then you need procedure to solve
system

a_0v_0 + a_1v_1 + ... + a_dv_d

where v_i are vectors of power series. In case of vectors
of dimenion 1 (which correspond to factor of order 1) you
can use 'guessAlgDep' with degree bound of 1 to solve
such problem. Vectors of arbitrary dimension can be handled
by underlying routines, but I need to check what is the
best to call them.

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

[Abhinav Baid]

I've changed the file according to the above comments.
 

4) In 'compute_bound' you seem to ignore ramified exponential
   parts.  This looks wrong.

Sorry, I don't get what I'm missing here. First, I check that the degree of 'ram' = 1 and then only the constant coefficient of 'expart' matters in further computation. So, I think it should work?
 

5) Infinity is used differently in bounds: other singularities
   bound denomiantors of a_i.  Infinity bounds differences
   between degrees of mumerators and denominators of a_i (so
   given bound on denominators infinity gives bound for degree
   of numarators).
 

Again, I apologize because I think I don't understand this. What are the a_i in question?
 

Concerning solving Hermite-Pade problem: if you have a solution
at your disposal (which is easy to obtain from a factor), then
'guessHolo' should be quite efficient.  More precisely, given
first order factor delta - r, the solution is exponential of
the integral of r/x which is an easy power series computations.
Using substitution S_e you can do this with Taylor series.
Given Taylor series is solution 'guessHolo' will produce the
factor (of course then you need to do S_{-e} to undo effect of S_e).

So, can guessHolo can be used only with first order factors? If so, then how do I handle the case of m (as defined in section 3.6) > 1, I think the same goes for guessAlgDep as well.
 

The above is a bit different from what van Hoej wrote.
If you want to follow van Hoej then you need procedure to solve
system

a_0v_0 + a_1v_1 + ... + a_dv_d

where v_i are vectors of power series.  In case of vectors
of dimenion 1 (which correspond to factor of order 1) you
can use 'guessAlgDep' with degree bound of 1 to solve
such problem.  Vectors of arbitrary dimension can be handled
by underlying routines, but I need to check what is the
best to call them.

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

Thanks,
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> > Concerning solving Hermite-Pade problem: if you have a solution=20
> > at your disposal (which is easy to obtain from a factor), then=20
> > 'guessHolo' should be quite efficient. More precisely, given=20
> > first order factor delta - r, the solution is exponential of=20
> > the integral of r/x which is an easy power series computations.=20
> > Using substitution S_e you can do this with Taylor series.=20
> > Given Taylor series is solution 'guessHolo' will produce the=20
> > factor (of course then you need to do S_{-e} to undo effect of S_e).=20

> >
>
> So, can guessHolo can be used only with first order factors? If so, then
> how do I handle the case of m (as defined in section 3.6) > 1, I think the
> same goes for guessAlgDep as well.

guessHolo can be used when you have a solution for a factor. For
first order factor finding solutions is easier, but you can find
solutions also for higher order factors.

For solving Hermite-Pade problem the following routine
should work. Solutions are given by columns of the
result. Note: this routine produces so called order
basis, so some columns are not solutions -- you
need to check which columns give you actual
solutions.


The routine uses a trick to encode vector Hermite-Pade
problem

\sum_i p_iv_i = 0

as

\sum_i p_i(z^m)f_i = 0


where f_i = pack(v_i) and pack(w) = \sum_{j=0}^{m-1} z^j(w_j(z^m))
that is we encode vector of polynomials w of dimension m as a single
polynomial

\sum_{j=0}^{m-1} z^j(w_j(z^m))

Then it uses general routine in FractionFreeFastGaussian.

----------------<cut here>----------------------------

)abbrev package VHPSOLV VectorHermitePadeSolver
VectorHermitePadeSolver : Exports == Implementation where
F ==> Expression(Integer)
SUP ==> SparseUnivariatePolynomial(F)
Exports ==> with
hp_solve : (List Vector(SUP), List(NonNegativeInteger)) -> Matrix SUP
+++ hp_solve(lv, eta) solves Hermite-Pade problem with degree
+++ bound eta
Implementation ==> add

FFFG ==> FractionFreeFastGaussian(F, SUP)

power_action(m : NonNegativeInteger
) : ((NonNegativeInteger, NonNegativeInteger, SUP) -> F) ==
(k, l, g) +-> DiffAction(k, m*l, g)$FFFG

hp_solve(lv, eta) ==
m := #first(lv)
lpp : List(SUP) := []
for v in lv repeat
#v ~= m =>
error "hp_solve: vectors must be of the same length"
pp : SUP := 0
for i in 1..m repeat
pp1 := multiplyExponents(v(i), m)
pp := pp + monomial(1, (i-1)::NonNegativeInteger)$SUP*pp1
lpp := cons(pp, lpp)
lpp := reverse!(lpp)
sumEta := reduce(_+, eta)
C : List(F) := [0 for i in 1..sumEta]
generalInterpolation(C, power_action(m), vector(lpp), eta)$FFFG


-----------------------<cut here>----------------------------

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

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On Wednesday, July 15, 2015 at 8:23:45 AM UTC+5:30, Waldek Hebisch wrote:
> >
> > 4) In 'compute_bound' you seem to ignore ramified exponential
> > parts. This looks wrong.
> >
>
> Sorry, I don't get what I'm missing here. First, I check that the degree of
>
> 'ram' = 1 and then only the constant coefficient of 'expart' matters in
> further computation. So, I think it should work?

When 'ram' = 1 this looks OK. But AFAICS ramified singularity
still may have generalized exponent contributing to constant term.
As an example consider operator in the example in the section
2.7. The operator has only ramified generalized exponents, but
the residue is nonzero. So in Lemma 27 in van Hoej thesis you
need to take into account all generalized exponents.

>
> > 5) Infinity is used differently in bounds: other singularities
> > bound denomiantors of a_i. Infinity bounds differences
> > between degrees of mumerators and denominators of a_i (so
> > given bound on denominators infinity gives bound for degree
> > of numarators).
> >
>
>
> Again, I apologize because I think I don't understand this. What are the
> a_i in question?

We are trying to find a factor R = \sum_i a_i\partial_x^i

so a_i are (unknown) coefficients of the factor. To prove
that operator is irreducible (which is one of possible results
of van Hoej algorithm) we need bounds on degree of denominators
and numerators of a_i. Only having such bounds we can conclude
that no solution to Hermite-Pade problem means that operator
is irreducible.

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

[Abhinav Baid]

On 07/17/2015 07:27 PM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>>> Concerning solving Hermite-Pade problem: if you have a solution=20
>>> at your disposal (which is easy to obtain from a factor), then=20
>>> 'guessHolo' should be quite efficient. More precisely, given=20
>>> first order factor delta - r, the solution is exponential of=20
>>> the integral of r/x which is an easy power series computations.=20
>>> Using substitution S_e you can do this with Taylor series.=20
>>> Given Taylor series is solution 'guessHolo' will produce the=20
>>> factor (of course then you need to do S_{-e} to undo effect of S_e).=20
>>>
>> So, can guessHolo can be used only with first order factors? If so, then
>> how do I handle the case of m (as defined in section 3.6) > 1, I think the
>> same goes for guessAlgDep as well.
> guessHolo can be used when you have a solution for a factor. For
> first order factor finding solutions is easier, but you can find
> solutions also for higher order factors.
>
> For solving Hermite-Pade problem the following routine
> should work. Solutions are given by columns of the
> result. Note: this routine produces so called order
> basis, so some columns are not solutions -- you
> need to check which columns give you actual
> solutions.

Thanks for the routine. I've added it to the file and changed
try_factorization accordingly.

[Abhinav Baid]

On 07/18/2015 12:08 AM, Waldek Hebisch wrote:
> Abhinav Baid wrote:
>> On Wednesday, July 15, 2015 at 8:23:45 AM UTC+5:30, Waldek Hebisch wrote:
>>> 4) In 'compute_bound' you seem to ignore ramified exponential
>>> parts. This looks wrong.
>>>
>> Sorry, I don't get what I'm missing here. First, I check that the degree of
>>
>> 'ram' = 1 and then only the constant coefficient of 'expart' matters in
>> further computation. So, I think it should work?
> When 'ram' = 1 this looks OK. But AFAICS ramified singularity
> still may have generalized exponent contributing to constant term.
> As an example consider operator in the example in the section
> 2.7. The operator has only ramified generalized exponents, but
> the residue is nonzero. So in Lemma 27 in van Hoej thesis you
> need to take into account all generalized exponents.

Okay. Is the current code fine, then?

>>> 5) Infinity is used differently in bounds: other singularities
>>> bound denomiantors of a_i. Infinity bounds differences
>>> between degrees of mumerators and denominators of a_i (so
>>> given bound on denominators infinity gives bound for degree
>>> of numarators).
>>>
>>
>> Again, I apologize because I think I don't understand this. What are the
>> a_i in question?
> We are trying to find a factor R = \sum_i a_i\partial_x^i
>
> so a_i are (unknown) coefficients of the factor. To prove
> that operator is irreducible (which is one of possible results
> of van Hoej algorithm) we need bounds on degree of denominators
> and numerators of a_i. Only having such bounds we can conclude
> that no solution to Hermite-Pade problem means that operator
> is irreducible.
>

Oh, I get the meaning of a_i now. But, in section 3.6, we are trying to
find a_i in k[x] and so, there'll be no x in the denominator. So, is
this about some other computation that I'm (not) doing in another
function? Also, speaking of bounds, could you please see if I'm using
'bound' and 'eb' correctly in try_factorization2?

Thanks,
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 07/18/2015 12:08 AM, Waldek Hebisch wrote:
> > Abhinav Baid wrote:
> >> On Wednesday, July 15, 2015 at 8:23:45 AM UTC+5:30, Waldek Hebisch wrote:
> >>> 4) In 'compute_bound' you seem to ignore ramified exponential
> >>> parts. This looks wrong.
> >>>
> >> Sorry, I don't get what I'm missing here. First, I check that the degree of
> >>
> >> 'ram' = 1 and then only the constant coefficient of 'expart' matters in
> >> further computation. So, I think it should work?
> > When 'ram' = 1 this looks OK. But AFAICS ramified singularity
> > still may have generalized exponent contributing to constant term.
> > As an example consider operator in the example in the section
> > 2.7. The operator has only ramified generalized exponents, but
> > the residue is nonzero. So in Lemma 27 in van Hoej thesis you
> > need to take into account all generalized exponents.
> Okay. Is the current code fine, then?

There is problem with computing trace. If is better to avoid
going via AlgebraicNumber here (in general we would like to
allow non-algebraic constants here) and the formula you use
look wrong. A simple routine to comput trace may look like:

get_trace(f : F, k : K) : F ==
min_pol := minPoly(k)
Sae := SimpleAlgebraicExtension(F, UP, min_pol)
fa := univariate(f, k, min_pol)
trace(reduce(fa)$Sae)$Sae

Note it computes trace with respect to one kernel passed
as the argument k. In general you need to call 'kernels'
to get list of kernels and compute trace for each algebraic
one. If there are nonalgebraic kernels you need to replace
them by rational numbers (almost all values are OK, you
only need to avoid division by 0).

Note: you need to do replacement consistently for all
singularities.

Also, van Hoej has "v'(e_i - e_j)" term in his formula. I
see nothing corresponding to this in your code.

> >>> 5) Infinity is used differently in bounds: other singularities
> >>> bound denomiantors of a_i. Infinity bounds differences
> >>> between degrees of mumerators and denominators of a_i (so
> >>> given bound on denominators infinity gives bound for degree
> >>> of numarators).
> >>>
> >>
> >> Again, I apologize because I think I don't understand this. What are the
> >> a_i in question?
> > We are trying to find a factor R = \sum_i a_i\partial_x^i
> >
> > so a_i are (unknown) coefficients of the factor. To prove
> > that operator is irreducible (which is one of possible results
> > of van Hoej algorithm) we need bounds on degree of denominators
> > and numerators of a_i. Only having such bounds we can conclude
> > that no solution to Hermite-Pade problem means that operator
> > is irreducible.
> >
> Oh, I get the meaning of a_i now. But, in section 3.6, we are trying to
> find a_i in k[x] and so, there'll be no x in the denominator. So, is
> this about some other computation that I'm (not) doing in another
> function?

We can write operator in monic form:

L = \partial_x + r

with rational r, or we can clear denomonators and write:

L = a_1\partial_x + a_0

with polynomial a_1 and a_0. In other places van Hoej assumes
monic operators, in particular bound are done for monic operators.
But in section 3.6 the Hermite-Pade part is done with a_i in k[x].
So a_m represents lcm of denominators of coefficients of monic
form.

BTW, I think you should now export function to compute generalized
exponents. One reason is that you need generalized exponents
in factor_global, so it makes sense to have functions that
computes them and returns in convenient form. Another reason
is that parts dealing with Chapter 2 of van Hoej thesis are
done now, so it would be good to do more testing. Function
giving generalized exponents would help in testing.
And of course generalized exponents are useful for other
purposes.


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

[Abhinav Baid]

How about the current code?

I've added an exported function to compute generalized exponents. It
takes an operator and a point p as inputs (if p ~= 0, then it is moved
to 0 and the result has x-p or 1/x if p = infinity in place of x). The
output is a list of equivalence class of generalized exponents up to
conjugation over Q((x)) [ecs, ecr, ect] where ecs is written in terms of
ecr=ect.

[Abhinav Baid]

I've found a bug in the existing code but am unable to pinpoint the
exact cause. Consider

qF := Expression(Integer)
uP := UP('x, qF)
rF := Fraction(uP)
l1 := LinearOrdinaryDifferentialOperator1(rF)
lF := LinearOrdinaryDifferentialOperatorFactorizerNew('x,0)
d1 := D()$l1
xx := x::uP::rF
x1 := xx^(-1)
fP := factorPolynomial$ExpressionFactorPolynomial(Integer,qF)
op1 := d1^2*(d1^3+xx^2*d1^2-x1^3*d1+(xx+1)/(xx-1))

The old factorizer returns the complete factorization

(11) -> factor(op1)$LODOF(qF,uP)

1 1 3 2 2 1 x + 1
(11) [- D - -,- D + -,D + x D - -- D + -----]
x x 3 x - 1
x
but the new one doesn't

(12) -> factor(op1,fP)$lF

(12)
[
4 5 4 5 2
5 2 4 4x - 1 3 3x - x + 6x - 6 2 - 4x - 12x +
24x - 12
D + x D + ------- D + ----------------- D +
----------------------- D
3 5 4 7 6 5
x x - x x - 2x + x
+
4
-----------------
3 2
x - 3x + 3x - 1
]

The problem here is that for the above, factor_minmult1 is called which
either produces a factorization and if not, then the operator is
supposed to be irreducible, i.e. it provably either factorizes or
returns the operator itself. However, for the above case, the result is
wrong as op1 clearly has non-trivial factors. Any pointers to what could
be the reason for this?

Thanks,
Abhinav.

[Waldek Hebisch]

AFAICS you should get a linear factor for adjoint corresponding
to exponential part 0 at x_0 = 0. Things to check:

- if you get correct factor in K((x))
- that you pass correct singularity and correct factor in K((x))
to factor_minmult1
- if bounds are OK

Basically the method I use is to put printouts in critical
places and check if computed values agree with expectation.
The ')trace' command may be useful, but freqently it
produces too much data.

BTW: You pass operators of order 2 to the old factorizer.
van Hoej method is supposed to handle them, but it seems
that your code ATM can not handle them. I would suggest
to firt make sure that order 2 works (the values involved
are simpler, so it is easier to check them for correctness).

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

[Waldek Hebisch]

In section 3.5 van Hoej gives conditions when his method is
guarenteed to find a factor. Basically, if there is a
factor of order d and a place p with d + 1 exponential
parts, then using his method at p will find a factor.
Your code tries to use place at infinity. 'gen_exp'
says that there are four exponential parts, so
in principle the method should work. But in general,
if you have reducible operator of order d, then there
must be factor of order at most d/2, so you need
[d/2 + 1] distinct exponential parts to conclude
that operator is irreducible.

So it seems that current problem is in factor_minmult1,
but you also need to take into account fact that
van Hoej may miss a factor if there is too few
exponential parts.

BTW: Factoring operators of order 2 seem to work. Please
disregard what I wrote about this before: I made a silly
mistake in input data.

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

[Abhinav Baid]

Thanks for the checklist. I've been trying to debug this and though this
has led to some bug fixes, the main problem still remains, namely that
some operators are not factored even though they are reducible. Here's
another example, similar to the previous one:

op2 := d1^2*(d1^3+xx^2*d1^2-x1^3*d1+1)

Here, the bug's even more visible because d1 is clearly a factor of the
above. Also, one common thing I've noted in the cases I've tested is
that factor_minmult1 fails when the factor in k((x)) passed to it has
degree > 1. So, for instance, for the operator mentioned in the previous
post op1, here are the singularity and factor passed:

[point= "infinity",

lpf =
5 - 3 4 - 3 1 3
1D + (- x + 10)D + (- 2x + 35 + O(x ))D
+
- 3 - 2 - 1 1 2
(- 2x - 2x - 2x + 48 + O(x ))D
+
- 3 - 2 - 1 1 - 2 - 1
(- x - 6x - 10x + 10 + O(x ))D - 4x - 12x - 24
- 40x
+
2
O(x )
,
dxt= 1]
2 3 4 2 3 4
D + (1 + x + O(x ))D + 1 + 2x + 2x + 4x + O(x )

Based on the above, the only possible factor that can now be generated
is the degree 3 one i.e. (d1^3+xx^2*d1^2-x1^3*d1+(xx+1)/(xx-1)).
However, that doesn't happen. This leads me to suspect whether the
correct singularity and factor are actually being passed. Do the above
arguments look reasonable? I'd like to remove the parts dealing with
skipped_factors. Is there some way in which I can change the loop in
lines 1208-1219 while retaining what it's supposed to do?

Thanks,
Abhinav.

[Abhinav Baid]

Okay, I figured out that I was passing the wrong singularity and factor
to factor_minmult1 and hence I changed the code so that factor_minmult1
is no longer used [1]. This allows both the operators mentioned
previously to factor. I now believe that the factorizer should work as
far as the methods from chapters 2 and 3 are concerned. Any comments or
bug reports are welcome.

Here's an example of an operator that the new factorizer works for but
the old one doesn't:

op := (d1^2+xx^3+x1^3)*(d1^2+xx^2-x1^3)

[1]
https://github.com/fandango-/fricas/commit/c721a400d6ea039a0691ee94ade3c772c3f402f0

[Abhinav Baid]

Hi all,

Just thought that I'd make a small post. I've not been able to work on
this much for the past week as college reopened. Waldek, have you gone
through the current code? Is it fine, or could you please point out
the errors?

Thanks,
Abhinav.

[Waldek Hebisch]

>
> Just thought that I'd make a small post. I've not been able to work on
> this much for the past week as college reopened. Waldek, have you gone
> through the current code? Is it fine, or could you please point out
> the errors?

In general it look OK. However, I noticed that you use
the same bounds for operator and its adjoint. Generalized
exponents of adjoint of L are related, but usually different
than generalized exponents of L. So, I would expect different
bound for adjoint. Using too small bound may lead to missing
a solution. Also, it is not clear if you can just pick
a basis element for a Hermite-Pade problem: columns of solution
matrix span solutions space. To prove that operator is
irreducible you need to check that no element of solution
space leads to a factor. Unfortunetely, currently 'fffg'
routine will always return several solutions. More precisely,
'fffg' only accepts bounds that lead to several solutions.
I will try to modify 'fffg' so that it accepts bounds leading
to smaller number of solutions (hopefully a unique solution).

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

[Abhinav Baid]

On 8/5/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
>>
>> Just thought that I'd make a small post. I've not been able to work on
>> this much for the past week as college reopened. Waldek, have you gone
>> through the current code? Is it fine, or could you please point out
>> the errors?
>
> In general it look OK. However, I noticed that you use
> the same bounds for operator and its adjoint. Generalized
> exponents of adjoint of L are related, but usually different
> than generalized exponents of L. So, I would expect different
> bound for adjoint. Using too small bound may lead to missing
> a solution.

Okay, I use a different bound for adjoint now. [1]

> Also, it is not clear if you can just pick
> a basis element for a Hermite-Pade problem: columns of solution
> matrix span solutions space. To prove that operator is
> irreducible you need to check that no element of solution
> space leads to a factor. Unfortunetely, currently 'fffg'
> routine will always return several solutions. More precisely,
> 'fffg' only accepts bounds that lead to several solutions.
> I will try to modify 'fffg' so that it accepts bounds leading
> to smaller number of solutions (hopefully a unique solution).
>

So, as the solutions are spanned by the columns of the matrix, I need
to find coefficients, not all zero, for each of them such that the
corresponding linear combination gives a dot product of zero with lv.
Is there a relevant function for this? I looked at expressIdealMember,
but it returns the trivial list of zeros.

> --
> 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/3MYCZGOyIOI/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.

[1] https://github.com/fandango-/fricas/commit/68a42e8464577563972b6ff6e95357ce6758a693

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 8/5/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> >>
> >> Just thought that I'd make a small post. I've not been able to work on
> >> this much for the past week as college reopened. Waldek, have you gone
> >> through the current code? Is it fine, or could you please point out
> >> the errors?
> >
> > In general it look OK. However, I noticed that you use
> > the same bounds for operator and its adjoint. Generalized
> > exponents of adjoint of L are related, but usually different
> > than generalized exponents of L. So, I would expect different
> > bound for adjoint. Using too small bound may lead to missing
> > a solution.
>
> Okay, I use a different bound for adjoint now. [1]

Good. AFAICS there is still serious problem: you pass too small
number of coefficients to Hermite-Pade solver. The rule is:
if you want to get back n coefficients you need to have at
least n equations. However, the number of equations is
related to order of approximation. In scalar Hermite-Pade
problem number of equations is the same as order of approximation.
So to get n coefficients you need approximation of order n,
which means that you need n terms for _each_ coefficient of
local factor. For example, if we seek factor of order 1
and degree bound is 10, that we need to determine 22 numbers
(11 for each coefficient of the factor). For this we need
22 terms of series. Note that missing coefficients are
effectively replaced by 0 (since we pass (list of vectors of)
polynomials to Hermite-Pade solver). So we are solving
wrong Hermite-Pade problem and conseqently we should
expect wrong solutions.

Additional remarks:
- If there is singularity with exponential
part of multiplicity 1, then failure of van Hoej algorithm
for this exponential part applied to operator and its
adjount means that operator is irreducible. In such
case there is no need to look at other singularities
or even other exponnetial parts. And if our implementation
fails to find factor for reducible operator this is bug
which should be fixed.
- In try_factorization2 you perform too many calls to
Hermite-Pade solver. Normally a single call per order
of factor should be enough (for possible exceptions
see remarks at the end)
- You have some redundant code, for example you call
'truncate' and then 'convertUTStoUP' which in turn
callse 'UTS2UP'. But 'UTS2UP' ignores unneeded terms
of the series, so call to 'truncate' is not needed.
- Could you try to use better variable names. 'i', 'j', 'k'
usually mean integer variables. Unless there is strong
reason (like trying to follow notation in a paper) please
do not use then for non-integer variables.

> > Also, it is not clear if you can just pick
> > a basis element for a Hermite-Pade problem: columns of solution
> > matrix span solutions space. To prove that operator is
> > irreducible you need to check that no element of solution
> > space leads to a factor. Unfortunetely, currently 'fffg'
> > routine will always return several solutions. More precisely,
> > 'fffg' only accepts bounds that lead to several solutions.
> > I will try to modify 'fffg' so that it accepts bounds leading
> > to smaller number of solutions (hopefully a unique solution).
> >
>
> So, as the solutions are spanned by the columns of the matrix, I need
> to find coefficients, not all zero, for each of them such that the
> corresponding linear combination gives a dot product of zero with lv.
> Is there a relevant function for this? I looked at expressIdealMember,
> but it returns the trivial list of zeros.

I have commited a new version of fffg which accepts requested
order of solution as an argument. Normally if requested order
is large enoung you will get single solution. If you look
for factor starting from factors of low order, then if you
take sufficiently many terms (that is order of Hermite-Pade
problem is large enough) you will get single solution. ATM
I do not know if there is a simple way to find out which
order is large enough. One way could be to increase number
of terms in Hermite-Pade when there are multiple solutions
and stop once we get only one.

Below is a modified version of 'hp_solve' which accepts
requested order of solution as an argument.

Note: New version of 'generalInterpolation' computes
defects of solutions, that is margin between degree
of solution and degree bound. If defect is negative
then given "solution" is really not a solution because
it does not satisfy degree bounds. 'hp_solve' rejects
them, so result of 'hp_solve' consists of true solutions
(in case there are no solution return value is matrix
with 0 columns).


----------------<cut here>-------------------

)abbrev package VHPSOLV VectorHermitePadeSolver
VectorHermitePadeSolver : Exports == Implementation where
F ==> Expression(Integer)
SUP ==> SparseUnivariatePolynomial(F)
Exports ==> with

hp_solve : (List Vector(SUP), List(NonNegativeInteger),
NonNegativeInteger) -> Matrix SUP
+++ hp_solve(lv, eta, K) solves Hermite-Pade problem with degree
+++ bound eta up to order K.

Implementation ==> add

FFFG ==> FractionFreeFastGaussian(F, SUP)

power_action(m : NonNegativeInteger
) : ((NonNegativeInteger, NonNegativeInteger, SUP) -> F) ==
(k, l, g) +-> DiffAction(k, m*l, g)$FFFG

hp_solve(lv, eta, K) ==
print("hp_solve"::OutputForm)
print(eta::OutputForm)

m := #first(lv)
lpp : List(SUP) := []
for v in lv repeat
#v ~= m =>
error "hp_solve: vectors must be of the same length"
pp : SUP := 0
for i in 1..m repeat
pp1 := multiplyExponents(v(i), m)
pp := pp + monomial(1, (i-1)::NonNegativeInteger)$SUP*pp1
lpp := cons(pp, lpp)
lpp := reverse!(lpp)

vd : Vector Integer := vector([ei::Integer for ei in eta])
C : List(F) := [0 for i in 1..K]
res1 := generalInterpolation(C, power_action(m), vector(lpp),
vd, K)$FFFG
print(vd::OutputForm)
n := #vd
n1 : NonNegativeInteger := 0
for i in 1..n | vd(i) >= 0 repeat
n1 := n1 + 1
print(n1::OutputForm)
res := new(n, n1, 0)$Matrix(SUP)
i1 : NonNegativeInteger := 0
for i in 1..n | vd(i) >= 0 repeat
i1 := i1 + 1
for j in 1..n repeat
res(j, i1) := res1(j, i)
res

------------------<cut here>-------------------
--
Waldek Hebisch
heb...@math.uni.wroc.pl

[Waldek Hebisch]

One more remark: ATM you should not call the old factorizer
from the new one. It only makes harder to find errors.
Just print message saying that given case is not handled
by the new factorizer.

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

[Abhinav Baid]

Yes, but how do I check this and appropriately call factor_minmult1?
(which is removed for now). The earlier check was using
skipped_factors but it was wrong. I think there should be a condition
in the loop after bounds are computed, but am not sure. Also, I have
fixed the code according to the other remarks.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 8/6/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> > Abhinav Baid wrote:
> >>
> > Additional remarks:
> > - If there is singularity with exponential
> > part of multiplicity 1, then failure of van Hoej algorithm
> > for this exponential part applied to operator and its
> > adjount means that operator is irreducible. In such
> > case there is no need to look at other singularities
> > or even other exponnetial parts. And if our implementation
> > fails to find factor for reducible operator this is bug
> > which should be fixed.
> Yes, but how do I check this and appropriately call factor_minmult1?
> (which is removed for now). The earlier check was using
> skipped_factors but it was wrong. I think there should be a condition
> in the loop after bounds are computed, but am not sure.

1) If you get no solution from Hermite-Pade solver, or if there
is single solution, but has no common factor with f, then
van Hoej method failed for given expomential part. In fact,
if you seek for factors increasing order by one, then
once you get a single reult from Hermite-Pade solver it
will be a factor of f or relatively prime to f. If there
are multiple solution then increasing number of terms in
expansion will eventually lead to single or no solution.
2) Checking for exponential part of multiplicity 1 is easy:
you need to check that there is generalized exponents
which is not equivalent to another one.

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

[Abhinav Baid]

Okay, I've added factor_minmult1 again now [1]. Is it fine?

> --
> 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/3MYCZGOyIOI/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.
>

[1] https://github.com/fandango-/fricas/commit/df5fb3e4ddc31a1678f94276957ffc64dc41dd15

--
Abhinav.

[Waldek Hebisch]

I am affraid that there are still problems. Quick remarks:

1) Why do you have so much code in 'factor_minmult1'? Two
calls to 'try_factorization' (one with operator and one
with adjoint) should be enough. If exponential part have
multiplicity 1 and both fail we know that operator is
irreducible. Note: you need to use 'degree(f)-1' as maximal
possible order of the factor, because we do not know
which factor will be found.
2) Code in 'try_factorization2' looks wrong. More precisely,
'sumEta' as order almost always will be wrong. Needed
order of solutions is related to value of 'acc':
if 'acc' is small than you have small number of terms
in input and normally you can not get solution accurate
to high order. OTOH, if 'acc' is large but requested
order is small, then you effectively waste time computing
many input terms, because Hermite-Pade solver will not
use them (and will generate _the same_ solution as if
you gave it smaller number of terms).
Also, you have a loop that looks like an attempt to
check solutions. After change that I made solutions
returned by Hermite-Pade solver will be accurate up
to requested order, so no need to check this. You
seem to expect exact solution (you check of "mr = 0"),
but that is rather unusual -- normally lv are truncations
of unfinite series which means that we will get nonzero
value around point of truncation.
AFAICS effect of two problems above is that
'try_factorization2' will frequently miss solutions.

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

[Abhinav Baid]

Thanks for the remarks.

On 8/12/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
>
> I am affraid that there are still problems. Quick remarks:
>
> 1) Why do you have so much code in 'factor_minmult1'? Two
> calls to 'try_factorization' (one with operator and one
> with adjoint) should be enough. If exponential part have
> multiplicity 1 and both fail we know that operator is
> irreducible. Note: you need to use 'degree(f)-1' as maximal
> possible order of the factor, because we do not know
> which factor will be found.

I've removed the last part of the code now. I actually make 4 calls to
try_factorization with maximal possible order as floor(degree(f) / 2)
and degree(f)-1 (once each with operator and adjoint), so that a lower
order factor is found first, if it exists. I've found this way to be
faster. Also, there's no repeated computation involved as I change the
min order to floor(degree(f) / 2) + 1 for the latter case.

> 2) Code in 'try_factorization2' looks wrong. More precisely,
> 'sumEta' as order almost always will be wrong. Needed
> order of solutions is related to value of 'acc':
> if 'acc' is small than you have small number of terms
> in input and normally you can not get solution accurate
> to high order. OTOH, if 'acc' is large but requested
> order is small, then you effectively waste time computing
> many input terms, because Hermite-Pade solver will not
> use them (and will generate _the same_ solution as if
> you gave it smaller number of terms).
> Also, you have a loop that looks like an attempt to
> check solutions. After change that I made solutions
> returned by Hermite-Pade solver will be accurate up
> to requested order, so no need to check this. You
> seem to expect exact solution (you check of "mr = 0"),
> but that is rather unusual -- normally lv are truncations
> of unfinite series which means that we will get nonzero
> value around point of truncation.
> AFAICS effect of two problems above is that
> 'try_factorization2' will frequently miss solutions.
>

Okay, I've changed the code again. Is this better?

--
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> On 8/12/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> >
> > I am affraid that there are still problems. Quick remarks:
> >
> > 1) Why do you have so much code in 'factor_minmult1'? Two
> > calls to 'try_factorization' (one with operator and one
> > with adjoint) should be enough. If exponential part have
> > multiplicity 1 and both fail we know that operator is
> > irreducible. Note: you need to use 'degree(f)-1' as maximal
> > possible order of the factor, because we do not know
> > which factor will be found.
> I've removed the last part of the code now. I actually make 4 calls to
> try_factorization with maximal possible order as floor(degree(f) / 2)
> and degree(f)-1 (once each with operator and adjoint), so that a lower
> order factor is found first, if it exists. I've found this way to be
> faster. Also, there's no repeated computation involved as I change the
> min order to floor(degree(f) / 2) + 1 for the latter case.

Looks OK.

> > 2) Code in 'try_factorization2' looks wrong. More precisely,
> > 'sumEta' as order almost always will be wrong. Needed
> > order of solutions is related to value of 'acc':
> > if 'acc' is small than you have small number of terms
> > in input and normally you can not get solution accurate
> > to high order. OTOH, if 'acc' is large but requested
> > order is small, then you effectively waste time computing
> > many input terms, because Hermite-Pade solver will not
> > use them (and will generate _the same_ solution as if
> > you gave it smaller number of terms).
> > Also, you have a loop that looks like an attempt to
> > check solutions. After change that I made solutions
> > returned by Hermite-Pade solver will be accurate up
> > to requested order, so no need to check this. You
> > seem to expect exact solution (you check of "mr = 0"),
> > but that is rather unusual -- normally lv are truncations
> > of unfinite series which means that we will get nonzero
> > value around point of truncation.
> > AFAICS effect of two problems above is that
> > 'try_factorization2' will frequently miss solutions.
> >
> Okay, I've changed the code again. Is this better?

Better, but I think there are still problems:

1) AFAICS, bound on degree of coefficients of factor is

cb := ceiling(bound(#lv-1))@Integer + eb)

cnt := #fl*(cb + 1) gives number of needed coefficients. So
to get single solution we need order [cnt/#first(fl)] where
[...] means rounding up. However, it makes sense use slightly
larger order, to decrease probability of wrong solutions. So
code to comute needed order may look like:

vdim := #first(fl)
acc := (cnt + vdim + 3) quo vdim

2) If ncols(hps) = 1 you call 'hp_solve' again, that is not
needed

Also, the following structure (pseudo code) would be simpler:

try_factorization2(fl, n, bound, eb) ==
....
acc := ...
repeat
lv := ...
hps := hp_solve(...
ncols(hps) = 1 => return members(column(hps, 1))
ncols(hps) = 0 => returm "failed"
-- increas acc
acc := ...

> --
> Abhinav.
>
> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, 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.
>

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

[Abhinav Baid]

Okay, I really hope the current code is alright [1].

[1] https://github.com/fandango-/fricas/commit/27f45f41a12fc5d30790e635a325262fdedaa94c

> 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/3MYCZGOyIOI/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:
>
> Okay, I really hope the current code is alright [1].
>
> [1] https://github.com/fandango-/fricas/commit/27f45f41a12fc5d30790e635a325262fdedaa94c

Looks OK. I have tried a few factorizarion examples. Most
works fine, but given:

qF := Expression(Integer)
uP := UP('x, qF)
rF := Fraction(uP)
l1 := LinearOrdinaryDifferentialOperator1(rF)
lF := LinearOrdinaryDifferentialOperatorFactorizerNew('x,0)
d1 := D()$l1
xx := x::uP::rF
x1 := xx^(-1)
fP := factorPolynomial$ExpressionFactorPolynomial(Integer,qF)
op1 := d1^2*(d1^3+xx^2*d1^2-x1^3*d1+(xx+1)/(xx-1))

I get:

(22) -> factor((d1^2+x1+xx)*op1, fP)$lF

>> Error detected within library code:
Denominator not equal to 1

The problem seem to be in lift_newton

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

[Abhinav Baid]

Yes, I was wrongly calling retract with shift instead of d*shift. I've
fixed it now [1].

[1] https://github.com/fandango-/fricas/commit/4f0457c967122b30577b2e77832f2735bde0f67c

On 8/13/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:

[Abhinav Baid]

Hi Waldek,

Are there any more problems as far as the code is concerned? If not,
then I can probably add documentation and make it merge-ready.

Thanks,
Abhinav.

--
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> Are there any more problems as far as the code is concerned? If not,
> then I can probably add documentation and make it merge-ready.

Examples that I tried work, including some which probably are
beyond capabilities of old factorizer. AFAICS there are
still few unimplemented things needed for "complete"
algorithm:

- polynomial factorizer should get list of algebraic kernels as
an argument. This is somewhat mechanical change, but may
force change of signatures of some (several??) routines,
- when there is only one exponential part we need to call
routine searching for exponential solutions (in principle
calling old factorizer here will work, but is likely to
be suboptimal).
- it would be nice to have eigenring method because it works
and should be efficient exactly in cases which are hard
for other approach.

However, since today is suggested 'pencils down' date all
of the above probably can wait past GSOC: the factorizer
you wrote is already quite capable.

Concerning merge, I think that it should be done in two
stages. First the new factorizer should be available
via explicit call to its routines. In the second stage
the new factorizer should be called as default factorizer,
whenever there is need to factor LODO-s. The second stage
will require filling the remaing cases (possibly via call
to old factorizer), profiling and comparizon of timings.
Clearly this will take time, so will happen after GSOC.
For the first stage I think we need some little cleanup.
First, AFAICS currently the only sensible value for the 'cen'
parameter is 0. It would be good to get rid of 'cen' as
parameter (this may be less trivial than it looks, because
there are some quirks in Spad compiler when handling constants
as type parameters). Second, it seems that you compute
several things (including generalized exponents) multiple
times, it would be good to compute things once and then
use computed result. In particular it would be nice
to have routine which computes all singularities and
generalized exponents for each singularity ('gen_exp'
works just for single singularity). Given list of
singularities and generalized exponents 'compute_bound'
could just use it (without calling 'factor_op'). And
'factor' could use the same list to decide at which singularity
to try van Hoej method.


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

[Abhinav Baid]

On 8/17/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> Are there any more problems as far as the code is concerned? If not,
>> then I can probably add documentation and make it merge-ready.
>
> Examples that I tried work, including some which probably are
> beyond capabilities of old factorizer. AFAICS there are
> still few unimplemented things needed for "complete"
> algorithm:
>
> - polynomial factorizer should get list of algebraic kernels as
> an argument. This is somewhat mechanical change, but may
> force change of signatures of some (several??) routines,
> - when there is only one exponential part we need to call
> routine searching for exponential solutions (in principle
> calling old factorizer here will work, but is likely to
> be suboptimal).
> - it would be nice to have eigenring method because it works
> and should be efficient exactly in cases which are hard
> for other approach.

Yes, I'll surely work on adding this later.

I've removed the cen parameter but this seems to create a warning
which doesn't seem to affect execution however:
convertL3toLL: not known that (LeftModule (UnivariateTaylorSeries
(Expression (Integer)) var (Zero))) is of mode (CATEGORY domain
(SIGNATURE coerce ($ (Variable var))) (SIGNATURE differentiate ($ $
(Variable var))) (IF (has (Expression (Integer)) (Algebra (Fraction
(Integer)))) (SIGNATURE integrate ($ $ (Variable var))) noBranch))
Is it fine to ignore this?
Also, I've added a function to compute all singularities and their
corresponding generalized exponents and use the return value in
factor_global and compute_bound. Anything more for the first stage?

--
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
>
> I've removed the cen parameter but this seems to create a warning
> which doesn't seem to affect execution however:
> convertL3toLL: not known that (LeftModule (UnivariateTaylorSeries
> (Expression (Integer)) var (Zero))) is of mode (CATEGORY domain
> (SIGNATURE coerce ($ (Variable var))) (SIGNATURE differentiate ($ $
> (Variable var))) (IF (has (Expression (Integer)) (Algebra (Fraction
> (Integer)))) (SIGNATURE integrate ($ $ (Variable var))) noBranch))
> Is it fine to ignore this?

Yes, it can be ignored.

> Also, I've added a function to compute all singularities and their
> corresponding generalized exponents and use the return value in
> factor_global and compute_bound. Anything more for the first stage?
>

It makes sense to export 'ge_minimal'. There is commented out
version of 'newtonpolygonPoints' -- it is better to remove it.
It would be good have a test file. Finally, the 'New' suffix
is likely to get obsolete sooner or later. I would suggest
name like LinearOrdinaryDifferentialOperatorFactorizer2 with
abbreviation LODOF2.

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

[Abhinav Baid]

On 8/21/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>>
>> I've removed the cen parameter but this seems to create a warning
>> which doesn't seem to affect execution however:
>> convertL3toLL: not known that (LeftModule (UnivariateTaylorSeries
>> (Expression (Integer)) var (Zero))) is of mode (CATEGORY domain
>> (SIGNATURE coerce ($ (Variable var))) (SIGNATURE differentiate ($ $
>> (Variable var))) (IF (has (Expression (Integer)) (Algebra (Fraction
>> (Integer)))) (SIGNATURE integrate ($ $ (Variable var))) noBranch))
>> Is it fine to ignore this?
>
> Yes, it can be ignored.
>

Well, seems like I spoke too soon. It actually affects execution.
Here's an example:

qF := Expression(Integer)
uP := UnivariatePolynomial('x, qF)
rF := Fraction(uP)
L1 := LinearOrdinaryDifferentialOperator1(rF)
lF := LinearOrdinaryDifferentialOperatorFactorizer2('x)
d := D()$L1

xx := x::uP::rF
x1 := xx^(-1)

fP := factor$ExpressionFactorPolynomial(Integer,qF)
op:=d^5+2*xx^3*d^3+6*xx^2*d^2+6*xx*d+xx^6*d+1
inf := "infinity" :: "infinity"
gen_exp(op,inf,fP,[])$lF

gives
>> System error:
The value 0 is not of type CONS.
The interesting thing here is that the result is actually getting
computed correctly. If I add a print statement before returning res in
gen_exp, I get
1 - 5 11 - 5 9 2 1
[[ecs= [0],ecr= x,ect= -],[ecs= [x + --,x + -],ecr= - x ,ect= -]]
x 4 4 x
as the result which is correct. Somehow, the interpreter signals an
error even when there is none. Any pointers on how to fix this?

>> Also, I've added a function to compute all singularities and their
>> corresponding generalized exponents and use the return value in
>> factor_global and compute_bound. Anything more for the first stage?
>>
>
> It makes sense to export 'ge_minimal'. There is commented out
> version of 'newtonpolygonPoints' -- it is better to remove it.
> It would be good have a test file. Finally, the 'New' suffix
> is likely to get obsolete sooner or later. I would suggest
> name like LinearOrdinaryDifferentialOperatorFactorizer2 with
> abbreviation LODOF2.
>

I've made changes according to the above. The test file has been added
to the src/input folder as lodof2.input. Also, I changed the
signatures of factor and gen_exp to allow the polynomial factorizer to
get a list of algebraic kernels as argument.

--
Abhinav.

[Waldek Hebisch]

This is rather old problem with Spad compiler: constant parameters
to types may get wrong values. The .spad compiles fine because
during compilation parameters are only compared and provided
value is consitently wrong (so values compare to equal where
it should). Most thing work OK because parameter is unused.
But printing routines use value of parameter and crash.

More preciely, due to bug power series domains get integer
zero as parameter, while the correct value is zero from
Expression(Integer).

Also printing type crashes. One can call 'gen_exp' like:

)set message type off
ge := gen_exp(op,inf,fP,[])$lF;

but attempt to print 'ge' crashes.

Concerning fixing the problem: proper fix would be to fix Spad
compiler, but this is a substantial change and will not happen
quickly. In existing code the MrvLimitPackage uses power
series with 0 parameter. To avoid the bug zero is assigned
to a variable ('zeroFE' in this case) and the variable is
passed as parameter. This is possible because power series
are only used inside MrvLimitPackage, no exported routine
has power series as (part of) argument or return value.

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

[Abhinav Baid]

Oh, but it makes sense to keep gen_exp exported right? Also, is there
a reason why something like UTS ==> UnivariateTaylorSeries(F, var,
0$F) is not supported. Given that you can take typed parameters and
then use them, doesn't it also make sense to allow specifying the type
of constant parameters you are using instead of the compiler doing
this work for you. Given how almost everywhere else, you have to
define the type of a variable or it's type should be clear from the
usage (function called, etc.), it seems to me that this is an
unnecessary anomaly. Are there cases when you want the compiler to
compute the type for you (even in this case, I think the usage I
mentioned above should be valid: many times you want the type to be
more general or specific than the one inferred, like when you use
retract or coerce for example), or is it a bug in the compiler? As
this problem seems to be unresolvable at present, is there anything
else that needs improvement?

--
Abhinav.

[Waldek Hebisch]

Yes. AFAICS generalized exponents can be expressed without
using series.

> Also, is there
> a reason why something like UTS ==> UnivariateTaylorSeries(F, var,
> 0$F) is not supported. Given that you can take typed parameters and
> then use them, doesn't it also make sense to allow specifying the type
> of constant parameters you are using instead of the compiler doing
> this work for you. Given how almost everywhere else, you have to
> define the type of a variable or it's type should be clear from the
> usage (function called, etc.), it seems to me that this is an
> unnecessary anomaly. Are there cases when you want the compiler to
> compute the type for you (even in this case, I think the usage I
> mentioned above should be valid: many times you want the type to be
> more general or specific than the one inferred, like when you use
> retract or coerce for example), or is it a bug in the compiler?

I consider this problem a compiler bug. But this is rather deep
one. As I wrote, fixing it probably will require substantial
changes to the compiler.

> As
> this problem seems to be unresolvable at present, is there anything
> else that needs improvement?

Well, crashing during some operations is unacceptable. So
the possibilities are:

- revert to having center as parameter

- or eliminate use of power series in exported signatures

P.S. Sorry for long silence, real life kept me busy with
other things.

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

[Abhinav Baid]

I

I went with option 1 as ge_minimal is used within unexported functions
as well [1]. Is that fine?

> P.S. Sorry for long silence, real life kept me busy with
> other things.
>

I'm in no hurry :)

[1] https://github.com/fandango-/fricas/commit/6476a153ab5afc55cfb3a4f6ba19de01e7d38d5e

--
Abhinav.

[Waldek Hebisch]

Abhinav Baid wrote:
> I went with option 1 as ge_minimal is used within unexported functions
> as well [1]. Is that fine?

Yes.

I have now commited the new factorizer -- I did a little
cosmetic changes to lodof2.spad and lodof2.input and adjusted
Makefile-s.

Thanks for your work. I hope that you will stay with us.

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

[Bill Page]

On 7 September 2015 at 20:33, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>> I went with option 1 as ge_minimal is used within unexported functions
>> as well [1]. Is that fine?
>
> Yes.

> ...

> Thanks for your work. I hope that you will stay with us.
>

+1

[Abhinav Baid]

Hi Waldek,

On 9/8/15, Waldek Hebisch <heb...@math.uni.wroc.pl> wrote:
> Abhinav Baid wrote:
>> I went with option 1 as ge_minimal is used within unexported functions
>> as well [1]. Is that fine?
>
> Yes.
>
> I have now commited the new factorizer -- I did a little
> cosmetic changes to lodof2.spad and lodof2.input and adjusted
> Makefile-s.

Awesome!

>
> Thanks for your work. I hope that you will stay with us.

Thank you for your valuable time, support and patience. This project
wouldn't have been possible without your help. I want to make the new
factorizer the default one, and hence, implement the eigenring method
as well. Also, the problem with the compiler not handling constants
properly is something I want to look at when I have time. So, I'll
definitely stay.

--
Abhinav.