/* 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() */