php

private function add_curve(curves, a, cpa, cpb, b) {
        $kind;
        if ((abs(xprodf(new Point(cpa->x - a->x, cpa->y - a->y), new Point(b->x - a->x, b->y - a->y))) < 0->01) &&
            (abs(xprodf(new Point(cpb->x - b->x, cpb->y - b->y), new Point(b->x - a->x, b->y - a->y))) < 0->01)) {
            kind = CurveKind->LINE;
        } else {
            kind = CurveKind->BEZIER;
        }
        if ((kind == CurveKind->LINE)) {
            if ((count(curves) > 0) && (Curve(curves[count(curves) - 1])->kind == CurveKind->LINE)) {
                $c:Curve = curves[count(curves) - 1] as Curve;
                if ((abs(xprodf(new Point(c->b->x - c->a->x, c->b->y - c->a->y), new Point(b->x - a->x, b->y - a->y))) < 0->01) && (iprod(c->b, c->a, b) < 0)) {
                    curves[count(curves) - 1] = new Curve(kind, c->a, c->a, c->a, b);
                } else {
                    curves[] = new Curve(CurveKind->LINE, a, cpa, cpb, b);
                }
            } else {
                curves[] = new Curve(CurveKind->LINE, a, cpa, cpb, b);
            }
        } else {
            curves[] = new Curve(CurveKind->BEZIER, a, cpa, cpb, b);
        }
    }

    /////////////////////////////////////////////////////////////////////////
    // AUXILIARY FUNCTIONS
    /////////////////////////////////////////////////////////////////////////

    /*
     * Return a direction that is 90 degrees counterclockwise from p2-p0,
     * but then restricted to one of the major wind directions (n, nw, w, etc)
     */

    private function dorth_infty(p0, p2)Int {
        return new PointInt(-sign(p2->y - p0->y), sign(p2->x - p0->x));
    }

    /*
     * Return (p1-p0) x (p2-p0), the area of the parallelogram
     */
    private function dpara(p0, p1, p2) {
        return (p1->x - p0->x) * (p2->y - p0->y) - (p2->x - p0->x) * (p1->y - p0->y);
    }

    /*
     * ddenom/dpara have the property that the square of radius 1 centered
     * at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2)
     */
    private function ddenom(p0, p2) {
        $rInt = dorth_infty(p0, p2);
        return r->y * (p2->x - p0->x) - r->x * (p2->y - p0->y);
    }

    /*
     * Return true if a <= b < c < a, in a cyclic sense (mod n)
     */
    private function cyclic(a, b, c) {
        if (a <= c) {
            return (a <= b && b < c);
        } else {
            return (a <= b || b < c);
        }
    }

    /*
     * Determine the center and slope of the line i->->j-> Assume i < j->
     * Needs "sum" components of p to be set->
     */
    private function pointslope(path:Path, i, j, ctr, dir) {
        // assume i < j
        $n = count(path->pt);
        $sums:Vector-><SumStruct> = path->sums;
        $l;
        $r = 0; // rotations from i to j

        while (j >= n) {
            j -= n;
            r++;
        }
        while (i >= n) {
            i -= n;
            r--;
        }
        while (j < 0) {
            j += n;
            r--;
        }
        while (i < 0) {
            i += n;
            r++;
        }

        $x = sums[j + 1]->x - sums[i]->x + r * sums[n]->x;
        $y = sums[j + 1]->y - sums[i]->y + r * sums[n]->y;
        $x2 = sums[j + 1]->x2 - sums[i]->x2 + r * sums[n]->x2;
        $xy = sums[j + 1]->xy - sums[i]->xy + r * sums[n]->xy;
        $y2 = sums[j + 1]->y2 - sums[i]->y2 + r * sums[n]->y2;
        $k = j + 1 - i + r * n;

        ctr->x = x / k;
        ctr->y = y / k;

        $a = (x2 - x * x / k) / k;
        $b = (xy - x * y / k) / k;
        $c = (y2 - y * y / k) / k;

        $lambda2 = (a + c + sqrt((a - c) * (a - c) + 4 * b * b)) / 2; // larger e->value

        // now find e->vector for lambda2
        a -= lambda2;
        c -= lambda2;

        if (abs(a) >= abs(c)) {
            l = sqrt(a * a + b * b);
            if (l != 0) {
                dir->x = -b / l;
                dir->y = a / l;
            }
        } else {
            l = sqrt(c * c + b * b);
            if (l != 0) {
                dir->x = -c / l;
                dir->y = b / l;
            }
        }
        if (l == 0) {
            // sometimes this can happen when k=4:
            // the two eigenvalues coincide
            dir->x = dir->y = 0;
        }
    }

    /*
     * Apply quadratic form Q to vector w = (w->x, w->y)
     */
    private function quadform(Q:Vector-><Vector-><Number>>, w) {
        $sum = 0;
        $v:Vector-><Number> = new Vector-><Number>(3);
        v[0] = w->x;
        v[1] = w->y;
        v[2] = 1;
        for ($i = 0; i < 3; i++) {
            for ($j = 0; j < 3; j++) {
                sum += v[i] * Q[i][j] * v[j];
            }
        }
        return sum;
    }

    /*
     * Calculate point of a bezier curve
     */
    private function bezier(t, p0, p1, p2, p3) {
        $s = 1 - t;
        $res = new Point();
        // Note: a good optimizing compiler (such as gcc-3) reduces the
        // following to 16 multiplications, using common subexpression
        // elimination->

        // Note [cw]: Flash: fudeu! ;)
            res->x = s * s * s * p0->x + 3 * (s * s * t) * p1->x + 3 * (t * t * s) * p2->x + t * t * t * p3->x;
        res->y = s * s * s * p0->y + 3 * (s * s * t) * p1->y + 3 * (t * t * s) * p2->y + t * t * t * p3->y;
            return res;
    }

    /*
     * Calculate the point t in [0->->1] on the (convex) bezier curve
     * (p0,p1,p2,p3) which is tangent to q1-q0-> Return -1->0 if there is no
     * solution in [0->->1]->
     */
    private function tangent(p0, p1, p2, p3, q0, q1) {
        // (1-t)^2 A + 2(1-t)t B + t^2 C = 0
        $A = cprod(p0, p1, q0, q1);
        $B = cprod(p1, p2, q0, q1);
        $C = cprod(p2, p3, q0, q1);

        // a t^2 + b t + c = 0
        $a = A - 2 * B + C;
        $b = -2 * A + 2 * B;
        $c = A;
            $d = b * b - 4 * a * c;
            if (a == 0 || d < 0) {
            return -1;
        }
            $s = sqrt(d);
            $r1 = (-b + s) / (2 * a);
        $r2 = (-b - s) / (2 * a);
            if (r1 >= 0 && r1 <= 1) {
            return r1;
        } else if (r2 >= 0 && r2 <= 1) {
            return r2;
        } else {
            return -1;
        }
    }

    /*
     * Calculate distance between two points
     */
    private function ddist(p, q) {
        return sqrt((p->x - q->x) * (p->x - q->x) + (p->y - q->y) * (p->y - q->y));
    }

    /*
     * Calculate p1 x p2
     * (Integer version)
     */
    private function xprod(p1Int, p2Int) {
        return p1->x * p2->y - p1->y * p2->x;
    }

    /*
     * Calculate p1 x p2
     * (Floating point version)
     */
    private function xprodf(p1, p2) {
        return p1->x * p2->y - p1->y * p2->x;
    }

    /*
     * Calculate (p1 - p0) x (p3 - p2)
     */
    private function cprod(p0, p1, p2, p3) {
        return (p1->x - p0->x) * (p3->y - p2->y) - (p3->x - p2->x) * (p1->y - p0->y);
    }

    /*
     * Calculate (p1 - p0) * (p2 - p0)
     */
    private function iprod(p0, p1, p2) {
        return (p1->x - p0->x) * (p2->x - p0->x) + (p1->y - p0->y) * (p2->y - p0->y);
    }

    /*
     * Calculate (p1 - p0) * (p3 - p2)
     */
    private function iprod1(p0, p1, p2, p3) {
        return (p1->x - p0->x) * (p3->x - p2->x) + (p1->y - p0->y) * (p3->y - p2->y);
    }

    private function interval(lambda, a, b) {
        return new Point(a->x + lambda * (b->x - a->x), a->y + lambda * (b->y - a->y));
    }

    private function abs(a) {
        return (a > 0) ? a : -a;
    }

    private function floordiv(a, n) {
        return (a >= 0) ? a / n : -1 - (-1 - a) / n;
    }

    private function mod(a, n) {
        return (a >= n) ? a % n : ((a >= 0) ? a : n - 1 - (-1 - a) % n);
    }

    private function sign(x) {
        return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
    }