# Curve that touches three points

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Curve that touches three points is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Draw a curve that touches 3 points (1 starting point, 2 medium, 3 final point)

1.  Do not use functions of a library, implement the curve() function yourself
2.  coordinates:(x,y) starting point (10,10) medium point (100,200) final point (200,10)

## Action!

```INCLUDE "H6:REALMATH.ACT"

TYPE Point=[INT x,y]

PROC QuadraticCurve(Point POINTER p1,p2,p3 REAL POINTER a,b,c)
REAL x1,y1,x2,y2,x3,y3,x11,x22,x33,m,n,tmp1,tmp2,tmp3,tmp4,r1

IntToRealForNeg(-1,r1)
IntToRealForNeg(p1.x,x1) IntToRealForNeg(p1.y,y1)
IntToRealForNeg(p2.x,x2) IntToRealForNeg(p2.y,y2)
IntToRealForNeg(p3.x,x3) IntToRealForNeg(p3.y,y3)

RealMult(x1,x1,x11) ;x11=x1^2
RealMult(x2,x2,x22) ;x22=x2^2
RealMult(x3,x3,x33) ;x33=x3^2

RealSub(x1,x2,m)   ;m=x1-x2
RealSub(x3,x2,n)   ;n=x3-x2
RealMult(m,n,tmp1) ;tmp1=m*n

IF IsNegative(tmp1) THEN
RealMult(m,r1,tmp1)
RealAssign(tmp1,m)  ;m=-m
FI

RealSub(y1,y2,tmp1)     ;tmp1=y1-y2
RealMult(n,tmp1,tmp2)   ;tmp2=n*(y1-y2)
RealSub(y3,y2,tmp1)     ;tmp1=y3-y2
RealMult(m,tmp1,tmp3)   ;tmp3=m*(y3-y2)
RealAdd(tmp2,tmp3,tmp1) ;tmp1=n*(y1-y2)+m*(y3-y2)

RealSub(x11,x22,tmp2)   ;tmp2=x1^2-x2^2
RealMult(n,tmp2,tmp3)   ;tmp3=n*(x1^2-x2^2)
RealSub(x33,x22,tmp2)   ;tmp2=x3^2-x2^2
RealMult(m,tmp2,tmp4)   ;tmp4=m*(x3^2-x2^2)
RealAdd(tmp3,tmp4,tmp2) ;tmp2=n*(x1^2-x2^2)+m*(x3^2-x2^2)

RealDiv(tmp1,tmp2,a)    ;a=(n*(y1-y2)+m*(y3-y2)) / (n*(x1^2-x2^2)+m*(x3^2-x2^2))

RealSub(x33,x22,tmp1)   ;tmp1=x3^2-x2^2
RealMult(tmp1,a,tmp2)   ;tmp2=(x3^2-x2^2)*a
RealSub(y3,y2,tmp1)     ;tmp1=y3-y2
RealSub(tmp1,tmp2,tmp3) ;tmp3=(y3-y2)-(x3^2-x2^2)*a

RealSub(x3,x2,tmp1)     ;tmp1=x3-x2
RealDiv(tmp3,tmp1,b)    ;b=((y3-y2)-(x3^2-x2^2)*a) / (x3-x2)

RealMult(a,x11,tmp1)    ;tmp1=a*x1^2
RealMult(b,x1,tmp2)     ;tmp2=b*x1
RealSub(y1,tmp1,tmp3)   ;tmp3=y1-a*x1^2
RealSub(tmp3,tmp2,c)    ;c=y1-a*x1^2-b*x1
RETURN

PROC DrawPoint(INT x,y)
Plot(x-2,y-2) DrawTo(x+2,y-2)
DrawTo(x+2,y+2) DrawTo(x-2,y+2)
DrawTo(x-2,y-2)
RETURN

INT FUNC Min(INT a,b)
IF a<b THEN RETURN (a) FI
RETURN (b)

INT FUNC Max(INT a,b)
IF a>b THEN RETURN (a) FI
RETURN (b)

INT FUNC CalcY(REAL POINTER a,b,c INT xi)
REAL xr,xr2,yr,tmp1,tmp2,tmp3
INT yi

IntToRealForNeg(xi,xr)  ;xr=x
RealMult(xr,xr,xr2)     ;xr2=x^2
RealMult(a,xr2,tmp1)    ;tmp1=a*x^2
RealMult(b,xr,tmp2)     ;tmp2=b*x
RealAdd(tmp1,tmp2,tmp3) ;tmp3=a*x^2+b*x
RealAdd(tmp3,c,yr)      ;y3=a*x^2+b*x+c
yi=Round(yr)
RETURN (yi)

PROC DrawCurve(Point POINTER p1,p2,p3)
REAL a,b,c
INT xi,yi,minX,maxX

QuadraticCurve(p1,p2,p3,a,b,c)

DrawPoint(p1.x,p1.y)
DrawPoint(p2.x,p2.y)
DrawPoint(p3.x,p3.y)

minX=Min(p1.x,p2.x)
minX=Min(minX,p3.x)
maxX=Max(p1.x,p2.x)
maxX=Max(maxX,p3.x)

yi=CalcY(a,b,c,minX)
Plot(minX,yi)
FOR xi=minX TO maxX
DO
yi=CalcY(a,b,c,xi)
DrawTo(xi,yi)
OD
RETURN

PROC Main()
BYTE CH=\$02FC,COLOR1=\$02C5,COLOR2=\$02C6
Point p1,p2,p3

Graphics(8+16)
Color=1
COLOR1=\$0C
COLOR2=\$02

p1.x=10 p1.y=10
p2.x=100 p2.y=180
p3.x=200 p3.y=10
DrawCurve(p1,p2,p3)

DO UNTIL CH#\$FF OD
CH=\$FF
RETURN```
Output:

## Ada

Find center and radius of circle that touches the 3 points. Solve with simple linear algebra. In this case no division by zero.

```with Ada.Text_Io;
with Ada.Numerics.Generic_Elementary_Functions;

procedure Three_Point_Circle is

type Real is new Float;

package Real_Math is
new Ada.Numerics.Generic_Elementary_Functions (Real);

package Real_Io is
new Ada.Text_Io.Float_Io (Real);

use Real_Io, Ada.Text_Io;

-- Point P1
X1 : constant Real := 10.0;
Y1 : constant Real := 10.0;

-- Point P2
X2 : constant Real := 100.0;
Y2 : constant Real := 200.0;

-- Point P3
X3 : constant Real := 200.0;
Y3 : constant Real :=  10.0;

-- Point P4 - midpoint between P1 and P2
X4 : constant Real := (X1 + X2) / 2.0;
Y4 : constant Real := (Y1 + Y2) / 2.0;
S4 : constant Real := (Y2 - Y1) / (X2 - X1); -- Slope P1-P2
A4 : constant Real := -1.0 / S4;             -- Slope P4-Center
-- Y4 = A4 * X4 + B4  <=>  B4 = Y4 - A4 * X4
B4 : constant Real := Y4 - A4 * X4;

-- Point P5 - midpoint between P2 and P3
X5 : constant Real := (X2 + X3) / 2.0;
Y5 : constant Real := (Y2 + Y3) / 2.0;
S5 : constant Real := (Y3 - Y2) / (X3 - X2); -- Slope P2-P3
A5 : constant Real := -1.0 / S5;             -- Slope P5-Center
-- Y5 = A5 * X5 + B5  <=>  B5 = Y5 - A5 * X5
B5 : constant Real := Y5 - A5 * X5;

-- Find center
-- Y = A4 * X + B4     -- Line 1
-- Y = A5 * X + B5     -- Line 2
-- Solve for X:
-- A4 * X + B4 = A5 * X + B5
-- A4 * X - A5 * X = B5 - B4
-- X * (A4 - A5) = B5 - B4
-- X = (B5 - B4) / (A4 - A5)
Xc : constant Real := (B5 - B4) / (A4 - A5);
Yc : constant Real := A4 * Xc + B4;
-- Radius
R  : constant Real := Real_Math.Sqrt ((X1 - Xc) ** 2 + (Y1 - Yc) ** 2);
begin
Real_Io.Default_Exp := 0;
Real_Io.Default_Aft := 1;

Put ("Center : "); Put ("("); Put (Xc);  Put (", ");  Put (Yc);  Put (")"); New_Line;
Put ("Radius : "); Put (R);  New_Line;
end Three_Point_Circle;
```
Output:
```Center : (105.0, 81.3)
Radius : 118.8
```

