.module/boot=3/boot=4        approximate_func;

{ POLY.DSP - Calculates the polynomial approximation to
             a function given by the coefficients.
             Uses 32 bit; y = f(x);

  INPUT:
     my0 = x lsb
     my1 = x msb
     ax0 = POLY_ORDER - 1
     i6  = -> coeffs (in pm) (Ci lsb, Ci msb, Ci-1 lsb .....)

  
  OUTPUT:
     mx0 = y msb
     ar  = y lsb

   f(x) is approximated by a polynomial:

   f(x) = C[0] + C[1]*X^1 .... + C[(POLY_ORDER-1)]*X^(POLY_ORDER-1)
  
}
#include "lpc.h"

.entry    poly_approx;

poly_approx:

  mx0 = pm(i6,m5); {c lsb}
  ar  = pm(i6,m5); {c msb}
  cntr = ax0;
  do approx_loop until ce;
    mr = mx0 * my1 (us)    , ay0 = pm(i6,m5); {c[i]lsb*xmsb, c[i-1] lsb}
    mr = mr + ar * my0 (su), ay1 = pm(i6,m5); {c[i]msb*xlsb, c[i-1] msb}
    mr0 = mr1;
    mr1 = mr2;                 {shift down by 16 bits}
    mr = mr + ar * my1 (ss);   {c[i]msb*xmsb}
    ar = mr0 + ay0;            {c[i]*x lsb + c[i-1] lsb}
    mx0 = ar;
  approx_loop: ar = mr1 + ay1 + c; {c[i]*x msb + c[i-1] msb}

rts;
     
.endmod;