home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / poly.cal < prev    next >
Encoding:
Text File  |  1995-07-17  |  17.6 KB  |  729 lines
  1. /* 
  2.  * A collection of functions designed for calculations involving
  3.  *     polynomials in one variable (by Ernest W. Bowen).
  4.  *
  5.  * On starting the program the independent variable has identifier x
  6.  *    and name "x", i.e. the user can refer to it as x, the
  7.  *    computer displays it as "x".  The name of the independent
  8.  *    variable is stored as varname, so, for example, varname = "alpha"
  9.  *    will change its name to "alpha".  At any time, the independent
  10.  *    variable has only one name.  For some purposes, a name like
  11.  *    "sin(t)" or "(a + b)" or "\lambda" might be useful;
  12.  *    names like "*" or "-27" are legal but might give expressions
  13.  *    that are difficult to intepret.
  14.  *
  15.  * Polynomial expressions may be constructed from numbers and the
  16.  *    independent variable and other polynomials by the algebraic
  17.  *    operations +, -, *, ^, and if the result is a polynomial /.
  18.  *    The operations // and % are defined to have the quotient and
  19.  *    remainder meanings as usually defined for polynomials.
  20.  *
  21.  * When polynomials are assigned to idenfifiers, it is convenient to
  22.  *    think of the polynomials as values.  For example, p = (x - 1)^2
  23.  *    assigns to p a polynomial value in the same way as q = (7 - 1)^2
  24.  *    would assign to q a number value.  As with number expressions
  25.  *    involving operations, the expression used to define the
  26.  *    polynomial is usually lost; in the above example, the normal
  27.  *    computer display for p will be  x^2 - 2x + 1.  Different
  28.  *    identifiers may of course have the same polynomial value.
  29.  * 
  30.  * The polynomial we think of as a_0 + a_1 * x + ... + a_n * x^n,
  31.  *    for number coefficients a_0, a_1, ... a_n may also be
  32.  *    constructed as pol(a_0, a_1, ..., a_n).  Note that here the
  33.  *    coefficients are to be in ascending power order.  The independent
  34.  *    variable is pol(0,1), so to use t, say, as an identifier for
  35.  *    this, one may assign  t = pol(0,1).  To simultaneously specify
  36.  *    an identifier and a name for the independent variable, there is
  37.  *    the instruction var, used as in identifier = var(name).  For
  38.  *    example, to use "t" in the way "x" is initially, one may give
  39.  *    the instruction  t = var("t").
  40.  *
  41.  * There are four parameters pmode, order, iod and ims for controlling
  42.  *    the format in which polynomials are displayed.
  43.  *    The parameter pmode may have values "alg" or "list": the
  44.  *    former gives a display as an algebraic formula, while
  45.  *    the latter only lists the coefficients.  Whether the terms or
  46.  *    coefficients are in ascending or descending power order is
  47.  *    controlled by order being "up" or "down".  If the
  48.  *    parameter iod (for integer-only display), the polynomial
  49.  *    is expressed in terms of a polynomial whose coefficients are
  50.  *    integers with gcd = 1, the leading coefficient having positive
  51.  *    real part, with where necessary a leading multiplying integer,
  52.  *    a Gaussian integer multiplier if the coefficients are complex
  53.  *    with a common complex factor, and a trailing divisor integer.
  54.  *    If a non-zero value is assigned to the parameter ims,
  55.  *    multiplication signs will be inserted where appropriate;
  56.  *    this may be useful if the expression is to be copied to a
  57.  *    program or a string to be used with eval. 
  58.  *
  59.  * For evaluation of polynomials the standard function is ev(p, t).
  60.  *    If p is a polynomial and t anything for which the relevant
  61.  *    operations can be performed, this returns the value of p
  62.  *    at t.  The function ev(p, t) also accepts lists or matrices
  63.  *    as possible values for p; each element of p is then evaluated
  64.  *    at t.  For other p, t is ignored and the value of p is returned.
  65.  *    If an identifier, a, say, is used for the polynomial, list or
  66.  *    matrix p, the definition
  67.  *            define a(t) = ev(a, t);
  68.  *    permits a(t) to be used for the value of a at t as if the
  69.  *    polynomial, list or matrix were a function.  For example,
  70.  *    if a = 1 + x^2, a(2) will return the value 5, just as if
  71.  *            define a(t) = 1 + t^2;
  72.  *    had been used.   However, when the polynomial definition is
  73.  *    used, changing the polynomial a will change a(t) to the value
  74.  *    of the new polynomial at t.  For example,
  75.  *    after 
  76.  *        L = list(x, x^2, x^3, x^4);
  77.         define a(t) = ev(a, t);
  78.  *    the loop
  79.  *        for (i = 0; i < 4; i++)
  80.  *            print ev(L[[i]], 5);
  81.  *    may be replaced by
  82.  *        for (i = 0; i < 4; i++) {
  83.  *            a = L[[i]];
  84.  *            print a(5);
  85.  *        } 
  86.  *      
  87.  * Matrices with polynomial elements may be added, subtracted and
  88.  *    multiplied as long as the usual rules for compatibility are
  89.  *    observed.  Also, matrices may be multiplied by polynomials,
  90.  *    i.e. if p is a     polynomial and A a matrix whose elements
  91.  *    may be numbers or polynomials, p * A returns the matrix of
  92.  *    the same shape as A with each element multiplied by p.
  93.  *    Square matrices may also be 'substituted for the variable' in
  94.  *    polynomials, e.g. if A is an m x m matrix, and
  95.  *    p = x^2 + 3 * x + 2, ev(p, A) returns the same as
  96.  *    A^2 + 3 * A + 2 * I, where I is the unit m x m matrix.  
  97.  *    
  98.  * On starting this program, three demonstration polynomials a, b, c
  99.  *    have been defined.  The functions a(t), b(t), c(t) corresponding
  100.  *    to a, b, c, and x(t) corresponding to x, have also been
  101.  *    defined, so the usual function notation can be used for
  102.  *    evaluations of a, b, c and x.  For x, as long as x identifies
  103.  *    the independent variable, x(t) should return the value of t,
  104.  *    i.e. it acts as an identity function.
  105.  *    
  106.  * Functions defined include:
  107.  *
  108.  *    monic(a) returns the monic multiple of a, i.e., if a != 0,
  109.  *        the multiple of    a with leading coefficient 1    
  110.  *    conj(a) returns the complex conjugate of a
  111.  *    ispmult(a,b) returns 1 or 0 according as a is or is not
  112.  *        a polynomial multiple of b
  113.  *    pgcd(a,b) returns the monic gcd of a and b 
  114.  *    pfgcd(a,b) returns a list of three polynomials (g, u, v)
  115.  *        where g = pgcd(a,b) and g = u * a + v * b.
  116.  *    plcm(a,b) returns the monic lcm of a and b
  117.  *
  118.  *    interp(X,Y,t) returns the value at t of the polynomial given
  119.  *        by Newtonian divided difference interpolation, where
  120.  *        X is a list of x-values, Y a list of corresponding
  121.  *        y-values.  If t is omitted, the interpolating
  122.  *        polynomial is returned.  A y-value may be replaced by
  123.  *        list (y, y_1, y_2, ...), where y_1, y_2, ... are
  124.  *        the reduced derivatives at the corresponding x;
  125.  *        i.e. y_r is the r-th derivative divided by fact(r).
  126.  *    mdet(A) returns the determinant of the square matrix A,
  127.  *        computed by an algorithm that does not require
  128.  *        inverses;  the built-in det function usually fails
  129.  *        for matrices with polynomial elements.  
  130.  *    D(a,n) returns the n-th derivative of a; if n is omitted,
  131.  *        the first derivative is returned.
  132.  *
  133.  * A first-time user can see what the initially defined polynomials
  134.  *    a, b and c are, and experiment with the algebraic operations
  135.  *    and other functions that have been defined by giving
  136.  *    instructions like:
  137.  *            a
  138.  *            b
  139.  *            c
  140.  *            (x^2 + 1) * a
  141.  *            a^27
  142.  *            a * b    
  143.  *            a % b
  144.  *            a // b
  145.  *            a(1 + x)
  146.  *            a(b)
  147.  *            conj(c)
  148.  *            g = pgcd(a, b)
  149.  *            g
  150.  *            a / g
  151.  *            D(a)
  152.  *            mat A[2,2] = {1 + x, x^2, 3, 4*x}
  153.  *            mdet(A)
  154.  *            D(A)
  155.  *            A^2
  156.  *            define A(t) = ev(A, t)
  157.  *            A(2)
  158.  *            A(1 + x)
  159.  *            define L(t) = ev(L, t)
  160.  *            L = list(x, x^2, x^3, x^4)
  161.  *            L(5)
  162.  *            a(L)
  163.  *            interp(list(0,1,2,3), list(2,3,5,7))
  164.  *            interp(list(0,1,2), list(0,list(1,0),2))
  165.  *
  166.  * One check on some of the functions is provided by the Cayley-Hamilton
  167.  *    theorem:  if A is any m x m matrix and I the m x m unit matrix,
  168.  *    and x is pol(0,1),
  169.  *            ev(mdet(x * I - A), A)
  170.  *    should return the zero m x m matrix.
  171.  */
  172.  
  173. obj poly {p};
  174.  
  175. define pol() {
  176.     local u,i,s;
  177.     obj poly u;
  178.     s = list();
  179.     for (i=1; i<= param(0); i++) append (s,param(i));
  180.     i=size(s) -1;
  181.     while (i>=0 && s[[i]]==0) {i--; remove(s)}
  182.     u.p = s;
  183.     return u;
  184. }
  185.  
  186. define ispoly(a) {
  187.     local y;
  188.     obj poly y;
  189.     return istype(a,y);
  190. }
  191.  
  192. define findlist(a) {
  193.     if (ispoly(a)) return a.p;
  194.     if (a) return list(a);
  195.     return list();
  196. }
  197.  
  198. pmode = "alg";    /* The other acceptable pmode is "list" */
  199. ims = 0;    /* To be non-zero if multiplication signs to be inserted */
  200. iod = 0;    /* To be non-zero for integer-only display */
  201. order = "down"    /* Determines order in which coefficients displayed */
  202.  
  203. define poly_print(a) {
  204.     local f, g, t;
  205.     if (size(a.p) == 0) {
  206.         print 0:;
  207.         return;
  208.     } 
  209.     if (iod) {
  210.         g = gcdcoeffs(a);
  211.         t = a.p[[size(a.p) - 1]] / g;
  212.         if (re(t) < 0) { t = -t; g = -g;}
  213.         if (g != 1) {
  214.             if (!isreal(t)) {
  215.                 if (im(t) > re(t)) g *= 1i;
  216.                 else if (im(t) <= -re(t)) g *= -1i;
  217.             }
  218.             if (isreal(g)) f = g;
  219.             else f = gcd(re(g), im(g));
  220.             if (num(f) != 1) {
  221.                 print num(f):;
  222.                 if (ims) print"*":;
  223.             }
  224.             if (!isreal(g)) {
  225.                 printf("(%d)", g/f);
  226.                 if (ims) print"*":;
  227.             }
  228.             if (pmode == "alg") print"(":;
  229.             polyprint(1/g * a);
  230.             if (pmode == "alg") print")":;
  231.             if (den(f) > 1) print "/":den(f):;
  232.             return;
  233.         }
  234.     }
  235.     polyprint(a);
  236. }
  237.  
  238. define polyprint(a) {
  239.     local s,n,i,c;
  240.     s = a.p;
  241.     n=size(s) - 1;
  242.     if (pmode=="alg") {
  243.         if (order == "up") {
  244.             i = 0;
  245.             while (!s[[i]]) i++;
  246.             pterm (s[[i]], i);
  247.             for (i++ ; i <= n; i++) {
  248.                 c = s[[i]];
  249.                 if (c) {
  250.                     if (isreal(c)) {
  251.                         if (c > 0) print" + ":;
  252.                         else {
  253.                             print" - ":;
  254.                             c = -c;
  255.                         }
  256.                     } 
  257.                     else print " + ":;
  258.                     pterm(c,i);
  259.                 }
  260.             }
  261.             return;
  262.         }
  263.         if (order == "down") {
  264.             pterm(s[[n]],n);
  265.             for (i=n-1; i>=0; i--) {
  266.                 c = s[[i]];
  267.                 if (c) {
  268.                     if (isreal(c)) {
  269.                         if (c > 0) print" + ":;
  270.                         else {
  271.                             print" - ":;
  272.                             c = -c;
  273.                         }
  274.                     } 
  275.                     else print " + ":;
  276.                     pterm(c,i);
  277.                 }
  278.             }
  279.             return;
  280.         }
  281.         quit "order to be up or down";
  282.     }
  283.     if (pmode=="list") {
  284.         plist(s);
  285.         return;
  286.     }
  287.     print pmode,:"is unknown mode";
  288. }
  289.         
  290.  
  291. define poly_neg(a) {
  292.     local s,i,y;
  293.     obj poly y;
  294.     s = a.p;
  295.     for (i=0; i< size(s); i++) s[[i]] = -s[[i]];
  296.     y.p = s;
  297.     return y;
  298. }
  299.  
  300. define poly_conj(a) {
  301.     local s,i,y;
  302.     obj poly y;
  303.     s = a.p;
  304.     for (i=0; i < size(s); i++) s[[i]] = conj(s[[i]]);
  305.     y.p = s;
  306.     return y;
  307. }
  308.  
  309. define poly_inv(a) = pol(1)/a;    /* This exists only for a of zero degree */
  310.  
  311. define poly_add(a,b) {
  312.     local sa, sb, i, y;
  313.     obj poly y;
  314.     sa=findlist(a); sb=findlist(b);
  315.     if (size(sa) > size(sb)) swap(sa,sb);
  316.     for (i=0; i< size(sa); i++) sa[[i]] += sb[[i]];    
  317.     while (i < size(sb)) append (sa, sb[[i++]]);
  318.     while (i > 0 && sa[[--i]]==0) remove (sa);
  319.     y.p = sa;
  320.     return y;
  321. }
  322.  
  323. define poly_sub(a,b) {
  324.      return a + (-b);
  325. }
  326.  
  327. define poly_cmp(a,b) {
  328.     local sa, sb;
  329.     sa = findlist(a);
  330.     sb=findlist(b);
  331.     return  (sa != sb);
  332. }
  333.  
  334. define poly_mul(a,b) {
  335.     local sa,sb,i, j, y;
  336.     if (ismat(a)) swap(a,b);
  337.     if (ismat(b)) {
  338.         y = b;
  339.         for (i=matmin(b,1); i <= matmax(b,1); i++)
  340.             for (j = matmin(b,2); j<= matmax(b,2); j++)
  341.                 y[i,j] = a * b[i,j];
  342.         return y;    
  343.     }
  344.     obj poly y;
  345.     sa=findlist(a); sb=findlist(b);
  346.     y.p = listmul(sa,sb);
  347.     return y;
  348. } 
  349.  
  350. define listmul(a,b) {
  351.     local da,db, s, i, j, u;
  352.     da=size(a)-1; db=size(b)-1;
  353.     s=list();
  354.     if (da >= 0 && db >= 0) {
  355.         for (i=0; i<= da+db; i++) { u=0;
  356.             for (j = max(0,i-db); j <= min(i, da); j++)
  357.             u += a[[j]]*b[[i-j]]; append (s,u);}}
  358.     return s;
  359. }
  360.  
  361. define ev(a,t) {
  362.     local v, i, j;
  363.     if (ismat(a)) {
  364.         v = a;
  365.         for (i = matmin(a,1); i <= matmax(a,1); i++)
  366.             for (j = matmin(a,2); j <= matmax(a,2); j++)
  367.                 v[i,j] = ev(a[i,j], t);
  368.         return v;
  369.     }
  370.     if (islist(a)) {
  371.         v = list();
  372.         for (i = 0; i < size(a); i++)
  373.             append(v, ev(a[[i]], t));
  374.         return v;
  375.     }
  376.     if (!ispoly(a)) return a;
  377.     if (islist(t)) {
  378.         v = list();
  379.         for (i = 0; i < size(t); i++)
  380.             append(v, ev(a, t[[i]]));
  381.         return v;
  382.     }    
  383.     if (ismat(t)) return evpm(a.p, t);
  384.     return evp(a.p, t); 
  385. }
  386.  
  387. define evp(s,t) {
  388.     local n,v,i;
  389.     n = size(s);
  390.     if (!n) return 0;
  391.     v = s[[n-1]];
  392.     for (i = n - 2; i >= 0; i--) v=t * v +s[[i]];
  393.     return v;
  394. }
  395.  
  396. define evpm(s,t) {
  397.     local m, n, V, i, I;
  398.     n = size(s);
  399.     m = matmax(t,1) - matmin(t,1);
  400.     if (matmax(t,2) - matmin(t,2) != m) quit "Non-square matrix";
  401.     mat V[m+1, m+1];
  402.     if (!n) return V;
  403.     mat I[m+1, m+1];
  404.     matfill(I, 0, 1);
  405.     V = s[[n-1]] * I;
  406.     for (i = n - 2; i >= 0; i--) V = t * V + s[[i]] * I;
  407.     return V;
  408. }
  409. pzero = pol(0);
  410. x = pol(0,1); 
  411. varname = "x";
  412. define x(t) = ev(x, t);
  413.  
  414. define iszero(a) {
  415.     if (ispoly(a))
  416.         return !size(a.p);
  417.     return a == 0;
  418. }
  419.  
  420. define isstring(a) = istype(a, " ");
  421.  
  422. define var(name) {
  423.     if (!isstring(name)) quit "Argument of var is to be a string";
  424.     varname = name;
  425.     return pol(0,1);
  426. }
  427.  
  428. define pcoeff(a) {
  429.         if (isreal(a)) print a:;
  430.         else print "(":a:")":;
  431. }
  432.  
  433. define pterm(a,n) {
  434.     if (n==0) {
  435.         pcoeff(a);
  436.         return;
  437.     }
  438.     if (n==1) {
  439.         if (a!=1) {
  440.             pcoeff(a);
  441.             if (ims) print"*":;
  442.         }
  443.         print varname:;
  444.         return;
  445.     }
  446.     if (a!=1) {
  447.         pcoeff(a);
  448.         if (ims) print"*":;
  449.     }
  450.     print varname:"^":n:;
  451. } 
  452.  
  453. define plist(s) {
  454.     local i, n;
  455.     n = size(s);
  456.     print "( ":;
  457.     if (order == "up") {
  458.         for (i=0; i< n-1 ; i++)
  459.             print s[[i]]:",",:;
  460.         if (n) print s[[i]],")":;
  461.         else print "0 )":;
  462.     }
  463.     else {
  464.         if (n) print s[[n-1]]:;
  465.         for (i = n - 2; i >= 0; i--)
  466.             print ", ":s[[i]]:;
  467.         print " )":;
  468.     }
  469. }
  470.  
  471. define deg(a) = size(a.p) - 1;
  472.  
  473. define polydiv(a,b) {
  474.     local q, r, d, u, i, m, n, sa, sb, sq;
  475.     obj poly q, r;
  476.     sa=findlist(a); sb = findlist(b); sq = list();
  477.     m=size(sa)-1; n=size(sb)-1;
  478.     if (n<0) quit "Zero divisor";
  479.     if (m<n) return list(pzero, a);
  480.     d = sb[[n]]; 
  481.     while ( m >= n) { u = sa[[m]]/d;
  482.         for (i = 0; i< n; i++) sa[[m-n+i]] -= u*sb[[i]];
  483.         push(sq,u); remove(sa); m--;
  484.         while (m>=n && sa[[m]]==0) { m--; remove(sa); push(sq,0)}}
  485.     while (m>=0 && sa[[m]]==0) { m--; remove(sa);}
  486.     q.p = sq;  r.p = sa;
  487.     return list(q, r);}
  488.         
  489. define poly_mod(a,b)  {
  490.     local u;
  491.     u=polydiv(a,b);
  492.     return u[[1]];
  493. }
  494.  
  495. define poly_quo(a,b) {
  496.     local p;
  497.     p = polydiv(a,b);
  498.     return p[[0]];
  499. }
  500.  
  501. define ispmult(a,b) = iszero(a % b);
  502.  
  503. define poly_div(a,b) {
  504.     if (!ispmult(a,b)) quit "Result not a polynomial";
  505.     return poly_quo(a,b);
  506. }
  507.  
  508. define pgcd(a,b) {
  509.     local r;
  510.     if (iszero(a) && iszero(b)) return pzero;
  511.     while (!iszero(b)) {
  512.         r = a % b;
  513.         a = b;
  514.         b = r;
  515.     }
  516.     return monic(a);
  517. }
  518.  
  519. define plcm(a,b) = monic( a * b // pgcd(a,b));
  520.   
  521. define pfgcd(a,b) {
  522.     local u, v, u1, v1, s, q, r, d, w;
  523.     u = v1 = pol(1); v = u1 = pol(0);
  524.     while (size(b.p) > 0) {s = polydiv(a,b);
  525.         q = s[[0]];
  526.         a = b; b = s[[1]]; u -= q*u1; v -= -q*v1;
  527.         swap(u,u1); swap(v,v1);}
  528.     d=size(a.p)-1; if (d>=0 && (w= 1/a.p[[d]]) !=1)
  529.          { a *= w; u *= w; v *= w;}
  530.     return list(a,u,v);
  531. }
  532.   
  533. define monic(a) {
  534.     local s, c, i, d, y;
  535.     if (iszero(a)) return pzero;
  536.     obj poly y;
  537.     s = findlist(a);
  538.     d = size(s)-1;
  539.     for (i=0; i<=d; i++) s[[i]] /= s[[d]];
  540.     y.p = s;
  541.     return y;
  542. }
  543.  
  544. define coefficient(a,n) = (n < size(a.p)) ? a.p[[n]] : 0;
  545.  
  546. define D(a, n) {
  547.     local i,j,v;
  548.     if (isnull(n)) n = 1;
  549.     if (!isint(n) || n < 1) quit "Bad order for derivative";
  550.      if (ismat(a)) {
  551.         v = a;
  552.         for (i = matmin(a,1); i <= matmax(a,1); i++)
  553.             for (j = matmin(a,2); j <= matmax(a,2); j++)
  554.                 v[i,j] = D(a[i,j], n);
  555.         return v;
  556.     }
  557.     if (!ispoly(a)) return 0;
  558.     return Dp(a,n);
  559. }
  560.  
  561. define Dp(a,n) {
  562.     local i, v;
  563.     if (n > 1) return Dp(Dp(a, n-1), 1);
  564.      obj poly v;
  565.     v.p=list();    
  566.     for (i=1; i<size(a.p); i++) append (v.p, i*a.p[[i]]);
  567.     return v;
  568. }
  569.  
  570.  
  571. define cgcd(a,b) {
  572.     if (isreal(a) && isreal(b)) return gcd(a,b);
  573.     while (a) {
  574.         b -= bround(b/a) * a;
  575.         swap(a,b);
  576.     }
  577.     if (re(b) < 0) b = -b; 
  578.     if (im(b) > re(b)) b *= -1i;
  579.     else if (im(b) <= -re(b)) b *= 1i;
  580.     return b;
  581. }
  582.  
  583. define gcdcoeffs(a) {
  584.     local s,i,g, c;
  585.     s = a.p;
  586.     g=0;
  587.     for (i=0; i < size(s) && g != 1; i++)
  588.         if (c = s[[i]]) g = cgcd(g, c);
  589.     return g;
  590. }
  591.  
  592. define interp(X, Y, t) = evalfd(makediffs(X,Y), t);
  593.  
  594. define makediffs(X,Y) {
  595.     local U, D, d, x, y, i, j, k, m, n, s;
  596.     U = D = list();
  597.     n = size(X);
  598.     if (size(Y) != n) quit"Arguments to be lists of same size";
  599.     for (i = n-1; i >= 0; i--) {
  600.         x = X[[i]];
  601.         y = Y[[i]];
  602.         m = size(U);
  603.         if (isnum(y)) {
  604.             d = y;
  605.             for (j = 0; j < m; j++) {
  606.                 d = D[[j]] = (D[[j]]-d)/(U[[j]] - x);
  607.             }
  608.             push(U, x);
  609.             push(D, y);
  610.         }
  611.         else {
  612.             s = size(y);
  613.             for (k = 0; k < s ; k++) {
  614.                 d = y[[k]];
  615.                 for (j = 0; j < m; j++) {
  616.                     d = D[[j]] = (D[[j]] - d)/(U[[j]] - x);
  617.                 }
  618.             }
  619.             for (j=s-1; j >=0; j--) {
  620.                 push(U,x);
  621.                 push(D, y[[j]]);
  622.             }
  623.         }
  624.     }
  625.     return list(U, D);
  626. }
  627.     
  628. define evalfd(T, t) {
  629.     local U, D, n, i, v;
  630.     if (isnull(t)) t = pol(0,1);
  631.     U = T[[0]];
  632.     D = T[[1]];
  633.     n = size(U);
  634.     v = D[[n-1]];
  635.     for (i = n-2; i >= 0; i--) 
  636.         v = v * (t - U[[i]]) + D[[i]];
  637.     return v;
  638. }
  639.  
  640.  
  641. define mdet(A) {
  642.     local n, i, j, k, I, J;
  643.     n = matmax(A,1) - (i = matmin(A,1));
  644.     if (matmax(A,2) - (j = matmin(A,2)) != n)
  645.         quit "Non-square matrix for mdet";
  646.     I = J = list();
  647.     k = n + 1;
  648.     while (k--) {
  649.         append(I,i++);
  650.         append(J,j++);
  651.     }
  652.     return M(A, n+1, I, J);
  653. }
  654.  
  655. define M(A, n, I, J) {
  656.     local v, J0, i, j, j1;
  657.     if (n == 1) return A[ I[[0]], J[[0]] ];
  658.     v = 0;
  659.     i = remove(I);
  660.     for (j = 0; j < n; j++) {
  661.         J0 = J;
  662.         j1 = delete(J0, j);
  663.         v += (-1)^(n-1+j) * A[i, j1] * M(A, n-1, I, J0);
  664.     }
  665.     return v;
  666. }
  667.  
  668. define mprint(A) {
  669.     local i,j;
  670.     if (!ismat(A)) quit "Argument to be a matrix";
  671.     for (i = matmin(A,1); i <= matmax(A,1); i++) {
  672.         for (j = matmin(A,2); j <= matmax(A,2); j++)
  673.             printf("%8.4d ", A[i,j]);
  674.         printf("\n");
  675.     }
  676. }
  677.     
  678. obj poly a;
  679. obj poly b;
  680. obj poly c;
  681.  
  682. define a(t) = ev(a,t);
  683. define b(t) = ev(b,t);
  684. define c(t) = ev(c,t);
  685.  
  686. a=pol(1,4,4,2,3,1);
  687. b=pol(5,16,8,1);
  688. c=pol(1+2i,3+4i,5+6i);
  689.  
  690. global lib_debug;
  691. if (lib_debug >= 0) {
  692.     print "obj poly {p} defined";
  693.     print "pol() defined";
  694.     print "poly_print(a) defined";
  695.     print "poly_add(a, b) defined";
  696.     print "poly_sub(a, b) defined";
  697.     print "poly_mul(a, b) defined";
  698.     print "poly_div(a, b) defined";
  699.     print "poly_quo(a,b) defined";
  700.     print "poly_mod(a,b) defined";
  701.     print "poly_neg(a) defined";
  702.     print "poly_conj(a) defined";
  703.     print "poly_cmp(a,b) defined";
  704.     print "iszero(a) defined";
  705.     print "plist(a) defined";
  706.     print "listmul(a,b) defined";
  707.     print "ev(a,t) defined";
  708.     print "evp(s,t) defined";
  709.     print "ispoly(a) defined";
  710.     print "isstring(a) defined";
  711.     print "var(name) defined";
  712.     print "pcoeff(a) defined";
  713.     print "pterm(a,n) defined";
  714.     print "deg(a) defined";
  715.     print "polydiv(a,b) defined";
  716.     print "D(a,n) defined";
  717.     print "Dp(a,n) defined";
  718.     print "pgcd(a,b) defined";
  719.     print "plcm(a,b) defined";
  720.     print "monic(a) defined";
  721.     print "pfgcd(a,b) defined";
  722.     print "interp(X,Y,x) defined";
  723.     print "makediffs(X,Y) defined";
  724.     print "evalfd(T,x) defined";
  725.     print "mdet(A) defined";
  726.     print "M(A,n,I,J) defined";
  727.     print "mprint(A) defined";
  728. }
  729.