sidef-scripts/Math/gaussian_divisors.sf

#!/usr/bin/ruby

# Author: Trizen
# Date: 12 June 2022
# https://github.com/trizen

# Find the factors and divisors of a Gaussian integer.

# See also:
#   https://oeis.org/A125271
#   https://oeis.org/A078930
#   https://oeis.org/A078910
#   https://oeis.org/A078911
#   https://projecteuler.net/problem=153
#   https://www.alpertron.com.ar/GAUSSIAN.HTM
#   https://en.wikipedia.org/wiki/Table_of_Gaussian_integer_factorizations

func gaussian_factors(Gauss z) {

    var n = z.norm
    var factors = []

    for p,e in (n.factor_exp) {
        if (p == 2) {
            var t = Gauss(1,1)
            while (z % t == 0) {
                factors << t
                z /= t
            }
        }
        elsif (p % 4 == 3) {
            var t = Gauss(p)
            while (z % t == 0) {
                factors << t
                z /= t
            }
        }
        else {
            for a,b in (sum_of_squares(p)) {

                var t1 = Gauss(a,b)
                var t2 = Gauss(b,a)

                while (z % t1 == 0) {
                    factors << t1
                    z /= t1
                }

                while (z % t2 == 0) {
                    factors << t2
                    z /= t2
                }
            }
        }
    }

    if (z != Gauss(1)) {
        factors << z
    }

    return factors.sort.run_length
}

func gaussian_divisors(Gauss n) {

    var d  = [Gauss(1), Gauss(-1), Gauss(0,1), Gauss(0,-1)]

    for p,e in (gaussian_factors(n)) {
        if ((p.imag == 0) && (p.real.abs == 1)) {
            next
        }
        if ((p.real == 0) && (p.imag.abs == 1)) {
            next
        }
        var r = Gauss(1)
        d << gather {
            e.times {
                r *= p
                d.each { |u|
                    take(u*r)
                }
            }
        }...
    }

    return d
}

for n in (1..5) {
    say "GaussDivisors(#{n}) = #{gaussian_divisors(Gauss(n))}"
}

say gaussian_divisors(Gauss(440,-55)).len     #=> 96
say 30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos}.len }      #=> OEIS: A125271

#
## Run some tests
#

assert_eq(
    1..5 -> map{gaussian_divisors(Gauss(_)).grep{.real.is_pos}.sum},
    [Gauss(1), Gauss(5), Gauss(4), Gauss(13), Gauss(12)]
)

assert_eq(
    1..48 -> map{gaussian_divisors(Gauss(_)).grep{.real.is_pos}.sum}.sum,
    Gauss(3654)
)

assert_eq(
    1..5 -> map{gaussian_divisors(Gauss(0, _)).grep{.real.is_pos}.sum},
    [Gauss(1), Gauss(5), Gauss(4), Gauss(13), Gauss(12)]
)

assert_eq(  # OEIS: A125271
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos}.len },
    %n[1, 1, 4, 2, 7, 6, 8, 2, 10, 3, 20, 2, 14, 6, 8, 12, 13, 6, 12, 2, 34, 4, 8, 2, 20, 15, 20, 4, 14, 6]
)

assert_eq(  # OEIS: A078930
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos}.sum.real },
    %n[1, 1, 5, 4, 13, 12, 20, 8, 29, 13, 56, 12, 52, 24, 40, 48, 61, 28, 65, 20, 144, 32, 60, 24, 116, 81, 112, 40, 104, 44]
)

assert_eq(  # OEIS: A078930
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos}.map{.real}.sum },
    %n[1, 1, 5, 4, 13, 12, 20, 8, 29, 13, 56, 12, 52, 24, 40, 48, 61, 28, 65, 20, 144, 32, 60, 24, 116, 81, 112, 40, 104, 44]
)

assert_eq(  # OEIS: A078930
    30.of { gaussian_divisors(Gauss(_)).grep{.imag.is_pos}.map{.imag}.sum },
    %n[1, 1, 5, 4, 13, 12, 20, 8, 29, 13, 56, 12, 52, 24, 40, 48, 61, 28, 65, 20, 144, 32, 60, 24, 116, 81, 112, 40, 104, 44]
)

assert_eq(  # OEIS: A078911
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos && .imag.is_pos}.map{.imag}.sum },
    %n[0, 0, 1, 0, 3, 3, 4, 0, 7, 0, 19, 0, 12, 5, 8, 12, 15, 5, 13, 0, 51, 0, 12, 0, 28, 25, 35, 0, 24, 7]
)

