INSTALL @lib$+"ARRAYLIB"

      Max% = 10000
      DIM vector(5), matrix(5,5)
      DIM x(Max%), x2(Max%), x3(Max%), x4(Max%), x5(Max%)
      DIM x6(Max%), x7(Max%), x8(Max%), x9(Max%), x10(Max%)
      DIM y(Max%), xy(Max%), x2y(Max%), x3y(Max%), x4y(Max%), x5y(Max%)

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

      sum_x = SUM(x())
      x2()  = x() * x()   : sum_x2  = SUM(x2())
      x3()  = x() * x2()  : sum_x3  = SUM(x3())
      x4()  = x2() * x2() : sum_x4  = SUM(x4())
      x5()  = x2() * x3() : sum_x5  = SUM(x5())
      x6()  = x3() * x3() : sum_x6  = SUM(x6())
      x7()  = x3() * x4() : sum_x7  = SUM(x7())
      x8()  = x4() * x4() : sum_x8  = SUM(x8())
      x9()  = x4() * x5() : sum_x9  = SUM(x9())
      x10() = x5() * x5() : sum_x10 = SUM(x10())

      sum_y = SUM(y())
      xy()  = x() * y()   : sum_xy  = SUM(xy())
      x2y() = x2() * y()  : sum_x2y = SUM(x2y())
      x3y() = x3() * y()  : sum_x3y = SUM(x3y())
      x4y() = x4() * y()  : sum_x4y = SUM(x4y())
      x5y() = x5() * y()  : sum_x5y = SUM(x5y())

      matrix() = \
      \ npts%,  sum_x,   sum_x2,  sum_x3,  sum_x4,  sum_x5, \
      \ sum_x,  sum_x2,  sum_x3,  sum_x4,  sum_x5,  sum_x6, \
      \ sum_x2, sum_x3,  sum_x4,  sum_x5,  sum_x6,  sum_x7, \
      \ sum_x3, sum_x4,  sum_x5,  sum_x6,  sum_x7,  sum_x8, \
      \ sum_x4, sum_x5,  sum_x6,  sum_x7,  sum_x8,  sum_x9, \
      \ sum_x5, sum_x6,  sum_x7,  sum_x8,  sum_x9,  sum_x10

      vector() = \
      \ sum_y,  sum_xy,  sum_x2y, sum_x3y, sum_x4y, sum_x5y

      PROC_invert(matrix())
      vector() = matrix().vector()

      @% = &2040A
      PRINT "Polynomial coefficients = "
      FOR term% = 5 TO 0 STEP -1
        PRINT ;vector(term%) " * x^" STR$(term%)
      NEXT