Check-ins Using Download Hex

Artifact bb90f8c0a7558186708b2d5fd0ba5d41700e76728ec6861c54ffe03a7a6b4662:

module interfac;

% Authors: A. C. Norman and P. M. A. Moore, 1981.

% Modifications by: Anthony C. Hearn.

fluid '(m!-image!-variable
        poly!-vector
        polyzero
        unknowns!-list
        varlist);


%**********************************************************************;
%
% Routines that are specific to REDUCE.
%  These are either routines that are not needed in the HASH system
%  (which is the other algebra system that this factorizer
%  can be plugged into) or routines that are specifically
%  redefined in the HASH system.


%---------------------------------------------------------------------;
% The following would normally live in section:  ALPHAS
%---------------------------------------------------------------------;

symbolic procedure assoc!-alpha(poly,alist);  assoc(poly,alist);


%---------------------------------------------------------------------;
% The following would normally live in section:  COEFFTS
%---------------------------------------------------------------------;

symbolic procedure termvector2sf v;
  begin scalar r,w;
    for i:=car getv(v,0) step -1 until 1 do <<
      w:=getv(v,i);
            % degree . coefft;
      r:=if car w=0 then cdr w else
        (mksp(m!-image!-variable,car w) .* cdr w) .+ r
    >>;
    return r
  end;

symbolic procedure force!-lc(a,n);
% force polynomial a to have leading coefficient as specified;
    (lpow a .* n) .+ red a;

symbolic procedure merge!-terms(u,v);
  merge!-terms1(1,u,v,car getv(v,0));

