begin
const string (1) snl = "
"

routine tobern2(long real a0, long real a1, long real a2, long real name b0, long real name b1, long real name b2)
     b0 = a0
     b1 = a0 + a1 / 2
     b2 = a0 + a1 + a2
end

{ uses de Casteljau’s algorithm }
long real fn evalbern2(long real b0, long real b1, long real b2, long real t)
    long real s = 1.0 - t
    long real b01 = s * b0 + t * b1
    long real b12 = s * b1 + t * b2
    result = s * b01 + t * b12
end

routine tobern3(long real a0, long real a1, long real a2, long real a3, long real name b0, long real name b1, long real name b2, long real name b3)
    b0 = a0
    b1 = a0 + a1 / 3
    b2 = a0 + a1 * 2/3 + a2 / 3
    b3 = a0 + a1 + a2 + a3
end

{ uses de Casteljau’s algorithm }
long real fn evalbern3(long real b0, long real b1, long real b2, long real b3, long real t)
    long real s = 1 - t
    long real b01  = s * b0  + t * b1
    long real b12  = s * b1  + t * b2
    long real b23  = s * b2  + t * b3
    long real b012 = s * b01 + t * b12
    long real b123 = s * b12 + t * b23
    result = s * b012 + t * b123
end

routine bern2to3(long real q0, long real q1, long real q2, long real name c0, long real name c1, long real name c2, long real name c3)
    c0 = q0
    c1 = q0 / 3   + q1 * 2/3
    c2 = q1 * 2/3 + q2 / 3
    c3 = q2
end

{ uses Horner's rule }
long real fn evalmono2(long real a0, long real a1, long real a2, long real t)
    result = a0 + (t * (a1 + (t * a2)))
end

{ uses Horner's rule }
long real fn evalmono3(long real a0, long real a1, long real a2, long real a3, long real t)
    result = a0 + (t * (a1 + (t * (a2 + (t * a3)))))
end

    long real p0m, p1m, p2m, q0m, q1m, q2m, r0m, r1m, r2m, r3m
    long real p0b2, p1b2, p2b2, q0b2, q1b2, q2b2
    long real p0b3, p1b3, p2b3, p3b3, q0b3, q1b3, q2b3, q3b3, r0b3, r1b3, r2b3, r3b3
    long real pc0, pc1, pc2, pc3, qc0, qc1, qc2, qc3

    long real x, y, m

    p0m = 1; p1m = 0; p2m = 0
    q0m = 1; q1m = 2; q2m = 3
    r0m = 1; r1m = 2; r2m = 3; r3m = 4

    tobern2 (p0m, p1m, p2m, p0b2, p1b2, p2b2)
    tobern2 (q0m, q1m, q2m, q0b2, q1b2, q2b2)
    print string ("Subprogram (1) examples:".snl)
    print string ("  mono ( ")
    print fl (p0m, 4); print string (", ")
    print fl (p1m, 4); print string (", ")
    print fl (p2m, 4); print string (") --> bern ( ")
    print fl (p0b2, 4); print string (", ")
    print fl (p1b2, 4); print string (", ")
    print fl (p2b2, 4); print string (")".snl)
    print string ("  mono ( ")
    print fl (q0m, 4); print string (", ")
    print fl (q1m, 4); print string (", ")
    print fl (q2m, 4); print string (") --> bern ( ")
    print fl (q0b2, 4); print string (", ")
    print fl (q1b2, 4); print string (", ")
    print fl (q2b2, 4); print string (")".snl)

    print string ("Subprogram (2) examples:".snl)
    x = 0.25
    y = evalbern2 (p0b2, p1b2, p2b2, x)
    print string ("  p ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 7.50
    y = evalbern2 (p0b2, p1b2, p2b2, x)
    print string ("  p ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 0.25
    y = evalbern2 (q0b2, q1b2, q2b2, x)
    print string ("  q ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 7.50
    y = evalbern2 (q0b2, q1b2, q2b2, x)
    print string ("  q ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline

    tobern3 (p0m, p1m, p2m, 0, p0b3, p1b3, p2b3, p3b3)
    tobern3 (q0m, q1m, q2m, 0, q0b3, q1b3, q2b3, q3b3)
    tobern3 (r0m, r1m, r2m, r3m, r0b3, r1b3, r2b3, r3b3)
    print string ("Subprogram (3) examples:".snl)
    print string ("  mono ( ")
    print fl (p0m, 4); print string (", ")
    print fl (p1m, 4); print string (", ")
    print fl (p2m, 4); print string (", ")
    print fl (0, 4); print string (") --> bern ( ")
    print fl (p0b3, 4); print string (", ")
    print fl (p1b3, 4); print string (", ")
    print fl (p2b3, 4); print string (", ")
    print fl (p3b3, 4); print string (")".snl)
    print string ("  mono ( ")
    print fl (q0m, 4); print string (", ")
    print fl (q1m, 4); print string (", ")
    print fl (q2m, 4); print string (", ")
    print fl (0, 4); print string (") --> bern ( ")
    print fl (q0b3, 4); print string (", ")
    print fl (q1b3, 4); print string (", ")
    print fl (q2b3, 4); print string (", ")
    print fl (q3b3, 4); print string (")".snl)
    print string ("  mono ( ")
    print fl (r0m, 4); print string (", ")
    print fl (r1m, 4); print string (", ")
    print fl (r2m, 4); print string (", ")
    print fl (r3m, 4); print string (") --> bern ( ")
    print fl (r0b3, 4); print string (", ")
    print fl (r1b3, 4); print string (", ")
    print fl (r2b3, 4); print string (", ")
    print fl (r3b3, 4); print string (")".snl)

    print string ("Subprogram (4) examples:".snl)
    x = 0.25
    y = evalbern3 (p0b3, p1b3, p2b3, p3b3, x)
    print string ("  p ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 7.50
    y = evalbern3 (p0b3, p1b3, p2b3, p3b3, x)
    print string ("  p ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 0.25
    y = evalbern3 (q0b3, q1b3, q2b3, q3b3, x)
    print string ("  q ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 7.50
    y = evalbern3 (q0b3, q1b3, q2b3, q3b3, x)
    print string ("  q ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 0.25
    y = evalbern3 (r0b3, r1b3, r2b3, r3b3, x)
    print string ("  r ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline
    x = 7.50
    y = evalbern3 (r0b3, r1b3, r2b3, r3b3, x)
    print string ("  r ( "); print fl (x, 4)
    print string (") = "); print fl (y, 4)
    newline

    bern2to3 (p0b2, p1b2, p2b2, pc0, pc1, pc2, pc3)
    bern2to3 (q0b2, q1b2, q2b2, qc0, qc1, qc2, qc3)
    print string ("Subprogram (5) examples:".snl)
    print string ("  bern ( ")
    print fl (p0b2, 4); print string (", ")
    print fl (p1b2, 4); print string (", ")
    print fl (p2b2, 4); print string (") --> bern ( ")
    print fl (pc0, 4); print string (", ")
    print fl (pc1, 4); print string (", ")
    print fl (pc2, 4); print string (", ")
    print fl (pc3, 4); print string (")".snl)
    print string ("  bern ( ")
    print fl (q0b2, 4); print string (", ")
    print fl (q1b2, 4); print string (", ")
    print fl (q2b2, 4); print string (") --> bern ( ")
    print fl (qc0, 4); print string (", ")
    print fl (qc1, 4); print string (", ")
    print fl (qc2, 4); print string (", ")
    print fl (qc3, 4); print string (")".snl)
endofprogram