Gaussian primes: Difference between revisions
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
|||
Line 170: | Line 170: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="julia">using LinearAlgebra |
||
using Plots |
using Plots |
||
using Primes |
using Primes |
Revision as of 19:22, 30 August 2022
A Gaussian Integer is a complex number such that its real and imaginary parts are both integers.
a + bi where a and b are integers and i is √-1.
The norm of a Gaussian integer is its product with its conjugate.
N(a + bi) = (a + bi)(a − bi) = a² + b²
A Gaussian integer is a Gaussian prime if and only if either: both a and b are non-zero and its norm is a prime number, or, one of a or b is zero and it is the product of a unit (±1, ±i) and a prime integer of the form 4n + 3.
Prime integers that are not of the form 4n + 3 can be factored into a Gaussian integer and its complex conjugate so are not a Gaussian prime.
E.G. 5 = (2 + i)(2 − i) So 5 is not a Gaussian prime
Gaussian primes are octogonally symmetrical on a real / imaginary Cartesian field. If a particular complex norm a² + b² is prime, since addition is commutative, b² + a² is also prime, as are the complex conjugates and negated instances of both.
- Task
Find and show, here on this page, the Gaussian primes with a norm of less than 100, (within a radius of 10 from the origin 0 + 0i on a complex plane.)
Plot the points corresponding to the Gaussian primes on a Cartesian real / imaginary plane at least up to a radius of 50.
- See also
F#
This task uses Extensible Prime Generator (F#)
// Gaussian primes. Nigel Galloway: July 29th., 2022
let isGP=function (n,0)|(0,n)->let n=abs n in n%4=3 && isPrime n |(n,g)->isPrime(n*n+g*g)
Seq.allPairs [-9..9] [-9..9]|>Seq.filter isGP|>Seq.iter(fun(n,g)->printf $"""%d{n}%s{match g with 0->"" |g->sprintf $"%+d{g}i"} """); printfn ""
- Output:
-9-4i -9+4i -8-7i -8-5i -8-3i -8+3i -8+5i -8+7i -7-8i -7-2i -7 -7+2i -7+8i -6-5i -6-1i -6+1i -6+5i -5-8i -5-6i -5-4i -5-2i -5+2i -5+4i -5+6i -5+8i -4-9i -4-5i -4-1i -4+1i -4+5i -4+9i -3-8i -3-2i -3 -3+2i -3+8i -2-7i -2-5i -2-3i -2-1i -2+1i -2+3i -2+5i -2+7i -1-6i -1-4i -1-2i -1-1i -1+1i -1+2i -1+4i -1+6i 0-7i 0-3i 0+3i 0+7i 1-6i 1-4i 1-2i 1-1i 1+1i 1+2i 1+4i 1+6i 2-7i 2-5i 2-3i 2-1i 2+1i 2+3i 2+5i 2+7i 3-8i 3-2i 3 3+2i 3+8i 4-9i 4-5i 4-1i 4+1i 4+5i 4+9i 5-8i 5-6i 5-4i 5-2i 5+2i 5+4i 5+6i 5+8i 6-5i 6-1i 6+1i 6+5i 7-8i 7-2i 7 7+2i 7+8i 8-7i 8-5i 8-3i 8+3i 8+5i 8+7i 9-4i 9+4i
J
Implementation:
isgpri=: {{
if. 1 p: (*+) y do. 1 return. end.
int=. |(+.y)-.0
if. 1=#int do. {.(1 p: int) * 3=4|int else. 0 end.
}}"0
Online plot of gaussian primes up to radius 100. (Hit "Run" in the upper right-hand corner.)
Plot of gaussian primes up to radius 50:
1j1#"1'#' (<"1]50++.(#~ isgpri * 50>:|) ,j./~i:100)} '+' (<50 50)} '|' 50}"1 '-' 50} 100 100$' '
| # # | # # # | # # # # # # # # # # # # | # # # # # # # # | # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # # | # # # # # # # # # # | # # # # # # # # | # # # # # # # # # | # # # # # # # # # # # # # | # # # # # # # # # # | # # # # # # # # # # # # | # # # # # # # # # # # # # # # | # # # # # # # # # # # # | # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # - - - # - - - # - - - - - - - - - - - # - - - - - - - # - - - # - - - - - - - # - - - # - - - # - - + - - # - - - # - - - # - - - - - - - # - - - # - - - - - - - # - - - - - - - - - - - # - - - # - - # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # | # # # # # # # # # # # # # # # # # # # # # # # # # # # # # | # # # # # # # # # # # # # | # # # # # # # # # # # # | # # # # # # # # # # # # # # # | # # # # # # # # # # # # | # # # # # # # # # # | # # # # # # # # # # # # # | # # # # # # # # # | # # # # # # # # | # # # # # # # # # # | # # # # # # # # # | # # # # # # # # # # # # # # # # # # | # # # # # # # # | # # # # # # # # | # # # # # # # # # # # # | # # # | # #
Gaussian primes less than radius 10 (sorted by radius):
10 10$(/: |)(#~ isgpri * 10>|) ,j./~i:10
_1j_1 _1j1 1j_1 1j1 _2j_1 _2j1 _1j_2 _1j2 1j_2 1j2
2j_1 2j1 _3 0j_3 0j3 3 _3j_2 _3j2 _2j_3 _2j3
2j_3 2j3 3j_2 3j2 _4j_1 _4j1 _1j_4 _1j4 1j_4 1j4
4j_1 4j1 _5j_2 _5j2 _2j_5 _2j5 2j_5 2j5 5j_2 5j2
_6j_1 _6j1 _1j_6 _1j6 1j_6 1j6 6j_1 6j1 _5j_4 _5j4
_4j_5 _4j5 4j_5 4j5 5j_4 5j4 _7 0j_7 0j7 7
_7j_2 _7j2 _2j_7 _2j7 2j_7 2j7 7j_2 7j2 _6j_5 _6j5
_5j_6 _5j6 5j_6 5j6 6j_5 6j5 _8j_3 _8j3 _3j_8 _3j8
3j_8 3j8 8j_3 8j3 _8j_5 _8j5 _5j_8 _5j8 5j_8 5j8
8j_5 8j5 _9j_4 _9j4 _4j_9 _4j9 4j_9 4j9 9j_4 9j4
Julia
using LinearAlgebra
using Plots
using Primes
"""
function isGaussianprime(n::Complex{T}) where T <: Integer
A Gaussian prime is a non-unit Gaussian integer m + ni divisible only by its associates and by the units
1, i, -1, -i and by no other Gaussian integers.
The Gaussian primes fall into one of three categories:
Gaussian integers with imaginary part zero and a prime real part m with |m| a real prime satisfying |m| = 3 mod 4
Gaussian integers with real part zero and an imaginary part n with |n| real prime satisfying |n| = 3 mod 4
Gaussian integers having both real and imaginary parts, and its complex norm (square of algebraic norm) is a real prime number
"""
function isGaussianprime(n::Complex{T}) where T <: Integer
r, c = abs(real(n)), abs(imag(n))
return isprime(r * r + c * c) || c == 0 && isprime(r) && (r - 3) % 4 == 0 || r == 0 && isprime(c) && (c - 3) % 4 == 0
end
function testgaussprimes(lim = 10)
testvals = map(c -> c[1] + im * c[2], collect(Iterators.product(-lim:lim, -lim:lim)))
gprimes = sort!(filter(c -> isGaussianprime(c) && norm(c) < lim, testvals), by = norm)
println("Gaussian primes within $lim of the origin on the complex plane:")
foreach(p -> print(lpad(p[2], 10), p[1] % 10 == 0 ? "\n" : ""), enumerate(gprimes)) # print
scatter(gprimes) # plot
end
testgaussprimes()
- Output:
Gaussian primes within 10 of the origin on the complex plane: 1 + 1im 1 - 1im -1 - 1im -1 + 1im 1 + 2im -2 + 1im 2 + 1im 2 - 1im -2 - 1im 1 - 2im -1 - 2im -1 + 2im 3 + 0im -3 + 0im 0 - 3im 0 + 3im -3 - 2im -2 + 3im 3 + 2im 3 - 2im -2 - 3im 2 + 3im 2 - 3im -3 + 2im 4 + 1im 4 - 1im -1 + 4im -4 - 1im -4 + 1im -1 - 4im 1 - 4im 1 + 4im 5 - 2im 2 + 5im -5 + 2im -5 - 2im 5 + 2im -2 + 5im 2 - 5im -2 - 5im 1 - 6im -6 + 1im 6 + 1im -6 - 1im -1 - 6im -1 + 6im 1 + 6im 6 - 1im -4 + 5im 5 + 4im -5 + 4im 4 + 5im 5 - 4im -5 - 4im 4 - 5im -4 - 5im 0 + 7im -7 + 0im 0 - 7im 7 + 0im 7 + 2im -2 + 7im -2 - 7im 2 - 7im 2 + 7im 7 - 2im -7 - 2im -7 + 2im 6 - 5im -6 - 5im 5 + 6im -5 - 6im 5 - 6im -6 + 5im -5 + 6im 6 + 5im 3 + 8im -8 + 3im 8 + 3im -3 + 8im -8 - 3im 8 - 3im 3 - 8im -3 - 8im 8 + 5im -5 - 8im -5 + 8im 5 - 8im -8 + 5im -8 - 5im 8 - 5im 5 + 8im -4 + 9im -4 - 9im 9 + 4im -9 + 4im 9 - 4im -9 - 4im 4 - 9im 4 + 9im
Mathematica/Wolfram Language
n = 100;
digs = Reap@Do[If[Norm[i + I j]^2 < n, If[PrimeQ[i + I j, GaussianIntegers -> True], Sow[i + I j]]],
{i,Floor[-Sqrt[n]], Ceiling[Sqrt[n]]}, {j, Floor[-Sqrt[n]], Ceiling[Sqrt[n]]}
];
Multicolumn[digs[[2, 1]], Appearance -> "Horizontal"]
n = 50^2;
digs = Table[If[Norm[i + I j]^2 < n, If[PrimeQ[i + I j, GaussianIntegers -> True], "*", " "], " "],
{i,Floor[-Sqrt[n]], Ceiling[Sqrt[n]]}, {j, Floor[-Sqrt[n]], Ceiling[Sqrt[n]]}
];
digs //= Map[StringJoin];
digs //= StringRiffle[#, "\n"] &;
digs
- Output:
-9-4 I -9+4 I -8-5 I -8-3 I -8+3 I -8+5 I -7-2 I -7 -7+2 I -6-5 I -6-I -6+I -6+5 I -5-8 I -5-6 I -5-4 I -5-2 I -5+2 I -5+4 I -5+6 I -5+8 I -4-9 I -4-5 I -4-I -4+I -4+5 I -4+9 I -3-8 I -3-2 I -3 -3+2 I -3+8 I -2-7 I -2-5 I -2-3 I -2-I -2+I -2+3 I -2+5 I -2+7 I -1-6 I -1-4 I -1-2 I -1-I -1+I -1+2 I -1+4 I -1+6 I -7 I -3 I 3 I 7 I 1-6 I 1-4 I 1-2 I 1-I 1+I 1+2 I 1+4 I 1+6 I 2-7 I 2-5 I 2-3 I 2-I 2+I 2+3 I 2+5 I 2+7 I 3-8 I 3-2 I 3 3+2 I 3+8 I 4-9 I 4-5 I 4-I 4+I 4+5 I 4+9 I 5-8 I 5-6 I 5-4 I 5-2 I 5+2 I 5+4 I 5+6 I 5+8 I 6-5 I 6-I 6+I 6+5 I 7-2 I 7 7+2 I 8-5 I 8-3 I 8+3 I 8+5 I 9-4 I 9+4 I * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Perl
#!/usr/bin/perl
use strict; # https://rosettacode.org/wiki/Gaussian_primes
use warnings;
use ntheory qw( is_prime );
my ($plot, @primes) = gaussianprimes(10);
print "Primes within 10\n", join(', ', @primes) =~ s/.{94}\K /\n/gr;
($plot, @primes) = gaussianprimes(50);
print "\n\nPlot within 50\n$plot";
sub gaussianprimes
{
my $size = shift;
my $plot = ( ' ' x (2 * $size + 1) . "\n" ) x (2 * $size + 1);
my @primes;
for my $A ( -$size .. $size )
{
my $limit = int sqrt $size**2 - $A**2;
for my $B ( -$limit .. $limit )
{
my $norm = $A**2 + $B**2;
if ( is_prime( $norm )
or ( $A==0 && is_prime(abs $B) && (abs($B)-3)%4 == 0)
or ( $B==0 && is_prime(abs $A) && (abs($A)-3)%4 == 0) )
{
push @primes, sprintf("%2d%2di", $A, $B) =~ s/ (\di)/+$1/r;
substr $plot, ($B + $size + 1) * (2 * $size + 2) + $A + $size + 1, 1, 'X';
}
}
}
return $plot, @primes;
}
- Output:
Primes within 10 -9-4i, -9+4i, -8-5i, -8-3i, -8+3i, -8+5i, -7-2i, -7+0i, -7+2i, -6-5i, -6-1i, -6+1i, -6+5i, -5-8i, -5-6i, -5-4i, -5-2i, -5+2i, -5+4i, -5+6i, -5+8i, -4-9i, -4-5i, -4-1i, -4+1i, -4+5i, -4+9i, -3-8i, -3-2i, -3+0i, -3+2i, -3+8i, -2-7i, -2-5i, -2-3i, -2-1i, -2+1i, -2+3i, -2+5i, -2+7i, -1-6i, -1-4i, -1-2i, -1-1i, -1+1i, -1+2i, -1+4i, -1+6i, 0-7i, 0-3i, 0+3i, 0+7i, 1-6i, 1-4i, 1-2i, 1-1i, 1+1i, 1+2i, 1+4i, 1+6i, 2-7i, 2-5i, 2-3i, 2-1i, 2+1i, 2+3i, 2+5i, 2+7i, 3-8i, 3-2i, 3+0i, 3+2i, 3+8i, 4-9i, 4-5i, 4-1i, 4+1i, 4+5i, 4+9i, 5-8i, 5-6i, 5-4i, 5-2i, 5+2i, 5+4i, 5+6i, 5+8i, 6-5i, 6-1i, 6+1i, 6+5i, 7-2i, 7+0i, 7+2i, 8-5i, 8-3i, 8+3i, 8+5i, 9-4i, 9+4i
Plot within 50 X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X XX XX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X XX XX X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
Phix
You can run this online here.
-- -- demo/rosetta/Gaussian_primes.exw -- ================================ -- with javascript_semantics function gaussian_primes(integer radius) integer sq_radius = radius*radius sequence res = {} for i=1 to radius do if remainder(i,4)=3 then res = append(res,{i*i,i,0}) end if integer i2 = i*i for j=i to radius do integer r = i2+j*j if r>sq_radius then exit end if if is_prime(r) then res = append(res,{r,j,i}) end if end for end for res = sort(res) return res end function include builtins\complex.e function gpp(integer i, j) return pad_head(complex_sprint({i,j}),6) end function function g4(integer i,j) sequence res = {gpp(i,j)} if i!=0 then res = append(res,gpp(-i,j)) if j!=0 then res = append(res,gpp(-i,-j)) end if end if if j!=0 then res = append(res,gpp(i,-j)) end if return res end function function reflect(sequence g) sequence res = {} for p in g do integer {n,i,j} = p res &= g4(i,j) if i!=j then res &= g4(j,i) end if end for return res end function sequence g = gaussian_primes(10) printf(1,"Gaussian primes with a norm less than 100 sorted by norm:\n%s\n", {join_by(reflect(g),1,10," ")}) --g = gaussian_primes(50) -- (radius of 50) g = gaussian_primes(150) -- (radius of 150) constant title = "Gaussian primes" include pGUI.e Ihandle dlg, canvas cdCanvas cddbuffer, cdcanvas integer cx, cy procedure plot4(integer i,j) for im=+1 to -1 by -2 do for jm=+1 to -1 by -2 do cdCanvasPixel(cddbuffer, cx+im*i, cy+jm*j, CD_YELLOW) end for end for end procedure function redraw_cb(Ihandle /*ih*/) integer {width, height} = IupGetIntInt(canvas, "DRAWSIZE") cx = floor(width/2) cy = floor(height/2) cdCanvasActivate(cddbuffer) cdCanvasClear(cddbuffer) for p in g do integer {n,i,j} = p plot4(i,j) plot4(j,i) end for cdCanvasFlush(cddbuffer) return IUP_DEFAULT end function function map_cb(Ihandle ih) cdcanvas = cdCreateCanvas(CD_IUP, ih) cddbuffer = cdCreateCanvas(CD_DBUFFER, cdcanvas) cdCanvasSetBackground(cddbuffer, CD_BLACK) return IUP_DEFAULT end function IupOpen() --canvas = IupCanvas("RASTERSIZE=320x320") canvas = IupCanvas("RASTERSIZE=340x340") IupSetCallbacks(canvas, {"MAP_CB", Icallback("map_cb"), "ACTION", Icallback("redraw_cb")}) dlg = IupDialog(canvas, `TITLE="%s",RESIZE=NO`,{title}) IupShow(dlg) if platform()!=JS then IupMainLoop() IupClose() end if
Output same as Raku
Python
''' python example for task rosettacode.org/wiki/Gaussian_primes '''
from matplotlib.pyplot import scatter
from sympy import isprime
from math import isqrt
def norm(c):
''' Task complex norm function '''
return c.real * c.real + c.imag * c.imag
def is_gaussian_prime(n):
'''
is_gaussian_prime(n)
A Gaussian prime is a non-unit Gaussian integer m + ni divisible only by its associates and by the units
1, i, -1, -i and by no other Gaussian integers.
The Gaussian primes fall into one of three categories:
Gaussian integers with imaginary part zero and a prime real part m with |m| a real prime satisfying |m| = 3 mod 4
Gaussian integers with real part zero and an imaginary part n with |n| real prime satisfying |n| = 3 mod 4
Gaussian integers having both real and imaginary parts, and its complex norm (square of algebraic norm) is a real prime number
'''
r, c = int(abs(n.real)), int(abs(n.imag))
return isprime(r * r + c * c) or c == 0 and isprime(r) and (r - 3) % 4 == 0 or r == 0 and isprime(c) and (c - 3) % 4 == 0
if __name__ == '__main__':
limitsquared = 100
lim = isqrt(limitsquared)
testvals = [complex(r, c) for r in range(-lim, lim) for c in range(-lim, lim)]
gprimes = sorted(filter(lambda c : is_gaussian_prime(c) and norm(c) < limitsquared, testvals), key=norm)
print(f'Gaussian primes within {isqrt(limitsquared)} of the origin on the complex plane:')
for i, c in enumerate(gprimes):
print(str(c).ljust(9), end='\n' if (i +1) % 10 == 0 else '')
scatter([c.real for c in gprimes], [c.imag for c in gprimes])
- Output:
Gaussian primes within 10 of the origin on the complex plane: (-1-1j) (-1+1j) (1-1j) (1+1j) (-2-1j) (-2+1j) (-1-2j) (-1+2j) (1-2j) (1+2j) (2-1j) (2+1j) (-3+0j) -3j 3j (3+0j) (-3-2j) (-3+2j) (-2-3j) (-2+3j) (2-3j) (2+3j) (3-2j) (3+2j) (-4-1j) (-4+1j) (-1-4j) (-1+4j) (1-4j) (1+4j) (4-1j) (4+1j) (-5-2j) (-5+2j) (-2-5j) (-2+5j) (2-5j) (2+5j) (5-2j) (5+2j) (-6-1j) (-6+1j) (-1-6j) (-1+6j) (1-6j) (1+6j) (6-1j) (6+1j) (-5-4j) (-5+4j) (-4-5j) (-4+5j) (4-5j) (4+5j) (5-4j) (5+4j) (-7+0j) -7j 7j (7+0j) (-7-2j) (-7+2j) (-2-7j) (-2+7j) (2-7j) (2+7j) (7-2j) (7+2j) (-6-5j) (-6+5j) (-5-6j) (-5+6j) (5-6j) (5+6j) (6-5j) (6+5j) (-8-3j) (-8+3j) (-3-8j) (-3+8j) (3-8j) (3+8j) (8-3j) (8+3j) (-8-5j) (-8+5j) (-5-8j) (-5+8j) (5-8j) (5+8j) (8-5j) (8+5j) (-9-4j) (-9+4j) (-4-9j) (-4+9j) (4-9j) (4+9j) (9-4j) (9+4j)
Raku
Plotting the points up to a radius of 150.
use List::Divvy;
my @next = { :1x, :1y, :2n },;
sub next-interval (Int $int) {
@next.append: (^$int).map: { %( :x($int), :y($_), :n($int² + .²) ) };
@next = |@next.sort: *.<n>;
}
my @gaussian = lazy gather {
my $interval = 1;
loop {
my @this = @next.shift;
@this.push: @next.shift while @next and @next[0]<n> == @this[0]<n>;
for @this {
.take if .<n>.is-prime || (!.<y> && .<x>.is-prime && (.<x> - 3) %% 4);
next-interval(++$interval) if $interval == .<x>
}
}
}
# Primes within a radius of 10 from origin
say "Gaussian primes with a norm less than 100 sorted by norm:";
say @gaussian.&before(*.<n> > 10²).map( {
my (\i, \j) = .<x y>;
flat ((i,j),(-i,j),(-i,-j),(i,-j),(j,i),(-j,i),(-j,-i),(j,-i)).map: {
.[0] ?? .[1] ?? (sprintf "%d%s%di", .[0], (.[1] ≥ 0 ?? '+' !! ''), .[1]) !! .[0] !! "{.[1]}i"
}} )».subst('1i', 'i', :g)».fmt("%6s")».unique.flat.batch(10).join: "\n" ;
# Plot points within a 150 radius
use SVG;
my @points = unique flat @gaussian.&before(*.<n> > 150²).map: {
my (\i, \j) = .<x y>;
do for (i,j),(-i,j),(-i,-j),(i,-j),(j,i),(-j,i),(-j,-i),(j,-i) {
:use['xlink:href'=>'#point', 'transform'=>"translate({500 + 3 × .[0]},{500 + 3 × .[1]})"]
}
}
'gaussian-primes-raku.svg'.IO.spurt: SVG.serialize(
svg => [
:width<1000>, :height<1000>,
:rect[:width<100%>, :height<100%>, :style<fill:black;>],
:defs[:g[:id<point>, :circle[:0cx, :0cy, :2r, :fill('gold')]]],
|@points
],
);
- Output:
Gaussian primes with a norm less than 100 sorted by norm: 1+i -1+i -1-i 1-i 2+i -2+i -2-i 2-i 1+2i -1+2i -1-2i 1-2i 3 -3 3i -3i 3+2i -3+2i -3-2i 3-2i 2+3i -2+3i -2-3i 2-3i 4+i -4+i -4-i 4-i 1+4i -1+4i -1-4i 1-4i 5+2i -5+2i -5-2i 5-2i 2+5i -2+5i -2-5i 2-5i 6+i -6+i -6-i 6-i 1+6i -1+6i -1-6i 1-6i 5+4i -5+4i -5-4i 5-4i 4+5i -4+5i -4-5i 4-5i 7 -7 7i -7i 7+2i -7+2i -7-2i 7-2i 2+7i -2+7i -2-7i 2-7i 6+5i -6+5i -6-5i 6-5i 5+6i -5+6i -5-6i 5-6i 8+3i -8+3i -8-3i 8-3i 3+8i -3+8i -3-8i 3-8i 8+5i -8+5i -8-5i 8-5i 5+8i -5+8i -5-8i 5-8i 9+4i -9+4i -9-4i 9-4i 4+9i -4+9i -4-9i 4-9i
Off-site SVG image: gaussian-primes-raku.svg
Wren
Plots the points up to a radius of 150 to produce a similar image to the Raku example.
import "dome" for Window
import "graphics" for Canvas, Color
import "./plot" for Axes
import "./complex" for Complex
import "./math2" for Int
import "./fmt" for Fmt
var norm = Fn.new { |c| c.real * c.real + c.imag * c.imag }
var GPrimes = []
var Radius = 150
for (r in -Radius+1...Radius) {
for (i in -Radius+1...Radius) {
if (i == 0) {
var m = r.abs
if (Int.isPrime(m) && (m - 3) % 4 == 0) GPrimes.add(Complex.new(r))
} else if (r == 0) {
var m = i.abs
if (Int.isPrime(m) && (m - 3) % 4 == 0) GPrimes.add(Complex.new(0, i))
} else {
var n = r * r + i * i
if (n < Radius * Radius && Int.isPrime(n)) GPrimes.add(Complex.new(r, i))
}
}
}
var gp10 = GPrimes.where { |p| norm.call(p) < 100 }.toList
gp10.sort { |i, j|
var ni = norm.call(i)
var nj = norm.call(j)
if (ni != nj) return ni < nj
if (i.real != j.real) return i.real > j.real
return i.imag > j.imag
}
System.print("Gaussian primes with a norm less than 100 sorted by norm:")
Fmt.tprint("($2.0z) ", gp10, 5)
GPrimes = GPrimes.map { |c| c.toPair }.toList
class Main {
construct new() {
Window.title = "Gaussian primes"
Canvas.resize(1000, 1000)
Window.resize(1000, 1000)
Canvas.cls(Color.black)
var axes = Axes.new(100, 900, 800, 800, -Radius..Radius, -Radius..Radius)
axes.plot(GPrimes, Color.yellow, "·")
}
init() {}
update() {}
draw(alpha) {}
}
var Game = Main.new()
- Output:
Terminal output:
Gaussian primes with a norm less than 100 sorted by norm: ( 1 + 1i) ( 1 - 1i) (-1 + 1i) (-1 - 1i) ( 2 + 1i) ( 2 - 1i) ( 1 + 2i) ( 1 - 2i) (-1 + 2i) (-1 - 2i) (-2 + 1i) (-2 - 1i) ( 3 + 0i) ( 0 + 3i) ( 0 - 3i) (-3 + 0i) ( 3 + 2i) ( 3 - 2i) ( 2 + 3i) ( 2 - 3i) (-2 + 3i) (-2 - 3i) (-3 + 2i) (-3 - 2i) ( 4 + 1i) ( 4 - 1i) ( 1 + 4i) ( 1 - 4i) (-1 + 4i) (-1 - 4i) (-4 + 1i) (-4 - 1i) ( 5 + 2i) ( 5 - 2i) ( 2 + 5i) ( 2 - 5i) (-2 + 5i) (-2 - 5i) (-5 + 2i) (-5 - 2i) ( 6 + 1i) ( 6 - 1i) ( 1 + 6i) ( 1 - 6i) (-1 + 6i) (-1 - 6i) (-6 + 1i) (-6 - 1i) ( 5 + 4i) ( 5 - 4i) ( 4 + 5i) ( 4 - 5i) (-4 + 5i) (-4 - 5i) (-5 + 4i) (-5 - 4i) ( 7 + 0i) ( 0 + 7i) ( 0 - 7i) (-7 + 0i) ( 7 + 2i) ( 7 - 2i) ( 2 + 7i) ( 2 - 7i) (-2 + 7i) (-2 - 7i) (-7 + 2i) (-7 - 2i) ( 6 + 5i) ( 6 - 5i) ( 5 + 6i) ( 5 - 6i) (-5 + 6i) (-5 - 6i) (-6 + 5i) (-6 - 5i) ( 8 + 3i) ( 8 - 3i) ( 3 + 8i) ( 3 - 8i) (-3 + 8i) (-3 - 8i) (-8 + 3i) (-8 - 3i) ( 8 + 5i) ( 8 - 5i) ( 5 + 8i) ( 5 - 8i) (-5 + 8i) (-5 - 8i) (-8 + 5i) (-8 - 5i) ( 9 + 4i) ( 9 - 4i) ( 4 + 9i) ( 4 - 9i) (-4 + 9i) (-4 - 9i) (-9 + 4i) (-9 - 4i)