symbolic procedure merge!-terms1(i,u,v,n);
  if i#>n then u
  else begin scalar a,b;
    a:=getv(v,i);
    if domainp u or not(mvar u=m!-image!-variable) then
      if not(car a=0) then errorf list("MERGING COEFFTS FAILED",u,a)
      else if cdr a then return cdr a
      else return u;
    b:=lt u;
    if tdeg b=car a then return
      (if cdr a then tpow b .* cdr a else b) .+
        merge!-terms1(i #+ 1,red u,v,n)
    else if tdeg b #> car a then return b .+ merge!-terms1(i,red u,v,n)
    else errorf list("MERGING COEFFTS FAILED ",u,a)
  end;

symbolic procedure list!-terms!-in!-factor u;
% ...;
  if domainp u then list (0 . nil)
  else (ldeg u . nil) . list!-terms!-in!-factor red u;

symbolic procedure try!-other!-coeffts(r,unknowns!-list,uv);
  begin scalar ldeg!-r,lc!-r,w;
    while not domainp r and (r:=red r) and not(w='complete) do <<
      if not depends!-on!-var(r,m!-image!-variable) then
        << ldeg!-r:=0; lc!-r:=r >>
      else << ldeg!-r:=ldeg r; lc!-r:=lc r >>;
      w:=solve!-next!-coefft(ldeg!-r,lc!-r,unknowns!-list,uv) >>
  end;


%---------------------------------------------------------------------;
% The following would normally live in section:  FACMISC
%---------------------------------------------------------------------;

symbolic procedure derivative!-wrt!-main!-variable(p,var);
% partial derivative of the polynomial p with respect to
% its main variable, var;
    if domainp p or (mvar p neq var) then nil
    else
     begin
      scalar degree;
      degree:=ldeg p;
      if degree=1 then return lc p; %degree one term is special;
      return (mksp(mvar p,degree-1) .* multf(degree,lc p)) .+
        derivative!-wrt!-main!-variable(red p,var)
     end;

symbolic procedure fac!-univariatep u;
% tests to see if u is univariate;
  domainp u or not multivariatep(u,mvar u);

symbolic procedure variables!.in!.form(a,sofar);
    if domainp a then sofar
    else <<
      if not memq(mvar a,sofar) then
        sofar:=mvar a . sofar;
      variables!.in!.form(red a,
        variables!.in!.form(lc a,sofar)) >>;


symbolic procedure degree!-in!-variable(p,v);
% returns the degree of the polynomial p in the
% variable v;
    if domainp p then 0
    else if lc p=0
     then errorf "Polynomial with a zero coefficient found"
    else if v=mvar p then ldeg p
    else max(degree!-in!-variable(lc p,v),
      degree!-in!-variable(red p,v));

symbolic procedure get!-height poly;
% find height (max coefft) of given poly;
  if null poly then 0
  else if numberp poly then abs poly
  else max(get!-height lc poly,get!-height red poly);


symbolic procedure poly!-minusp a;
    if a=nil then nil
    else if domainp a then minusp a
    else poly!-minusp lc a;

symbolic procedure poly!-abs a;
    if poly!-minusp a then negf a
    else a;

symbolic procedure fac!-printfactors l;
% procedure to print the result of factorize!-form;
% ie. l is of the form: (c . f)
%  where c is the numeric content (may be 1)
%  and f is of the form: ( (f1 . e1) (f2 . e2) ... (fn . en) )
%    where the fi's are s.f.s and ei's are numbers;
<< terpri();
  if not (car l = 1) then printsf car l;
  for each item in cdr l do
    printsf !*p2f mksp(prepf car item,cdr item) >>;


%---------------------------------------------------------------------;
% The following would normally live in section:  FACPRIM
%---------------------------------------------------------------------;

symbolic procedure invert!.poly(u,var);
% u is a non-trivial primitive square free multivariate polynomial.
% assuming var is the top-level variable in u, this effectively
% reverses the position of the coeffts: ie
%   a(n)*var**n + a(n-1)*var**(n-1) + ... + a(0)
% becomes:
%   a(0)*var**n + a(1)*var**(n-1) + ... + a(n) .               ;
  begin scalar w,invert!-sign;
    w:=invert!.poly1(red u,ldeg u,lc u,var);
    if poly!-minusp lc w then <<
      w:=negf w;
      invert!-sign:=-1 >>
    else invert!-sign:=1;
    return invert!-sign . w
  end;

symbolic procedure invert!.poly1(u,d,v,var);
% d is the degree of the poly we wish to invert.
% assume d > ldeg u always, and that v is never nil;
  if (domainp u) or not (mvar u=var) then
    (var to d) .* u .+ v
  else invert!.poly1(red u,d,(var to (d-ldeg u)) .* (lc u) .+ v,var);


symbolic procedure trailing!.coefft(u,var);
% u is multivariate poly with var as the top-level variable. we find
% the trailing coefft - ie the constant wrt var in u;
  if domainp u then u
  else if mvar u=var then trailing!.coefft(red u,var)
  else u;


%---------------------------------------------------------------------;
% The following would normally live in section:  IMAGESET
%---------------------------------------------------------------------;

symbolic procedure make!-image!-lc!-list(u,imset);
  reversip make!-image!-lc!-list1(u,imset,
    for each x in imset collect car x);

symbolic procedure make!-image!-lc!-list1(u,imset,varlist);
% If IMSET=((x1 . a1, x2 . a2, ... , xn . an)) (ordered) where xj is
% the variable and aj its value, then this fn creates n images of U wrt
% sets S(i) where S(i)= ((x1 . a1), ... , (xi . ai)). The result is an
% ordered list of pairs: (u(i) . X(i+1)) where u(i)= U wrt S(i) and
% X(i) = (xi, ... , xn) and X(n+1) = NIL.  VARLIST = X(1).
% (Note. the variables tagged to u(i) should be all those
% appearing in u(i) unless it is degenerate). The returned list is
% ordered with u(1) first and ending with the number u(n);
  if null imset then nil
  else if domainp u then list(!*d2n u . cdr varlist)
  else if mvar u=caar imset then
    begin scalar w;
      w:=horner!-rule!-for!-one!-var(
        u,caar imset,cdar imset,polyzero,ldeg u) . cdr varlist;
      return
        if polyzerop car w then list (0 . cdr w)
        else (w . make!-image!-lc!-list1(car w,cdr imset,cdr varlist))
    end
  else make!-image!-lc!-list1(u,cdr imset,cdr varlist);

symbolic procedure horner!-rule!-for!-one!-var(u,x,val,c,degg);
  if domainp u or not(mvar u=x)
    then if zerop val then u else addf(u,multf(c,!*n2f(val**degg)))
  else begin scalar newdeg;
    newdeg:=ldeg u;
    return horner!-rule!-for!-one!-var(red u,x,val,
       if zerop val then lc u
        else addf(lc u,
                  multf(c,!*n2f(val**(idifference(degg,newdeg))))),
                            newdeg)
  end;

symbolic procedure make!-image(u,imset);
% finds image of u wrt image set, imset, (=association list);
  if domainp u then u
  else if mvar u=m!-image!-variable then
    adjoin!-term(lpow u,!*n2f evaluate!-in!-order(lc u,imset),
                        make!-image(red u,imset))
  else !*n2f evaluate!-in!-order(u,imset);

symbolic procedure evaluate!-in!-order(u,imset);
% makes an image of u wrt imageset, imset, using horner's rule. result
% should be purely numeric;
  if domainp u then !*d2n u
  else if mvar u=caar imset then
    horner!-rule(evaluate!-in!-order(lc u,cdr imset),
      ldeg u,red u,imset)
  else evaluate!-in!-order(u,cdr imset);

symbolic procedure horner!-rule(c,degg,a,vset);
% c is running total and a is what is left;
  if domainp a
    then if zerop cdar vset then !*d2n a
          else (!*d2n a)+c*((cdar vset)**degg)
  else if not(mvar a=caar vset)
   then if zerop cdar vset then evaluate!-in!-order(a,cdr vset)
         else evaluate!-in!-order(a,cdr vset)+c*((cdar vset)**degg)
  else begin scalar newdeg;
    newdeg:=ldeg a;
    return horner!-rule(if zerop cdar vset
                          then evaluate!-in!-order(lc a,cdr vset)
                         else evaluate!-in!-order(lc a,cdr vset)
      +c*((cdar vset)**(idifference(degg,newdeg))),newdeg,red a,vset)
  end;


%---------------------------------------------------------------------;
% The following would normally live in section:  MHENSFNS
%---------------------------------------------------------------------;

symbolic procedure max!-degree(u,n);
% finds maximum degree of any single variable in U (n is max so far);
  if domainp u then n
  else if igreaterp(n,ldeg u) then
    max!-degree(red u,max!-degree(lc u,n))
  else max!-degree(red u,max!-degree(lc u,ldeg u));

symbolic procedure diff!-over!-k!-mod!-p(u,k,v);
% derivative of u wrt v divided by k (=number);
  if domainp u then nil
  else if mvar u = v then
    if ldeg u = 1 then quotient!-mod!-p(lc u,modular!-number k)
    else adjoin!-term(mksp(v,isub1 ldeg u),
      quotient!-mod!-p(
        times!-mod!-p(modular!-number ldeg u,lc u),
        modular!-number k),
      diff!-over!-k!-mod!-p(red u,k,v))
  else adjoin!-term(lpow u,
    diff!-over!-k!-mod!-p(lc u,k,v),
    diff!-over!-k!-mod!-p(red u,k,v));

symbolic procedure diff!-k!-times!-mod!-p(u,k,v);
% differentiates u k times wrt v and divides by (k!) ie. for each term
% a*v**n we get [n k]*a*v**(n-k) if n>=k and nil if n<k where
% [n k] is the binomial coefficient;
  if domainp u then nil
  else if mvar u = v then
    if ldeg u < k then nil
    else if ldeg u = k then lc u
    else adjoin!-term(mksp(v,ldeg u - k),
      times!-mod!-p(binomial!-coefft!-mod!-p(ldeg u,k),lc u),
      diff!-k!-times!-mod!-p(red u,k,v))
  else adjoin!-term(lpow u,
    diff!-k!-times!-mod!-p(lc u,k,v),
    diff!-k!-times!-mod!-p(red u,k,v));

symbolic procedure spreadvar(u,v,slist);
% find all the powers of V in U and merge their degrees into SLIST.
% We ignore the constant term wrt V;
  if domainp u then slist
  else <<
    if mvar u=v and not member(ldeg u,slist) then slist:=ldeg u . slist;
    spreadvar(red u,v,spreadvar(lc u,v,slist)) >>;


%---------------------------------------------------------------------;
% The following would normally live in section:  UNIHENS
%---------------------------------------------------------------------;

symbolic procedure root!-squares(u,sofar);
  if null u then pmam!-sqrt sofar
  else if domainp u then pmam!-sqrt(sofar+(u*u))
  else root!-squares(red u,sofar+(lc u * lc u));

%---------------------------------------------------------------------;
% The following would normally live in section:  VECPOLY
%---------------------------------------------------------------------;

symbolic procedure poly!-to!-vector p;
% spread the given univariate polynomial out into POLY-VECTOR;
    if domainp p then putv(poly!-vector,0,!*d2n p)
    else <<
      putv(poly!-vector,ldeg p,lc p);
      poly!-to!-vector red p >>;

symbolic procedure vector!-to!-poly(p,d,v);
% Convert the vector P into a polynomial of degree D in variable V;
  begin
    scalar r;
    if d#<0 then return nil;
    r:=!*n2f getv(p,0);
    for i:=1:d do
      if getv(p,i) neq 0 then r:=((v to i) .* getv(p,i)) .+ r;
    return r
  end;



endmodule;


end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]