LEX    'HORNER'      * HORNER'S SCHEME FOR THE EVALUATION OF
        ID     #EC           * POLYNOMIALS UP TO DEGREE 13.
        MSG    0
        POLL   0
        ENTRY  hORNER
        CHAR   #F
        KEY    'INTERP'      * FOR SYNTAX SEE BELOW
        TOKEN  13
FUNCD1  EQU    #2F8C0        * A SAFE PLACE TO STORE D1
POP1R   EQU    #0E8FD        * POP 1 REAL, ERROR ON COMPLEX
SPLITA  EQU    #0C6BF        * SPLIT A INTO A 15-FORM IN (A,B)
MP2-15  EQU    #0C43A        * MULTIPLY TWO 15-FORMS
AD2-15  EQU    #0C363        * ADD TWO 15-FORMS
uRES12  EQU    #0C994        * 15-FORM IN (A,B) TO 12-FORM IN C(W)
FNRTN1  EQU    #0F216        * FNCTION EXIT
        ENDTXT
        NIBHEX 888888888888888  * ALL NUMERIC PARAMETERS
        NIBHEX 3F            * AT LEAST 3, AT MOST 15 OF THEM!
hORNER  R0=C                 * SAVE C IN R0 - C(S) HOLDS NUMBER OF PARAMETERS
        CD1EX                * ACTUALLY PASSED, THEN SAVE D1
        D1=(5) FUNCD1        * IN FUNCD1.
        DAT1=C A
        CD1EX
        CD0EX
        R2=C                 * SAVE D0 IN R2
        CD0EX
        C=R0
        C=C-1  S             * DECREMENT C(S) BY 3 BECAUSE DEGREE IS 3 LESS
        C=C-1  S             * THAN NUMBER OF PARAMETERS
        C=C-1  S
        R3=C                 * SAVE DEGREE IN R3
        GOSBVL POP1R         * GET THE ARGUMENT OF THE POLYNOMIAL
        GOSBVL SPLITA
        R0=A                 * SAVE IN R0 AND R1
        A=B    W             * REGISTER B WON'T TALK TO SCRATCH REGISTERS
        R1=A
        D1=D1+ 16            * SHIFT ALONG STACK FOR NEXT PARAMETER
        GOSBVL POP1R         * COLLECT IT (THIS ONE IS THE LEADING 
        GOSBVL SPLITA        * COEFFICIENT)
loop    C=R1                 * RECOVER
        D=C    W             * THE POLYNOMIAL
        C=R0                 * ARGUMENT.
        GOSBVL MP2-15        * MULTIPLY
        C=B    W             * COPY (A,B) INTO (C,D)
        D=C    W
        C=A    W
        D1=D1+ 16            * MOVE D1 TO PICK UP NEXT COEFFICIENT
        GOSBVL POP1R
        GOSBVL SPLITA
        GOSBVL AD2-15        * ADD
        CR3EX                * R3 CONTAINS THE COUNTER
        C=C-1  S             * WHICH IS NOW DECREMENTED
        CR3EX                * AND REPLACED
        GONC   loop          * CARRY SET IF NO MORE COEFFICIENTS
        GOSBVL uRES12        * RESULT IS PUT INTO C(W) FOR FNRTN1
        D1=(5) FUNCD1        * RESTORE D1 FOR OUTPUT
        A=DAT1 A             * PUT IT IN A(A)
        D1=A
        A=R2
        GOSBVL FNRTN1        * EXIT
        END


SYNTAX: INTERP(a0,a1[,a2[,a3[,a4[,a5[,a6[,a7.....[,a13]]]...],x), where the
subscripted a's are the coefficients and x is the polynomial argument.
In standard polynomial notation, we want to evaluate:
f(x)=a0+a1x+a2x^2+a3x^3+....+a13x^13 .
This is a very poor method numerically as it leads to severe rounding errors
through exponentiation. Horner's scheme, however, is not subject to that
since no exponentiation is done:
f(x)=a0+x(a1+x(a2+x(a3+x(....(a12+a13x))...).