%%HP: T(3)A(D)F(.);
\<<
@ Display version while working
"     SPLINE V2.0
  " 1 DISP

  6 CF 7 CF 8 CF 9 CF           @ Prepare flags 6-9

  IF DUP TYPE 6 SAME            @ Check if 'D1' or 'D2' are specified
  THEN
    CASE DUP 'D1' SAME
      THEN 6 SF                 @ Set flag 6 if 'D1' specified
      END DUP 'D2' SAME
      THEN 7 SF                 @ Set flag 7 if 'D2' specified
      END DUP \->STR ": Unknown flag" + DOERR           
    END
    DROP                        @ Drop 'D1' or 'D2' from stack
  END

  IF DUP TYPE 5 \=/             @ See if the end-derivatives are specified
  THEN { 0 0 }                  @ If not, insert {0 0} (as a place
holder)
  ELSE 8 SF                     @ If yes, then set flag 8
  END
  
  IF OVER TYPE 0 SAME           @ If we have a number n in level 2 then we 
  THEN OVER 2 + ROLLD           @ expect n pairs of coordinates above it
  ELSE SWAP ROT                 @ Otherwise we expect two arrays in levels 2
  END                           @ and 3.  In either case, we move the { s1 s2
}
                                @ list form level 1 to the top of the stack.

  IF DUP TYPE 0 \=/             @ If we don't have a number in level 1
  THEN 9 SF                     @ Then the coordinates are given as arrays
  DUP SIZE EVAL                 @ Determine the size of array
  END

                                @ Set up the local variables
     DUP 1 - { } { } 0  0   0  0 0 0 0 0   0    0   0   0
  \->   n  k  x   y  h \Gl \Gm s d m a b \Gl0 \Gmn s11 s1n

  @ Begin the main program
  \<<

    IF 9 FC?C
    THEN 1 n                    @ Flag 9 is clear so data is in coord. pairs
      FOR j                     @ Convert coordinate pairs into lists x and y
          C\->R 'y' STO+ 'x' STO+
      NEXT
    ELSE
      OBJ\-> EVAL \->LIST 'y' STO       @ Store array of y_i into the list y
      OBJ\-> DROP OVER - k /            @ Compute the mesh size
      0 k
      FOR j
        DUP j * 3 PICK + 3 ROLLD        @ Generate the x mesh
      NEXT
      DROP2 n \->LIST 'x' STO           @ Store the x values into the list x
    END

    EVAL 's1n' STO 's11' STO            @ Read the end-derivative values

    x EVAL                              @ Compute the list of
    1 k                                 @ interval lengths h_j = x_{j+1} - x_j
    FOR j OVER - n ROLLD
    NEXT
    DROP
    k \->LIST 'h' STO

    1 k                                 @ Compute the list of slopes s_j
    FOR j                               @ s_j = ( y_{j+1} - y_j ) / h_j
      y j GETI 3 ROLLD GET - NEG h j GET /
    NEXT
    k \->LIST 's' STO

    @ Compute the elements d_j of the list d:
    IF 8 FS?
    THEN                                @ End-derivatives are specified
      1 '\Gl0' STO
      1 '\Gmn' STO
      s 1 GET s11 - h 1 GET / 6 *
    ELSE 0
    END

    @ Still computing d:
    1 n 2 -
    FOR j
      s j GETI 3 ROLLD GET - NEG h j GETI 3 ROLLD GET + / 6 *
    NEXT

    @ End of computation of d:
    IF 8 FS?C
    THEN
      s1n s k GET - h k GET / 6 *
    ELSE 0
    END
    n \->LIST 'd' STO


    @ Compute lambda_j:
    h OBJ\-> 1 - 1 SWAP
    FOR j
      DUP 3 PICK + / k ROLLD
    NEXT
    DROP
    n 2 - \->LIST '\Gl' STO

    @ Compute gamma_j:
    \Gl OBJ\-> 1 SWAP
    FOR j
      NEG 1 +
      n 2 - ROLL
    NEXT
    n 2 - \->LIST '\Gm' STO

    @ Compute the moments m_j:
    n IDN 2 *
    2 k
    FOR j
      j DUP 1 - 2 \->LIST
      \Gm j 1 - GET
      PUT
      j DUP 1 + 2 \->LIST
      \Gl j 1 - GET
      PUT
    NEXT
    2 \Gl0 PUT
    n SQ 1 - \Gmn PUT
    INV
    d OBJ\-> \->ARRY *
    'm' STO

    @ Compute a_j:
    1 k
    FOR j
      m j GETI 3 ROLLD GET - h j GET * 6 / s j GET +
    NEXT
    k \->LIST 'a' STO

    @ Compute b_j:
    1 k
    FOR j
      y j GET m j GET h j GET SQ * 6 / -
    NEXT
    k \->LIST 'b' STO

    @ Now we compute the individual arcs of the spline:
    CASE 6 FS?C
      THEN                              @ Will compute y'(x)
        1 k
        FOR j
          m j 1 + GET 'X' x j GET - SQ *
          m j GET x j 1 + GET 'X' - SQ * -
          h j GET / 2 / a j GET +
        NEXT
      END 7 FS?C
      THEN                              @ Will compute y''(x)
      1 k
      FOR j
        m j GET x j 1 + GET 'X' - *
        m j 1 + GET 'X' x j GET - * +
        h j GET / NEXT
      END
      1 k
      FOR j                             @ Will compute y(x)
      m j GET x j 1 + GET 'X' - 3 ^ *
      m j 1 + GET 'X' x j GET - 3 ^ * +
      h j GET / 6 /
      a j GET 'X' x j GET - * +
      b j GET +
      NEXT
    END

    @ Create the output program:
    "\<<\-> X\<<CASE " n ROLLD
    k 2
    FOR j
      'X' " " + x j GET + " \>= THEN " + SWAP + " END " +
      j ROLLD
      -1
    STEP
   " END EVAL\>>\>>"

   @ Concatenate all parts:
   1 n
   FOR j +
   NEXT
   OBJ\->
  \>>
\>>