Curve that touches three points

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!

Library: Action! Tool Kit

Library: Action! Real Math

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:

Screenshot from Atari 8-bit computer

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

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:

File:C3p.png

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()

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

HP-48G emulator screenshot

« → 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

File:XPL0 Curve.gif

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