Re: vector version of gcdext?

[John Cremona <john.cremona@gmail.com>]

Thanks Bill, that is a clever trick to divide into two rather than take the factors one at a time.    It turned out that in my application, it was sufficient, and faster, to use the 2-termgcdext n times, with each irreducuble factor against the rest.

John

On Mon, 8 Oct 2018 at 20:21, Bill Allombert <Bill.Allombert@math.u-bordeaux.fr> wrote:

On Thu, Oct 04, 2018 at 10:31:58AM +0100, John Cremona wrote:
> On Wed, 3 Oct 2018 at 20:04, John Cremona <john.cremona@gmail.com> wrote:
>
> > Thanks for the replies.  I was planning to use this with polynomials in
> > one variable over a cyclotomic field.  I'll use mathnf.
>
> To answer Bill in a bit more detail.  I have a square matrix M with
> square-free char poly, f(x), factored as f=f1*f2*...*fn, for simplicity
> take n=2 so f=f1*f2*f3.  I need the extended gcd of h1=f2f3, h2=f1f3,
> h3=f1f2 and then writing 1=g1h1+g2h2+g3h3, substiituting M for x we get the
> identity as a sum of three orthgonal idempotent matrices, which we can then
> use to block-diagonalise M as a sirect sum of 3 blocks, each with char poly
> one of the fi.

So you want the dual version (you give [f1,f2,f3] and you get
[g1,g2,g3]).

dualmultigcdext(V)=
{
  my(n=#V);
  if(n==1,[[1],V[1]],
    my([A1,P1]=dualmultigcdext(V[1..n\2])
      ,[A2,P2]=dualmultigcdext(V[n\2+1..n])
      ,[u,v,d]=bezout(P1,P2));
    [concat(v*A1,u*A2),P1*P2/d])
}

Let V a vector, return [W,P] such that
sum W[i]/V[i] = 1/P and P = lcm V[i]

? V=primes(100);
? [W,P]=dualmultigcdext(V);
? W*[1/n|n<-V]~==1/P
%33 = 1

Cheers,
Bill.