From 0895026182f24d15a745d33c1a3c3745ea45bc17 Mon Sep 17 00:00:00 2001 From: trizen Date: Fri, 12 Jun 2020 16:00:10 +0300 Subject: [PATCH] - 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 | 1 + lib/Sidef/Types/Number/Number.pm | 125 +++++++++++++++++++++++++++-- lib/Sidef/Types/Number/Number.pod | 60 ++++++++++++++ scripts/Tests/gaussian_integers.sf | 51 ++++++++++++ 4 files changed, 230 insertions(+), 7 deletions(-) create mode 100644 scripts/Tests/gaussian_integers.sf diff --git a/MANIFEST b/MANIFEST index 1c1d5fb6d..aca314707 100644 --- a/MANIFEST +++ b/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 diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm index 5a14c49cb..5434eba5f 100644 --- a/lib/Sidef/Types/Number/Number.pm +++ b/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); diff --git a/lib/Sidef/Types/Number/Number.pod b/lib/Sidef/Types/Number/Number.pod index d7952fcb0..edbb8b333 100644 --- a/lib/Sidef/Types/Number/Number.pod +++ b/lib/Sidef/Types/Number/Number.pod @@ -1021,6 +1021,16 @@ Aliases: I =cut +=head2 cadd + +Number.cadd() -> I + +Return the + +Aliases: I + +=cut + =head2 catalan n.catalan @@ -1040,6 +1050,16 @@ Cube root of C. =cut +=head2 cdiv + +Number.cdiv() -> I + +Return the + +Aliases: I + +=cut + =head2 ceil n.ceil @@ -1116,6 +1136,16 @@ Example: =cut +=head2 cinvmod + +Number.cinvmod() -> I + +Return the + +Aliases: I + +=cut + =head2 circular_permutations Number.circular_permutations() -> I @@ -1140,6 +1170,26 @@ Return the =cut +=head2 cmod + +Number.cmod() -> I + +Return the + +Aliases: I + +=cut + +=head2 cmul + +Number.cmul() -> I + +Return the + +Aliases: I + +=cut + =head2 combinations Number.combinations() -> I @@ -1286,6 +1336,16 @@ Return the =cut +=head2 csub + +Number.csub() -> I + +Return the + +Aliases: I + +=cut + =head2 cyclotomic Number.cyclotomic() -> I diff --git a/scripts/Tests/gaussian_integers.sf b/scripts/Tests/gaussian_integers.sf new file mode 100644 index 000000000..a9f7beba7 --- /dev/null +++ b/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!"