```txt

"3. x^2 + 1.99999999996 x + 1.00000000006"

Fortran

module fitting
contains

  function polyfit(vx, vy, d)
    implicit none
    integer, intent(in)                   :: d
    integer, parameter                    :: dp = selected_real_kind(15, 307)
    real(dp), dimension(d+1)              :: polyfit
    real(dp), dimension(:), intent(in)    :: vx, vy

    real(dp), dimension(:,:), allocatable :: X
    real(dp), dimension(:,:), allocatable :: XT
    real(dp), dimension(:,:), allocatable :: XTX

    integer :: i, j

    integer     :: n, lda, lwork
    integer :: info
    integer, dimension(:), allocatable :: ipiv
    real(dp), dimension(:), allocatable :: work

    n = d+1
    lda = n
    lwork = n

    allocate(ipiv(n))
    allocate(work(lwork))
    allocate(XT(n, size(vx)))
    allocate(X(size(vx), n))
    allocate(XTX(n, n))

    ! prepare the matrix
    do i = 0, d
       do j = 1, size(vx)
          X(j, i+1) = vx(j)**i
       end do
    end do

    XT  = transpose(X)
    XTX = matmul(XT, X)

    ! calls to LAPACK subs DGETRF and DGETRI
    call DGETRF(n, n, XTX, lda, ipiv, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if
    call DGETRI(n, XTX, lda, ipiv, work, lwork, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if

    polyfit = matmul( matmul(XTX, XT), vy)

    deallocate(ipiv)
    deallocate(work)
    deallocate(X)
    deallocate(XT)
    deallocate(XTX)

  end function

end module

Example

program PolynomalFitting
  use fitting
  implicit none

  ! let us test it
  integer, parameter      :: degree = 2
  integer, parameter      :: dp = selected_real_kind(15, 307)
  integer                 :: i
  real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
  real(dp), dimension(11) :: y = (/ 1,   6,  17,  34, &
                                   57,  86, 121, 162, &
                                   209, 262, 321 /)
  real(dp), dimension(degree+1) :: a

  a = polyfit(x, y, degree)

  write (*, '(F9.4)') a

end program

{{out}} (lower powers first, so this seems the opposite of the Python output):

   1.0000
   2.0000
   3.0000

FreeBASIC

Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
    Dim As Integer n=Ubound(matrix,1)
    Redim ans(0):Redim ans(1 To n)
    Dim As Double b(1 To n,1 To n),r(1 To n)
    For c As Integer=1 To n 'take copies
        r(c)=rhs(c)
        For d As Integer=1 To n
            b(c,d)=matrix(c,d)
        Next d
    Next c
    #macro pivot(num)
    For p1 As Integer  = num To n - 1
        For p2 As Integer  = p1 + 1 To n
            If Abs(b(p1,num))<Abs(b(p2,num)) Then
                Swap r(p1),r(p2)
                For g As Integer=1 To n
                    Swap b(p1,g),b(p2,g)
                Next g
            End If
        Next p2
    Next p1
    #endmacro
    For k As Integer=1 To n-1
        pivot(k)              'full pivoting
        For row As Integer =k To n-1
            If b(row+1,k)=0 Then Exit For
            Var f=b(k,k)/b(row+1,k)
            r(row+1)=r(row+1)*f-r(k)
            For g As Integer=1 To n
                b((row+1),g)=b((row+1),g)*f-b(k,g)
            Next g
        Next row
    Next k
    'back substitute
    For z As Integer=n To 1 Step -1
        ans(z)=r(z)/b(z,z)
        For j As Integer = n To z+1 Step -1
            ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
        Next j
        Next    z
    End Sub

    'Interpolate through points.
    Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
        Var n=Ubound(x_values)
        Redim p(0):Redim p(1 To n)
        Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
        For a As Integer=1 To n
            rhs(a)=y_values(a)
            For b As Integer=1 To n
                matrix(a,b)=x_values(a)^(b-1)
            Next b
        Next a
        'Solve the linear equations
        GaussJordan(matrix(),rhs(),p())
    End Sub

 '
### ===================== SET UP THE POINTS ============


    Dim As Double x(1 To ...)={0,1,2,3,4,5,6,7,8,9,10}
    Dim As Double y(1 To ...)={1,6,17,34,57,86,121,162,209,262,321}

    Redim As Double Poly(0)
    'Get the polynomial Poly()
    Interpolate(x(),y(),Poly())

    'print coefficients to console
    print "Polynomial Coefficients:"
    print
    For z As Integer=1 To Ubound(Poly)
        If z=1 Then
            Print "constant term  ";tab(20);Poly(z)
        Else
            Print tab(8); "x^";z-1;" =  ";tab(20);Poly(z)
        End If
    Next z

    sleep

Polynomial Coefficients:

constant term       1
       x^ 1 =       2
       x^ 2 =       3
       x^ 3 =       0
       x^ 4 =       0
       x^ 5 =       0
       x^ 6 =       0
       x^ 7 =       0
       x^ 8 =       0
       x^ 9 =       0
       x^ 10 =      0

GAP

PolynomialRegression := function(x, y, n)
	local a;
	a := List([0 .. n], i -> List(x, s -> s^i));
	return TransposedMat((a * TransposedMat(a))^-1 * a * TransposedMat([y]))[1];
end;

x := [0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10];
y := [1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321];

# Return coefficients in ascending degree order
PolynomialRegression(x, y, 2);
# [ 1, 2, 3 ]

gnuplot

# The polynomial approximation
f(x) = a*x**2 + b*x + c

# Initial values for parameters
a = 0.1
b = 0.1
c = 0.1

# Fit f to the following data by modifying the variables a, b, c
fit f(x) '-' via a, b, c
   0   1
   1   6
   2  17
   3  34
   4  57
   5  86
   6 121
   7 162
   8 209
   9 262
  10 321
e

print sprintf("\n --- \n Polynomial fit: %.4f x^2 + %.4f x + %.4f\n", a, b, c)

Go

Library gonum/matrix

package main

import (
    "fmt"

    "github.com/gonum/matrix/mat64"
)

var (
    x = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
    y = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}

    degree = 2
)

func main() {
    a := Vandermonde(x, degree)
    b := mat64.NewDense(len(y), 1, y)
    c := mat64.NewDense(degree+1, 1, nil)

    qr := new(mat64.QR)
    qr.Factorize(a)

    err := c.SolveQR(qr, false, b)
    if err != nil {
        fmt.Println(err)
    } else {
        fmt.Printf("%.3f\n", mat64.Formatted(c))
    }
}

func Vandermonde(a []float64, degree int) *mat64.Dense {
    x := mat64.NewDense(len(a), degree+1, nil)
    for i := range a {
        for j, p := 0, 1.; j <= degree; j, p = j+1, p*a[i] {
            x.Set(i, j, p)
        }
    }
    return x
}

⎡1.000⎤
⎢2.000⎥
⎣3.000⎦

Library go.matrix

Least squares solution using QR decomposition and package [http://github.com/skelterjohn/go.matrix go.matrix].

package main

import (
    "fmt"

    "github.com/skelterjohn/go.matrix"
)

var xGiven = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
var yGiven = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
var degree = 2

func main() {
    m := len(yGiven)
    n := degree + 1
    y := matrix.MakeDenseMatrix(yGiven, m, 1)
    x := matrix.Zeros(m, n)
    for i := 0; i < m; i++ {
        ip := float64(1)
        for j := 0; j < n; j++ {
            x.Set(i, j, ip)
            ip *= xGiven[i]
        }
    }

    q, r := x.QR()
    qty, err := q.Transpose().Times(y)
    if err != nil {
        fmt.Println(err)
        return
    }
    c := make([]float64, n)
    for i := n - 1; i >= 0; i-- {
        c[i] = qty.Get(i, 0)
        for j := i + 1; j < n; j++ {
            c[i] -= c[j] * r.Get(i, j)
        }
        c[i] /= r.Get(i, i)
    }
    fmt.Println(c)
}

{{out}} (lowest order coefficient first)

[0.9999999999999758 2.000000000000015 2.999999999999999]

Haskell

Uses module Matrix.LU from [http://hackage.haskell.org/package/dsp hackageDB DSP]

import Data.List
import Data.Array
import Control.Monad
import Control.Arrow
import Matrix.LU

ppoly p x = map (x**) p

polyfit d ry = elems $ solve mat vec  where
   mat = listArray ((1,1), (d,d)) $ liftM2 concatMap ppoly id [0..fromIntegral $ pred d]
   vec = listArray (1,d) $ take d ry

{{out}} in GHCi:

*Main> polyfit 3 [1,6,17,34,57,86,121,162,209,262,321]
[1.0,2.0,3.0]

HicEst

REAL :: n=10, x(n), y(n), m=3, p(m)

   x = (0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10)
   y = (1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321)

   p = 2 ! initial guess for the polynom's coefficients

   SOLVE(NUL=Theory()-y(nr), Unknown=p, DataIdx=nr, Iters=iterations)

   WRITE(ClipBoard, Name) p, iterations

FUNCTION Theory()
   ! called by the solver of the SOLVE function. All variables are global
   Theory = p(1)*x(nr)^2 + p(2)*x(nr) + p(3)
 END

SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt):
p(1)=2.997135145; p(2)=2.011348347; p(3)=0.9906627242; iterations=19;

Hy

(import [numpy [polyfit]])

(setv x (range 11))
(setv y [1 6 17 34 57 86 121 162 209 262 321])

(print (polyfit x y 2))

J

   Y=:1 6 17 34 57 86 121 162 209 262 321
   (%. ^/~@x:@i.@#) Y
1 2 3 0 0 0 0 0 0 0 0

Note that this implementation does not use floating point numbers, so we do not introduce floating point errors. Using exact arithmetic has a speed penalty, but for small problems like this it is inconsequential.

The above solution fits a polynomial of order 11. If the order of the polynomial is known to be 3 (as is implied in the task description) then the following solution is probably preferable:

   Y %. (i.3) ^/~ i.#Y
1 2 3

(note that this time we used floating point numbers, so that result is approximate rather than exact - it only looks exact because of how J displays floating point numbers (by default, J assumes six digits of accuracy) - changing (i.3) to (x:i.3) would give us an exact result, if that mattered.)

Java

import java.util.Arrays;
import java.util.function.IntToDoubleFunction;
import java.util.stream.IntStream;

public class PolynomialRegression {
    private static void polyRegression(int[] x, int[] y) {
        int n = x.length;
        int[] r = IntStream.range(0, n).toArray();
        double xm = Arrays.stream(x).average().orElse(Double.NaN);
        double ym = Arrays.stream(y).average().orElse(Double.NaN);
        double x2m = Arrays.stream(r).map(a -> a * a).average().orElse(Double.NaN);
        double x3m = Arrays.stream(r).map(a -> a * a * a).average().orElse(Double.NaN);
        double x4m = Arrays.stream(r).map(a -> a * a * a * a).average().orElse(Double.NaN);
        double xym = 0.0;
        for (int i = 0; i < x.length && i < y.length; ++i) {
            xym += x[i] * y[i];
        }
        xym /= Math.min(x.length, y.length);
        double x2ym = 0.0;
        for (int i = 0; i < x.length && i < y.length; ++i) {
            x2ym += x[i] * x[i] * y[i];
        }
        x2ym /= Math.min(x.length, y.length);

        double sxx = x2m - xm * xm;
        double sxy = xym - xm * ym;
        double sxx2 = x3m - xm * x2m;
        double sx2x2 = x4m - x2m * x2m;
        double sx2y = x2ym - x2m * ym;

        double b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
        double c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
        double a = ym - b * xm - c * x2m;

        IntToDoubleFunction abc = (int xx) -> a + b * xx + c * xx * xx;

        System.out.println("y = " + a + " + " + b + "x + " + c + "x^2");
        System.out.println(" Input  Approximation");
        System.out.println(" x   y     y1");
        for (int i = 0; i < n; ++i) {
            System.out.printf("%2d %3d  %5.1f\n", x[i], y[i], abc.applyAsDouble(x[i]));
        }
    }

    public static void main(String[] args) {
        int[] x = IntStream.range(0, 11).toArray();
        int[] y = new int[]{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
        polyRegression(x, y);
    }
}

y = 1.0 + 2.0x + 3.0x^2
 Input  Approximation
 x   y     y1
 0   1    1.0
 1   6    6.0
 2  17   17.0
 3  34   34.0
 4  57   57.0
 5  86   86.0
 6 121  121.0
 7 162  162.0
 8 209  209.0
 9 262  262.0
10 321  321.0

Julia

The least-squares fit problem for a degree n can be solved with the built-in backslash operator (coefficients in increasing order of degree):

polyfit(x::Vector, y::Vector, deg::Int) = collect(v ^ p for v in x, p in 0:deg) \ y

x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
@show polyfit(x, y, 2)

polyfit(x, y, 2) = [1.0, 2.0, 3.0]

Kotlin

// version 1.1.51

fun polyRegression(x: IntArray, y: IntArray) {
    val n = x.size
    val r = 0 until n
    val xm = x.average()
    val ym = y.average()
    val x2m = r.map { it * it }.average()
    val x3m = r.map { it * it * it }.average()
    val x4m = r.map { it * it * it * it }.average()
    val xym = x.zip(y).map { it.first * it.second }.average()
    val x2ym = x.zip(y).map { it.first * it.first * it.second }.average()

    val sxx = x2m - xm * xm
    val sxy = xym - xm * ym
    val sxx2 = x3m - xm * x2m
    val sx2x2 = x4m - x2m * x2m
    val sx2y = x2ym - x2m * ym

    val b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    val c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    val a = ym - b * xm - c * x2m

    fun abc(xx: Int) = a + b * xx + c * xx * xx

    println("y = $a + ${b}x + ${c}x^2\n")
    println(" Input  Approximation")
    println(" x   y     y1")
    for (i in 0 until n) {
        System.out.printf("%2d %3d  %5.1f\n", x[i], y[i], abc(x[i]))
    }
}

fun main(args: Array<String>) {
    val x = IntArray(11) { it }
    val y = intArrayOf(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
    polyRegression(x, y)
}

y = 1.0 + 2.0x + 3.0x^2

 Input  Approximation
 x   y     y1
 0   1    1.0
 1   6    6.0
 2  17   17.0
 3  34   34.0
 4  57   57.0
 5  86   86.0
 6 121  121.0
 7 162  162.0
 8 209  209.0
 9 262  262.0
10 321  321.0

Lua

function eval(a,b,c,x)
    return a + (b + c * x) * x
end

function regression(xa,ya)
    local n = #xa

    local xm = 0.0
    local ym = 0.0
    local x2m = 0.0
    local x3m = 0.0
    local x4m = 0.0
    local xym = 0.0
    local x2ym = 0.0

    for i=1,n do
        xm = xm + xa[i]
        ym = ym + ya[i]
        x2m = x2m + xa[i] * xa[i]
        x3m = x3m + xa[i] * xa[i] * xa[i]
        x4m = x4m + xa[i] * xa[i] * xa[i] * xa[i]
        xym = xym + xa[i] * ya[i]
        x2ym = x2ym + xa[i] * xa[i] * ya[i]
    end
    xm = xm / n
    ym = ym / n
    x2m = x2m / n
    x3m = x3m / n
    x4m = x4m / n
    xym = xym / n
    x2ym = x2ym / n

    local sxx = x2m - xm * xm
    local sxy = xym - xm * ym
    local sxx2 = x3m - xm * x2m
    local sx2x2 = x4m - x2m * x2m
    local sx2y = x2ym - x2m * ym

    local b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    local c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    local a = ym - b * xm - c * x2m

    print("y = "..a.." + "..b.."x + "..c.."x^2")

    for i=1,n do
        print(string.format("%2d %3d  %3d", xa[i], ya[i], eval(a, b, c, xa[i])))
    end
end

local xa = {0, 1,  2,  3,  4,  5,   6,   7,   8,   9,  10}
local ya = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
regression(xa, ya)

y = 1 + 2x + 3x^2
 0   1    1
 1   6    6
 2  17   17
 3  34   34
 4  57   57
 5  86   86
 6 121  121
 7 162  162
 8 209  209
 9 262  262
10 321  321

Maple

with(CurveFitting);
PolynomialInterpolation([[0, 1], [1, 6], [2, 17], [3, 34], [4, 57], [5, 86], [6, 121], [7, 162], [8, 209], [9, 262], [10, 321]], 'x');

Result:

3*x^2+2*x+1

Mathematica

Using the built-in "Fit" function.

data = Transpose@{Range[0, 10], {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}};
Fit[data, {1, x, x^2}, x]

Second version: using built-in "InterpolatingPolynomial" function.

Simplify@InterpolatingPolynomial[{{0, 1}, {1, 6}, {2, 17}, {3, 34}, {4, 57}, {5, 86}, {6, 121}, {7, 162}, {8, 209}, {9, 262}, {10, 321}}, x]

Result:

1 + 2x + 3x^2

MATLAB

Matlab has a built-in function "polyfit(x,y,n)" which performs this task. The arguments x and y are vectors which are parametrized by the index suck that $point\_\{i\}\; =\; (x\_\{i\},y\_\{i\})$ and the argument n is the order of the polynomial you want to fit. The output of this function is the coefficients of the polynomial which best fit these x,y value pairs.

 x = [0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10];
>> y = [1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321];
>> polyfit(x,y,2)

ans =

   2.999999999999998   2.000000000000019   0.999999999999956

=={{header|MK-61/52}}== Part 1: ПC С/П ПD ИП9 + П9 ИПC ИП5 + П5 ИПC x^2 П2 ИП6 + П6 ИП2 ИПC * ИП7

• П7 ИП2 x^2 ИП8 + П8 ИПC ИПD * ИПA + ПA ИП2 ИПD * ИПB + ПB ИПD КИП4 С/П БП 00

''Input'': В/О x<sub>1</sub> С/П y<sub>1</sub> С/П x<sub>2</sub> С/П y<sub>2</sub> С/П ...

Part 2:
<lang>ИП5	ПC	ИП6	ПD	П2	ИП7	П3	ИП4	ИПD	*
ИПC	ИП5	*	-	ПD	ИП4	ИП7	*	ИПC	ИП6
*	-	П7	ИП4	ИПA	*	ИПC	ИП9	*	-
ПA	ИП4	ИП3	*	ИП2	ИП5	*	-	П3	ИП4
ИП8	*	ИП2	ИП6	*	-	П8	ИП4	ИПB	*
ИП2	ИП9	*	-	ИПD	*	ИП3	ИПA	*	-
ИПD	ИП8	*	ИП7	ИП3	*	-	/	ПB	ИПA
ИПB	ИП7	*	-	ИПD	/	ПA	ИП9	ИПB	ИП6
*	-	ИПA	ИП5	*	-	ИП4	/	П9	С/П

''Result'': Р9 = a₀, РA = a₁, РB = a₂.

=={{header|Modula-2}}==

MODULE PolynomialRegression;
FROM FormatString IMPORT FormatString;
FROM RealStr IMPORT RealToStr;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;

PROCEDURE Eval(a,b,c,x : REAL) : REAL;
BEGIN
    RETURN a + b*x + c*x*x;
END Eval;

PROCEDURE Regression(x,y : ARRAY OF INTEGER);
VAR
    n,i : INTEGER;
    xm,x2m,x3m,x4m : REAL;
    ym : REAL;
    xym,x2ym : REAL;
    sxx,sxy,sxx2,sx2x2,sx2y : REAL;
    a,b,c : REAL;
    buf : ARRAY[0..63] OF CHAR;
BEGIN
    n := SIZE(x)/SIZE(INTEGER);

    xm := 0.0;
    ym := 0.0;
    x2m := 0.0;
    x3m := 0.0;
    x4m := 0.0;
    xym := 0.0;
    x2ym := 0.0;
    FOR i:=0 TO n-1 DO
        xm := xm + FLOAT(x[i]);
        ym := ym + FLOAT(y[i]);
        x2m := x2m + FLOAT(x[i]) * FLOAT(x[i]);
        x3m := x3m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]);
        x4m := x4m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]);
        xym := xym + FLOAT(x[i]) * FLOAT(y[i]);
        x2ym := x2ym + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(y[i]);
    END;
    xm := xm / FLOAT(n);
    ym := ym / FLOAT(n);
    x2m := x2m / FLOAT(n);
    x3m := x3m / FLOAT(n);
    x4m := x4m / FLOAT(n);
    xym := xym / FLOAT(n);
    x2ym := x2ym / FLOAT(n);

    sxx := x2m - xm * xm;
    sxy := xym - xm * ym;
    sxx2 := x3m - xm * x2m;
    sx2x2 := x4m - x2m * x2m;
    sx2y := x2ym - x2m * ym;

    b := (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
    c := (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
    a := ym - b * xm - c * x2m;

    WriteString("y = ");
    RealToStr(a, buf);
    WriteString(buf);
    WriteString(" + ");
    RealToStr(b, buf);
    WriteString(buf);
    WriteString("x + ");
    RealToStr(c, buf);
    WriteString(buf);
    WriteString("x^2");
    WriteLn;

    FOR i:=0 TO n-1 DO
        FormatString("%2i %3i  ", buf, x[i], y[i]);
        WriteString(buf);
        RealToStr(Eval(a,b,c,FLOAT(x[i])), buf);
        WriteString(buf);
        WriteLn;
    END;
END Regression;

TYPE R = ARRAY[0..10] OF INTEGER;
VAR
    x,y : R;
BEGIN
    x := R{0,1,2,3,4,5,6,7,8,9,10};
    y := R{1,6,17,34,57,86,121,162,209,262,321};
    Regression(x,y);

    ReadChar;
END PolynomialRegression.

Octave

x = [0:10];
y = [1,   6,  17,  34,  57,  86, 121, 162, 209, 262, 321];
coeffs = polyfit(x, y, 2)

PARI/GP

Lagrange interpolating polynomial:

polinterpolate([0,1,2,3,4,5,6,7,8,9,10],[1,6,17,34,57,86,121,162,209,262,321])

In newer versions, this can be abbreviated:

polinterpolate([0..10],[1,6,17,34,57,86,121,162,209,262,321])

3*x^2 + 2*x + 1

Least-squares fit:

V=[1,6,17,34,57,86,121,162,209,262,321]~;
M=matrix(#V,3,i,j,(i-1)^(j-1));Polrev(matsolve(M~*M,M~*V))

Code thanks to [http://pari.math.u-bordeaux.fr/archives/pari-users-1105/msg00006.html Bill Allombert]

3*x^2 + 2*x + 1

Least-squares polynomial fit in its own function:

lsf(X,Y,n)=my(M=matrix(#X,n+1,i,j,X[i]^(j-1))); Polrev(matsolve(M~*M,M~*Y~))
lsf([0..10], [1,6,17,34,57,86,121,162,209,262,321], 2)

Perl

This script depends on the Math::MatrixReal CPAN module to compute matrix determinants.

use strict;
use warnings;
use feature 'say';

#This is a script to calculate an equation for a given set of coordinates.
#Input will be taken in sets of x and y. It can handle a grand total of 26 pairs.
#For matrix functions, we depend on the Math::MatrixReal package.
use Math::MatrixReal;

=pod

Step 1: Get each x coordinate all at once (delimited by " ") and each for y at once
on the next prompt in the same format (delimited by " ").
=cut

sub getPairs() {
    my $buffer = <STDIN>;
    chomp($buffer);
    return split(" ", $buffer);
}
say("Please enter the values for the x coordinates, each delimited by a space. \(Ex: 0 1 2 3\)");
my @x = getPairs();
say("Please enter the values for the y coordinates, each delimited by a space. \(Ex: 0 1 2 3\)");
my @y = getPairs();
#This whole thing depends on the number of x's being the same as the number of y's
my $pairs = scalar(@x);

=pod

Step 2: Devise the base equation of our polynomial using the following idea
There is some polynomial of degree n (n == number of pairs - 1) such that
f(x)=ax^n + bx^(n-1) + ... yx + z
=cut

#Create an array of coefficients and their degrees with the format ("coefficent degree")
my @alphabet;
my @degrees;
for(my $alpha = "a", my $degree = $pairs - 1; $degree >= 0; $degree--, $alpha++) {
    push(@alphabet, "$alpha");
    push(@degrees, "$degree");
}


=pod

Step 3: Using the array of coeffs and their degrees, set up individual equations solving for
each coordinate pair. Why put it in this format? It interfaces witht he Math::MatrixReal package better this way.
=cut

my @coeffs;
for(my $count = 0; $count < $pairs; $count++) {
    my $buffer = "[ ";
    foreach (@degrees) {
        $buffer .= (($x[$count] ** $_) . " ");
    }
    push(@coeffs, ($buffer . "]"));
}
my $row;
foreach (@coeffs) {
    $row .= ("$_\n");
}

=pod

Step 4: We now have rows of x's raised to powers. With this in mind, we create a coefficient matrix.
=cut

my $matrix = Math::MatrixReal->new_from_string($row);
my $buffMatrix = $matrix->new_from_string($row);

=pod

Step 5: Now that we've gotten the matrix to do what we want it to do, we need to calculate the various determinants of the matrices
=cut

my $coeffDet = $matrix->det();

=pod

Step 6: Now that we have the determinant of the coefficient matrix, we need to find the determinants of the coefficient matrix with each column (1 at a time) replaced with the y values.
=cut

#NOTE: Unlike in Perl, matrix indices start at 1, not 0.
for(my $rows = my $column = 1; $column <= $pairs; $column++) {
    #Reassign the values in the current column to the y values
    foreach (@y) {
        $buffMatrix->assign($rows, $column, $_);
        $rows++;
    }
    #Find the values for the variables a, b, ... y, z in the original polynomial
    #To round the difference of the determinants, I had to get creative
    my $buffDet = $buffMatrix->det() / $coeffDet;
    my $tempDet = int(abs($buffDet) + .5);
    $alphabet[$column - 1] = $buffDet >= 0 ? $tempDet : 0 - $tempDet;
    #Reset the buffer matrix and the row counter
    $buffMatrix = $matrix->new_from_string($row);
    $rows = 1;
}


=pod

Step 7: Now that we've found the values of a, b, ... y, z of the original polynomial, it's time to form our polynomial!
=cut

my $polynomial;
for(my $i = 0; $i < $pairs-1; $i++) {
    if($alphabet[$i] == 0) {
        next;
    }
    if($alphabet[$i] == 1) {
        $polynomial .= ($degrees[$i] . " + ");
    }
    if($degrees[$i] == 1) {
        $polynomial .= ($alphabet[$i] . "x" . " + ");
    }
    else {
        $polynomial .= ($alphabet[$i] . "x^" . $degrees[$i] . " + ");
    }
}
#Now for the last piece of the poly: the y-intercept.
$polynomial .= $alphabet[scalar(@alphabet)-1];

print("An approximating polynomial for your dataset is $polynomial.\n");

Please enter the values for the x coordinates, each delimited by a space. (Ex: 0 1 2 3)
0 1 2 3 4 5 6 7 8 9 10
Please enter the values for the y coordinates, each delimited by a space. (Ex: 0 1 2 3)
1 6 17 34 57 86 121 162 209 262 321
An approximating polynomial for your dataset is 3x^2 + 2x + 1.

Perl 6

We'll use a Clifford algebra library.

use Clifford;

constant @x1 = <0 1 2 3 4 5 6 7 8 9 10>;
constant @y = <1 6 17 34 57 86 121 162 209 262 321>;

constant $x0 = [+] @e[^@x1];
constant $x1 = [+] @x1 Z* @e;
constant $x2 = [+] @x1 »**» 2  Z* @e;

constant $y  = [+] @y Z* @e;

my $J = $x1 ∧ $x2;
my $I = $x0 ∧ $J;

my $I2 = ($I·$I.reversion).Real;

.say for
(($y ∧ $J)·$I.reversion)/$I2,
(($y ∧ ($x2 ∧ $x0))·$I.reversion)/$I2,
(($y ∧ ($x0 ∧ $x1))·$I.reversion)/$I2;

1
2
3

Phix

constant x = {0,1,2,3,4,5,6,7,8,9,10}
constant y = {1,6,17,34,57,86,121,162,209,262,321}
constant n = length(x)

function regression()
atom {xm, ym, x2m, x3m, x4m, xym, x2ym} @= 0
    for i=1 to n do
        atom xi = x[i],
             yi = y[i]
        xm += xi
        ym += yi
        x2m += power(xi,2)
        x3m += power(xi,3)
        x4m += power(xi,4)
        xym += xi*yi
        x2ym += power(xi,2)*yi
    end for
    xm /= n
    ym /= n
    x2m /= n
    x3m /= n
    x4m /= n
    xym /= n
    x2ym /= n
    atom Sxx = x2m-power(xm,2),
         Sxy = xym-xm*ym,
         Sxx2 = x3m-xm*x2m,
         Sx2x2 = x4m-power(x2m,2),
         Sx2y = x2ym-x2m*ym,
         B = (Sxy*Sx2x2-Sx2y*Sxx2)/(Sxx*Sx2x2-power(Sxx2,2)),
         C = (Sx2y*Sxx-Sxy*Sxx2)/(Sxx*Sx2x2-power(Sxx2,2)),
         A = ym-B*xm-C*x2m
    return {C,B,A}
end function

atom {a,b,c} = regression()

function f(atom x)
    return a*x*x+b*x+c
end function

printf(1,"y=%gx^2+%gx+%g\n",{a,b,c})
printf(1,"\n  x   y  f(x)\n")
for i=1 to n do
  printf(1," %2d %3d   %3g\n",{x[i],y[i],f(x[i])})
end for

y=3x^2+2x+1

  x   y  f(x)
  0   1     1
  1   6     6
  2  17    17
  3  34    34
  4  57    57
  5  86    86
  6 121   121
  7 162   162
  8 209   209
  9 262   262
 10 321   321

Alternatively, a simple plot, (as per [[Polynomial_regression#Racket|Racket]]):

include pGUI.e

constant x = {0,1,2,3,4,5,6,7,8,9,10}
constant y = {1,6,17,34,57,86,121,162,209,262,321}

IupOpen()

Ihandle plot = IupPlot("GRID=YES, MARGINLEFT=50, MARGINBOTTOM=40")
             -- (just add ", AXS_YSCALE=LOG10" for a nice log scale)
IupPlotBegin(plot, 0)
for i=1 to length(x) do
    IupPlotAdd(plot, x[i], y[i])
end for
{} = IupPlotEnd(plot)

Ihandle dlg = IupDialog(plot)
IupSetAttributes(dlg, "RASTERSIZE=%dx%d", {640, 480})
IupSetAttribute(dlg, "TITLE", "simple plot")
IupShow(dlg)

IupMainLoop()
IupClose()

PowerShell

function qr([double[][]]$A) {
    $m,$n = $A.count, $A[0].count
    $pm,$pn = ($m-1), ($n-1)
    [double[][]]$Q = 0..($m-1) | foreach{$row = @(0) * $m; $row[$_] = 1; ,$row}
    [double[][]]$R = $A | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
    foreach ($h in 0..$pn) {
        [double[]]$u = $R[$h..$pm] | foreach{$_[$h]}
        [double]$nu = $u | foreach {[double]$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
        $u[0] -= if ($u[0] -lt 1) {$nu} else {-$nu}
        [double]$nu = $u | foreach {$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
        [double[]]$u = $u | foreach { $_/$nu}
        [double[][]]$v = 0..($u.Count - 1) | foreach{$i = $_; ,($u | foreach{2*$u[$i]*$_})}
        [double[][]]$CR = $R | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
        [double[][]]$CQ = $Q | foreach{$row = $_; ,@(0..$pm | foreach{$row[$_]})}
        foreach ($i in  $h..$pm) {
            foreach ($j in  $h..$pn) {
                $R[$i][$j] -=  $h..$pm | foreach {[double]$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CR[$_][$j]} {$sum}
            }
        }
        if (0 -eq $h)  {
            foreach ($i in  $h..$pm) {
                foreach ($j in  $h..$pm) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i][$_]*$CQ[$_][$j]} {$sum}
                }
            }
        } else  {
            $p = $h-1
            foreach ($i in  $h..$pm) {
                foreach ($j in  0..$p) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
                }
                foreach ($j in  $h..$pm) {
                    $Q[$i][$j] -=  $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
                }
            }
        }
    }
    foreach ($i in  0..$pm) {
        foreach ($j in  $i..$pm) {$Q[$i][$j],$Q[$j][$i] = $Q[$j][$i],$Q[$i][$j]}
    }
    [PSCustomObject]@{"Q" = $Q; "R" = $R}
}

function leastsquares([Double[][]]$A,[Double[]]$y) {
    $QR = qr $A
    [Double[][]]$Q = $QR.Q
    [Double[][]]$R = $QR.R
    $m,$n = $A.count, $A[0].count
    [Double[]]$z = foreach ($j in  0..($m-1)) {
            0..($m-1) | foreach {$sum = 0} {$sum += $Q[$_][$j]*$y[$_]} {$sum}
    }
    [Double[]]$x = @(0)*$n
    for ($i = $n-1; $i -ge 0; $i--) {
        for ($j = $i+1; $j -lt $n; $j++) {
            $z[$i] -= $x[$j]*$R[$i][$j]
        }
        $x[$i] = $z[$i]/$R[$i][$i]
    }
    $x
}

function polyfit([Double[]]$x,[Double[]]$y,$n) {
    $m = $x.Count
    [Double[][]]$A = 0..($m-1) | foreach{$row = @(1) * ($n+1); ,$row}
    for ($i = 0; $i -lt $m; $i++) {
        for ($j = $n-1; 0 -le $j; $j--) {
            $A[$i][$j] = $A[$i][$j+1]*$x[$i]
        }
    }
    leastsquares $A $y
}

function show($m) {$m | foreach {write-host "$_"}}

$A = @(@(12,-51,4), @(6,167,-68), @(-4,24,-41))
$x = @(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
$y = @(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
"polyfit "
"X^2 X constant"
"$(polyfit $x $y 2)"

polyfit
X^2 X constant
3 1.99999999999998 1.00000000000005

Python

 x = [0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10]
>>> y = [1,   6,  17,  34,  57,  86, 121, 162, 209, 262, 321]
>>> coeffs = numpy.polyfit(x,y,deg=2)
>>> coeffs
array([ 3.,  2.,  1.])

Substitute back received coefficients.

 yf = numpy.polyval(numpy.poly1d(coeffs), x)
>>> yf
array([   1.,    6.,   17.,   34.,   57.,   86.,  121.,  162.,  209., 262.,  321.])

Find max absolute error:

 '%.1g' % max(y-yf)
'1e-013'

Example

For input arrays x' and y':

 x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]

 p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N')
>>> print p
       2
1.085 N + 10.36 N - 0.6164

Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see [[Query Performance|'''Example''' section for Python language on ''Query Performance'' page]]).

R

The easiest (and most robust) approach to solve this in R is to use the base package's ''lm'' function which will find the least squares solution via a QR decomposition:

x <- c(0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10)
y <- c(1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321)
coef(lm(y ~ x + I(x^2)))

(Intercept)           x      I(x^2)
          1           2           3

Alternately, use poly:

coef(lm(y ~ poly(x, 2, raw=T)))

         (Intercept) poly(x, 2, raw = T)1 poly(x, 2, raw = T)2
                   1                    2                    3

Racket

#lang racket
(require math plot)

(define xs '(0 1  2  3  4  5   6   7   8   9  10))
(define ys '(1 6 17 34 57 86 121 162 209 262 321))

(define (fit x y n)
  (define Y (->col-matrix y))
  (define V (vandermonde-matrix x (+ n 1)))
  (define VT (matrix-transpose V))
  (matrix->vector (matrix-solve (matrix* VT V) (matrix* VT Y))))

(define ((poly v) x)
  (for/sum ([c v] [i (in-naturals)])
    (* c (expt x i))))

(plot (list (points   (map vector xs ys))
            (function (poly (fit xs ys 2)))))

[[File:polyreg-racket.png]]

REXX

/* REXX ---------------------------------------------------------------
* Implementation of http://keisan.casio.com/exec/system/14059932254941
*--------------------------------------------------------------------*/
xl='0 1  2  3  4  5   6   7   8   9  10'
yl='1 6 17 34 57 86 121 162 209 262 321'
n=11
Do i=1 To n
  Parse Var xl x.i xl
  Parse Var yl y.i yl
  End
xm=0
ym=0
x2m=0
x3m=0
x4m=0
xym=0
x2ym=0
Do i=1 To n
  xm=xm+x.i
  ym=ym+y.i
  x2m=x2m+x.i**2
  x3m=x3m+x.i**3
  x4m=x4m+x.i**4
  xym=xym+x.i*y.i
  x2ym=x2ym+(x.i**2)*y.i
  End
xm =xm /n
ym =ym /n
x2m=x2m/n
x3m=x3m/n
x4m=x4m/n
xym=xym/n
x2ym=x2ym/n
Sxx=x2m-xm**2
Sxy=xym-xm*ym
Sxx2=x3m-xm*x2m
Sx2x2=x4m-x2m**2
Sx2y=x2ym-x2m*ym
B=(Sxy*Sx2x2-Sx2y*Sxx2)/(Sxx*Sx2x2-Sxx2**2)
C=(Sx2y*Sxx-Sxy*Sxx2)/(Sxx*Sx2x2-Sxx2**2)
A=ym-B*xm-C*x2m
Say 'y='a'+'||b'*x+'c'*x**2'
Say ' Input  "Approximation"'
Say ' x   y     y1'
Do i=1 To 11
  Say right(x.i,2) right(y.i,3) format(fun(x.i),5,3)
  End
Exit
fun:
  Parse Arg x
  Return a+b*x+c*x**2

y=1+2*x+3*x**2
 Input  "Approximation"
 x   y     y1
 0   1     1.000
 1   6     6.000
 2  17    17.000
 3  34    34.000
 4  57    57.000
 5  86    86.000
 6 121   121.000
 7 162   162.000
 8 209   209.000
 9 262   262.000
10 321   321.000

Ruby

require 'matrix'

def regress x, y, degree
  x_data = x.map { |xi| (0..degree).map { |pow| (xi**pow).to_r } }

  mx = Matrix[*x_data]
  my = Matrix.column_vector(y)

  ((mx.t * mx).inv * mx.t * my).transpose.to_a[0].map(&:to_f)
end

'''Testing:'''

p regress([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
          [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321],
          2)

[1.0, 2.0, 3.0]

Scala

{{Out}}See it yourself by running in your browser [https://scastie.scala-lang.org/NklZH2LlScCpfsN4NSfFvA Scastie (remote JVM)].

object PolynomialRegression extends App {
  private def xy = Seq(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321).zipWithIndex.map(_.swap)

  private def polyRegression(xy: Seq[(Int, Int)]): Unit = {
    val r = xy.indices

    def average[U](ts: Iterable[U])(implicit num: Numeric[U]) = num.toDouble(ts.sum) / ts.size

    def x3m: Double = average(r.map(a => a * a * a))
    def x4m: Double = average(r.map(a => a * a * a * a))
    def x2ym = xy.reduce((a, x) => (a._1 + x._1 * x._1 * x._2, 0))._1.toDouble / xy.size
    def xym = xy.reduce((a, x) => (a._1 + x._1 * x._2, 0))._1.toDouble / xy.size

    val x2m: Double = average(r.map(a => a * a))
    val (xm, ym) = (average(xy.map(_._1)), average(xy.map(_._2)))
    val (sxx, sxy) = (x2m - xm * xm, xym - xm * ym)
    val sxx2: Double = x3m - xm * x2m
    val sx2x2: Double = x4m - x2m * x2m
    val sx2y: Double = x2ym - x2m * ym
    val c: Double = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    val b: Double = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
    val a: Double = ym - b * xm - c * x2m

    def abc(xx: Int) = a + b * xx + c * xx * xx

    println(s"y = $a + ${b}x + ${c}x^2")
    println(" Input  Approximation")
    println(" x   y     y1")
    xy.foreach {el => println(f"${el._1}%2d ${el._2}%3d  ${abc(el._1)}%5.1f")}
  }

  polyRegression(xy)

}

Sidef

func regress(x, y, degree) {
    var A = Matrix.build(x.len, degree+1, {|i,j|
        x[i]**j
    })

    var B = Matrix.column_vector(y...)
    ((A.transpose * A)**(-1) * A.transpose * B).transpose[0]
}

func poly(x) {
    3*x**2 + 2*x + 1
}

var coeff = regress(
    10.of { _ },
    10.of { poly(_) },
    2
)

say coeff

[1, 2, 3]

Stata

See '''[http://www.stata.com/help.cgi?fvvarlist Factor variables]''' in Stata help for explanations on the ''c.x##c.x'' syntax.

. clear
. input x y
0 1
1 6
2 17
3 34
4 57
5 86
6 121
7 162
8 209
9 262
10 321
end

. regress y c.x##c.x

      Source |       SS           df       MS      Number of obs   =        11
-------------+----------------------------------   F(2, 8)         =         .
       Model |      120362         2       60181   Prob > F        =         .
    Residual |           0         8           0   R-squared       =    1.0000
-------------+----------------------------------   Adj R-squared   =    1.0000
       Total |      120362        10     12036.2   Root MSE        =         0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |          2          .        .       .            .           .
             |
     c.x#c.x |          3          .        .       .            .           .
             |
       _cons |          1          .        .       .            .           .
------------------------------------------------------------------------------

Tcl

package require math::linearalgebra

proc build.matrix {xvec degree} {
    set sums [llength $xvec]
    for {set i 1} {$i <= 2*$degree} {incr i} {
        set sum 0
        foreach x $xvec {
            set sum [expr {$sum + pow($x,$i)}]
        }
        lappend sums $sum
    }

    set order [expr {$degree + 1}]
    set A [math::linearalgebra::mkMatrix $order $order 0]
    for {set i 0} {$i <= $degree} {incr i} {
        set A [math::linearalgebra::setrow A $i [lrange $sums $i $i+$degree]]
    }
    return $A
}

proc build.vector {xvec yvec degree} {
    set sums [list]
    for {set i 0} {$i <= $degree} {incr i} {
        set sum 0
        foreach x $xvec y $yvec {
            set sum [expr {$sum + $y * pow($x,$i)}]
        }
        lappend sums $sum
    }

    set x [math::linearalgebra::mkVector [expr {$degree + 1}] 0]
    for {set i 0} {$i <= $degree} {incr i} {
        set x [math::linearalgebra::setelem x $i [lindex $sums $i]]
    }
    return $x
}

# Now, to solve the example from the top of this page
set x {0   1   2   3   4   5   6   7   8   9  10}
set y {1   6  17  34  57  86 121 162 209 262 321}

# build the system A.x=b
set degree 2
set A [build.matrix $x $degree]
set b [build.vector $x $y $degree]
# solve it
set coeffs [math::linearalgebra::solveGauss $A $b]
# show results
puts $coeffs

This will print: 1.0000000000000207 1.9999999999999958 3.0 which is a close approximation to the correct solution.

=={{header|TI-89 BASIC}}==

DelVar x
seq(x,x,0,10) → xs
{1,6,17,34,57,86,121,162,209,262,321} → ys
QuadReg xs,ys
Disp regeq(x)

seq(''expr'',''var'',''low'',''high'') evaluates ''expr'' with ''var'' bound to integers from ''low'' to ''high'' and returns a list of the results. → is the assignment operator. QuadReg, "quadratic regression", does the fit and stores the details in a number of standard variables, including regeq, which receives the fitted quadratic (polynomial) function itself. We then apply that function to the (undefined as ensured by DelVar) variable x to obtain the expression in terms of x, and display it.

3.·x2 + 2.·x + 1.

Ursala

The fit function defined below returns the coefficients of an nth-degree polynomial in order of descending degree fitting the lists of inputs x and outputs y. The real work is done by the dgelsd function from the lapack library. Ursala provides a simplified interface to this library whereby the data can be passed as lists rather than arrays, and all memory management is handled automatically.

#import std
#import nat
#import flo

(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"

test program:

x = <0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.>
y = <1.,6.,17.,34.,57.,86.,121.,162.,209.,262.,321.>

#cast %eL

example = fit2(x,y)

<3.000000e+00,2.000000e+00,1.000000e+00>

VBA

Excel VBA has built in capability for line estimation.

Option Base 1
Private Function polynomial_regression(y As Variant, x As Variant, degree As Integer) As Variant
    Dim a() As Double
    ReDim a(UBound(x), 2)
    For i = 1 To UBound(x)
        For j = 1 To degree
            a(i, j) = x(i) ^ j
        Next j
    Next i
    polynomial_regression = WorksheetFunction.LinEst(WorksheetFunction.Transpose(y), a, True, True)
End Function
Public Sub main()
    x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]
    y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}]
    result = polynomial_regression(y, x, 2)
    Debug.Print "coefficients   : ";
    For i = UBound(result, 2) To 1 Step -1
        Debug.Print Format(result(1, i), "0.#####"),
    Next i
    Debug.Print
    Debug.Print "standard errors: ";
    For i = UBound(result, 2) To 1 Step -1
        Debug.Print Format(result(2, i), "0.#####"),
    Next i
    Debug.Print vbCrLf
    Debug.Print "R^2 ="; result(3, 1)
    Debug.Print "F   ="; result(4, 1)
    Debug.Print "Degrees of freedom:"; result(4, 2)
    Debug.Print "Standard error of y estimate:"; result(3, 2)
End Sub

coefficients   : 1,         2,            3,
standard errors: 0,         0,            0,

R^2 = 1
F   = 7,70461300500498E+31
Degrees of freedom: 8
Standard error of y estimate: 2,79482284961344E-14

zkl

Using the GNU Scientific Library

var [const] GSL=Import("zklGSL");	// libGSL (GNU Scientific Library)
xs:=GSL.VectorFromData(0,  1,  2,  3,  4,  5,   6,   7,   8,   9,  10);
ys:=GSL.VectorFromData(1,  6, 17, 34, 57, 86, 121, 162, 209, 262, 321);
v :=GSL.polyFit(xs,ys,2);
v.format().println();
GSL.Helpers.polyString(v).println();
GSL.Helpers.polyEval(v,xs).format().println();

1.00,2.00,3.00
 1 + 2x + 3x^2
1.00,6.00,17.00,34.00,57.00,86.00,121.00,162.00,209.00,262.00,321.00

Or, using lists: Uses the code from [[Multiple regression#zkl]].

Example:

polyfit(T(T(0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0)),
        T(T(1.0,6.0,17.0,34.0,57.0,86.0,121.0,162.0,209.0,262.0,321.0)), 2)
.flatten().println();

L(1,2,3)