NB. Functions for factoring gaussian integers into primes.
NB. Described in "Prime Factorization in the Gaussian integers"
NB. Submitted to Journal of J, June 2013

NB. Norm
N=: *+

NB. -1 is a Quadratic residue witness

qrw=: 3 : 0
p=.x: y
e=.(p-1)%2 4
y=.0
while. y~:p-1 do.
  n=.x: ?&.<:p
  'y x'=.p&|@:(n&^)e
end.
x
)

NB. split p cong to 1 mod 4
splitm1=: (,0j1*+)@(+. qrw j. 1:)"0

NB. integer query
iq=: =<.

NB. Gaussian integer query
giq=:*./@:iq@:(**|)@:+."0

NB. gaussian integer fuzz removal
giclean=:(**|)&.+.

NB. main factorization function
gifactor=:3 : 0
pn=.q: N y
'p1 p2 p3'=.}.&.>(4&| </. ]) 1 2 3,pn
pf=.(#p2)#1j1
pf=.pf,;<@(-:@#{.])/.~p3
sp1=.,splitm1 ~.p1
y=.giclean y%*/pf
for_sp. sp1 do.
   while. giq q=.y%sp do.
   y=.giclean q
   pf=.pf,sp
   end.
end.
y,(/:|)pf
)

NB. gaussian prime query
giprimeq=:2 = #@gifactor

NB. gaussian integer prime query extended to 0
gipq=:0:`giprimeq@.(0&~:)"0