Thanks for the review. While working on this, I came across the
following error:
Function: ?*? : (UnivariatePolynomial(x,Fraction(Integer)),%) -> % is
missing from domain: UnivariateLaurentSeries(Fraction(Integer),x,0)
Internal Error
The function * with signature $
(UnivariatePolynomial x (Fraction (Integer)))$ is missing from
domain
UnivariateLaurentSeries(Fraction (Integer))x(0 . 1)
I get this while multiplying li and ri (both of type LL) whose values are:
li := x^-1*(xD) - x^-2
ri := (xD) + 1
Any reason for this? I think this is related to the previous warnings:
not known that (LeftModule UP) is of mode (CATEGORY domain (SIGNATURE
coerce ($ (Variable var))) (SIGNATURE differentiate ($ $ (Variable
var))) (IF (has F (Algebra (Fraction (Integer)))) (SIGNATURE integrate
($ $ (Variable var))) noBranch))
Thanks,
Abhinav.
On 06/08/2015 03:26 AM, Waldek Hebisch wrote:
Abhinav Baid wrote:
Could you please look through the lift_newton and associated helper
functions that I've implemented? [1]
The implementation doesn't seem to be bug-free though as I get a "Heap
exhausted during garbage collection" error.
Also, I'm not sure if I'm calling lift_newton correctly and whether all
the required parameters are being passed. Could you please report any
errors you find so I can correct them?
ATM 'lift_newton' looks strange. What is the purpose of
first ('llaur') and last ('acc') parameter to 'lift_newton'?
And why do you have loop is 'factor_newton2' passing
coefficients of 'res1' and 'res2' to 'lift_newton'?
Lifting process is supposed to work on all coefficients
of the operator simultanously, so I do not see why you
pass single coefficient to 'lift_newton'.
Concerning 'nm_block': the 'round_off' stuff there looks
wrong. Also, why do you return a differential operator
from 'nm_block'?
In general 'lift_newton' looks more complicated than needed.
A single lifting step should look like:
-- compute error term
ei := f - li*ri
-- extract lower order part of error term to polynomial
pi := nm_block(ei, ...)
-- find correction terms
....
-- add correction terms to li and ri
lc := coefs_operator(...)
rc := coefs_operator(...)
li := li + lc
ri := ri + rc
As I wrote the 'ei := f - li*ri' part is going to take most of
execution time so later it will make sense to optimize it. But
ATM I strongly suggest to use the simplest form. You can test
code adding printouts, so print 'ei', 'pi', computed polynomial
correction terms, lc and rc. Run this in a loop say 10 times.
You will see if valuations of 'ei' increase and if all variables
have expected values. Once that works you can change factor_newton2
to return operators. More precisely, lifting step should
return stream of [lc, rc] pairs. You will need extra routine
to convert such stream into pair of operators. In fact, from
stream of [lc, rc] pairs you will get streams of coefficients
and from them you cab buid Lauren series and operators.
But at first write a loop that just prints results of few
lifting steps.
Concerning helper routines: I do not see why would you need
both 'coefs_operator' and 'coefs_operator1'. Namely,
given polynomial, slope and a starting point there is
essentially one way to make operator. More precisely,
natural convention is that with slope = n/d increase
in degree of polynomial by 1 corresponds to increase of
order of operator by d (and increase of degree of coefficient
by n). You seem to use different convention, where
increase of degree of polynomial by 1 corresponds to
increase of order by 1. Using this convention you must
ensure that degrees of all nonzero terms of polynomial
are divisible be d. More precisely, you may allow some
shift so differences of degrees of terms must be
divisible by d. But I see no reason to allow extra
shifts. AFAICS the first convention is better, because
you get polynomials of lower degree and there is no
risk of getting extra terms of wrong degree.
Regardless of convention, 'floor' or 'ceiling' in
'coefs_operator' looks wrong. Also I do not see why
you want to pass differential operator to 'coefs_operator'.
To put it differently, you have two basic types, namely
differential operators and polynomials. You need a function
which extracts part of given valuation from operator and
produces a polynomial (AFAICS your 'nm_block' should do
this) and a function which takes a polynomial and
produces a differential operator with given valuation.
But functions needs extra argument which indicates which
term of differential operator corresponds to constant
polynomial (sometimes this is clear, but sometimes you
will need an extra computation for this). Note that
'nm_block' and 'coefs_operator' are supposed to be
inverses, namely 'nm_block(coefs_operator(p, ...). ...)'
should give you back 'p'. 'coefs_operator(nm_block(op, ...),...)'
should give you back part of 'op' of given valuation.
So you can test them independently of other routines.
Currently you need no more functions to "convert" between
operators and polynomial. Trying to add some extra
similar function is likely lead to confusion and bugs.
Also, keeping convention that 'coefs_operator' takes
polynomial as first argument and produces operator
and 'nm_block' is doing the opposite means that
type checker will find errors for you: with such
convention logic error have good chance to lead to
type error.
Some extra remarks: for testing it makes sense to export
functions like 'coefs_operator' and 'nm_block' which
in final version will be local. Also, try to use better
variable names: 'res1' and 'res2' tell nothing to reader
while 'l' and 'r' explaim their roles. Or maybe 'lfac'
or 'lop' or 'li' to distingish from polynomials.
One extra remark: in 'factor_newton2' you recompute
Neton polygon. Passing 'k' to 'factor_newton2' would
look clearer.
--
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 [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
