/* PARI/GP program gfact for finding decomposition */
/* of Gaussian integer into Gaussian primes and units */
/* Michael Somos, Jan 08, 2003 */

/* vector from matrix and back */
vecmat(m)=vector(matsize(m)[1],n,m[n,])
{matvec(v)=local(m); m=if(v==[],[;],matrix(length(v),length(v[1])));
   for(n=1,length(v),m[n,]=v[n]);m}

/* return a+I*b where p=a^2+b^2 */
{g4n1(p)=
  if(p%4!=1,if(p==2,return(1+I),print("p="p);error("g4n1: not 4n+1")));
  forstep(pa=1,sqrtint(p),2,if(issquare(p-pa^2),return(pa+I*sqrtint(p-pa^2))))
} /* end g4n1 */

/* is z a normalized complex number? */
normal(z)=real(z)>0&imag(z)>=0    

/* Gaussian integer or rational factorization of z=a+b*I */
/* returns matrix of factorization ala factor() */
/* if flg zero also splits into Gaussian primes */
{gfact(z,flg=0)= local(t,e1,e2,nm,p,gp,fact);
  if(z==0|z==1|((flg!=0)&(type(z)!=type(I))),
     return(factor(simplify(z))));   
  fact=[]; until(1,
  if(flg!=0,
     if((t=gcd(real(z),imag(z)))!=1,
     fact=vecmat(factor(t));z/=t); 
  if(type(real(z))!=type(0)|type(imag(z))!=type(0),
     error("gfact: not gaussian integer")); 
  if(z==1,break);
  ); /* end if flg!=0 */
  nm=factor(norm(z));    
  for(n=1,matsize(nm)[1],
     p=nm[n,1];e2=nm[n,2];
     if(p%4==3,e2/=2;z/=p^e2;fact=concat(fact,[[p,e2]]) ,
       gp=g4n1(p); e1=valuation(z,p); t=z/gp^e1;
       if(1!=gcd(t,gp), e1=e2-e1); e2=e2-e1;
       z/=gp^e1*conj(gp)^e2;
       if(e1,fact=concat(fact,[[gp,e1]]));
       if(e2,fact=concat(fact,[[conj(gp),e2]]));
     ); /* end if p%4==3 */
  ); /* end for n */
  t=0;while(!normal(z),z/=I;t++);
  if(t,fact=concat(fact,[[I,t]]));
  if(z==1,break);
  print("residual factor="z);error("gfact: residual factor");
  ); /* end until(1,) */
  return(matvec(fact));
} /* end gfact() */