assert_eq(  # OEIS: A078911
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos && .imag.is_pos}.map{.real}.sum },
    %n[0, 0, 1, 0, 3, 3, 4, 0, 7, 0, 19, 0, 12, 5, 8, 12, 15, 5, 13, 0, 51, 0, 12, 0, 28, 25, 35, 0, 24, 7]
)

assert_eq(  # OEIS: A078911
    30.of { gaussian_divisors(Gauss(_)).grep{.real>=0 || .imag>=0}.sum.real },
    %n[0, 0, 1, 0, 3, 3, 4, 0, 7, 0, 19, 0, 12, 5, 8, 12, 15, 5, 13, 0, 51, 0, 12, 0, 28, 25, 35, 0, 24, 7]
)

assert_eq(  # OEIS: A078911
    30.of { gaussian_divisors(Gauss(_)).grep{.real>=0 && .imag>0}.sum.real },
    %n[0, 0, 1, 0, 3, 3, 4, 0, 7, 0, 19, 0, 12, 5, 8, 12, 15, 5, 13, 0, 51, 0, 12, 0, 28, 25, 35, 0, 24, 7]
)

assert_eq(  # OEIS: A078911
    30.of { gaussian_divisors(Gauss(_)).grep{.real>=0 || .imag>=0}.sum.imag },
    %n[0, 0, 1, 0, 3, 3, 4, 0, 7, 0, 19, 0, 12, 5, 8, 12, 15, 5, 13, 0, 51, 0, 12, 0, 28, 25, 35, 0, 24, 7]
)

assert_eq(  # OEIS: A078910
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos || .imag.is_pos}.map{.real}.sum },
    %n[1, 1, 4, 4, 10, 9, 16, 8, 22, 13, 37, 12, 40, 19, 32, 36, 46, 23, 52, 20, 93, 32, 48, 24, 88, 56, 77, 40, 80, 37]
)

assert_eq(  # OEIS: A078910
    30.of { gaussian_divisors(Gauss(_)).grep{.real.is_pos || .imag.is_pos}.map{.imag}.sum },
    %n[1, 1, 4, 4, 10, 9, 16, 8, 22, 13, 37, 12, 40, 19, 32, 36, 46, 23, 52, 20, 93, 32, 48, 24, 88, 56, 77, 40, 80, 37]
)

assert_eq(  # OEIS: A078910
    30.of { gaussian_divisors(Gauss(_)).grep{.real>0 && .imag>=0}.sum.real },
    %n[1, 1, 4, 4, 10, 9, 16, 8, 22, 13, 37, 12, 40, 19, 32, 36, 46, 23, 52, 20, 93, 32, 48, 24, 88, 56, 77, 40, 80, 37]
)

assert_eq(
    30.of { gaussian_divisors(Gauss(_)).sum },
    30.of { Gauss(0,0) }
)

__END__
GaussDivisors(1) = [Gauss(-1, 0), Gauss(0, -1), Gauss(0, 1), Gauss(1, 0)]
GaussDivisors(2) = [Gauss(-2, 0), Gauss(-1, -1), Gauss(-1, 0), Gauss(-1, 1), Gauss(0, -2), Gauss(0, -1), Gauss(0, 1), Gauss(0, 2), Gauss(1, -1), Gauss(1, 0), Gauss(1, 1), Gauss(2, 0)]
GaussDivisors(3) = [Gauss(-3, 0), Gauss(-1, 0), Gauss(0, -3), Gauss(0, -1), Gauss(0, 1), Gauss(0, 3), Gauss(1, 0), Gauss(3, 0)]
GaussDivisors(4) = [Gauss(-4, 0), Gauss(-2, -2), Gauss(-2, 0), Gauss(-2, 2), Gauss(-1, -1), Gauss(-1, 0), Gauss(-1, 1), Gauss(0, -4), Gauss(0, -2), Gauss(0, -1), Gauss(0, 1), Gauss(0, 2), Gauss(0, 4), Gauss(1, -1), Gauss(1, 0), Gauss(1, 1), Gauss(2, -2), Gauss(2, 0), Gauss(2, 2), Gauss(4, 0)]
GaussDivisors(5) = [Gauss(-5, 0), Gauss(-2, -1), Gauss(-2, 1), Gauss(-1, -2), Gauss(-1, 0), Gauss(-1, 2), Gauss(0, -5), Gauss(0, -1), Gauss(0, 1), Gauss(0, 5), Gauss(1, -2), Gauss(1, 0), Gauss(1, 2), Gauss(2, -1), Gauss(2, 1), Gauss(5, 0)]