Numerically evaluate sine integral $ Si(x) $

Si.py

#!/usr/bin/env python3

from itertools import count

# This is stable up to about x=34.
# Above that, reimplement with rational arithmetic.

def Si(x):
    """Sine integral.  Numerically stable up to about x=34."""
    sign = +1
    fact = 1
    numerator = x
    acc = 0
    for n in count():
        tnp1 = 2 * n + 1
        if n:
            fact *= 2 * n * tnp1
        abs_term = numerator / (tnp1 * fact)
        if acc - abs_term == acc:
            # further terms won't affect the result.
            break
        acc += sign * abs_term
        numerator *= x**2
        sign = -sign
    return acc

if __name__ == '__main__':
    import numpy as np
    import matplotlib.pyplot as plt

    x = np.linspace(0, 15, 500)
    y = np.array([Si(xx) for xx in x])

    plt.plot(x, y)
    plt.show()