%********************************************************************
module reduce_patches$
%********************************************************************
%  Routines for finding leading derivatives and others
%  Author: Tony Hearn, Winfried Neun, Arthur Norman, since 2003

symbolic fluid '(!*fjwflag)$ % for diffp() below

%>>>>>>>> A 2001 patch, now in the development version:

    symbolic procedure noncomp u$
       !*ncmp and noncomp1 u$
 
    symbolic procedure noncomp1 u$
      if null pairp u then nil
       else if pairp car u then noncomfp u
       else flagp(car u,'noncom) or noncomlistp cdr u$
 
    symbolic procedure noncomlistp u$
       pairp u and (noncomp1 car u or noncomlistp cdr u)$
 
%>>>>>>>> The new fix which becomes necessary due to the above 2001 patch:

symbolic procedure mchcomb(u,v,op)$
   begin integer n;
      n := length u - length v +1;
      if n<1 then return nil
       else if n=1 then return mchsarg(u,v,op)
       else if not smemqlp(frlis!*,v) then return nil;
      return for each x in comb(u,n) join
           if null ncmp!* then mchsarg((op . x) . setdiff(u,x),v,op)
            else (if null y then nil
                   else mchsarg((op . x) . car y,
                                if cdr y then reverse v else v,op))
            where y = mchcomb2(x,u,nil,nil,nil)
   end$

symbolic procedure mchcomb2(u,v,w,bool1,bool2)$
   if null u then nconc(reversip w,v) . bool2
    else if car u = car v
           then if noncomp car u then mchcomb2(cdr u,cdr v,w,t,bool2)
                 else mchcomb2(cdr u,cdr v,w,bool1,bool2)
    else if noncomp car u
     then if bool1 then nil
           else mchcomb2(u,cdr v,car v . w,t,if bool2 then bool2 else t)
    else mchcomb2(u,cdr v,car v . w,bool1,bool2)$

%>>>>>>>>>>  the next fixes some noncom bug

symbolic procedure diffp(u,v)$
   begin scalar n,w,x,y,z; integer m;
	n := cdr u;
	u := car u;
	if n>1 and noncomp u
	  then return addsq(multsq(simpdf {u,v},simpexpt {u,n - 1}),
			    multpq(u .** 1,diffp(u . (n - 1),v)))
	 else if u eq v and (w := 1 ./ 1) then go to e
	 else if atom u then go to f
	 else if (not atom car u and (w:= difff(u,v)))
		  or (car u eq '!*sq and (w:= diffsq(cadr u,v)))
	  then go to c
	 else if x := get(car u,'dfform) then return apply3(x,u,v,n)
	 else if x:= get(car u,dfn_prop u) then nil
	 else if car u eq 'plus and (w := diffsq(simp u,v))
	  then go to c
	 else go to h;
	y := x;
	z := cdr u;
    a:  w := diffsq(simp car z,v) . w;
	if caar w and null car y then go to h;
	y := cdr y;
	z := cdr z;
	if z and y then go to a
	 else if z or y then go to h;
	y := reverse w;
	z := cdr u;
	w := nil ./ 1;
    b:  if caar y
	  then w := addsq(multsq(car y,simp subla(pair(caar x,z),
						   cdar x)),
			  w);
	x := cdr x;
	y := cdr y;
	if y then go to b;
    c:
    e: if (x := atsoc(u,wtl!*))
	  then w := multpq('k!* .** (-cdr x),w);
	m := n-1;
	return rationalizesq if n=1 then w
		else if flagp(dmode!*,'convert)
		     and null(n := int!-equiv!-chk
					   apply1(get(dmode!*,'i2d),n))
		 then nil ./ 1
		else multsq(!*t2q((u .** m) .* n),w);
    f:  if not depends(u,v)
	   and (not (x:= atsoc(u,powlis!*))
		 or not depends(cadddr x,v))
	   and null !*depend
	  then return nil ./ 1;
	w := list('df,u,v);
	w := if x := opmtch w then simp x else mksq(w,1);
	go to e;
    h:  if car u eq 'df
	  then if depends(cadr u,v)
			and not(cadr u eq v)
		 then
	  <<if !*fjwflag and eqcar(cadr u, 'int) then
	      if caddr cadr u eq v then
		 << w := 'df . cadr cadr u . cddr u;  go to j >>
	      else if !*allowdfint and
		 not_df_p(w := diffsq(simp!* cadr cadr u, v))
	      then <<
		 w := 'df . reval{'int, mk!*sq w, caddr cadr u} . cddr u;
		 go to j >>;
	   if (x := find_sub_df(w:= cadr u . derad(v,cddr u),
					   get('df,'kvalue)))
			  then <<w := simp car x;
				 for each el in cdr x do
				    for i := 1:cdr el do
					w := diffsq(w,car el);
				 go to e>>
		       else w := 'df . w
		>>
		else if null !*depend then return nil ./ 1
		else w := {'df,u,v}
	 else w := {'df,u,v};
   j:   if (x := opmtch w) then w := simp x
	 else if not depends(u,v) and null !*depend then return nil ./ 1
	 else w := mksq(w,1);
      go to e
   end$


%>>>>>>>> A fix by Tony Hearn, 17. Sep 2003

symbolic procedure lterm(u,kern)$
   begin scalar x,y;
      u := simp!* u;
      y := denr u;
      tstpolyarg(y,u);
      u := numr u;
      kern := !*a2k kern;
      if domainp u then return if null u then 0 else !*ff2a(u,y)
       else if mvar u eq kern then return !*ff2a(lt u .+ nil,y);
      x := updkorder kern;
      u := reorder u;
      if mvar u eq kern then u := lt u .+ nil;
      setkorder x;
      u := reorder u;
      return !*ff2a(u,y)
   end$

%>>>>>>>> A factorizer fix by Arthur Norman, 7. Dec 2003

symbolic fluid '(alphalist hensel!-growth!-size number!-of!-factors)$

symbolic procedure divide!-all!-alphas n$
% Multiply the factors by n mod p and alter the alphas accordingly.
  begin scalar om,m,nn;
    om:=set!-modulus hensel!-growth!-size;
    nn:=modular!-number n;
    m:=modular!-expt(
          modular!-reciprocal nn,
          number!-of!-factors #- 1);
    alphalist:=for each a in alphalist collect
      (times!-mod!-p(nn,car a) . times!-mod!-p(m,cdr a));
    set!-modulus om
  end$

%>>>>>>>> A gcd catastrophic error fix by Winfried Neun, 13. Sep 2006

symbolic procedure gcdf(u,v)$
%  U and V are standard forms.
%  Value is the gcd of U and V, complete only if *GCD is true.
begin scalar !*exp,!*rounded; 
%  The next line was to prevent numerators moving to denominators
%  as in weight x=1,y=2$ wtlevel 4$ wtest:=(z^4-z^3*y-z^3*x+z^2*y^2
%  +2*z^2*x*y+z^2*x^2-3*z*x^2*y-z*x^3+x^4)/z^5; wtest where z=>a;
%  However, the results are formally correct without it, and it
%  causes other problems.
%  if wtl!* then return 1;
!*exp := t;
u := if domainp u or domainp v or not !*ezgcd
%       or dmode!* memq '(!:rn!: !:rd!:)  % Should be generalized.
or dmode!* % I don't know what to do in this case.
or free!-powerp u or free!-powerp v
then gcdf1(u,v)
     else ezgcdf(u,v);
  return if minusf u then negf u else u
end$

%>>>>>>>> A fix by Winfried Neun, 20. Sep 2006 
%>>>>>>>> to simplify df(int(f,y),x,y) --> df(f,x)
%>>>>>>>> An alternative is
%>>>>>>>> on allowdfint$
%>>>>>>>> on dfint$

copyd('oldsimpdf,'simpdf)$

symbolic procedure simpdf (li)$
begin scalar intvar,intvar2,vars,restvars;

 if (pairp car li) and (caar li = 'int)
   then << intvar := caddar li;
           vars := append(li,nil)>> % will be destroyed
   else return oldsimpdf(li);

 if idp intvar and ( restvars := mymemq(intvar,vars,nil) )
   then
     if (pairp cdr restvars) and numberp (intvar2 := cadr restvars)
      then  << vars := car vars .
                 (intvar . ( intvar2 . append( cdr vars, cddr
   restvars)));
               return oldsimpdf (vars) >>
     else  <<  vars := car vars . (intvar . append( cdr vars, cdr
   restvars));
               return oldsimpdf (vars) >>;

 return oldsimpdf(li);
end$

symbolic procedure mymemq (u , v, v1)$
% EQ version of Member
% hard truncating the list in front of the item found
if not pairp v  then nil
   else if eq( u ,car v) then << if v1 then rplacd(v1,nil) ; v>>
   else mymemq(u ,cdr v, v)$

%>>>>>>>>>> The normal REDUCE algebraic mode function cons
% converts standard quotient lists into prefix form which to
% convert back into standard quotient form would take very very long
% for large expressions. The new function sqcons returns
% a list of standard quotient expressions.

symbolic procedure sq!*cons(x)$
<< 'list .  cons (aeval car x, cdr aeval cadr x)>>$

put('sqcons,'psopfn,'sq!*cons)$

%--------
symbolic procedure sq!*rest(x)$
<< 'list . cddr aeval car x>>$

put('sqrest,'psopfn,'sq!*rest)$

%--------
symbolic procedure sq!*first(x)$
cadr aeval car x$

put('sqfirst,'psopfn,'sq!*first)$

%--------
symbolic procedure sq!*second(x)$
caddr aeval car x$

put('sqsecond,'psopfn,'sq!*second)$

%--------
symbolic procedure sq!*third(x)$
cadddr aeval car x$

put('sqthird,'psopfn,'sq!*third)$

%--------
symbolic procedure sq!*part(x)$
% This procedure is only equivalent to part(a,b) if the first
% argument to sqpart is an algebraic list and it the second
% argument is not 0.
begin scalar c1,c2$
 c1:=aeval car x$
 c2:=aeval cadr x$
 return
 if (c2=0) and not pairp c1 then -1 
                            else nth(c1,add1 c2)$
end$

put('sqpart,'psopfn,'sq!*part)$

%--------
symbolic procedure sq!*reverse(x)$
<< 'list . reverse cdr aeval car x>>$

put('sqreverse,'psopfn,'sq!*reverse)$

%--------
symbolic procedure sq!*append(x)$
<< 'list . append(cdr aeval car x,cdr aeval cadr x)>>$

put('sqappend,'psopfn,'sq!*append)$


%>>>>>>>>>> To avoid printing warnings of the compiler that
% PSL constants are non-local variables

symbolic procedure symbid(u,vars)$
<<if fname!* and
     null(ftype!* memq '(macro smacro)) and
     not(atsoc(u,vars) or fluidp u or globalp u or null u or u eq t or
         flagp(u,'share) or !*comp or get(u,'constant!?))
  then lprim list("nonlocal use of undeclared variable",u,
                  "in procedure",fname!*);
  u
>>$


%>>>>>>>>>> To avoid "***** READ Buffer overflow, Truncating"
% when reading from a file too long numbers or strings
% One needs an extra modified REDUCE and then also

% If the corresponding patch has not yet been used and REDUCE been
% re-compiled then uncomment the following:

%load token!-decls$

%fluid '(bigtokenbuffersize)$
%lisp <<tokenbuffer := gtwarray(25000)$
%       bigtokenbuffersize := 200000; >>$

%symbolic procedure readinbuf ()$
%<<tokch:= channelreadchar tokchannel;
%  putstrbyt(inf tokenbuffer, currentchar, tokch);
%  chtokentype := tokentypeofchar tokch;
%  if wlessp(currentchar,bigtokenbuffersize) then
%  currentchar := currentchar #+ 1           else
%  if currentchar #= bigtokenbuffersize then <<
%    errorprintf( "***** READ Buffer overflow, Truncating");
%    currentchar:=bigtokenbuffersize #+ 1
%  >>                                   else
%  currentchar :=  bigtokenbuffersize #+ 1
%>>$

%>>>>>>>>>> To avoid a crash of gcdlist involving rational numbers
% 4 Oct 2007 by Winfried Neun

symbolic procedure gcdlist_aux(l);
if null dmode!* then gcdlist l
else if length l = 1 then car l
else if car l = 1 then 1
else if length l = 2 then  gcdf(car l, cadr l)
else gcdlist_aux (gcdf(car l, cadr l) . cddr l);

%>>>>>>>>>> To speed up substitutions
% by Tony Hearn

symbolic procedure equalreval u;
  (if !*evallhseqp or not atom car u and flagp(caar u,'immediate)
     then list('equal,reval car u,x)
    else list('equal,car u,x))
   where x= (if eqcar(y,'!*sq) then aeval y else reval y)        % <-- was: where x = reval y
             where y=cadr u;
 
symbolic procedure subeval0 u;
   % This is the general evaluator for SUB forms.  All but the last
   % argument are assumed to be substitutions.  These can either be
   % an explicit rule with a lhs and rhs separated by an equal sign,
   % a list of such rules, or something that evaluates to this.
   begin scalar x,y,z,ns,caddrx;                                  % <-- caddrx added
   % Check for spurious substitutions.
   while cdr u do <<if not eqcar(car u,'equal) then x := car u . x
                     else if not(cadar u = (y := aeval caddar u)) % <-- was reval istead of aeval
                      then x := {caar u,cadar u,y} . x;
                    u := cdr u>>;
   if null x then return car u else u := nconc(reversip x,u);
   % Separate assignments from expression.
      if u member sublist!* then return mk!*sq !*p2q mksp('sub . u,1)
       else sublist!* := u . sublist!*;
      if null(u and cdr u)
       then rederr "SUB requires at least 2 arguments"; % F.J. Wright.
      (while cdr u do
          <<x := reval car u;
            if getrtype x eq 'list then u := append(cdr x,cdr u)
             else <<if not eqexpr x then errpri2(car u,t);
                    y := cadr x;
                    if null getrtype y then y := !*a2kwoweight y;
                    caddrx:=reval caddr x;                          % <-- inserted
                    if getrtype caddrx then ns := (y . caddrx) . ns % <-- now 2x with caddrx
                     else z := (y . caddrx) . z;                    % <-- now with caddrx
                    u := cdr u>>>>) where !*evallhseqp=nil;
      x := aeval car u;
      return subeval1(append(ns,z),x)
end;

%>>>>>>>> Fix of an infinite substitution loop generated in odesolve()
%         when sin^2 => 1-cos^2 or cos^2=>1-sin^2 was issued.

load_package trigsimp$ % load_package instead of load to load recursively all

symbolic procedure get_trig_arguments(term, args)$
%% Return a SET of all the arguments of the trig functions in the
%% expression.  (Note that trig functions are unary!)  The
%% arguments may themselves be general expressions -- they need not
%% be kernels!
if atom term then args else
begin scalar f, r;
f := car term;                    % function or operator
if f = '!*sq then << term := reval term; f := car term>>;
          % because SQs in equations may not be transformed to prefix
r := cdr term;                    % arguments or operands
if f memq '(sin cos sinh cosh) then return
  if not member(r := car r, args) then r . args else args;
for each j in r do args := get_trig_arguments(j, args);
return args
end$

%>>>>>>>> Fix of error message ".. invalid as distributive polynomial exponent"
% 3 June 2008 by Eberhard Schruefer 

% The extension of a2dip allows exponentials in *non*-distributive variables.
% Remaining problem: if a variable is u^p then in the expression
% a3*u^(3*p)+ a1*u^p the power u^(3*p) is not represented as (u^p)^3 .
% So exponents are still to be avoided in crineq.red .

symbolic procedure a2dip u;
%     Converts the algebraic (prefix) form u to a distributive poly.
%     We assume that all variables used have been previously
%     defined in dipvars!*, but a check is also made for this
if atom u then a2dipatom u
    else if not atom car u or not idp car u
     then typerr(car u,"dipoly operator")
       % Handling expt separately because the exponents should
       % not be simplified as domain elements.
    else if car u='expt
            then if vevzero!? car a2dip cadr u and vevzero!? car a2dip caddr u
                    then dipfmon(simp!* u,evzero())
                  else dipfnpow(a2dip cadr u,caddr u)
    else (if x then apply(x,list for each y in cdr u collect a2dip y)
         else a2dipatom u)
          where x=get(car u,'dipfn)$

%>>>>>>>> Fix of error message "***** Non-numeric argument in arithmetic"
% resulting under on complex from factorization but being caused by the
% procedure solvequartic
% 22 July 2008 by Eberhard Schruefer

symbolic procedure solvequartic(a4,a3,a2,a1,a0)$
   % Solve the quartic equation a4*x**4+a3*x**3+a2*x**2+a1*x+a0 = 0,
   % where the ai are standard quotients, using technique described in
   % Section 3.8.3 of Abramowitz and Stegun;
   begin scalar x,y,yy,cx,z,s,l,zz1,zz2,dmode,neg,!*msg,!*numval,
                a1cr,a2cr,a3cr,xcr,ycr,zcr;
      % Convert equation to monomial form.
      dmode := dmode!*;
      a3 := quotsq(a3,a4);
      a2 := quotsq(a2,a4);
      a1 := quotsq(a1,a4);
      a0 := quotsq(a0,a4);
      % Build and solve the resultant cubic equation.  We select the
      % real root if there is only one; or if there are three, we choose
      % one that yields real coefficients for the quadratics.  If no
      % roots are known to be real, we use an arbitrary one.
      yy := subtrsq(exptsq(a3,2),multfq(4,a2));
      x := solvecubic(!*f2q 1,
                      negsq a2,
                      subs2!* subtrsq(multsq(a1,a3),multfq(4,a0)),
                      subs2!* negsq addsq(exptsq(a1,2),
                                          multsq(a0,yy)));
      cx := car x;
      % Now check for real roots of the cubic.
      for each rr in x do if one_real rr then s := append(s,list rr);
      x := if (l := length s)=1 then car s else cx;
      % Now solve the two equivalent quadratic equations.
      a3 := quotsqf(a3,2); yy := quotsqf(yy,4);
      % select real coefficient for quadratic if possible.
      y := addsq(yy,x);
      if l<2 then go to zz;
loop: if not pos_num negsq y then go to zz else if l=1 then
        <<x := cx; y := addsq(yy,x); go to zz>>;
      l := l-1; s := cdr s; x := car s;
      y := addsq(yy,x); go to loop;
  zz: y := sqrtq y;
      x := quotsqf(x,2);
      z := sqrtq subtrsq(exptsq(x,2),a0);
     % the following test is needed, according to some editions of
     % Abramowitz and Stegun, to select the correct signs
     % (for the terms z) in the quadratics to produce correct roots.
     % Unfortunately, this test may fail for coefficients which are not
     % numeric because of the inability to recognize zero.
      !*numval := t;
      on rounded,complex;
      a1cr := resimp a1; a2cr := resimp a2; a3cr := resimp a3;
      xcr := resimp x; ycr := resimp y; zcr := resimp z;
      if null numr
         (zz1 :=
           resimp subtrsq(a1cr,addsq(multsq(subtrsq(a3cr,ycr),addsq(xcr,zcr)),
                 multsq(addsq(a3cr,ycr),subtrsq(xcr,zcr))))) then go to rst;
      if null numr
         (zz2 :=
           resimp subtrsq(a1cr,addsq(multsq(subtrsq(a3cr,ycr),subtrsq(xcr,zcr)),
                 multsq(addsq(a3cr,ycr),addsq(xcr,zcr)))))
         then <<neg := t; go to rst>>;
      if domainp numr zz1 and domainp numr zz2
             and numberp denr zz1 and numberp denr zz2 and
         numr simp list('sign,list('difference,list('norm,mk!*sq zz1),
           list('norm,mk!*sq zz2)))=1 then neg := t;
 rst: off rounded,complex;
      if dmode then onoff(get(dmode,'dname),t);
      if neg then z := negsq z;
      return append(solvequadratic(!*f2q 1,subtrsq(a3,y),subtrsq(x,z)),
                    solvequadratic(!*f2q 1,addsq(a3,y),addsq(x,z)))
   end$

endmodule$

end$