- Optimized Number ipow2(n) and ipow10(n) when the result is a na…

…tive integer.
trizen · Dec 7, 2023 · b787581 · b787581
1 parent d8d9ca3
commit b787581
Show file tree

Hide file tree

Showing 2 changed files with 94 additions and 5 deletions.
diff --git a/lib/Sidef/Types/Number/Number.pm b/lib/Sidef/Types/Number/Number.pm
@@ -37,6 +37,7 @@ package Sidef::Types::Number::Number {
         ONE   => bless(\(my $one   = 1)),
         TWO   => bless(\(my $two   = 2)),
         THREE => bless(\(my $three = 3)),
+        TEN   => bless(\(my $three = 10)),
         ZERO  => bless(\(my $zero  = 0)),
         MONE  => bless(\(my $mone  = -1)),
 
@@ -4327,9 +4328,20 @@ package Sidef::Types::Number::Number {
     sub ipow2 {
         my ($n) = @_;
 
-        $n = _any2si($$n) // goto &nan;
+        $n = $$n;
+
+        if (ref($n)) {
+            $n = _any2si($n) // goto &nan;
+        }
 
-        return ZERO if $n < 0;
+        ($n == 0) and return ONE;
+        ($n == 1) and return TWO;
+        ($n < 0)  and return ZERO;
+
+        if (CORE::log(2) * $n < CORE::log(ULONG_MAX)) {
+            my $r = (HAS_NEW_PRIME_UTIL ? Math::Prime::Util::powint(2, $n) : Math::Prime::Util::GMP::powint(2, $n));
+            return bless \$r;
+        }
 
         my $r = Math::GMPz::Rmpz_init();
         Math::GMPz::Rmpz_setbit($r, $n);
@@ -4339,9 +4351,20 @@ package Sidef::Types::Number::Number {
     sub ipow10 {
         my ($n) = @_;
 
-        $n = _any2si($$n) // goto &nan;
+        $n = $$n;
+
+        if (ref($n)) {
+            $n = _any2si($n) // goto &nan;
+        }
 
-        return ZERO if $n < 0;
+        ($n == 0) and return ONE;
+        ($n == 1) and return TEN;
+        ($n < 0)  and return ZERO;
+
+        if (CORE::log(10) * $n < CORE::log(ULONG_MAX)) {
+            my $r = (HAS_NEW_PRIME_UTIL ? Math::Prime::Util::powint(10, $n) : Math::Prime::Util::GMP::powint(10, $n));
+            return bless \$r;
+        }
 
         my $r = Math::GMPz::Rmpz_init();
         Math::GMPz::Rmpz_ui_pow_ui($r, 10, $n);
@@ -13178,7 +13201,7 @@ package Sidef::Types::Number::Number {
 
             # say "Small k: ($n, $k, $m)";
 
-            if ($k <= 1e6) {
+            if ($k <= 1e7) {
                 my $bin = Math::GMPz::Rmpz_init();
                 if (Math::GMPz::Rmpz_fits_ulong_p($n) and Math::GMPz::Rmpz_cmp_ui($n, 1e5) <= 0) {
                     Math::GMPz::Rmpz_bin_uiui($bin, Math::GMPz::Rmpz_get_ui($n), $k);
@@ -13191,6 +13214,8 @@ package Sidef::Types::Number::Number {
                 return $bin;
             }
 
+            # TODO: find a faster algorithm
+
             my $t   = Math::GMPz::Rmpz_init();
             my $u   = Math::GMPz::Rmpz_init();
             my $bin = Math::GMPz::Rmpz_init_set_ui(1);

diff --git a/scripts/Extended tests/binomialmod.sf b/scripts/Extended tests/binomialmod.sf
@@ -50,4 +50,68 @@ for n in (-upto .. upto), k in (-upto .. upto), m in (-upto .. upto) {
     assert_eq(binomialmod(2*n, n, n), powmod(2, n, n))
 }
 
+# OEIS: A080469
+%n[36, 57, 121, 132, 552, 8397, 7000713, 9692541, 36294723, 564033861].each {|n|
+    assert_eq(binomialmod(3*n, n, n), powmod(3, n, n))
+}
+
+# OEIS: A260640
+%n[1, 3, 6, 12, 21, 35, 44, 55, 60, 70, 78, 88, 90, 99, 102, 110, 117, 119, 120, 133, 156, 171, 176, 180, 184, 204, 207, 220, 225, 230, 231, 234, 238, 240, 247, 252, 255, 285, 286, 300, 312, 341, 342, 348, 360, 368, 372, 391, 403, 408, 414, 425, 434, 460, 462, 465, 468, 481, 483, 494, 495, 504, 506, 510, 550, 555, 561, 572, 574, 585, 600].all {|n|
+    assert_eq(binomialmod(3*n, n, n), 0)
+}
+
+# OEIS: A109642
+%n[4, 15, 57, 765, 1025, 2097, 4947, 9189, 103599, 216927, 4346128, 1558269, 1977777].each_kv {|n,m|
+    assert_eq(binomialmod(3*m, m, m), powmod(3, n+1, m))
+}
+
+# OEIS: A109760
+%n[4, 365, 400, 685, 3200, 6400, 12550, 12800, 16525, 25600, 51200, 225125, 70463125, 271094125, 431434441].each {|n|
+    assert_eq(binomialmod(5*n, n, n), powmod(5, n, n))
+}
+
+# OEIS: A290040
+%n[260, 1056, 1060, 3460, 3905, 4428, 5000, 5060, 5512, 5860, 6372, 6596, 7460, 8200, 8908, 9612, 9860, 10660, 11556, 12260, 12625, 13060, 14600, 14660, 14744, 14796, 15460, 16260, 17060, 17800, 17860, 18425, 18496, 18660, 19396, 20260, 21717, 21860, 22168, 22248, 22660, 24260, 24616, 25164, 26660, 27108, 27400, 27460, 28872, 29060, 29128, 29860].each {|m|
+    assert(m.divisors.slice(1).any {|d| binomialmod(m+d, d, m) == 1 })
+}
+
+# OEIS: A260209
+%n[1, 3, 25, 245, 121, 169, 867, 3249, 6877, 9251, 961, 15059, 57154, 61017, 68479, 106742, 201898, 208376, 107736, 176435, 330398, 237158, 158447, 213867, 903264, 856884, 21218, 755634, 1259386, 944906, 161290, 531991, 150152, 656914, 1287658, 592826, 640874].each_kv {|n,A|
+    var p = prime(n+1)
+    assert_eq(binomialmod(2*p - 1, p-1, p**4), (A*p + 1) % (p**4))
+}
+
+# OEIS: A263429
+%n[2, 3, 5, 16843].each_kv{|n,p|
+    assert_eq(binomialmod(2*p - 1, p-1, p**(n+1)), 1)
+}
+
+assert_eq(binomialmod(2*4514260853041 - 1, 4514260853041 - 1, 4514260853041.isqrt), 1)
+assert_eq(binomialmod(2*283686649 - 1, 283686649, 283686649.isqrt), 1)
+
+#assert_eq(binomialmod(2*283686649 - 1, 283686649, 283686649), 1)                       # this takes 3.5s
+#assert_eq(binomialmod(2*283686649 - 1, 283686649, 283686649**2), 1)                    # this takes 7.5s
+#assert_eq(binomialmod(2*4514260853041 - 1, 4514260853041 - 1, 4514260853041), 1)       # FIXME: this takes a lot of time
+
+assert_eq(binomialmod(2*16843 - 1, 16843-1, 16843**4), 1)
+assert_eq(binomialmod(2*2124679 - 1, 2124679-1, 2124679**2), 1)
+
+assert_eq(binomialmod(2*2001341 - 1, 2001341-1, 2001341), 1)
+assert_eq(binomialmod(2*16024189487 - 1, 16024189487-1, 16024189487), 1)
+
+assert_eq(
+    100.of{|n| binomial(2*n, n) % binomial(2*n - 2, n-1) }.slice(1),
+    100.of{|n| binomialmod(2*n, n, binomial(2*n - 2, n-1)) }.slice(1),
+)
+
+# Morley’s congruence
+primes(5,100).each {|p|
+    assert((-1)**((p-1)/2) * binomialmod(p-1, (p-1)/2, p**3) -> is_congruent(powmod(4, p-1, p**3), p**3))
+}
+
+assert_eq(
+    100.of{|n| binomial(2*n, n) % ((n+1)*(n+2)) },
+    100.of{|n| binomialmod(2*n, n, (n+1)*(n+2)) },
+)
+
 say "** Test passed!"