Gaussian integers

The problem with Gaussian integers in Forth is the problem of overflow. It's a minor problem that multiplication of two gintegers might need double word length for the components of the result, but a major problem that the norm (the sum of the square of the components) of gintegers might need double word length, since the components of the result of division of two gintegers are results of division with the norm. 

In Forth there are no standard words for division where the denominator is a double. It's possible to define of course, but this need to be done in assembler to obtain a fair speed. Perhaps it can be done in floating point with 20 decimals precision, but that seems tricky according to portability. Or there might be some other trick?

This first approach to Gaussian integers will be of poor range: just 3 decimals in 32 bit systems and 6 decimals in 64-bit systems. I might be able to double the range in the future but that's uncertain, but I will implement Gaussian integers in big integers, even if there might be problems with speed for big Gaussian integers.

The input of a ginteger is simply two integers on the stack, re im, where the imaginary component is at the top of the stack. But it's nice with a typographical output. The standard output of integers, like the dot, also include spaces, which don't allow desirable outputs like 

10-7i 

In 'Starting Forth' the output formatting procedures are described and the word for printing signed integers without spaces could be like this

: .sign-w/o-space \ n --
  dup 0< swap abs 0 <# #s rot sign #> type ;

which can be used for typographic output of Gaussian integers:

: g. \ a b --
  2dup 0. d= if d. exit then
  swap dup                           \ b a a
  if over 
     if dup .sign-w/o-space
     else dup .
     then
  then swap dup 0< 0=                \ a b f
  if dup \ a b b
     if swap if ." +" then dup 1 >   \ b f
        if .sign-w/o-space else drop then ." i" space
     else 2drop
     then
  else nip dup                       \ b b 
     if dup -1 <                     \ b f
        if .sign-w/o-space else ." -" drop then ." i" space
     else drop
     then
  then ;

Now 1 -1 g. results in 1-i, 0 1 g. in i and 10 7 g. in 10+7i.

And so a word printing the numbers on the stack interpreted as gintegers:

: .gs depth 2 < if exit then 2>r recurse 2r> 2dup g. ;

To drop, swap, rot gintegers one have to use 2drop, 2swap, 2rot and so on.

Arithmetics

: gnorm \ a b -- n 
  dup * swap dup * + ;

: g< \ a b c d -- flag
  gnorm -rot gnorm u> ;

: g+ \ a b c d -- a+c b+d
  under+ under+ ;

: g- \ a b c d -- a-c b-d
  negate under+ negate under+ ;

: g* { a b c d -- ac-bd ad+bc }  \ (a+ib)(c+id)
  a c * b d * - a d * b c * + ;

: g/ { a b c d -- (ac+bd)/(c^2+d^2) (bc-ad)/(c^2+d^2) }
  c d gnorm
  a c * b d * + over /
  b c * a d * - rot / ;

There is a need for a more general interpretation of modulo for gintegers that guaranties that the rest of a division has a smaller norm than the norm of the denominator. This is well explained in this paper.

The main idea is to use a generalized modulo for integers: the rest of n/m should be the smallest non negative integer r such that n=mq±r. This correspond to a generalized modulo for Gaussian integers, where all division of components use "rounded" results instead of the floor:

: /' ( a b -- c ) dup 2/ under+ / ;   \ c=floor((a+b/2)/b)

: g/' { a b c d -- (ac+bd)/'(c^2+d^2) (bc-ad)/'(c^2+d^2) }
  c d gnorm
  a c * b d * + over /'
  b c * a d * - rot /' ;

Thinking about it I found a way to use floats to increase the range:

\ --- alternative for about 50 % greater range -----

: g* { a b c d -- e f } 
  a c m* b d m* d- d>s a d m* b c m* d+ d>s ;

: g/ { a b c d -- e f }
  c dup m* d dup m* d+ d>f
  a c m* b d m* d+ d>f fover f/
  b c m* a d m* d- d>f frot f/
  fswap f>d d>s f>d d>s ;

: g/' { a b c d -- e f }
  c dup m* d dup m* d+ d>f
  a c m* b d m* d+ d>f fover f/ fround
  b c m* a d m* d- d>f frot f/ fround
  fswap f>d d>s f>d d>s ;

\ --------------------------------------------------

