Inverse survival function of standard normal: pretty-print the polynomial approximation

inverse_survival_function.py

print("""What is the "inverse survival function" of the standard normal?

The cdf is based on erf(x), the survival function is based on erfc(x),
and the inverse survival is based on 1/erfc(x). The error function
erf(x) is not an elementary function, but it can be approximated as
1 - 1 / (p(x))^4 , where p(x) is a 4th degree polynomial in x, with
coefficients a = 0.278393, b = 0.230389, c = 0.000972, d = 0.078108.
[https://en.wikipedia.org/wiki/Error_function] This means that
1/erfc(x) can be approximated as (p(x))^4 .

Specifically, (1 + a x + b x^2 + c x^3 + d x^4)^4 .

Wolfram|Alpha will expand that for us, but it's not nicely organized
by exponent of x.

Hence :

""")

wolfram_output = "a^4 x^4 + 4 a^3 b x^5 + 4 a^3 c x^6 + 4 a^3 d x^7 + 4 a^3 x^3 + 6 a^2 b^2 x^6 + 12 a^2 b c x^7 + 12 a^2 b d x^8 + 12 a^2 b x^4 + 6 a^2 c^2 x^8 + 12 a^2 c d x^9 + 12 a^2 c x^5 + 6 a^2 d^2 x^10 + 12 a^2 d x^6 + 6 a^2 x^2 + 4 a b^3 x^7 + 12 a b^2 c x^8 + 12 a b^2 d x^9 + 12 a b^2 x^5 + 12 a b c^2 x^9 + 24 a b c d x^10 + 24 a b c x^6 + 12 a b d^2 x^11 + 24 a b d x^7 + 12 a b x^3 + 4 a c^3 x^10 + 12 a c^2 d x^11 + 12 a c^2 x^7 + 12 a c d^2 x^12 + 24 a c d x^8 + 12 a c x^4 + 4 a d^3 x^13 + 12 a d^2 x^9 + 12 a d x^5 + 4 a x + b^4 x^8 + 4 b^3 c x^9 + 4 b^3 d x^10 + 4 b^3 x^6 + 6 b^2 c^2 x^10 + 12 b^2 c d x^11 + 12 b^2 c x^7 + 6 b^2 d^2 x^12 + 12 b^2 d x^8 + 6 b^2 x^4 + 4 b c^3 x^11 + 12 b c^2 d x^12 + 12 b c^2 x^8 + 12 b c d^2 x^13 + 24 b c d x^9 + 12 b c x^5 + 4 b d^3 x^14 + 12 b d^2 x^10 + 12 b d x^6 + 4 b x^2 + c^4 x^12 + 4 c^3 d x^13 + 4 c^3 x^9 + 6 c^2 d^2 x^14 + 12 c^2 d x^10 + 6 c^2 x^6 + 4 c d^3 x^15 + 12 c d^2 x^11 + 12 c d x^7 + 4 c x^3 + d^4 x^16 + 4 d^3 x^12 + 6 d^2 x^8 + 4 d x^4 + 1"

terms = wolfram_output.split(" + ")

MAX_EXPONENT = 16

# Initialize empty list of lists.
polynomial = []
for i in range(MAX_EXPONENT + 1):
    polynomial += [[]]  # [[]] * 17 won't work right

# Populate the list of lists, carefully keeping things ordered by the
# degree of each term.
for t in terms:
    if not 'x' in t:
        coefficient_str = t
        degree_str = '0'
    elif not '^' in t:
        coefficient_str, degree_str = t.split(' x')
        degree_str = '1'  # reassign
    else:
        coefficient_str, degree_str = t.split(' x^')
    degree = int(degree_str)
    polynomial[degree] += [coefficient_str]

# Print the polynomial, in order of descending powers of x.
for degree in range(MAX_EXPONENT, -1, -1):
    string_list = polynomial[degree]
    coefficient_str = ' + '.join(string_list)
    if degree == 0:
        print(coefficient_str)  # Just say the constant, not x^0.
    elif degree == 1:
        print("(%s) x +" % (coefficient_str))  # Just say x, not x^1.
    else:
        print("(%s) x^%i +" % (coefficient_str, degree))