Re: High Precision Approximate GCD

• Previous message: Paul Burridge: "Re: Confused by symbol in equation"
• In reply to: Dave Rusin: "Re: High Precision Approximate GCD"
• Next in thread: Z. Zeng: "Most recent papers (Re: High Precision Approximate GCD)"
• Messages sorted by: [ date ] [ thread ]

Date: Fri, 23 Jul 2004 23:40:39 GMT

Hi, Dave:

You are in DeKalb, aren't you? By the way, I was. Please say hi to Biswa.

The linear equation you proposed:

For given (i) polynomial p and q, (ii) the degree of u; solve linear system
    f*p+g*q = u for (f,g,u)
is a possible way. However, my experience is that this system can often be
poorly conditioned.

I met with one of the authors of EpsilonGCD and QuasiGCD about the the
failure cases. His explanation was that they set a tight criteria. If the
code suspect that an accurate GCD is in doubt, the code would output the
word "fail".

Back to the GCD problem. Computing GCD is a two stage process. The first
stage is to compute the degree of the GCD and an initial approximation. The
second stage is to refine the GCD iteratively. The problem eventually comes
down to solving a quadratic system. Namely, for give p and q, find (u,v,w)
such that

                   u*v = p
                   u*w = q

In short, this overdetermined system is nonsingular. Therefore, an ill-posed
problem (i.e., finding GCD) is converted to a well-posed one, which is
usually well conditioned. Therefore, an accurate solution can be achieved.

You may want to check out my paper "Computing multiple roots of inexact
polynomials" that is now available at Math. Comp. website

http://www.ams.org/journals/mcom/0000-000-00/

and will be on a printed issue soon. Anyway, the problem of computing
approximate GCD is now well solved.

Zhonggang Zeng

"Dave Rusin" <rusin@vesuvius.math.niu.edu> wrote in message
news:cdrn0n$pn1$1@news.math.niu.edu...
> In article <EGaMc.16317$8_6.5658@attbi_s04>, Z. Zeng <zzeng@neiu.edu>
wrote:
>
> >Let
> >
> >p = (x-1)^i*(x-2)^j*(x-3)^k*(x-4)^l
> >q = p'
> >
> >GCD(p,q) = (x-1)^(i-1) * (x-2)^(j-1) * (x-3)^(k-1) * (x-4)^(l-1)
> >
> >For different sets of [i,j,k,l], the following table lists the error in
> >computed approximate GCD's
> >using hardware precision (16 digits)
> >
> > [i,j,k,l] EpsilonGCD QRGCD UvGCD*
> >------------------------------------------------------
> > 2,1,1,0 Fail 1e-13 7e-16
>
>
> Let me see if I understand this: you have taken the polynomial p,
> presumably in expanded form x^4-7*x^3+17*x^2-17*x+6, and its
> derivative 4*x^3-21*x^2+34*x-17, and one of these algorithms
> _failed_ to find a gcd from these data?
>
> As I wrote privately to the OP, one can view the search for the
> gcd as a set of linear equations to be solved, which is certainly
> a pretty well-understood area of numerical analysis. In the OP's
> case, there was a clear understanding that the gcd was to be of
> degree 1 but the technique is similar whenever the degree of the
> gcd is known (obviously not always an easy thing to discover or
> even to define). We simply want to find polynomials P and Q of
> suitable degrees so that p P + q Q is a monic polynomial of
> the specified degree. In this example, a linear gcd requires
> deg(Q)=3 and deg(P)=2. So here it is in Maple:
>
> P:=add( a||i * x^(2-i), i=0..2):
> Q:=add( b||i * x^(3-i), i=0..3):
> total := expand( p*P + q*Q - (x^1 + junk) ):
> eqs:={ seq( coeff( total, x, i ), i= 1 .. 6 ) }:
> answer:=solve(eqs):
> taylor( subs(answer,p*P+q*Q) , x, 7 );
>
> The result is -1 + x , which is the correct answer. It is simple to
> adapt this code to other degrees for p and q if the degree of the
> gcd is also known.
>
> I don't deny that the "solve" command masks quite a bit of complexity
> which would be present when the the coefficients are floating-point
> values or when the polynomials are only _approximately_ equal to
> multiples of the gcd, but it is rather astounding to me that there
> would be software which would fail with such simple initial data.
> Could you elaborate? Moreover, could you explain why the naive approach
> I have suggested is not a good basis for solving problems like the OP had?
>
>
> FWIW, I tried more subtle cases like the ones the OP proposed. I randomly
> picked one root, picked five more with the kinds of variation s/he
proposed,
> and held to the condition that the linear coefficient be zero.
> Starting with roots r1=0.2491072884 and s1 = 0.2491079210, respectively,
> for the two polynomials, and using random leading coefficients as well,
> I was led to these two quartics:
> p:=0.4689144526*x^4 -0.3119802476*x^3
+.05837848017*x^2 -.0006056607140:
> q:=.09329956241*x^4 -.06206752682*x^3
+.01161291099*x^2 -.0001204532483:
> whose roots are, respectively,
> -0.08316541550, 0.2491371046, 0.2494694059, 0.2498833494
> -0.08315593602, 0.2491118087, 0.2491901349, 0.2501039165
> that is, the round-off error while computing the coefficients of the
> expanded forms of p and q had moved the roots around nearly as
> much as the distance to the other roots the OP was _not_ interested in.
> At this point I believe it is very much a judgement call to decide
> what the degree of the gcd is ! Insisting, as the OP did, that the gcd
> be of degree 1, the code I showed gives a gcd of approximately
> x - 0.2491426331
> when I stick to Maple's default precision for the calculations and
> blithely ignore the extra degree of freedom in the "solve" step.
> (There are ways to make a concrete choice but I don't know which is
"best".)
>
> dave

• Previous message: Paul Burridge: "Re: Confused by symbol in equation"
• In reply to: Dave Rusin: "Re: High Precision Approximate GCD"
• Next in thread: Z. Zeng: "Most recent papers (Re: High Precision Approximate GCD)"
• Messages sorted by: [ date ] [ thread ]