Python library or algorithm to generate arc geometry from three coordinate pairs?

Shapely PyQGIS, and GeoDjango use the same API based on the GEOS library:

With Shapely:

with a list of points:

from shapely.geometry import Point, LineString, mapping
pt1 = Point(0,0)
pt2 = Point(20,20)
pt3 = Point (50,50)
line = LineString([pt1,pt2,pt3])
#GeoJSON format
print mapping(line)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (20.0, 20.0), (50.0, 50.0))}

or with a list of coordinates:

line = LineString([(0, 0), (20,20),(50,50)])
#GeoJSON format
print mapping(line)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (20.0, 20.0), (50.0, 50.0))}

with PyQGIS:

with a list of points:

line = QgsGeometry.fromPolyline([QgsPoint(0,0),QgsPoint(20,20),QgsPoint(50,50)])
#GeoJSON format
print line.exportToGeoJSON()
{ "type": "LineString", "coordinates": [ [0, 0], [20, 20], [50, 50] ] }

with GeoDjango (look at Django: GEOS API)

with a list of coordinates:

from django.contrib.gis.geos import LineString
line = LineString((0, 0), (20, 20), (50, 50))
# GeoJSON format
print line.json
{ "type": "LineString", "coordinates": [ [ 0.0, 0.0 ], [ 20.0, 20.0 ], [ 50.0, 50.0 ] ] }

The resulting line is made up of two segments: (0,0) to (20,20) and (20,20) to (50,50)

Building on @gene 's answer above if you're looking to approximate a curve as a sequence of points using spline interpolation. In Python you can do this through the scipy.interpolate library. Particularly 1d interpolation.

For example

import numpy as np
import scipy.interpolate

coords = np.array([[0, 0], [25, 10], [50, 50]])

#The curve fits as a quadratic equation on three points
f = scipy.interpolate.interp1d(coords[:, 0], coords[:, 1], kind='quadratic')

#New points will be evenly distributed along x
new_x = np.linspace(np.min(coords[:, 0]), np.max(coords[:, 0]), 10)
new_y = f(new_x)

new_coords = np.vstack([new_x, new_y]).T

With a bit more work (play with shapely.interpolate) you can get a curve in segments of equal like so:

from shapely.geometry import Point, LineString, mapping

fine_x = np.linspace(np.min(coords[:, 0]), np.max(coords[:, 0]), 1000)
fine_y = f(fine_x)

fine_coords = zip(fine_x, fine_y)
fine_line = LineString(fine_coords)

even_line = LineString([
    np.array(fine_line.interpolate(i))
    for i in np.arange(0, fine_line.length, 5) #points 5 units apart
])