File:Cubic graph special points repeated.svg

#!/usr/bin/env python
import re
def eval_at(x): return a * x ** 3 + b * x ** 2 + c * x + d
def is_integer(x): return abs(x - int(x)) < 1e-5
def print_keyval(dic, re_key):
 print('\t'.join(sorted(['%s=%s' % (key, dic[key])
                         for key in dic if re.search(re_key, key)])))
## http://wolframalpha.com/input?i=expand+%28x-r1%29%28x-r2%29%5E2
## http://wolframalpha.com/input?i=differentiate+(x-r1)(x-r2)(x-r2)
## http://wolframalpha.com/input?i=solve+2r1*r2-2r1*x%2Br2^2-4r2*x%2B3x^2%3D0+for+x
## http://wolframalpha.com/input?i=differentiate+(x-r2)(x-(2r1%2Br2)%2F3)
## http://wolframalpha.com/input?i=solve+2x-2r1%2F3-4r2%2F3%3D0+for+x
## http://wolframalpha.com/input?i=y%3D(x-4)(x-1)^2
## http://www.wolframalpha.com/input?i=factorise+3x^2-12x%2B9

## Cubic bezier control points from left to right, to be populated
pts = [[-1,None], [None,None], [None,None], [5,None]]
n   = 4 ## search range = [-n, n]

for  r2 in range(n, -1 - n, -1):
 for r1 in range(n, r2    , -1):
  if r1 * r2 == 0: continue
  r_1 = (2 * r1 + r2) / 3.0
  r_2 = r2
  if is_integer(r_1) and r_1 != 0: r_1 = int(r_1)
  else:                            continue
  r__ = (r1 + 2 * r2) / 3.0
  if is_integer(r__) and r__ != 0: r__ = int(r__)
  else:                            continue
  if len(set([r1, r_1, r__])) < 3: continue
  a = 1
  b = -2 * r2 - r1
  c = r2 * (2 * r1 + r2)
  d = -r1 * r2 * r2
  p = pts[0][0]
  q = pts[3][0]
  pts[0][1] = eval_at(p)
  pts[3][1] = eval_at(q)
  pts[1][0] = (2 * p + q) / 3.0
  pts[1][1] = a * p * p * q + b * (p * p + 2 * p * q) / 3.0 + c * (2 * p + q) / 3.0 + d
  pts[2][0] = (p + 2 * q) / 3.0
  pts[2][1] = a * p * q * q + b * (q * q + 2 * p * q) / 3.0 + c * (p + 2 * q) / 3.0 + d
  print_keyval(locals(), r'^(r_*\d?|[abcd]|pts)$')