123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186 |
- #!/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)]
|