square_product_subsets.sf

#!/usr/bin/ruby

# Find subsets of integers whose product is a square, using Gaussian elimination on a GF(2) matrix of vector exponents.

# Code inspired by:
#   https://github.com/martani/Quadratic-Sieve/blob/master/matrix.c

# See also:
#   https://btravers.weebly.com/uploads/6/7/2/9/6729909/quadratic_sieve_slides.pdf

func gaussian_elimination(A, n) {

    var m = A.end
    var I = (m+1).of {|i| 1 << i }

    var nrow = -1

    for col in (0 .. min(m, n)) {
        var npivot = -1

        for row in (nrow+1 .. m) {
            if (A[row].bit(col)) {
                npivot = row
                nrow++
                break
            }
        }

        next if (npivot < 0)

        if (npivot != nrow) {
            A.swap(npivot, nrow)
            I.swap(npivot, nrow)
        }

        for row in (nrow+1 .. m) {
            if (A[row].bit(col)) {
                A[row] ^= A[nrow]
                I[row] ^= I[nrow]
            }
        }
    }

    return I
}

func exponents_signature(factors) {
    var sig = 0

    for p,e in (factors) {
        sig.setbit!(p.prime_count-1) if e.is_odd
    }

    return sig
}

func find_square_subsets(Q) {

    var max_prime = 2

    var A = gather {
        Q.each { |n|
            var factors = n.factor_exp || next
            max_prime = (factors[-1][0] `max` max_prime)
            take(exponents_signature(factors))
        }
    }

    if (A.len < max_prime.prime_count) {
        A += (max_prime.prime_count - A.len -> of(0))
    }

    var I  = gaussian_elimination(A, max_prime.prime_count-1)
    var LR = (A.end - A.rindex_by { !.is_zero })

    var square_subsets = []

    for solution in (I.last(LR)) {

        var prod  = 1
        var terms = []

        for i in (^Q) {
            if (solution.bit(i)) {

                prod  *= Q[i]
                terms << Q[i]

                square_subsets << [terms...] if prod.is_square
            }
        }
    }

    return square_subsets
}

var Q = [10, 97, 24, 35, 75852, 54, 12, 13, 11,
         33, 37, 48, 57, 58, 63, 68, 377, 15,
         20, 26, 7, 3, 17, 29, 43, 41, 4171, 78]

#Q = [10, 24, 35, 52, 54, 78]

var S = find_square_subsets(Q)

S.each { |solution|
    say solution
}

__END__
[12, 48]
[10, 24, 35, 12, 63]
[24, 54]
[24, 12, 13, 58, 377]
[10, 24, 15]
[10, 24, 12, 20]
[24, 12, 13, 26]
[10, 24, 35, 12, 7]
[12, 3]
[68, 17]
[24, 12, 58, 29]
[75852, 43]
[12, 11, 33]
[97, 75852, 4171]
[24, 13, 78]