Fast quadratic polynomial fit from three points

poly2fit_from_3points.py

def poly2fit_from_3points(xs,ys):
    # p0,p1,p2
    # quadratic function:
    # y=a*x*x + b*x + c
    #
    # [1] p0y=a*p0x*p0x + b*p0x + c
    # [2] p1y=a*p1x*p1x + b*p1x + c
    # [3] p2y=a*p2x*p2x + b*p2x + c
    # [4] isolate c in [1]:    c = p0y - a*p0x*p0x - b*p0x
    # [5] insert [4]c in [2]: p1y = a*p1x*p1x + b*p1x + p0y - a*p0x*p0x - b*p0x
    # [6] swap eq sides of [5]: a*p1x*p1x + b*p1x + p0y - a*p0x*p0x - b*p0x = p1y
    # [7] isolate a terms in [6]: a*p1x*p1x - a*p0x*p0x = p1y - p0y + b*(p0x - p1x)
    # [8] isolate a in [7]: a = (p1y - p0y + b*(p0x - p1x)) / (p1x*p1x - p0x*p0x)
    # note: (p1x*p1x - p0x*p0x) != 0
    # note: p1x*p1x != p0x*p0x
    # [9] insert [4]c in [3]: p2y=a*p2x*p2x + b*p2x + p0y - a*p0x*p0x - b*p0x
    # [10] reduce occurences of a & b in [9]: p2y=a*(p2x*p2x - p0x*p0x) + b*(p2x-p0x) + p0y
    # [11] insert [8]a in [10]: p2y=((p1y - p0y + b*(p0x - p1x)) / (p1x*p1x - p0x*p0x))*(p2x*p2x - p0x*p0x) + b*(p2x-p0x) + p0y
    # [12] isolate b in [11]: 
    # p2y=((p1y - p0y + b*(p0x - p1x)) / (p1x*p1x - p0x*p0x))*(p2x*p2x - p0x*p0x) + b*(p2x-p0x) + p0y
    # p2y - p0y=(((p1y - p0y) + b*(p0x - p1x) )*(p2x*p2x - p0x*p0x)) / (p1x*p1x - p0x*p0x) + b*(p2x-p0x) 
    #           (((p1y - p0y) + b*(p0x - p1x) )*(p2x*p2x - p0x*p0x)) 
    #           = (p1y - p0y)*(p2x*p2x - p0x*p0x) + b*(p0x - p1x)*(p2x*p2x - p0x*p0x)
    # p2y - p0y=(((p1y - p0y)*(p2x*p2x - p0x*p0x) + b*(p0x - p1x)*(p2x*p2x - p0x*p0x)) / (p1x*p1x - p0x*p0x) + b*(p2x-p0x) 
    # p2y - p0y=(p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) + b*(p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) + b*(p2x-p0x) 
    # p2y - p0y - (p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) = + b*(p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) + b*(p2x-p0x) 
    # p2y - p0y - (p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) = b * ((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x))
    # b = (p2y - p0y - (p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x)) / ((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x))
    # note: ((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x)) != 0
    #
    # b = (p2y - p0y - (p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x)) / ((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x))
    # a = (p1y - p0y + b*(p0x - p1x)) / (p1x*p1x - p0x*p0x)
    # c = p0y - a*p0x*p0x - b*p0x
    
    p0x,p1x,p2x = xs
    p0y,p1y,p2y = ys
    #p0x,p0y = p0
    #p1x,p1y = p1
    #p2x,p2y = p2
    #print("p0",p0)
    #print("p1",p1)
    #print("p2",p2)
    assert(p1x*p1x != p0x*p0x);
    #print((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x))
    assert((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x) != 0)
    b = (p2y - p0y - (p1y - p0y)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x)) / ((p2x-p0x) + (p0x - p1x)*(p2x*p2x - p0x*p0x)/(p1x*p1x - p0x*p0x))
    a = (p1y - p0y + b*(p0x - p1x)) / (p1x*p1x - p0x*p0x)
    c = p0y - a*p0x*p0x - b*p0x
    #print("p0y - (a*p0x*p0x + b*p0x + c)", p0y - (a*p0x*p0x + b*p0x + c))
    #print("p1y - (a*p1x*p1x + b*p1x + c)", p1y - (a*p1x*p1x + b*p1x + c))
    #print("p2y - (a*p2x*p2x + b*p2x + c)", p2y - (a*p2x*p2x + b*p2x + c))
    
    return (a,b,c)
  
# foo = list(zip(Xs,Ys))
# %timeit poly2fit_from_3points(foo)
# 7.47 µs ± 77.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
#
# %timeit np.polyfit(Xs,Ys,2)
# 131 µs ± 1.67 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)