Using floats made a test running on my computer about 10 times slower, when I repeated testing gpollard2 (see below). But then I used numbers with random components in the range of 9 decimals, while in the first test a range of 6 decimals. The alternative g* g/ and g/' manage at least 4 decimals in 32-bit systems and 9 decimals in 64-bit systems.

: gmod ( a b c d -- e f )
  2over 2over g/' g* g- ;

: ggcd ( a b c d -- e f ) \ gcd for Gaussian integers
  0 0 loc{ a b c d s t }  \ Euclid's algorithm
  a b c d g<
  if a c to a to c
     b d to b to d
  then
  begin c d or \ while c+id<>0
  while c to s d to t
     a b c d gmod to d to c
     t to b s to a
  repeat a b ;

Of course the divisor is up to multiplication with unit: 1, i, -1, -i. To bring some order and simplicity to the matter I use a normal form: the normal form of a Gaussian integer is the representant where the real component is greater or equal then zero and greater or equal then the imaginary component:

: gnormal ( a b -- c d )  \ c>=d and c>=0
  2dup abs swap abs >
  if 0 1 g* then over 0< 
  if negate swap negate swap then ;

Due to Wolfram Mathworld a Gaussian integer a+ib is a prime if and only if it is satisfying one of the following properties:

1. If both  and  are nonzero then,  is a Gaussian prime iff  is an ordinary prime;
2. If , then  is a Gaussian prime iff  is an ordinary prime and ;
3. If , then  is a Gaussian prime iff  is an ordinary prime and .

: gprime1 \ n --flag
  abs dup prime swap 3 and 3 = and ;

: gprime \ a b -- flag 
  over 0= if nip gprime1 exit then
  dup 0= if drop gprime1 exit then
  gnorm prime ;

Since prime is only defined for single integers, gprime only works properly if a²+b² is a single cell integer.

To print the norm and the number a+ib in a table form:

: .normgp ( a b norm -- )
  cr 8 .r space g. ;

The next word print the norm and the normal form of all primes with a certain norm.

: .gprime loc{ norm -- }
  norm 2 = if 1 1 norm .normgp exit then
  norm sqrtf 1+ 0
  do i 0
     ?do i j gnorm norm =
        if i j gprime 
           if j i norm .normgp then
        then
     loop
  loop ;

: .gps ( n -- ) 1+ 2 do i .gprime loop ;

The word .gps lists all norms and primes on normal form up to numbers with the norm n. 

Note, that if a+ib is a prime then so is (a±ib)i^k, k=0,1,2,3.

Factorization

The Pollard rho factorization method seems to work for Gaussian integers with some minor changes:

: gfunc ( a b x y -- u v )
  2dup g* -1 0 g+ 2swap gmod ;

First I tested with the function z²+1 which don't seems to factorize all Gaussian integers, for example not 4+3i. Then I tried z²-1 and with that function I've tested millions of cases without problem.

: gpollard1 ( a b c d -- p q flag )  \ a+ib not having factor 2
  2dup loc{ a b alpha1 alpha2 beta1 beta2 } 0.
  begin 2drop a b alpha1 alpha2 gfunc to alpha2 to alpha1
     a b 2dup beta1 beta2 gfunc gfunc to beta2 to beta1
     alpha1 alpha2 beta1 beta2 g-
     a b ggcd 2dup a b d=
     if false exit
     then 2dup gnorm 1 = 0=
  until true ;

: gpollard2 ( a b -- p q )
  false loc{ flag }
  2dup gnorm 1 = if exit then
  2dup 2 0 gmod d0= if 2drop 1 1 exit then
  2dup gprime if exit then -1 1
  do i 0
     do 2dup j i gpollard1
        if true to flag leave then 2drop
     loop flag if leave then
  loop 2swap 2drop ;

: gfac ( a b -- p1 q1 ... pk qk )
  2dup gpollard2 2over 2over gnorm -rot gnorm =
  if 2drop exit
  then 2dup 2>r g/ recurse 2r> recurse ;

795 649 2dup g. 795+649i  ok
gfac .gs 8-5i 4+9i -6+5i -1-i  ok

g* g* g* g. 795+649i  ok

7891 9785 2dup g. 7891+9785i  ok
gfac .gs -47-10i 19+184i -1+i  ok
g* g* g. 7891+9785i  ok