Factorize a Gaussian integer

Ruby, 258 256 249 246+8 = 264 257 254 bytes

Uses the -rprime flag.

Geez, what a mess.

Uses this algorithm from stackoverflow.

->c{m=->x,y{x-y*eval("%d+%di"%(x/y).rect)};a=c.abs2.prime_division.flat_map{|b,e|b%4<2?(1..e).map{k=(2..d=b).find{|n|n**(~-b/2)%b==b-1}**(~-b/4)%b+1i;d,k=k,m[d,k]while k!=0;c/=d=m[c,d]==0?d:d.conj;d}:(c/=b<3?(b=1+1i)**e:b**e/=2;[b]*e)};a[0]*=c;a}

Try it online!

Python 2, 250 239 223 215 bytes

e,i,w=complex,int,abs
def f(*Z):
 if Z:
	z=Z[0];q=i(w(z));Q=4*q*q
	while Q>0:
 	 a=Q/q-q;b=Q%q-q;x=e(a,b)
 	 if w(x)>1:
		y=z/x
		if w(y)>1 and y==e(i(y.real),i(y.imag)):f(x,y);z=Q=0
 	 Q-=1
	if z:print z
	f(*Z[1:])

Try it online!

• -11 bytes when using Multiple Function Arguments
• -2²*² bytes when using one variable to parse couples (a,b)
• -2³ bytes when mixing tabs and spaces: thanks to ovs

Some explanation recursively decompose a complex to two complexes until no decomposition is possible...

Jelly, 61 55 50 bytes

ÆiḤp/-,1p`¤×€×1,ıFs2S€⁸÷ÆiḞƑ$ƇỊÐḟ1;Ṫð,÷@\ḟ1
Ç€F$ÐL

Try it online! (Header and Footer formats the output)

-6 bytes thanks to @EricTheOutgolfer

-5 bytes thanks to @caird coinheringaahing

How it Works

ÆiḤp/-,1p`¤×€×1,ıFs2S€⁸÷ÆiḞƑ$ƇỊÐḟ1;Ṫð,÷@\ḟ1 - helper: outputs a factor pair of the input
ÆiḤp/                    - creates a list of possible factors a+bi, a,b>=0
      -,1p`¤×€           - extend to the other three quadrants 
              ×1,ıFs2S€  - convert to  actual complex numbers 
⁸÷                       - get quotient with input complex number
  ÆiḞƑ$Ƈ                    - keep only Gaussian numbers (those unchanged when complex parts are floored)
     ỊÐḟ                 - remove units (i,-i,1,-1)
        1;               - append a 1 to deal with primes having no non-unit factors
          Ṫð,÷@\         - convert to a factor pair
                ḟ1       - remove 1s
Ç€F$ÐL
Ç€      - factor each number
   $    - and
  F     - flatten the list
    ÐL  - until factoring each number and flattening does not change the list