Spring 2018 PMATH340 BQF Reduction Algorithm

bqf-reduction.py

#!/usr/bin/env python3

# Joseph Winters
# Reduction Algorithm
# July 2018

# This script implements the reduction of primitive, positive definite binary
# quadratic forms, which are those of the form:
#     ax^2 + bxy + cy^2
# having:
# * gcd(a,b,c) = 1 (primitive)
# * a > 0 and D = b^2 - 4ac < 0 (positive definite)

# This implementation is based on the version of the algorithm presented by
# Nick Rollick in PMATH 340 at the University of Waterloo during the Spring
# 2018 term.

import sys
import math

USAGE_MESSAGE = "Usage: python bqf-reduction.py a b c"
ERROR_WRONG_NUMBER_OF_ARGUMENTS = "The wrong number of arguments was entered.\n" + USAGE_MESSAGE
ERROR_ARGUMENTS_NOT_INTEGERS = "The given arguments are not integers.\n" + USAGE_MESSAGE
ERROR_NOT_PRIMITIVE = "The given integers do not represent a primitive bqf.\n" + USAGE_MESSAGE
ERROR_NOT_POSITIVE_DEFINITE = "The given integers do not represent a positive definite bqf.\n" + USAGE_MESSAGE

def meets_condition1(a, b, c):
    return abs(b) <= a

def meets_condition2(a, b, c):
    return a <= c

def apply_a_m(old_a, old_b, old_c):
    b_not = old_b % (2 * old_a)
    if b_not > old_a:
        b_not -= 2 * old_a
    m = (b_not - old_b) // (2 * old_a)  # NB: This division is exact.
    return old_a, b_not, (old_a*m*m + old_b*m + old_c)

def apply_s(old_a, old_b, old_c):
    return old_c, -old_b, old_a

def main():
    # Need three integers as input.
    if len(sys.argv) != 4:
        print(ERROR_WRONG_NUMBER_OF_ARGUMENTS)
        sys.exit(1)
    try:
        a = int(sys.argv[1])
        b = int(sys.argv[2])
        c = int(sys.argv[3])
    except ValueError:
        print(ERROR_ARGUMENTS_NOT_INTEGERS)
        sys.exit(1)
    # The bqf represented by the input integers must be primitive.
    if math.gcd(a, math.gcd(b, c)) != 1:
        print(ERROR_NOT_PRIMITIVE)
        sys.exit(1)
    # The bqf represented by the input integers must be primitive.
    if (a <= 0) or (b*b - 4*a*c >= 0):
        print(ERROR_NOT_POSITIVE_DEFINITE)
        sys.exit(1)
    # Main loop of the algorithm.
    while not (meets_condition1(a,b,c) and meets_condition2(a,b,c)):
        if not meets_condition1(a,b,c):
            a, b, c = apply_a_m(a, b, c)
        if not meets_condition2(a,b,c):
            a, b, c = apply_s(a, b, c)
    # Cleaning up final conditions.
    if (a == c) and (b < 0):
        a, b, c = apply_s(a, b, c)
    elif (abs(b) == a) and (b < 0):
        a, b, c = a, a, c
    return a, b, c

if __name__ == "__main__":
    a, b, c = main()
    print("{}x^2 + {}xy + {}y^2".format(a, b, c))