new file: Math/achilles_numbers.sf

new file:   Math/square_form_factorization_method.sf

Author: trizen

Math/achilles_numbers.sf

@@ -0,0 +1,41 @@
+#!/usr/bin/ruby
+
+# Generate Achilles and strong Achilles numbers.
+
+# See also:
+#   https://rosettacode.org/wiki/Achilles_numbers
+
+var P = 2.powerful(1e6)
+var achilles = Set(P.grep{ !.is_power }...)
+var strong_achilles = achilles.grep { achilles.has(.phi) }
+
+say "First 50 Achilles numbers:"
+achilles.sort.first(50).slices(10).each { .map{'%4s'%_}.join(' ').say }
+
+say "\nFirst 30 strong Achilles numbers:"
+strong_achilles.sort.first(30).slices(10).each { .map{'%5s'%_}.join(' ').say }
+
+say "\nNumber of Achilles numbers with:"
+achilles.to_a.group_by{.len}.sort_by{|k| k.to_i }.each_2d{|a,b|
+    say "#{a} digits: #{b.len}"
+}
+
+__END__
+First 50 Achilles numbers:
+  72  108  200  288  392  432  500  648  675  800
+ 864  968  972 1125 1152 1323 1352 1372 1568 1800
+1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
+3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
+5324 5400 5408 5488 6075 6125 6272 6728 6912 7200
+
+First 30 strong Achilles numbers:
+  500   864  1944  2000  2592  3456  5000 10125 10368 12348
+12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
+49392 50000 52488 55296 61731 64827 67500 69984 78608 80000
+
+Number of Achilles numbers with:
+2 digits: 1
+3 digits: 12
+4 digits: 47
+5 digits: 192
+6 digits: 664

Math/square_form_factorization_method.sf

@@ -0,0 +1,76 @@
+#!/usr/bin/ruby
+
+# Daniel Shanks's Square Form Factorization (SquFoF).
+
+# See also:
+#   https://rosettacode.org/wiki/Square_form_factorization
+
+const multipliers = divisors(3*5*7*11).grep { _ > 1 }
+
+func sff(N) {
+
+    N.is_prime  && return 1         # n is prime
+    N.is_square && return N.isqrt   # n is square
+
+    multipliers.each {|k|
+
+        var P0 = isqrt(k*N)       # P[0]=floor(sqrt(N)
+        var Q0 = 1                # Q[0]=1
+        var Q  = (k*N - P0*P0)    # Q[1]=N-P[0]^2 & Q[i]
+        var P1 = P0               # P[i-1] = P[0]
+        var Q1 = Q0               # Q[i-1] = Q[0]
+        var P  = 0                # P[i]
+        var Qn = 0                # P[i+1]
+        var b  = 0                # b[i]
+
+        while (!Q.is_square) {              # until Q[i] is a perfect square
+            b = idiv(isqrt(k*N) + P1, Q)    # floor(floor(sqrt(N+P[i-1])/Q[i])
+            P = (b*Q - P1)                  # P[i]=b*Q[i]-P[i-1]
+            Qn = (Q1 + b*(P1 - P))          # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
+            (Q1, Q, P1) = (Q, Qn, P)
+        }
+
+        b  = idiv(isqrt(k*N) + P, Q)        # b=floor((floor(sqrt(N)+P[i])/Q[0])
+        P1 = (b*Q0 - P)                     # P[i-1]=b*Q[0]-P[i]
+        Q  = (k*N - P1*P1)/Q0               # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
+        Q1 = Q0                             # Q[i-1] = Q[0]
+
+        loop {
+            b  = idiv(isqrt(k*N) + P1, Q)   # b=floor(floor(sqrt(N)+P[i-1])/Q[i])
+            P  = (b*Q - P1)                 # P[i]=b*Q[i]-P[i-1]
+            Qn = (Q1 + b*(P1 - P))          # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
+            break if (P == P1)              # until P[i+1]=P[i]
+            (Q1, Q, P1) = (Q, Qn, P)
+        }
+        with (gcd(N,P)) {|g|
+            return g if g.is_ntf(N)
+        }
+    }
+
+    return 0
+}
+
+[  11111, 2501, 12851, 13289, 75301, 120787, 967009, 997417,  4558849,
+   7091569, 13290059, 42854447, 223553581, 2027651281,
+].each {|n|
+    var v = sff(n)
+    if    (v == 0) { say "The number #{n} is not factored." }
+    elsif (v == 1) { say "The number #{n} is a prime."      }
+    else           { say "#{n} = #{[n/v, v].sort.join(' * ')}" }
+}
+
+__END__
+11111 = 41 * 271
+2501 = 41 * 61
+12851 = 71 * 181
+13289 = 97 * 137
+75301 = 257 * 293
+120787 = 43 * 2809
+967009 = 601 * 1609
+997417 = 257 * 3881
+4558849 = 383 * 11903
+7091569 = 2663 * 2663
+13290059 = 3119 * 4261
+42854447 = 4423 * 9689
+223553581 = 11213 * 19937
+2027651281 = 44021 * 46061

README.md

@@ -124,6 +124,7 @@ A nice collection of day-to-day Sidef scripts.
     * [Sierpinski triangle](./Image/sierpinski_triangle.sf)
     * [Voronoi diagram](./Image/voronoi_diagram.sf)
 * Math
+    * [Achilles numbers](./Math/achilles_numbers.sf)
     * [AGM calculate pi](./Math/AGM_calculate_pi.sf)
     * [Aitken's array](./Math/aitken_s_array.sf)
     * [Akiyama-tanigawa numerators](./Math/akiyama-tanigawa_numerators.sf)
@@ -666,6 +667,7 @@ A nice collection of day-to-day Sidef scripts.
     * [Sqrt convergents](./Math/sqrt_convergents.sf)
     * [Square-full numbers](./Math/square-full_numbers.sf)
     * [Square congruence lookup factorization](./Math/square_congruence_lookup_factorization.sf)
+    * [Square form factorization method](./Math/square_form_factorization_method.sf)
     * [Square product subsets](./Math/square_product_subsets.sf)
     * [Square root arithmetic-harmonic mean](./Math/square_root_arithmetic-harmonic_mean.sf)
     * [Square root good rational approximations](./Math/square_root_good_rational_approximations.sf)