- Added support for rational complex number operations.

Provided by the following new Number methods:

	cmod(a,b,m)	#=> (a%m, b%m)
	cadd(a,b,x,y) 	#=> (a+x, b+y)
	csub(a,b,x,y)	#=> (a-x, b-y)
	cmul(a,b,x,y)	#=> (a*x - b*y, a*y + b*x)
	cdiv(a,b,x,y)	#=> ((a*x + b*y)/(x*x + y*y), (b*x - a*y)/(x*x + y*y))

Also added the `complex_invmod(a, b, m)` method, which returns a pair of integers (x,y) such that cmod(cmul(a,b,x,y), m) == (1, 0).

MANIFEST

@@ -752,6 +752,7 @@ scripts/Tests/function_composition.sf
 scripts/Tests/functional_modules.sf
 scripts/Tests/functional_style.sf
 scripts/Tests/gather_take.sf
+scripts/Tests/gaussian_integers.sf
 scripts/Tests/gcd.sf
 scripts/Tests/get_value_cyclic_references.sf
 scripts/Tests/getopt.sf

lib/Sidef/Types/Number/Number.pm

@@ -6973,14 +6973,125 @@ package Sidef::Types::Number::Number {
     *expmod = \&modpow;
     *powmod = \&modpow;
 
+    sub complex_mod {
+        my ($x, $y, $m) = @_;
+        _valid(\$y, \$m);
+        ((bless \__mod__($$x, $$m)), (bless \__mod__($$y, $$m)));
+    }
+
+    *cmod = \&complex_mod;
+
+    sub complex_add {
+        my ($x_re, $x_im, $y_re, $y_im) = @_;
+
+        $x_im //= ZERO;
+        $y_re //= ZERO;
+        $y_im //= ZERO;
+
+        _valid(\$x_im, \$y_re, \$y_im);
+        ((bless \__add__($$x_re, $$y_re)), (bless \__add__($$x_im, $$y_im)));
+    }
+
+    *cadd = \&complex_add;
+
+    sub complex_sub {
+        my ($x_re, $x_im, $y_re, $y_im) = @_;
+
+        $x_im //= ZERO;
+        $y_re //= ZERO;
+        $y_im //= ZERO;
+
+        _valid(\$x_im, \$y_re, \$y_im);
+        ((bless \__sub__($$x_re, $$y_re)), (bless \__sub__($$x_im, $$y_im)));
+    }
+
+    *csub = \&complex_sub;
+
+    sub complex_mul {
+        my ($x_re, $x_im, $y_re, $y_im) = @_;
+
+        # (a + b*i) * (x + y*i) = (a*x - b*y) + (a*y + b*x)*i
+
+        $x_im //= ZERO;
+        $y_re //= ZERO;
+        $y_im //= ZERO;
+
+        _valid(\$x_im, \$y_re, \$y_im);
+
+#<<<
+        (
+            (bless \__sub__(__mul__($$x_re, $$y_re), __mul__($$x_im, $$y_im))),
+            (bless \__add__(__mul__($$x_re, $$y_im), __mul__($$x_im, $$y_re))),
+        );
+#>>>
+    }
+
+    *cmul = \&complex_mul;
+
+    sub complex_div {
+        my ($x_re, $x_im, $y_re, $y_im) = @_;
+
+        # (a + b*i) / (x + y*i) = (a*x + b*y)/(x^2 + y^2) + (b*x - a*y)/(x^2 + y^2)*i
+
+        $x_im //= ZERO;
+        $y_re //= ZERO;
+        $y_im //= ZERO;
+
+        _valid(\$x_im, \$y_re, \$y_im);
+
+        my $den = __add__(__mul__($$y_re, $$y_re), __mul__($$y_im, $$y_im));
+
+#<<<
+        (
+            (bless \__div__(__add__(__mul__($$x_re, $$y_re), __mul__($$x_im, $$y_im)), $den)),
+            (bless \__div__(__sub__(__mul__($$x_im, $$y_re), __mul__($$x_re, $$y_im)), $den)),
+        );
+#>>>
+    }
+
+    *cdiv = \&complex_div;
+
+    sub complex_invmod {
+        my ($x, $y, $m) = @_;
+
+        _valid(\$y, \$m);
+
+        $x = _any2mpz($$x) // return (nan(), nan());
+        $y = _any2mpz($$y) // return (nan(), nan());
+        $m = _any2mpz($$m) // return (nan(), nan());
+
+        my $t = Math::GMPz::Rmpz_init();
+
+        Math::GMPz::Rmpz_mul($t, $x, $x);
+        Math::GMPz::Rmpz_addmul($t, $y, $y);
+
+        if (Math::GMPz::Rmpz_invert($t, $t, $m)) {
+
+            my $c0 = Math::GMPz::Rmpz_init();
+            my $c1 = Math::GMPz::Rmpz_init();
+
+            Math::GMPz::Rmpz_mul($c0, $x, $t);
+            Math::GMPz::Rmpz_mul($c1, $y, $t);
+            Math::GMPz::Rmpz_neg($c1, $c1);
+            Math::GMPz::Rmpz_mod($c0, $c0, $m);
+            Math::GMPz::Rmpz_mod($c1, $c1, $m);
+
+            return ((bless \$c0), (bless \$c1));
+        }
+
+        return (nan(), nan());    # no inverse
+    }
+
+    *cinvmod = \&complex_invmod;
+
     sub complex_pow {
         my ($x, $y, $n) = @_;
 
         _valid(\$y, \$n);
 
-        $x = _any2mpz($$x) // goto &nan;
-        $y = _any2mpz($$y) // goto &nan;
-        $n = _any2mpz($$n) // goto &nan;
+        $x = _any2mpz($$x) // return (nan(), nan());
+        $y = _any2mpz($$y) // return (nan(), nan());
+        $n = _any2mpz($$n) // return (nan(), nan());
 
         my $c0 = Math::GMPz::Rmpz_init_set_ui(1);
         my $c1 = Math::GMPz::Rmpz_init_set_ui(0);
@@ -7020,10 +7131,10 @@ package Sidef::Types::Number::Number {
 
         _valid(\$y, \$n, \$m);
 
-        $x = _any2mpz($$x) // goto &nan;
-        $y = _any2mpz($$y) // goto &nan;
-        $n = _any2mpz($$n) // goto &nan;
-        $m = _any2mpz($$m) // goto &nan;
+        $x = _any2mpz($$x) // return (nan(), nan());
+        $y = _any2mpz($$y) // return (nan(), nan());
+        $n = _any2mpz($$n) // return (nan(), nan());
+        $m = _any2mpz($$m) // return (nan(), nan());
 
         my $c0 = Math::GMPz::Rmpz_init_set_ui(1);
         my $c1 = Math::GMPz::Rmpz_init_set_ui(0);

lib/Sidef/Types/Number/Number.pod

@@ -1021,6 +1021,16 @@ Aliases: I<first>
 
 =cut
 
+=head2 cadd
+
+Number.cadd() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_add>
+
+=cut
+
 =head2 catalan
 
     n.catalan
@@ -1040,6 +1050,16 @@ Cube root of C<n>.
 
 =cut
 
+=head2 cdiv
+
+Number.cdiv() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_div>
+
+=cut
+
 =head2 ceil
 
     n.ceil
@@ -1116,6 +1136,16 @@ Example:
 
 =cut
 
+=head2 cinvmod
+
+Number.cinvmod() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_invmod>
+
+=cut
+
 =head2 circular_permutations
 
 Number.circular_permutations() -> I<Obj>
@@ -1140,6 +1170,26 @@ Return the
 
 =cut
 
+=head2 cmod
+
+Number.cmod() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_mod>
+
+=cut
+
+=head2 cmul
+
+Number.cmul() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_mul>
+
+=cut
+
 =head2 combinations
 
 Number.combinations() -> I<Obj>
@@ -1286,6 +1336,16 @@ Return the
 
 =cut
 
+=head2 csub
+
+Number.csub() -> I<Obj>
+
+Return the
+
+Aliases: I<complex_sub>
+
+=cut
+
 =head2 cyclotomic
 
 Number.cyclotomic() -> I<Obj>

scripts/Tests/gaussian_integers.sf

@@ -0,0 +1,51 @@
+#!/usr/bin/ruby
+
+var r = (-2 .. 2)
+
+for a in (r), b in (r), c in (r), d in (r) {
+
+    assert_eq([Complex(a,b) + Complex(c,d) -> reals], [cadd(a,b,c,d)])
+    assert_eq([Complex(a,b) - Complex(c,d) -> reals], [csub(a,b,c,d)])
+    assert_eq([Complex(a,b) * Complex(c,d) -> reals], [cmul(a,b,c,d)])
+
+    if (c*c + d*d != 0) {
+        assert_eq([Complex(a,b) / Complex(c,d) -> reals], [cdiv(a,b,c,d)].map{.float})
+    }
+}
+
+for a in (r), b in (r) {
+
+    var n = irand(0, 100)
+    var m = irand(100, 1000)
+
+    assert_eq([Complex(a,b)**n -> reals], [cpow(a,b,n)])
+    assert_eq([Complex(a,b)**n -> reals].map { .mod(m) } , [cpowmod(a,b,n,m)])
+}
+
+func gaussian_sum(n) {
+
+    var total = [0, 0]
+
+    for k in (1..n) {
+        total = [cadd(total..., cdiv(cpow(0, 1, k-1), k))]
+    }
+
+    [cmul(total..., n!)]
+}
+
+var arr = 10.of(gaussian_sum)
+
+assert_eq(arr.map{.head}, [0, 1, 2, 4, 16, 104, 624, 3648, 29184, 302976])
+assert_eq(arr.map{.tail}, [0, 0, 1, 3, 6, 30, 300, 2100, 11760, 105840])
+
+do {
+    var m = 10001
+    var a = 43
+    var b = 97
+
+    assert_eq([cdiv(1, 0, a, b)], [a / (a*a + b*b), -b / (a*a + b*b)])
+    assert_eq([a * invmod(a*a + b*b, m), -b * invmod(a*a + b*b, m)].map{.mod(m)}, [complex_invmod(a, b, m)])
+    assert_eq([cmod(cmul(a, b, complex_invmod(a, b, m)), m)], [1, 0])
+}
+
+say "** Test passed!"