new file: Math/partial_sums_of_core_function.sf

Author: trizen

Math/eisenstein_integers.sf

@@ -15,6 +15,10 @@ class Eisenstein(a, b) {   # represents: a + b*w, where w = (-1 + i*sqrt(3))/2
         (a == c.a) && (b == c.b)
     }
 
+    method add (Number z) {
+        Eisenstein(a+z, b)
+    }
+
     method add (Eisenstein z) {
         var (c,d) = (z.a, z.b)
         Eisenstein(a+c, b+d)
@@ -33,6 +37,10 @@ class Eisenstein(a, b) {   # represents: a + b*w, where w = (-1 + i*sqrt(3))/2
         a*a - a*b + b*b
     }
 
+    method mul (Number z) {
+        Eisenstein(a*z, b*z)
+    }
+
     method mul (Eisenstein z) {
         var (c,d) = (z.a, z.b)
         Eisenstein(a*c - b*d, b*c + a*d - b*d)

Math/partial_sums_of_core_function.sf

@@ -0,0 +1,54 @@
+#!/usr/bin/ruby
+
+# Sub-linear algorithm for computing partial sums of the core(n) function.
+
+# OEIS sequence:
+#   https://oeis.org/A069891
+
+# Sequence A069891(10^n):
+#   S(10^1)  = 38
+#   S(10^2)  = 3233
+#   S(10^3)  = 328322
+#   S(10^4)  = 32926441
+#   S(10^5)  = 3289873890
+#   S(10^6)  = 328984021545
+#   S(10^7)  = 32898872196712
+#   S(10^8)  = 3289866649713946
+#   S(10^9)  = 328986818422458525
+#   S(10^10) = 32898680588469254505
+#   S(10^11) = 3289868138800129869623
+
+func squarefree_sum(n) {    # A066779
+    sum(1..n.isqrt, {|k|
+        moebius(k) * k*k * faulhaber(idiv(n, k*k), 1)
+    })
+}
+
+func core_partial_sum(n) {  # A069891
+    n.dirichlet_sum(
+        {|k| k.is_square ? 1 : 0 },
+        {|k| abs(k*mu(k)) },
+        {|k| k.isqrt },
+        {|k| squarefree_sum(k) },
+    )
+}
+
+func f(n) {     # A046970
+    n.factor_prod {|p|
+        1 - p*p
+    }
+}
+
+func core_partial_sum_faster(n) {   # A069891
+    sum(1..n.isqrt, {|k|
+        f(k) * faulhaber(idiv(n, k*k), 1)
+    })
+}
+
+say core_partial_sum(10**5)             #=> 3289873890
+say core_partial_sum_faster(10**5)      #=> 3289873890
+
+assert_eq(
+    100.of(core_partial_sum),
+    100.of(core_partial_sum_faster)
+)

Math/quadratic_integers.sf

@@ -5,66 +5,66 @@
 # See also:
 #   https://en.wikipedia.org/wiki/Quadratic_integer
 
-class Quadratic(a, b, w = 2) {   # represents: a + b*sqrt(w)
+class QuadraticInteger(a, b, w = 2) {   # represents: a + b*sqrt(w)
 
     method to_s {
-        "Quadratic(#{a}, #{b}, #{w})"
+        "QuadraticInteger(#{a}, #{b}, #{w})"
     }
 
-    method ==(Quadratic c) {
+    method ==(QuadraticInteger c) {
         (a == c.a) && (b == c.b) && (w == c.w)
     }
 
     method conjugate {
-        Quadratic(a, -b, w)
+        QuadraticInteger(a, -b, w)
     }
 
     method norm {
         a*a - w*b*b
     }
 
     method add (Number c) {
-        Quadratic(a+c, b, w)
+        QuadraticInteger(a+c, b, w)
     }
 
-    method add (Quadratic z) {
+    method add (QuadraticInteger z) {
         var (c,d) = (z.a, z.b)
-        Quadratic(a+c, b+d, w)
+        QuadraticInteger(a+c, b+d, w)
     }
 
     __CLASS__.alias_method(:add, '+')
 
     method sub (Number c) {
-        Quadratic(a-c, b, w)
+        QuadraticInteger(a-c, b, w)
     }
 
-    method sub (Quadratic z) {
+    method sub (QuadraticInteger z) {
         var (c,d) = (z.a, z.b)
-        Quadratic(a-c, b-d, w)
+        QuadraticInteger(a-c, b-d, w)
     }
 
     __CLASS__.alias_method(:sub, '-')
 
     method mul (Number c) {
-        Quadratic(a*c, b*c, w)
+        QuadraticInteger(a*c, b*c, w)
     }
 
-    method mul (Quadratic z) {
+    method mul (QuadraticInteger z) {
         var (c,d) = (z.a, z.b)
-        Quadratic(a*c + b*d*w, a*d + b*c, w)
+        QuadraticInteger(a*c + b*d*w, a*d + b*c, w)
     }
 
     __CLASS__.alias_method(:mul, '*')
 
     method mod (Number m) {
-        Quadratic(a % m, b % m, w)
+        QuadraticInteger(a % m, b % m, w)
     }
 
     __CLASS__.alias_method(:mod, '%')
 
     method pow(Number n) {
         var x = self
-        var c = Quadratic(1, 0, w)
+        var c = QuadraticInteger(1, 0, w)
 
         for bit in (n.digits(2)) {
             c *= x if bit
@@ -79,7 +79,7 @@ class Quadratic(a, b, w = 2) {   # represents: a + b*sqrt(w)
     method powmod(Number n, Number m) {
 
         var x = self
-        var c = Quadratic(1, 0, w)
+        var c = QuadraticInteger(1, 0, w)
 
         for bit in (n.digits(2)) {
             (c *= x) %= m if bit        #=
@@ -100,11 +100,11 @@ func is_quadratic_pseudoprime (n, r=2) {
 
     return true if (r <= 0)
 
-    var x = Quadratic(r, 1, r+2).powmod(n, n)
+    var x = QuadraticInteger(r, 1, r+2).powmod(n, n)
 
     x.a == r || return false
 
-    var y = Quadratic(r, -1, r+2).powmod(n, n)
+    var y = QuadraticInteger(r, -1, r+2).powmod(n, n)
 
     y.a == r || return false
 
@@ -114,20 +114,20 @@ func is_quadratic_pseudoprime (n, r=2) {
 say is_quadratic_pseudoprime(43)    #=> true
 say is_quadratic_pseudoprime(97)    #=> true
 
-with (Quadratic(1, 1, 2)) {|q|
+with (QuadraticInteger(1, 1, 2)) {|q|
     say 15.of { q.pow(_).a }        #=> A001333
     say 15.of { q.pow(_).b }        #=> A000129
 }
 
-with (Quadratic(1, 1, 3)) {|q|
+with (QuadraticInteger(1, 1, 3)) {|q|
     say 15.of { q.pow(_).a }        #=> A026150
     say 15.of { q.pow(_).b }        #=> A002605
 }
 
 var n = (274177-1)
 var m = (2**64 + 1)
 
-with (Quadratic(3, 4, 2)) {|q|
+with (QuadraticInteger(3, 4, 2)) {|q|
     var r = q.powmod(n, m)
     say gcd(r.a-1, m)       #=> 2741177
     say gcd(r.b, m)         #=> 2741177
@@ -140,7 +140,7 @@ func is_gaussian_quadratic_pseudoprime(n) {
     return false if (n <= 1)
     return true  if (n <= 3)
 
-    static x = Quadratic(1, -1, -2)
+    static x = QuadraticInteger(1, -1, -2)
 
     given (n%8) {
         when ([1,3]) {

Math/quaternion_integer_primality_test.sf

@@ -28,15 +28,15 @@
 #   4765950001, 104510424961, 373523673601, 1240013784241, 4575844294561, 4866233601601
 
 # See also:
-#   https://en.wikipedia.org/wiki/Quaternion
+#   https://en.wikipedia.org/wiki/QuaternionInteger
 
-class Quaternion(a, b=0, c=0, d=0) {
+class QuaternionInteger(a, b=0, c=0, d=0) {
 
     method to_s {
-        "Quaternion(#{a}, #{b}, #{c}, #{d})"
+        "QuaternionInteger(#{a}, #{b}, #{c}, #{d})"
     }
 
-    method ==(Quaternion z) {
+    method ==(QuaternionInteger z) {
         (a == z.a) && (b == z.b) && (c == z.c) && (d == z.d)
     }
 
@@ -62,7 +62,7 @@ class Quaternion(a, b=0, c=0, d=0) {
             )
         }
 
-        return Quaternion(a2,b2,c2,d2)
+        return QuaternionInteger(a2,b2,c2,d2)
     }
 }
 
@@ -79,13 +79,13 @@ func quaternion_primality_test(n, tries = 3) {
         return n.is_prime
     }
 
-    var z = Quaternion(a,b,c,d)
+    var z = QuaternionInteger(a,b,c,d)
     var r = z.powmod(n,n)
 
     (
         (r == z) ||
-        (r == Quaternion(a)) ||
-        (r == Quaternion(a, n-b, n-c, n-d))
+        (r == QuaternionInteger(a)) ||
+        (r == QuaternionInteger(a, n-b, n-c, n-d))
     ) && (
         (tries > 0) ? __FUNC__(n, tries-1) : true
     )

Math/quaternion_integers.sf

@@ -1,4 +1,4 @@
-#!/usr/bin/ruby
+7#!/usr/bin/ruby
 
 # Simple implementation of quaternion integers.
 
@@ -88,7 +88,7 @@ var m = (2**64 + 1)
 
 with (Quaternion(2,3,4,5)) {|q|
     var r = q.powmod(n, m)
-    say gcd(r.b, m)         #=> 2741177
-    say gcd(r.c, m)         #=> 2741177
-    say gcd(r.d, m)         #=> 2741177
+    say gcd(r.b, m)         #=> 274177
+    say gcd(r.c, m)         #=> 274177
+    say gcd(r.d, m)         #=> 274177
 }

README.md

@@ -355,6 +355,7 @@ A nice collection of day-to-day Sidef scripts.
     * [Order factorization method](./Math/order_factorization_method.sf)
     * [Ore's harmonic numbers](./Math/ore_s_harmonic_numbers.sf)
     * [Partial sums of 2 to the bigomega of n](./Math/partial_sums_of_2_to_the_bigomega_of_n.sf)
+    * [Partial sums of core function](./Math/partial_sums_of_core_function.sf)
     * [Partial sums of dedekind psi function](./Math/partial_sums_of_dedekind_psi_function.sf)
     * [Partial sums of dedekind psi function recursive](./Math/partial_sums_of_dedekind_psi_function_recursive.sf)
     * [Partial sums of euler totient function](./Math/partial_sums_of_euler_totient_function.sf)