Source code in the jq programming language

### Complex numbers

# For clarity:
def real: first;
def imag: last;

# Complex or real
def norm:
  def sq: .*.;
  if type == "number" then sq
  else map(sq) | add
  end;

# Complex or real
def abs:
  if type == "array" then norm | sqrt
  elif . < 0 then - . else .
  end;

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

# Input should be an integer
def isPrime:
  . as $n
  | if ($n < 2)         then false
    elif ($n % 2 == 0)  then $n == 2
    elif ($n % 3 == 0)  then $n == 3
    else 5
    | until( . <= 0;
        if .*. > $n then -1
	elif ($n % . == 0) then 0
        else . + 2
        |  if ($n % . == 0) then 0
           else . + 4
           end
        end)
     | . == -1
     end;

# Given a stream of non-null values, 
# group by the values of `x|filter` that occur in a run.
def runs(stream; filter):
  foreach (stream, null) as $x ({emit: false, array: []};
    if $x == null 
    then .emit = .array
    elif .array == [] or ($x|filter) == (.array[0]|filter) 
    then .array += [$x] | .emit = false
    else {emit: .array, array: [$x]}
    end;
    select(.emit).emit)  ;

# emit a stream of Gaussian primes with real and imaginary parts within the given radius
def GPrimes($Radius):
  def check: abs | isPrime and ((. - 3) % 4 == 0);
  ($Radius | norm) as $R2
  | range(-$Radius; $Radius + 1) as $r
  | range(-$Radius; $Radius + 1) as $i
  | [$r, $i]
  | norm as $norm
  | select( $norm < $R2 )
  | if   $i == 0  then select($r|check)
    elif $r == 0  then select($i|check)
    else               select($norm|isPrime)
    end ;
    
def plot($Radius):
  "# X Y",
  (GPrimes($Radius) | "\(first) \(last)");

# For the first task:
"Gaussian primes with norm < 100, grouped by the real part:",
 (runs(GPrimes(10); real) | map(tostring) | join(" "))

# For the plot:
# plot(100)