## ALGOL 68

Translation of: Ada

As with the Ada and FreeBASIC translation of it, finds the centre and radius of a circle through the points.

```BEGIN # draw a curve that passes through 3 points - translated from the Ada  #
# sample - finds the centre and radius of a circle through the points  #

REAL x1 = 10;                    # Point P1                              #
REAL y1 = 10;

REAL x2 = 100                    # Point P2                              #;
REAL y2 = 200;

REAL x3 = 200;                   # Point P3                              #
REAL y3 =  10;

REAL x4 = ( x1 + x2 ) / 2;       # Point P4 - midpoint between P1 and P2 #
REAL y4 = ( y1 + y2 ) / 2;
REAL s4 = ( y2 - y1 ) / ( x2 - x1 );                   # Slope P1-P2     #
REAL a4 = -1 / s4;                                     # Slope P4-Centre #
# Y4 = A4 * X4 + B4  <=>  B4 = Y4 - A4 * X4                              #
REAL b4 = y4 - a4 * x4;

REAL x5 = ( x2 + x3 ) / 2;       # Point P5 - midpoint between P2 and P3 #
REAL y5 = ( y2 + y3 ) / 2;
REAL s5 = ( y3 - y2 ) / ( x3 - x2 );                   # Slope P2-P3     #
REAL a5 = -1.0 / s5;                                   # Slope P5-Centre #
# Y5 = A5 * X5 + B5  <=>  B5 = Y5 - A5 * X5                              #
REAL b5 = y5 - a5 * x5;

# Find centre                                                            #
# Y = A4 * X + B4     -- Line 1                                          #
# Y = A5 * X + B5     -- Line 2                                          #
# Solve for X:                                                           #
# A4 * X + B4 = A5 * X + B5                                              #
# A4 * X - A5 * X = B5 - B4                                              #
# X * (A4 - A5) = B5 - B4                                                #
# X = (B5 - B4) / (A4 - A5)                                              #
REAL xc = ( b5 - b4 ) / ( a4 - a5 );
REAL yc = a4 * xc + b4;
# Radius                                                                 #
REAL r = sqrt( ( x1 - xc ) ^ 2 + ( y1 - yc ) ^ 2 );

print( ( "Centre : (", fixed( xc, -6, 1 ), ", ", fixed( yc, -6, 1 ), " )", newline ) );
print( ( "Radius : ",  fixed( r, -6, 1 ), newline ) )

END```
Output:
```Centre : ( 105.0,   81.3 )
Radius :  118.8
```

## AutoHotkey

```QuadraticCurve(p1,p2,p3){ ; Y = aX^2 + bX + c
x1:=p1.1, y1:=p1.2, x2:=p2.1, y2:=p2.2, x3:=p3.1, y3:=p3.2
m:=x1-x2, n:=x3-x2, m:= ((m*n)<0?-1:1) * m
a:=(n*(y1-y2)+m*(y3-y2)) / (n*(x1**2 - x2**2) + m*(x3**2 - x2**2))
b:=((y3-y2) - (x3**2 - x2**2)*a) / (x3-x2)
c:=y1 - a*x1**2 - b*x1
return [a,b,c]
}
```

Examples:

```P1 := [10,10], P2 := [100,200], P3 := [200,10]
v := QuadraticCurve(p1,p2,p3)
a := v.1, b:= v.2, c:= v.3
for i, X in [10,100,200]{
Y := a*X**2 + b*X + c	; Y = aX^2 + bX + c
res .= "[" x ", " y "]`n"
}
MsgBox % "Y = " a " X^2 " (b>0?"+":"") b " X " (c>0?"+":"") c " `n" res
```
for plotting, use code from RosettaCode: Plot Coordinate Pairs

Outputs:

```Y = -0.021111 X^2 +4.433333 X -32.222222
[10, 10.000000]
[100, 200.000000]
[200, 10.000000]```

## F#

This task uses Lagrange_Interpolation#F#

```// Curve that touches three points. Nigel Galloway: September 13th., 2023
open Plotly.NET
let points=let a=LIF([10;100;200],[10;200;10]).Expression in [10.0..200.0]|>List.map(fun n->(n,(Evaluate.evaluate (Map.ofList ["x",n]) a).RealValue))
Chart.Point(points)|>Chart.show
```
Output:

## FreeBASIC

Translation of: Ada
```' Point P1
Dim As Double X1 = 10.0
Dim As Double Y1 = 10.0

' Point P2
Dim As Double X2 = 100.0
Dim As Double Y2 = 200.0

' Point P3
Dim As Double X3 = 200.0
Dim As Double Y3 =  10.0

' Point P4 - midpoint between P1 and P2
Dim As Double X4 = (X1 + X2) / 2.0
Dim As Double Y4 = (Y1 + Y2) / 2.0
Dim As Double S4 = (Y2 - Y1) / (X2 - X1) ' Slope P1-P2
Dim As Double A4 = -1.0 / S4             ' Slope P4-Center
' Y4 = A4 * X4 + B4  <=>  B4 = Y4 - A4 * X4
Dim As Double B4 = Y4 - A4 * X4

' Point P5 - midpoint between P2 and P3
Dim As Double X5 = (X2 + X3) / 2.0
Dim As Double Y5 = (Y2 + Y3) / 2.0
Dim As Double S5 = (Y3 - Y2) / (X3 - X2) ' Slope P2-P3
Dim As Double A5 = -1.0 / S5             ' Slope P5-Center
' Y5 = A5 * X5 + B5  <=>  B5 = Y5 - A5 * X5
Dim As Double B5 = Y5 - A5 * X5

' Find center
' Y = A4 * X + B4     ' Line 1
' Y = A5 * X + B5     ' Line 2
' Solve for X:
' A4 * X + B4 = A5 * X + B5
' A4 * X - A5 * X = B5 - B4
' X * (A4 - A5) = B5 - B4
' X = (B5 - B4) / (A4 - A5)
Dim As Double Xc = (B5 - B4) / (A4 - A5)
Dim As Double Yc = A4 * Xc + B4
' Radius
Dim As Double R  = Sqr((X1 - Xc) ^ 2 + (Y1 - Yc) ^ 2)

Print Using "Center : (###.#, ###.#)"; Xc; Yc
Print Using "Radius : ###.#"; R

Sleep
```

## Go

Library: Go Graphics

There are, of course, an infinity of curves which can be fitted to 3 points. The most obvious solution is to fit a quadratic curve (using Lagrange interpolation) and so that's what we do here.

As we're not allowed to use library functions to draw the curve, we instead divide the x-axis of the curve between successive points into equal segments and then join the resulting points with straight lines.

The resulting 'curve' is then saved to a .png file where it can be viewed with a utility such as EOG.

```package main

import "github.com/fogleman/gg"

var p = [3]gg.Point{{10, 10}, {100, 200}, {200, 10}}

func lagrange(x float64) float64 {
return (x-p[1].X)*(x-p[2].X)/(p[0].X-p[1].X)/(p[0].X-p[2].X)*p[0].Y +
(x-p[0].X)*(x-p[2].X)/(p[1].X-p[0].X)/(p[1].X-p[2].X)*p[1].Y +
(x-p[0].X)*(x-p[1].X)/(p[2].X-p[0].X)/(p[2].X-p[1].X)*p[2].Y
}

func getPoints(n int) []gg.Point {
pts := make([]gg.Point, 2*n+1)
dx := (p[1].X - p[0].X) / float64(n)
for i := 0; i < n; i++ {
x := p[0].X + dx*float64(i)
pts[i] = gg.Point{x, lagrange(x)}
}
dx = (p[2].X - p[1].X) / float64(n)
for i := n; i < 2*n+1; i++ {
x := p[1].X + dx*float64(i-n)
pts[i] = gg.Point{x, lagrange(x)}
}
return pts
}

