/*---------------------------------------------------------------------------
bspline.asm:    Cubic Bspline Polynomial Evaluation
-----------------------------------------------------------------------------
Description:

The Cubic Bspline polynomial is commonly used to represent curves and surfaces
as a compact list of control points, or "knots". 

This routine evaluates the Cubic Bspline polynomial:

              3         2                      3    2                   3
        -(u-1) p  + [ 3u (u-2) + 4 ] p  - [ 3(u  + u  + u) + 1) ] p  + u p 
                1                     2                            3      4
x(u) = ----------------------------------------------------------------------
                                           6.0

-----------------------------------------------------------------------------
Program Characteristics:

Calling Values:
    r0  = number of U values
    f14 = 1.0
    i8  = PM  pointer to P1 --> P4
    i0  = DM pointer to U values
    i1  = DM pointer to constants { 3.0, 4.0, 0.1666666666 }
    i2  = DM pointer to results x(U)

Return Value
    i2  = pointer to results x(U)

Computation Time 
    = 15N + 11 cycles per bspline evaluation
    = 600ns per bspline evaluation @ 40ns instruction cycle
    = 1.7 million bspline evaluations/sec @ 40ns instruction cycle

-----------------------------------------------------------------------------
Notes:

The register defines used in the code below refer to the various sections 
of the factored eqn. being assembled:

               3         2                      3    2                   3
         -(u-1) p  + [ 3u (u-2) + 4 ] p  - [ 3(u  + u  + u) + 1) ] p  + u p 
                 1                     2                            3      4
          --M1--       ---A6---                ----A5-----
                       -----M2-----          -------M3---------
         ---A1----    -------A2---------    ----------A3------------    -A4-
         ------------AX-----------------    ---------------AY---------------

This assembly code which this notation produces is a bit easier to follow
than using the actual registers, but both are somewhat cumbersome.

-----------------------------------------------------------------------------
Author:     Jim Donahue, Analog Devices DSP Division
Revised:    11-SEP-91
References: Foley, James, and VanDam, Andries, et.al. "Computer Graphics,
            Principles and Practice", Addison-Wesley 1990, p.491-495
----------------------------------------------------------------------------*/
.GLOBAL bspline;

/* Multiplier Input defines */

#define U       f10     /* U, the parametric value  */

#define X       f0      /* result x(U)              */

#define M1      f0
#define M2      f1
#define M3      f2

#define MUS     f5      /* U^2, MULT version        */
#define MUM1    f6      /* U - 1 */
#define MUM2    f3      /* U - 2 */

#define P1      f4      /* parametric coef #1       */
#define P2      f4      /* parametric coef #2       */
#define P3      f4      /* parametric coef #3       */
#define P4      f4      /* parametric coef #4       */

#define THREE   f7      /* 3.0                      */
#define ISIX    f15     /* 0.16666666666            */
#define ONE     f14


/* ALU Input defines */

#define AUC     f12     /* U^3, ALU version         */
#define AUM1    f9      /* U-1, ALU version         */
#define A1      f15
#define A2      f11
#define A3      f9
#define A4      f10
#define A5      f9
#define A6      f11
#define FOUR    f15     /* 4.0                      */
#define AU      f8      /* U, ALU version           */
#define AUS     f13     /* U^2, ALU version         */
#define AX      f8      /* sum of 1st 2 terms       */
#define AY      f12     /* sum of last 2 terms      */

.SEGMENT /pm    pm_code;
bspline:

    U  = dm(i0,m1);
    AU = dm(i0,m0);
    MUM1 = U - ONE;

    lcntr = r0, do bslp until lce;

        M1   = MUM1 * MUM1,                             AUM1 = MUM1;
        M1   = M1 * MUM1,       MUM2 = AUM1 - ONE;
        MUS  = U * U;   
        AUC  = U * MUS,         AUS  = MUS;
        M2   = MUM2 * MUS,      A5  = AU + AUS,         THREE = dm(i1,m0);
        A6   = M2 * THREE,      M3  = A5 + AUC,         FOUR = dm(i1,m0);
        A5   = M3 * THREE,      M2  = A6 + FOUR,        P1 = pm(i8,m8);
        A1   = M1 * P1,         M3  = A5 + ONE,         P2 = pm(i8,m8);
        A2   = M2 * P2,                                 P3 = pm(i8,m8);
        A3   = M3 * P3,         AX  = A2 - A1,          P4 = pm(i8,m8);
        A4   = AUC * P4;
                                AY  = A4 - A3,          U = dm(i0,m1);
                                X  = AX + AY,           ISIX = dm(i1,m0);
        X = X * ISIX,                                   AU = dm(i0,m0);
bslp:                           MUM1 = U - ONE,         dm(i2,m0) = X;

    rts;

.ENDSEG;