Deadinsidick:
double h = 2 * area(r1,r2, O1O2) / O1O2;
        double OE = sqrt(r1 * r1 - h * h);
        //cout « O1O2 « " " « h « " " « OE « "\n";
        point E = O1 + (O2-O1) * (OE / O1O2);
        point vec = rot90(O2-O1);
        point M1 = E + ((h / O1O2) * vec);
        point M2 = E - ((h / O1O2) * vec);

p0.x=a.x+(b.x-a.x)*m/l;
    p0.y=a.y+(b.y-a.y)*m/l;
    point p1,p2;
    p1.x=p0.x+(b.y-a.y)/l*h;
    p1.y=p0.y-(b.x-a.x)/l*h;
    p2.x=p0.x-(b.y-a.y)/l*h;
    p2.y=p0.y+(b.x-a.x)/l*h;

double m=(r1*r1-r2*r2+l*l)/2/l;   // m-O1H