func main() {
const n = 50 // more than enough for this
dc := gg.NewContext(210, 210)
dc.SetRGB(1, 1, 1) // White background
dc.Clear()
for _, pt := range getPoints(n) {
dc.LineTo(pt.X, pt.Y)
}
dc.SetRGB(0, 0, 0) // Black curve
dc.SetLineWidth(1)
dc.Stroke()
dc.SavePNG("quadratic_curve.png")
}
```

## J

```   NB. coordinates:(x,y) starting point (10,10) medium point (100,200) final point (200,10)

X=: 10 100 200
Y=: 10 200 10

NB. matrix division computes polynomial coefficients
NB. %. implements singular value decomposition
NB. in other words, we can also get best fit polynomials of lower order.

polynomial=: (Y %. (^/ ([: i. #)) X)&p.

assert 10 200 10 -: polynomial X  NB. test

Filter=: (#~`)(`:6)

Round=: adverb def '<.@:(1r2&+)&.:(%&m)'
assert 100 120 -: 100 8 Round 123  NB. test, round 123 to nearest multiple of 100 and of 8

NB. libraries not permitted, character cell graphics are used.

GRAPH=: 50 50 \$ ' '  NB. is an array of spaces

NB. place the axes
GRAPH=: '-' [`(([:<0; i.@:#)@:])`]} GRAPH
GRAPH=: '|' [`(([:<0;~i.@:#)@:])`]} GRAPH
GRAPH=: '+' [`((<0;0)"_)`]} GRAPH           NB. origin

NB. clip the domain.
EXES=: ((<:&(>./X) *. (<./X)&<:))Filter 5 * i. 200
WHYS=: polynomial EXES

NB. draw the curve
1j1 #"1 |. 'X' [`((<"1 WHYS ;&>&:([: 1 Round %&5) EXES)"_)`]} GRAPH

NB. were we to use a library:
load'plot'
'title 3 point fit' plot (j. polynomial) i.201
```

## Julia

To make things more specific, the example below finds the circle determined by the points. The curve is then the arc between the 3 points.

```using Plots

struct Point; x::Float64; y::Float64; end
# Find a circle passing through the 3 points
const p1 = Point(10, 10)
const p2 = Point(100, 200)
const p3 = Point(200, 10)
const allp = [p1, p2, p3]

# set up problem matrix and solve.
# if (x - a)^2 + (y - b)^2 = r^2 then for some D, E, F, x^2 + y^2 + Dx + Ey + F = 0
# therefore Dx + Ey + F = -x^2 - y^2
v = zeros(Int, 3)
m = zeros(Int, 3, 3)
for row in 1:3
m[row, 1:3] .= [allp[row].x, allp[row].y, 1]
v[row] = -(allp[row].x)^2 - (allp[row].y)^2
end
q = (m \ v)  # [-210.0, -162.632, 3526.32]
a, b, r = -q[1] / 2, -q[2] / 2, sqrt((q[1]^2/4) + q[2]^2/4 - q[3])

println("The circle with center at x = \$a, y = \$b and radius \$r.")

x = a-r:0.25:a+r
y0 = sqrt.(r^2 .- (x .- a).^2)
plt = plot(x, y0 .+ b, color = :red)
plot!(x, b .- y0, color = :red)
scatter!([p.x for p in allp], [p.y for p in allp],  markersize = r / 10)
```
Output:
```The circle with center at x = 105.0, y = 81.31578947368422 and radius 118.78948534384199.
```

## Lambdatalk

We find a curve interpolating three points using a bezier algorithm. A bezier curve built on 3 points, p0, p1, p2 doesn't interpolate p1. We compute a new point q symetric of the middle of p0, p2 with respect to p1. The curve built on p0, q, p2 interpolates p0, p1, p2.

### bezier interpolation of 3 points

```p(t) = 1*p0(1-t)2 + 2*p1(1-t)t + 1*p2t2

{def interpol
{lambda {:p0 :p1 :p2 :t :u}       // u =1-t
{+ {* 1 {A.get 0 :p0} :u :u}
{* 2 {A.get 0 :p1} :u :t}
{* 1 {A.get 0 :p2} :t :t}}
{+ {* 1 {A.get 1 :p0} :u :u}
{* 2 {A.get 1 :p1} :u :t}
{* 1 {A.get 1 :p2} :t :t}} }}
-> interpol
```

### two useful functions

```{def middle
{lambda {:p1 :p2}     // compute the middle point of p1 and p2
{A.new
{/ {+ {A.get 0 :p1} {A.get 0 :p2}} 2}
{/ {+ {A.get 1 :p1} {A.get 1 :p2}} 2} }}}
-> middle

{def symetric        // compute the symmetric point of p1 with respect to p2
{lambda {:p1 :p2}
{A.new
{- {* 2 {A.get 0 :p2}} {A.get 0 :p1} }
{- {* 2 {A.get 1 :p2}} {A.get 1 :p1} } }}}
-> symetric
```

### computing the curve

```{def curve
{lambda {:pol :n}
{S.map {{lambda {:p0 :p1 :p2 :n :i}
{interpol :p0 :p1 :p2 {/ :i :n} {- 1 {/ :i :n}}}
} {A.get 0 :pol} {A.get 1 :pol} {A.get 2 :pol} :n}
{S.serie -1 {+ :n 1}} }}}
-> curve
```

### drawing a point

```{def dot
{lambda {:pt}
{circle
{@ cx="{A.get 0 :pt}" cy="{A.get 1 :pt}" r="5"
stroke="#0ff" fill="transparent" stroke-width="2"}}}}
-> dot
```

### defining points

```{def P0 {A.new 150 180}}         -> P0
{def P1 {A.new 300 250}}         -> P1
{def P2 {A.new 150 330}}         -> P2

{def P02 {middle {P0} {P2}}}     -> P02
{def P20 {symetric {P02} {P1}}}  -> P20

{def P10 {middle {P1} {P0}}}     -> P10
{def P01 {symetric {P10} {P2}}}  -> P01

{def P21 {middle {P2} {P1}}}     -> P21
{def P12 {symetric {P21} {P0}}}  -> P12
```

### drawing points and curves

```{svg {@ width="500" height="500" style="background:#444;"}
{polyline {@ points="{curve {A.new {P0} {P20} {P2}} 20}"
stroke="#f00" fill="transparent" stroke-width="4"}}
{polyline {@ points="{curve {A.new {P1} {P01} {P0}} 20}"
stroke="#0f0" fill="transparent" stroke-width="4"}}
{polyline {@ points="{curve {A.new {P2} {P12} {P1}} 20}"
stroke="#00f" fill="transparent" stroke-width="4"}}

{dot {P0}} {dot {P1}} {dot {P2}}

{dot {P02}} {dot {P20}}
{dot {P10}} {dot {P01}}
{dot {P21}} {dot {P12}}
}
```

See the result in http://lambdaway.free.fr/lambdawalks/?view=bezier_3

## Mathematica /Wolfram Language

### Built-in

```pts = {{10, 10}, {100, 200}, {200, 10}};
cs = Circumsphere[pts]
Graphics[{PointSize[Large], Point[pts], cs}]
```
Output:

Outputs the circle:

`Sphere[{105, 1545/19}, (5 Sqrt[203762])/19]`

and a graphical representation of the input points and the circle.

### Alternate implementation

```pts = {{10, 10}, {100, 200}, {200, 10}};
createCircle[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] :=
With[{a = Det[({{x1, y1, 1}, {x2, y2, 1}, {x3, y3, 1}})],
d = -Det[({{x1^2 + y1^2, y1, 1}, {x2^2 + y2^2, y2,
1}, {x3^2 + y3^2, y3, 1}})],
e = Det[({{x1^2 + y1^2, x1, 1}, {x2^2 + y2^2, x2, 1}, {x3^2 + y3^2,
x3, 1}})],
f = -Det[({{x1^2 + y1^2, x1, y1}, {x2^2 + y2^2, x2,
y2}, {x3^2 + y3^2, x3, y3}})]},
Circle[{-(d/(2 a)), -(e/(2 a))}, Sqrt[(d^2 + e^2)/(4 a^2) - f/a]]]
cs = createCircle[pts]
Graphics[{PointSize[Large], Point[pts], cs}]
```
Output:

Outputs the circle:

`Circle[{105, 1545/19}, (5 Sqrt[203762])/19]`

and a graphical representation of the input points and the circle.

## Nim

Translation of: Go
Library: imageman
```import imageman

type
FPoint = tuple[x, y: float]
FPoints3 = array[3, FPoint]

func lagrange(p: FPoints3; x: float): float =
(x-p[1].x) * (x-p[2].x) / (p[0].x-p[1].x) / (p[0].x-p[2].x) * p[0].y +
(x-p[0].x) * (x-p[2].x) / (p[1].x-p[0].x) / (p[1].x-p[2].x) * p[1].y +
(x-p[0].x) * (x-p[1].x) / (p[2].x-p[0].x) / (p[2].x-p[1].x) * p[2].y

func points(p: FPoints3; n: int): seq[Point] =
result.setLen(2 * n + 1)
var dx = (p[1].x - p[0].x) / float(n)
for i in 0..<n:
let x = p[0].x + dx * float(i)
result[i] = (x.toInt, p.lagrange(x).toInt)
dx = (p[2].x - p[1].x) / float(n)
for i in n..2*n:
let x = p[1].x + dx * float(i - n)
result[i] = (x.toInt, p.lagrange(x).toInt)

const N = 50

const P: FPoints3 =[(10.0, 10.0), (100.0, 200.0), (200.0, 10.0)]

var img = initImage[ColorRGBF](210, 210)
img.fill(ColorRGBF([float32 1, 1, 1]))    # White background.
let color = ColorRGBF([float32 0, 0, 0])  # Black.
img.drawPolyline(closed = false, color, P.points(N))
img.savePNG("curve.png", compression = 9)
```

## ooRexx

### Version 1

```/* REXX ***************************************************************
* Compute the polynome satisfying  three given Points
**********************************************************************/
pl='(10,10) (100,200) (200,10)'
Do i=1 To 3
Parse Var pl '(' x.i ',' y.i ')' pl
s.i=x.i**2 x.i 1 y.i
End
Parse Value lingl() With a b c
If a<>0 Then
gl=a'*x**2'
Else
gl=''
If b>0 & gl<>'' Then b='+'||b
If b<>0 Then gl=gl||b'*x'
If c>0 & gl<>'' Then c='+'||c
If c<>0 Then gl=gl||c
Say 'y='gl
Say 'x / f(x) / y'
Do i=1 To 3
Say x.i '/' fun(x.i) '/' y.i
End
Exit

fun:
Parse Arg x
Return a*x**2+b*x+c

lingl: Procedure  Expose s.
/************************************************* Version: 25.11.1996 *
* Lösung eines linearen Gleichungssystems
* 22.11.1996 PA neu
***********************************************************************/
Numeric Digits 12
Do i=1 to 3
l=s.i
Do j=1 By 1 While l<>''
Parse Var l a.1.i.j l
End
m=j-1
End
n=i-1
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.1.i.j,6,2)
End
Call dbg s
End
Do ie=1 To i-1
u=ie
v=ie+1
Do kk=ie To n
If a.u.kk.ie<>0 Then Leave
End
Select
When kk=ie Then Nop
When kk>n Then Call ex 'eine Katastrophe'
Otherwise Do
Do jj=1 To m
temp=a.u.ie.jj
a.u.ie.jj=a.u.kk.jj
a.u.kk.jj=temp
End
Do ip=1 To n
s=''
Do jp=1 To m
s=s format(a.u.ip.jp,12,2)
End
Call dbg s
End
End
End

Do i=1 To n
Do j=1 To m
If i<=ie Then
a.v.i.j=a.u.i.j
Else
a.v.i.j=a.u.i.j*a.u.ie.ie-a.u.i.ie*a.u.ie.j
End
End

Call dbg copies('-',70)
Do i=1 To n
Do j=1 To m
If a.v.i.j<>0 Then Leave
End
Select
When j=m Then Call ex 'Widersprü�chliches Gleichungssystem'
When j>m Then Call ex 'Gleichungen sind linear abhängig'
Otherwise Nop
End
End
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.v.i.j,12,2)
End
Call dbg s
End
End
n1=n+1
Do i=n To 1 By -1
i1=i+1
x.i=a.v.i.n1/a.v.i.i
sub=0
Do j=i+1 To n
sub=sub+a.v.i.j*x.j
End
x.i=x.i-sub/a.v.i.i
End

sol=''
Do i=1 To n
sol=sol x.i
End
Return sol

ex:
Say arg(1)
Exit

dbg: Return
```
Output:
```y=-0.021111111111*x**2+4.43333333333*x-32.2222222222
x / f(x) / y
10 / 10.0000000 / 10
100 / 200.000000 / 200
200 / 10.0000008 / 10```

### Version 2 using fraction arithmetic

```/* REXX ***************************************************************
* Compute the polynome satisfying  three given Points
**********************************************************************/
Numeric Digits 20
pl='(10,10) (100,200) (200,10)'
Do i=1 To 3
Parse Var pl '(' x.i ',' y.i ')' pl
s.i=x.i**2 x.i 1 y.i
End
abc=lingl()
a=abc[1]
b=abc[2]
c=abc[3]
If a~numerator<>0 Then
gl=a'*x**2'
Else
gl=''
If b~numerator<>0 Then gl=gl'+'||b'*x'
If c~numerator<>0 Then gl=gl'+'||c
o='y='gl
o=replr(o,'-(','+(-')
o=replr(o,'=-(','=(-')
o=replr(o,'=','=+')
Say o
Say 'x / f(x) / y'
Do i=1 To 3
Say x.i '/' fun(x.i) '/' y.i
End
Exit

fun:
Parse Arg x
Return a*x**2+b*x+c

lingl: Procedure  Expose s.
/************************************************* Version: 25.11.1996 *
* Lösung eines linearen Gleichungssystems
* 22.11.1996 PA neu
***********************************************************************/
Numeric Digits 20
Do i=1 to 3
l=s.i
Do j=1 By 1 While l<>''
Parse Var l a.1.i.j l
fa.1.i.j=.fraction~new(a.1.i.j,1)
End
m=j-1
End
n=i-1
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.1.i.j,20)
End
Call dbg s
End
Do ie=1 To i-1
u=ie
v=ie+1
Do kk=ie To n
If a.u.kk.ie<>0 Then Leave
End
Select
When kk=ie Then Nop
When kk>n Then Call ex 'eine Katastrophe'
Otherwise Do
Do jj=1 To m
temp=a.u.ie.jj
a.u.ie.jj=a.u.kk.jj
a.u.kk.jj=temp
ftemp=fa.u.ie.jj
fa.u.ie.jj=fa.u.kk.jj
fa.u.kk.jj=ftemp
End
Do ip=1 To n
s=''
Do jp=1 To m
s=s format(a.u.ip.jp,20)
End
Call dbg s
End
End
End

Do i=1 To n
Do j=1 To m
If i<=ie Then Do
a.v.i.j=a.u.i.j
fa.v.i.j=fa.u.i.j
End
Else Do
a.v.i.j=a.u.i.j*a.u.ie.ie-a.u.i.ie*a.u.ie.j
fa.v.i.j=fa.u.i.j*fa.u.ie.ie-fa.u.i.ie*fa.u.ie.j
End
End
End

Call dbg copies('-',70)
Do i=1 To n
Do j=1 To m
If a.v.i.j<>0 Then Leave
End
Select
When j=m Then Call ex 'Widersprü�chliches Gleichungssystem'
When j>m Then Call ex 'Gleichungen sind linear abhängig'
Otherwise Nop
End
End
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.v.i.j,20)
End
Call dbg s
End
End
n1=n+1
Do i=n To 1 By -1
x.i=a.v.i.n1/a.v.i.i
fx.i=fa.v.i.n1/fa.v.i.i
sub=0
fsub=.fraction~new(0,1)
Do j=i+1 To n
sub=sub+a.v.i.j*x.j
fsub=fsub+fa.v.i.j*fx.j
End
x.i=x.i-sub/a.v.i.i
fx.i=fx.i-fsub/fa.v.i.i
End

Return .array~of(fx.1,fx.2,fx.3)

ex:
Say arg(1)
Exit

dbg: Return
--REQUIRES fraction.cls

::class fraction public inherit stringlike orderable comparable

::method init                                 /* initialize a fraction          */
expose numerator denominator                /* expose the state data          */
Numeric Digits 20
Use Strict Arg numerator = 0, denominator = 1 /* access the two numbers       */
numerator += 0                              /* force rounding                 */
denominator += 0

anum=abs(numerator)
aden=abs(denominator)
x=gcd2(anum,aden)
anum=anum/x
aden=aden/x
If sign(denominator)<>sign(numerator) Then
numerator=-anum
Else
numerator=anum
denominator=aden

::method '[]' class                           /* create a new fraction          */
forward message("NEW")                      /* just a synonym for NEW         */

-- read-only attributes for numerator and denominator
::attribute numerator GET
::attribute denominator GET

::method '+'                                  /* addition method                */
expose numerator denominator                /* access the state values        */
Numeric Digits 20
Use Strict Arg adder = .nil                 /* get the operand                */

if arg(1,'o') Then                          /* prefix plus operation?         */
Return self                               /* don't do anything with this    */

if adder~isa(.string) Then                  /* if just a simple number,       */
adder = self~class~new(adder)             /* convert to a fraction          */

rnum=self~numerator*adder~denominator+,
self~denominator*adder~numerator
rdenom=self~denominator*adder~denominator

Return self~class~new(rnum,rdenom)

::method '-'                                  /* subtraction method             */
expose numerator denominator                /* access the state values        */
Numeric Digits 20
Use Strict Arg adder = .nil                 /* get the operand                */

if arg(1,'o') Then do                       /* prefix minus operation?        */
rdenom=self~denominator
rnum=-self~numerator
End
Else Do
if adder~isa(.string) Then                /* if just a simple number,       */
adder = self~class~new(adder)           /* convert to a fraction          */

rnum=self~numerator*adder~denominator-,
self~denominator*adder~numerator
rdenom=self~denominator*adder~denominator
End

Return self~class~new(rnum,rdenom)

::method '*'                                  /* multiplication method          */
expose numerator denominator                /* access the state values        */
Numeric Digits 20
Use Strict Arg adder = .nil                 /* get the operand                */

if adder~isa(.string) Then                  /* if just a simple number,       */
adder = self~class~new(adder)             /* convert to a fraction          */

rnum=self~numerator*adder~numerator
rdenom=self~denominator*adder~denominator

Return self~class~new(rnum,rdenom)

::method '/'                                  /* division method                */
expose numerator denominator                /* access the state values        */
Numeric Digits 20
Use Strict Arg adder = .nil                 /* get the operand                */

if adder~isa(.string) Then                  /* if just a simple number,       */
adder = self~class~new(adder)             /* convert toa fraction           */

rnum=self~numerator*adder~denominator
rdenom=self~denominator*adder~numerator

Return self~class~new(rnum,rdenom)

::method 'value'                              /* the fraction' numeric Value    */
expose numerator denominator                /* access the state values        */
Return numerator/denominator

::method string                               /* format as a string value       */
If self~denominator=1 Then
Return '('self~numerator')'
Else
Return '('self~numerator'/'self~denominator')' /* format as '(a,b)'         */

::class "Stringlike" PUBLIC MIXINCLASS object

-- This unknown method forwards all method invocations to the object's string value,
-- effectively adding all of the string methods to the class
::method unknown UNGUARDED                    /* create an unknown method       */
Use Arg msgname, args                       /* get the message and arguments  */
/* just forward to the string val.*/
forward to(self~string) message(msgname) arguments(args)

::ROUTINE gcd2
/**********************************************************************
* Compute greatest common divider
**********************************************************************/
Numeric Digits 20
Parse Arg a,b
if b = 0 Then Return abs(a)
Return GCD2(b,a//b)
::ROUTINE replr
/* REXX ***************************************************************
* Replace,in s, occurrences of old by new and return the changed string
* ooRexx has the builtin function changestr
**********************************************************************/
Parse Arg s,new,old
Do i=1 To 2 Until p=0
p=pos(old,s)
If p>0 Then
s=left(s,p-1)||new||substr(s,p+length(old))
End
Return s
```
Output:
```y=-(19/900)*x**2+(133/30)*x-(290/9)
x / f(x) / y
10 / (10) / 10
100 / (200) / 200
```

### Version 3 computing the circumcircle (among many other things)

```/* REXX ****************************************************************
* Triangle computes data about a given triangle
* The circumcircle is what we need here
***********************************************************************/
call triangle 10 10 200 10 100 200
Exit
triangle:
/***********************************************************************
* Triangle Computations
* 940810 PA  new
* 220624 a mere 38 years later completed and anglisized
***********************************************************************/
Parse Arg ax ay bx by cx cy
If ax='?' Then Do
Say 'REXX Triangle ax ay bx by cx cy'
Say ' computes miscellaneous data about this triangle'
Exit
End
If ax='' Then Do
d='D 0 0 10 0 5 10'
Parse Var d . ax ay bx by cx cy .
End
Else
d='D' ax ay bx by cx cy .
Say ''
Say 'Triangle ABC:'
A='P' ax ay  ; Say 'A' rep(A)
B='P' bx by  ; Say 'B' rep(B)
C='P' cx cy  ; Say 'C' rep(C)

areal=a(. ax ay bx by cx cy)
If areal<1e-3 Then
Call ex 'This isn''t a Triangle!! area='areal
Say ''
Say 'Triangle''s sides:'
al=dist(B,C) ; Say 'B-C a='round(al)
bl=dist(C,A) ; Say 'C-A b='round(bl)
cl=dist(A,B) ; Say 'A-B c='round(cl)

/* c**2=a**2+b**2-2*a*b*cos(gamma) */
cnvf=180/rxcalcpi() -- 57.2957796
alpha=rxCalcarccos((bl**2+cl**2-al**2)/(2*bl*cl),,'R')*cnvf
beta =rxCalcarccos((al**2+cl**2-bl**2)/(2*al*cl),,'R')*cnvf
gamma=rxCalcarccos((al**2+bl**2-cl**2)/(2*al*bl),,'R')*cnvf
Say ''
Say 'Triangle''s angles:'
Say 'alpha='round(alpha)
Say 'beta ='round(beta)
Say 'gamma='round(gamma)
Say 'sum  ='round(alpha+beta+gamma)

Say ''
Say 'Angle-bisectors:'
wsa=ws(A,C,B); Say 'wsA' left(rep(wsA),20)
wsb=ws(B,A,C); Say 'wsB' left(rep(wsB),20)
wsc=ws(C,A,B); Say 'wsC' left(rep(wsC),20)

ha=normale(A,g(B,C))
Call dbg  'HA' rep(ha) ha
hb=normale(B,g(A,C))
Call dbg  'HB' rep(hb) hb
hc=normale(C,g(B,A))
Call dbg  'Hc' rep(hc) hc
HSP=sp(ha,hc)
If HSP='?' Then
HSP=sp(ha,hb)
Say ''
Say 'Orthocenter:' rep(HSP)

/***********************************************************************
* Perimeter and Area
***********************************************************************/
Say ''
Say 'Perimeter:' round(u(d))
Say 'Area:     ' round(a(d))

/***********************************************************************
* Circumcircle
***********************************************************************/
U=sp(ss(A,B),ss(B,C))
Call dbg  'ss(A,B)='ss(A,B)
Call dbg  'ss(B,c)='ss(B,c)
Say ''
Say 'Center of circumcircle    :' rep(U)
Say 'Radius                    :' round(dist(U,A))

/***********************************************************************
* Inscribed circle
***********************************************************************/
I=sp(wsa,wsb)
Say ''
Say 'Center of inscribed circle:' rep(I)
Say 'Radius                    :' round(rho(d))

/***********************************************************************
* Centroid
***********************************************************************/
Call dbg  MP(B,C)
Call dbg  MP(C,A)
sa=g(A,MP(B,C)); Call dbg  'sa='sa  rep(sa)
sb=g(B,MP(C,A)); Call dbg  'sb='sb  rep(sb)
S=sp(sa,sb)
Say ''
Say 'centroid:' rep(S)

/***********************************************************************
* Feuerbach Circle
***********************************************************************/
MAB='P' (ax+bx)/2 (ay+by)/2
MBC='P' (bx+cx)/2 (by+cy)/2
MCA='P' (cx+ax)/2 (cy+ay)/2
F=sp(ss(MAB,MBC),ss(MBC,MCA))
Say ''
Say 'Center of Feuerbach Circle:' rep(F)
Say 'Radius                    :' round(dist(F,MAB))

/***********************************************************************
* Euler's Line  contains the following points:
* Centroid
* Center of circumcircle
* Orthocenter
* Center of Feuerbach Circle
***********************************************************************/
Call dbg 'Centroid..................' rep(S)
Call dbg 'Center of circumcircle....' rep(U)
Call dbg 'Orthocenter...............' rep(HSP)
Call dbg 'Center of Feuerbach Circle' rep(F)

Say ''
If abs(al-bl)<1.e-5 & abs(bl-cl)<1.e-5 Then
Say 'Equilateral Triangle - no Eulersche Gerade'
Else Do
Say 'Euler''s Line:'
su=rep(g(S,U));   Say 'S-U' su
sh=rep(g(S,HSP)); Say 'S-H' sh
sf=rep(g(S,F));   Say 'S-F' sf
uh=rep(g(U,HSP)); Say 'U-H' uh
End
Exit

round: Procedure
Numeric Digits 6
Parse Arg z
Return z+0

rep: Procedure Expose sigl
/***********************************************************************
* Show representation of a point or a line
***********************************************************************/
Parse Arg type a
Select
When type='P' Then Do
Parse Var a x y
res='('||round(x)||'/'||round(y)||')'
End
When type='g' Then Do
Parse Var a x1 y1 x2 y2
Select
When x1=x2 Then
res='x='||round(x1)
When y1=y2 Then
res='y='round(y1)
Otherwise Do
k=(y2-y1)/(x2-x1)
d=round(y1-k*x1)
Select
When d>0 Then d='+'d
When d=0 Then d=''
Otherwise Nop
End
If k=1 Then
res='y=x'd
Else
res='y='round(k)'*x'd
End
End
End
Otherwise Do
Say 'sigl='sigl
Say 'Unsupported type' type
res='???'
End
End
Return res

a: Procedure
/***********************************************************************
* Area (Heron's formula)
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
c=dist('P' ax ay,'P' bx by)
a=dist('P' bx by,'P' cx cy)
b=dist('P' cx cy,'P' ax ay)
s=(a+b+c)/2
res=rxCalcsqrt(s*(s-a)*(s-b)*(s-c))
Return res

rho: Procedure Expose ax ay bx by cx cy
/***********************************************************************
* Radius of inscribed circle
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
c=dist('P' ax ay,'P' bx by)
a=dist('P' bx by,'P' cx cy)
b=dist('P' cx cy,'P' ax ay)
s=(a+b+c)/2
res=rxCalcsqrt((s-a)*(s-b)*(s-c)/s)
Return res

u: Procedure
/***********************************************************************
* Perimeter
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
Return dist('P' ax ay,'P' bx by)+,
dist('P' bx by,'P' cx cy)+,
dist('P' cx cy,'P' ax ay)

dist: Procedure
/***********************************************************************
* Distance between two points
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return rxCalcsqrt((x2-x1)**2+(y2-y1)**2)

g: Procedure
/***********************************************************************
* Intern representation of a line though two points
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return 'g' x1 y1 (x1+(x2-x1)) (y1+(y2-y1))

MP: Procedure
/***********************************************************************
* Center of a line segment
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return 'P' ((x1+x2)/2) ((y2+y1)/2)

sp: Procedure
/***********************************************************************
* Intersection of two lines
***********************************************************************/
Parse Arg . xa ya xb yb . , . xc yc xd yd .
z=(xc-xa)*(yd-yc) - (yc-ya)*(xd-xc)
n=(xb-xa)*(yd-yc) - (yb-ya)*(xd-xc)
If n=0 Then Do
If z=0 Then
Call dbg 'lines are identical' z'/'n xa ya xb yb xc yc xd yd
Else
Call dbg 'lines are paralllel' z'/'n xa ya xb yb xc yc xd yd
Return '?'
End
Else Do
t=z/n
x=xa+(xb-xa)*t
y=ya+(yb-ya)*t
Call dbg x y
Return 'P' x y
End

euler: Procedure Expose S U HSP
/***********************************************************************
* Schwerpunkt, Umkreismittelpunkt, Höhenschnittpunkt
***********************************************************************/
Parse Arg . sx sy . ux uy . hx hy
Say 'Euler:' sx sy ux uy hx hy
eg=g(S,U);  Say rep(eg)
eg2=g(S,HSP);  Say rep(eg2)
eg3=g(U,HSP);  Say rep(eg3)
Return

dist2:Procedure
/***********************************************************************
* Distance of a point C from a line AB
***********************************************************************/
Parse Arg ax ay bx by cx cy
Say 'A('ax'/'ay')' 'B('bx'/'by')' 'C('cx'/'cy')'
gx.1=ax
gx.2=bx-ax
gy.1=ay
gy.2=by-ay

Select
When gx.2=0 & gy.2=0 Then
Call ex 'g isn''t a line'
When gx.2=0 Then Do
xf=1
yf=0
c=-ax
End
When gy.2=0 Then Do
xf=0
yf=1
c=-ay
End
Otherwise Do
xf=1/gx.2
yf=-1/gy.2
c=-((ay/gy.2)+(ax/gx.2))
End
End
call dbg xf'*x+'yf'*y-'c'=0'

d=abs((xf*cx+yf*cy-c)/rxCalcsqrt(xf**2+yf**2))
Call dbg 'd='d
Return d

normale: Procedure
/***********************************************************************
* compute the line through point C that is normal to line A-B
***********************************************************************/
Parse Arg . ax ay . , . bx by cx cy .
vx=cx-bx
vy=cy-by
res='g' ax ay ax+vy ay-vx
Call dbg  res
Return res

ss: Procedure
/***********************************************************************
* compute the perpendicular bisector of a line segment
***********************************************************************/
Parse Arg . ax ay . , . bx by .
Call dbg 'A('ax'/'ay')' 'B('bx'/'by')'
If ax=bx & ay=by Then
Call ex 'AB isn''t a line segment'
mx=(ax+bx)/2
my=(ay+by)/2
vx=bx-ax
vy=by-ay
Select
When vx=0 Then Parse Value 1 0 With sx sy
When vy=0 Then Parse Value 0 1 With sx sy
Otherwise Do
sx=vy
sy=-vx
End
End
Call dbg    'g' mx my (mx+sx) (my+sy)
Return 'g' mx my (mx+sx) (my+sy)

ws: Procedure
/***********************************************************************
* compute the angular symmetric of point A
***********************************************************************/
Parse Arg . ax ay . , . bx by . , . cx cy .
ebl=rxCalcsqrt((bx-ax)**2+(by-ay)**2)
ecl=rxCalcsqrt((cx-ax)**2+(cy-ay)**2)
--Say 'AB   ' (bx-ax)/ebl (by-ay)/ebl
--Say 'AC   ' (cx-ax)/ecl (cy-ay)/ecl
--Say 'AB+AC' ((bx-ax)/ebl+(cx-ax)/ecl) ((by-ay)/ebl+(cy-ay)/ecl)
res='g' ax ay ax+((bx-ax)/ebl+(cx-ax)/ecl)*10,
ay+((by-ay)/ebl+(cy-ay)/ecl)*10
Return res

dbg:
Return
Say      arg(1)

ex:
Say arg(1)
Exit

::requires rxMath library
```
Output:
```Triangle ABC:
A (10/10)
B (200/10)
C (100/200)

Triangle's sides:
B-C a=214.709
C-A b=210.238
A-B c=190

Triangle's angles:
alpha=64.6538
beta =62.2415
gamma=53.1047
sum  =180.000

Angle-bisectors:
wsA y=0.632831*x+3.67169
wsB y=-0.603732*x+130.74
wsC y=-47.4947*x+4949.47

Orthocenter: (100.000/57.3684)

Perimeter: 614.947
Area:      18050

Center of circumcircle    : (105/81.3158)
Radius                    : 118.789

Center of inscribed circle: (102.764/68.7042)
Radius                    : 58.7042

centroid: (103.333/73.3333)

Center of Feuerbach Circle: (102.500/69.3421)
Radius                    : 59.3947

Euler's Line:
S-U y=4.78947*x-421.579
S-H y=4.78947*x-421.579
S-F y=4.78948*x-421.579
U-H y=4.78947*x-421.579```

## Perl

Hilbert curve task code repeated here, with the addition that the 3 task-required points are marked. Satisfies the letter-of-the-law of task specification while (all in good fun) subverting the spirit of the thing.

```use SVG;
use List::Util qw(max min);

use constant pi => 2 * atan2(1, 0);

# Compute the curve with a Lindemayer-system
%rules = (
A => '-BF+AFA+FB-',
B => '+AF-BFB-FA+'
);
\$hilbert = 'A';
\$hilbert =~ s/([AB])/\$rules{\$1}/eg for 1..6;

# Draw the curve in SVG
(\$x, \$y) = (0, 0);
\$theta   = pi/2;
\$r       = 5;

for (split //, \$hilbert) {
if (/F/) {
push @X, sprintf "%.0f", \$x;
push @Y, sprintf "%.0f", \$y;
\$x += \$r * cos(\$theta);
\$y += \$r * sin(\$theta);
}
elsif (/\+/) { \$theta += pi/2; }
elsif (/\-/) { \$theta -= pi/2; }
}

\$max =  max(@X,@Y);
\$xt  = -min(@X)+10;
\$yt  = -min(@Y)+10;
\$svg = SVG->new(width=>\$max+20, height=>\$max+20);
\$points = \$svg->get_path(x=>\@X, y=>\@Y, -type=>'polyline');
\$svg->rect(width=>"100%", height=>"100%", style=>{'fill'=>'black'});
\$svg->polyline(%\$points, style=>{'stroke'=>'orange', 'stroke-width'=>1}, transform=>"translate(\$xt,\$yt)");
my \$task = \$svg->group( id => 'task-points', style => { stroke => 'red', fill => 'red' },);
\$task->circle( cx =>  10, cy =>  10, r => 1, id => 'point1' );
\$task->circle( cx => 100, cy => 200, r => 1, id => 'point2' );
\$task->circle( cx => 200, cy =>  10, r => 1, id => 'point3' );

open  \$fh, '>', 'curve-3-points.svg';
print \$fh  \$svg->xmlify(-namespace=>'svg');
close \$fh;
```

## Phix

Translation of: zkl
Library: Phix/pGUI
Library: Phix/online

You can run this online here.

```with javascript_semantics
include pGUI.e

Ihandle dlg, canvas
cdCanvas cddbuffer, cdcanvas

enum X, Y
constant p = {{10,10},{100,200},{200,10}}

function lagrange(atom x)
return (x - p[2][X])*(x - p[3][X])/(p[1][X] - p[2][X])/(p[1][X] - p[3][X])*p[1][Y] +
(x - p[1][X])*(x - p[3][X])/(p[2][X] - p[1][X])/(p[2][X] - p[3][X])*p[2][Y] +
(x - p[1][X])*(x - p[2][X])/(p[3][X] - p[1][X])/(p[3][X] - p[2][X])*p[3][Y]
end function

function getPoints(integer n)
sequence pts = {}
atom {dx,pt,cnt} := {(p[2][X] - p[1][X])/n, p[1][X], n}
for j=1 to 2 do
for i=0 to cnt do
atom x := pt + dx*i;
pts = append(pts,{x,lagrange(x)});
end for
{dx,pt,cnt} = {(p[3][X] - p[2][X])/n, p[2][X], n+1};
end for
return pts
end function

procedure draw_cross(sequence xy)
integer {x,y} = xy
cdCanvasLine(cddbuffer, x-3, y, x+3, y)
cdCanvasLine(cddbuffer, x, y-3, x, y+3)
end procedure

function redraw_cb(Ihandle /*ih*/)
cdCanvasActivate(cddbuffer)
cdCanvasSetForeground(cddbuffer, CD_BLUE)
cdCanvasBegin(cddbuffer,CD_OPEN_LINES)
atom {x,y} = {p[1][X], p[1][Y]}; -- curve starting point
cdCanvasVertex(cddbuffer, x, y)
sequence pts = getPoints(50)
for i=1 to length(pts) do
{x,y} = pts[i]
cdCanvasVertex(cddbuffer, x, y)
end for
cdCanvasEnd(cddbuffer)
cdCanvasSetForeground(cddbuffer, CD_RED)
for i=1 to length(p) do draw_cross(p[i]) end for
cdCanvasFlush(cddbuffer)
return IUP_DEFAULT
end function

function map_cb(Ihandle ih)
cdcanvas = cdCreateCanvas(CD_IUP, ih)
cddbuffer = cdCreateCanvas(CD_DBUFFER, cdcanvas)
cdCanvasSetBackground(cddbuffer, CD_WHITE)
return IUP_DEFAULT
end function

procedure main()
IupOpen()

canvas = IupCanvas(NULL)
IupSetAttribute(canvas, "RASTERSIZE", "220x220")
IupSetCallback(canvas, "MAP_CB", Icallback("map_cb"))

dlg = IupDialog(canvas,"DIALOGFRAME=YES")
IupSetAttribute(dlg, "TITLE", "Quadratic curve")
IupSetCallback(canvas, "ACTION", Icallback("redraw_cb"))

IupMap(dlg)
IupSetAttribute(canvas, "RASTERSIZE", NULL)
IupShowXY(dlg,IUP_CENTER,IUP_CENTER)
if platform()!=JS then
IupMainLoop()
IupClose()
end if
end procedure
main()
```

## PHP

Finds a circle that passes through those three points. Adapted from AutoHotkey.

```function circle_3_points(\$x1, \$y1, \$x2, \$y2, \$x3, \$y3) {
\$x4 = (\$x1 + \$x2) / 2.0;
\$y4 = (\$y1 + \$y2) / 2.0;
\$s4 = (\$y2 - \$y1) / (\$x2 - \$x1);
\$a4 = -1.0 / \$s4;
\$b4 = \$y4 - \$a4 * \$x4;
\$x5 = (\$x2 + \$x3) / 2.0;
\$y5 = (\$y2 + \$y3) / 2.0;
\$s5 = (\$y3 - \$y2) / (\$x3 - \$x2);
\$a5 = -1.0 / \$s5;
\$b5 = \$y5 - \$a5 * \$x5;
\$xc = (\$b5 - \$b4) / (\$a4 - \$a5);
\$yc = \$a4 * \$xc + \$b4;
\$r  = √(²(\$x1 - \$xc) + ²(\$y1 - \$yc));
return [\$xc, \$yc, \$r];
}
```

## Raku

(formerly Perl 6)

Kind of bogus. There are an infinite number of curves that pass through those three points. I'll assume a quadratic curve. Lots of bits and pieces borrowed from other tasks to avoid relying on library functions.

Saved as a png for wide viewing support. Note that png coordinate systems have 0,0 in the upper left corner.

```use Image::PNG::Portable;

# the points of interest
my @points = (10,10), (100,200), (200,10);

# Solve for a quadratic line that passes through those points
my (\a, \b, \c) = (rref ([.[0]², .[0], 1, .[1]] for @points) )[*;*-1];

# Evaluate quadratic equation
sub f (\x) { a×x² + b×x + c }

my (\$w, \$h) = 500, 500;  # image size
my \$scale   = 2;         # scaling factor

my \$png = Image::PNG::Portable.new: :width(\$w), :height(\$h);

my (\$lastx, \$lasty) = 8, f(8).round;
(9 .. 202).map: -> \$x {
my \$f = f(\$x).round;
line(\$lastx, \$lasty, \$x, \$f, \$png, [0,255,127]);
(\$lastx, \$lasty) = \$x, \$f;
}

# Highlight the defining points
dot( | \$_, \$(255,0,0), \$png, 2) for @points;

\$png.write: 'Curve-3-points-perl6.png';

# Assorted helper routines
sub rref (@m) {
return unless @m;
my (\$lead, \$rows, \$cols) = 0, @m, @m[0];
for ^\$rows -> \$r {
\$lead < \$cols or return @m;
my \$i = \$r;
until @m[\$i;\$lead] {
++\$i == \$rows or next;
\$i = \$r;
++\$lead == \$cols and return @m;
}
@m[\$i, \$r] = @m[\$r, \$i] if \$r != \$i;
@m[\$r] »/=» \$ = @m[\$r;\$lead];
for ^\$rows -> \$n {
next if \$n == \$r;
@m[\$n] »-=» @m[\$r] »×» (@m[\$n;\$lead] // 0);
}
++\$lead;
}
@m
}

sub line(\$x0 is copy, \$y0 is copy, \$x1 is copy, \$y1 is copy, \$png, @rgb) {
my \$steep = abs(\$y1 - \$y0) > abs(\$x1 - \$x0);
(\$x0,\$y0,\$x1,\$y1) »×=» \$scale;
if \$steep {
(\$x0, \$y0) = (\$y0, \$x0);
(\$x1, \$y1) = (\$y1, \$x1);
}
if \$x0 > \$x1 {
(\$x0, \$x1) = (\$x1, \$x0);
(\$y0, \$y1) = (\$y1, \$y0);
}
my \$Δx = \$x1 - \$x0;
my \$Δy = abs(\$y1 - \$y0);
my \$error = 0;
my \$Δerror = \$Δy / \$Δx;
my \$y-step = \$y0 < \$y1 ?? 1 !! -1;
my \$y = \$y0;
next if \$y < 0;
for \$x0 .. \$x1 -> \$x {
next if \$x < 0;
if \$steep {
\$png.set(\$y, \$x, |@rgb);
} else {
\$png.set(\$x, \$y, |@rgb);
}
\$error += \$Δerror;
if \$error ≥ 0.5 {
\$y += \$y-step;
\$error -= 1.0;
}
}
}

sub dot (\$X is copy, \$Y is copy, @rgb, \$png, \$radius = 3) {
(\$X, \$Y) »×=» \$scale;
for (\$X X+ -\$radius .. \$radius) X (\$Y X+ -\$radius .. \$radius) -> (\$x, \$y) {
\$png.set(\$x, \$y, |@rgb) if ( \$X - \$x + (\$Y - \$y) × i ).abs <= \$radius;
}
}
```

See Curve-3-points-perl6.png (offsite .png image)

## RPL

HP-49+ RPL has a built-in function named `LAGRANGE` that returns the curve equation from the 3 points, but let's consider it as belonging to a library.

Works with: HP version 48G
```« → a b c
« { 0 0 0 }
c a - C→R SWAP / c b - RE /
b a - C→R SWAP / c b - RE / -
1 SWAP PUT
b a - C→R SWAP / OVER 1 GET b a + RE * -
2 SWAP PUT
a IM OVER 1 GET a RE SQ * - OVER 2 GET a RE * -
3 SWAP PUT
{ 'X^2' 'X' 1 } * ∑LIST
» » 'PAREQ' STO

« # 131d # 64d PDIM 0 210 DUP2 XRNG YRNG
FUNCTION
(10,10) (100,200) (200,10)
3 DUPN PAREQ STEQ
ERASE DRAW
1 3 START PIXOFF NEXT  @ switch off the 3 points on the curve
{ } PVIEW
{ PPAR EQ } PURGE
» 'TASK' STO
```

## Wren

Translation of: Go
Library: DOME
```import "graphics" for Canvas, Color, Point
import "dome" for Window

class Game {
static init() {
Window.title = "Quadratic curve"
var width = 210
var height = 210
Window.resize(width, height)
Canvas.resize(width, height)
Canvas.cls(Color.white)
var n = 50
var p = [Point.new(10, 10), Point.new(100, 200), Point.new(200, 10)]
var col = Color.black // black curve
quadratic(n, p, col)
}

static update() {}

static draw(alpha) {}

static lagrange(p, x) {
return (x-p[1].x)*(x-p[2].x)/(p[0].x-p[1].x)/(p[0].x-p[2].x)*p[0].y +
(x-p[0].x)*(x-p[2].x)/(p[1].x-p[0].x)/(p[1].x-p[2].x)*p[1].y +
(x-p[0].x)*(x-p[1].x)/(p[2].x-p[0].x)/(p[2].x-p[1].x)*p[2].y
}

static quadratic(n, p, col) {
var pts = List.filled(2*n+1, null)
var dx = (p[1].x - p[0].x) / n
for (i in 0...n) {
var x = p[0].x + dx*i
pts[i] = Point.new(x, lagrange(p, x))
}
dx = (p[2].x - p[1].x) / n
for (i in n...2*n+1) {
var x = p[1].x + dx*(i-n)
pts[i] = Point.new(x, lagrange(p, x))
}
var prev = pts[0]
for (pt in pts.skip(1)) {
Canvas.line(prev.x, prev.y, pt.x, pt.y, col)
prev = pt
}
}
}
```

## XPL0

```def     X0=10., Y0=10.,  X1=100., Y1=200.,  X2=200., Y2=10.;

func real Lagrange(X);
real X;
return  (X-X1) * (X-X2) / (X0-X1) / (X0-X2) * Y0 +
(X-X0) * (X-X2) / (X1-X0) / (X1-X2) * Y1 +
(X-X0) * (X-X1) / (X2-X0) / (X2-X1) * Y2;

def  Offset = 205;
real X;
[SetVid(\$13);   \VGA 320x200x8 graphics
X:= X0;
Move(fix(X), Offset-fix(Y0));
repeat  X:= X + 1.;
Line(fix(X), Offset-fix(Lagrange(X)), 3\cyan\);
until   X >= X2;
Point(fix(X0), Offset-fix(Y0), 14\yellow\);
Point(fix(X1), Offset-fix(Y1), 14);
Point(fix(X2), Offset-fix(Y2), 14);
]```

## zkl

Translation of: Go

Uses Image Magick and the PPM class from http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#zkl

```const X=0, Y=1;   // p.X == p[X]
var p=L(L(10.0, 10.0), L(100.0, 200.0), L(200.0, 10.0));  // (x,y)

fcn lagrange(x){  // float-->float
(x - p[1][X])*(x - p[2][X])/(p[0][X] - p[1][X])/(p[0][X] - p[2][X])*p[0][Y] +
(x - p[0][X])*(x - p[2][X])/(p[1][X] - p[0][X])/(p[1][X] - p[2][X])*p[1][Y] +
(x - p[0][X])*(x - p[1][X])/(p[2][X] - p[0][X])/(p[2][X] - p[1][X])*p[2][Y]
}

fcn getPoints(n){  // int-->( (x,y) ..)
pts:=List.createLong(2*n+1);
dx,pt,cnt := (p[1][X] - p[0][X])/n, p[0][X], n;
do(2){
foreach i in (cnt){
x:=pt + dx*i;
pts.append(L(x,lagrange(x)));
}
dx,pt,cnt = (p[2][X] - p[1][X])/n, p[1][X], n+1;
}
pts
}

fcn main{
var [const] n=50; // more than enough for this
img,color := PPM(210,210,0xffffff), 0;     // white background, black curve
foreach x,y in (p){ img.cross(x.toInt(),y.toInt(), 0xff0000) } // mark 3 pts

a,b := p[0][X].toInt(), p[0][Y].toInt(); // curve starting point
foreach x,y in (getPoints(n)){
x,y = x.toInt(),y.toInt();
img.line(a,b, x,y, color);	 // can only deal with ints
a,b = x,y;
}
img.writeJPGFile("quadraticCurve.zkl.jpg");
}();```
Output:

Image at quadratic curve