1module trigsmp2$  % TrigSimp executable code
2
3% small revision by Winfried Neun 3. Nov. 2008
4% need to take care of the dependencies (depl*) in case a df is in the form,
5% e.g. trigsimp(cos((f1 - f2)/4)**4*df(f1,x,y),sin);
6
7% Revised by Francis J. Wright <f.j.wright@maths.qmul.ac.uk>
8% Revision Time-stamp: <FJW, 07 November 2008>
9
10% (FJW) To do:
11%   check with non-integer number domains
12
13% These variables control rules in trigsmp1:
14share hyp_preference, trig_preference$
15
16fluid '(!*complex dmode!*)$
17
18
19%%%%%%%%%%%%
20% TrigSimp %
21%%%%%%%%%%%%
22
23fluid '(depl!*);
24
25symbolic procedure trigsimp!*(u);
26   (trigsimp(reval car u, revlis cdr u) where depl!* = depl!*);
27
28put('trigsimp, 'psopfn, 'trigsimp!*)$
29
30%% FJW: trigsimp is defined to autoload
31
32symbolic procedure trigsimp(f, options);
33   %% Map trigsimp1 over possible structures:
34   if atom f then f                     % nothing to simplify!
35   else if car f eq 'equal then         % equation
36      'equal . for each ff in cdr f collect trigsimp(ff, options)
37   else if car f eq 'list then          % list
38      'list . for each ff in cdr f collect trigsimp(ff, options)
39   else if car f eq 'mat then           % matrix
40      'mat . for each ff in cdr f collect
41         for each fff in ff collect trigsimp(fff, options)
42   else trigsimp1(f, options);          % scalar
43
44
45symbolic procedure trigsimp1(f, options);
46   %% The main TrigSimp driver.
47   begin scalar dname, trigpreference, hyppreference,
48         tanpreference, tanhpreference,
49         direction, mode, keepalltrig, onlytan, opt_args;
50
51      onlytan := not or(smember('sin,f), smember('cos,f),
52         smember('sinh,f), smember('cosh,f),
53         smember('csc,f), smember('sec,f),
54         smember('csch,f), smember('sech,f));
55
56      %% Return quickly if simplification not appropriate:
57      if onlytan and not or(smember('tan,f), smember('cot,f),
58         smember('tanh,f), smember('coth,f),
59         smember('exp,f), smember('e,f)) then
60         return f;
61
62      if (dname := get(dmode!*, 'dname)) then <<
63         %% Force integer domain mode:
64         off dname;
65         f := prepsq simp!* f
66      >>;
67
68      %% Process optional arguments:
69      for each u in options do
70         if u memq '(sin cos) then
71            if trigpreference then
72               (u eq trigpreference) or
73                  rederr "Incompatible options: use either sin or cos."
74            else trigpreference := u
75         else if u memq '(sinh cosh) then
76            if hyppreference then
77               (u eq hyppreference) or
78                  rederr "Incompatible options: use either sinh or cosh."
79            else hyppreference := u
80         else if u eq 'tan then tanpreference := t
81         else if u eq 'tanh then tanhpreference := t
82         else if u memq '(expand combine compact) then
83            if direction then
84               (u eq direction) or
85                  rederr "Incompatible options: use either expand or combine or compact."
86            else direction := u
87         else if u memq '(hyp trig expon) then
88            if mode then
89               (u eq mode) or
90                  rederr "Incompatible options: use either hyp or trig or expon."
91            else mode := u
92         else if u eq 'keepalltrig then keepalltrig := t
93         else if eqcar(u, 'quotient) and not(u member opt_args) then
94            %% optional trig arg of the form `x/2'
95            opt_args := u . opt_args
96         else
97            rederr {"Option", u, "invalid.", " Allowed options are",
98               "sin or cos, tan, cosh or sinh, tanh,",
99               "expand or combine or compact,",
100               "hyp or trig or expon, keepalltrig."};
101
102      %% Set defaults and globals:
103      if trigpreference then
104         (if tanpreference then         % reverse trig preference
105            trigpreference := if trigpreference eq 'sin then 'cos else 'sin)
106      else
107         trigpreference := 'sin;
108      trig_preference := trigpreference;
109
110      if hyppreference then
111         (if tanhpreference then        % reverse hyp preference
112            hyppreference := if hyppreference eq 'sinh then 'cosh else 'sinh)
113      else
114         hyppreference := 'sinh;
115      hyp_preference := hyppreference;
116
117      direction or (direction := 'expand);
118
119      %% Application:
120      %% algebraic let trig_normalize!*;
121      if trigpreference eq 'sin
122      then algebraic let trig_normalize2sin!*
123      else algebraic let trig_normalize2cos!*;
124      if hyppreference eq 'sinh
125      then algebraic let trig_normalize2sinh!*
126      else algebraic let trig_normalize2cosh!*;
127
128      %% f := algebraic f;
129
130      if not keepalltrig or direction memq '(combine compact)
131      then f := algebraic(f where trig_standardize!*);
132
133      if mode then f :=
134         if mode eq 'trig then
135            behandle algebraic(f where hyp2trig!*)
136         else if mode eq 'hyp then
137            << f := behandle(f);  algebraic(f where trig2hyp!*) >>
138         else if mode eq 'expon then
139            algebraic(f where trig2exp!*);
140
141      if direction eq 'expand then
142         algebraic(begin scalar u;
143            %% Handling of dependent variables
144            let trig_expand_addition!*;
145            %% f := f;
146            symbolic(u := subs_symbolic_multiples(reval f, opt_args));
147            symbolic(f := car u);       % substituted term
148            let trig_expand_multiplication!*;
149            f := sub(symbolic cadr u, f); % unsubstitution equations
150            clearrules trig_expand_addition!*,
151               trig_expand_multiplication!*
152         end)
153      else if direction eq 'combine then <<
154         f := algebraic(f where trig_combine!*);
155         if onlytan and keepalltrig then
156            f := algebraic(f where subtan!*)
157      >>;
158      %% algebraic clearrules(trig_normalize!*);
159      algebraic clearrules trig_normalize2sin!*, trig_normalize2cos!*,
160         trig_normalize2sinh!*, trig_normalize2cosh!*;
161
162      if direction eq 'compact then algebraic <<
163         load_package compact;
164         %% f := f where trig_expand!*;
165         f := (f where trig_expand_addition!*,
166            trig_expand_multiplication!*);
167         f := compact(f, {sin(x)**2+cos(x)**2=1})
168      >>;
169
170      if tanpreference then
171         f := if trigpreference eq 'sin then
172            algebraic(f where sin ~x => cos x * tan x)
173         else
174            algebraic(f where cos ~x => sin x / tan x);
175      if tanhpreference then
176         f := if hyppreference eq 'sinh then
177            algebraic(f where sinh ~x => cosh x * tanh x)
178         else
179            algebraic(f where cosh ~x => sinh x / tanh x);
180
181      if dname then <<
182         %% Resimplify using global domain mode:
183         on dname;
184         f := prepsq simp!* f
185      >>;
186
187      return f
188   end;
189
190
191symbolic procedure more_variables(a, b);
192   length find_indets(a, nil) > length find_indets(b, nil);
193
194symbolic procedure find_indets(term, vars);
195   % Watch out!!!  Expect to see the exponential function as "e" here
196   if numberp term then vars            % FJW number (integer)
197   else if atom term then               % FJW variable
198      (if not memq(term, vars) then term . vars)
199   else if cdr term then <<             % FJW examine function arguments only
200      term := cdr term;
201      vars := find_indets(car term, vars);
202      if cdr term then find_indets(cdr term, vars) else vars
203   >> else                              % FJW nullary function
204      find_indets(car term, vars);
205
206
207% auxiliary variables
208
209algebraic operator auxiliary_symbolic_var!*$
210
211symbolic procedure subs_symbolic_multiples(term, opt_args);
212   %% This procedure replaces trig arguments in `term' that differ
213   %% only by a (rational) numerical factor by their lowest common
214   %% denominator, e.g. x/3 and x/4 would be replaced by x' = x/12, so
215   %% that x/3 -> 4x' and x/4 -> 3x'.
216   %% Assumes `term' is a prefix expression.
217   %% `opt_args' is an initial list of user-specified trig arguments.
218   %% Returns a Lisp list:
219   %%   {substituted term, unsubstitution equation list}.
220   if term = 0 then '(0 (list)) else
221   begin scalar arg_list, unsubs, j;
222      opt_args := union(opt_args, nil); % make into set
223      arg_list := get_trig_arguments(term, opt_args);
224      arg_list := for each arg in arg_list collect simp!* arg;
225      j := 0;
226      while arg_list do
227      begin scalar x, x_nu, x_lcm;
228         j := j + 1;
229         x := car arg_list;
230         x_lcm := denr(x_nu := numberget x); % integer
231
232         %% Find args that differ only by a numerical factor, and find
233         %% the lcm of their denominators.  Delete args that have been
234         %% processed.
235         begin scalar tail;  tail := arg_list;
236            while cdr tail do
237            begin scalar y, q, y_den;
238               y := cadr tail;  q := quotsq(x, y);
239               if atom numr q and atom denr q then <<
240                  y_den := integer_content denr y;
241                  %% Integer arithmetic, division guaranteed:
242                  x_lcm := (x_lcm * y_den) / gcdn(x_lcm,y_den);
243                  %% Delete the argument:
244                  cdr tail :=  cddr tail
245               >> else
246                  tail := cdr tail
247            end
248         end;
249
250         arg_list := cdr arg_list;
251
252         if x_lcm neq 1 then <<
253            x := quotsq(x, x_nu); % primitive part
254	    if not domainp numr x then <<
255               x := !*q2a x;
256               depl!* := append(depl!*,
257                   sublis(list (reval x .
258                           list('auxiliary_symbolic_var!*,j)),depl!*));
259% in case of a df(x,...) in the term. This would be nullified. WN
260               term := algebraic
261                  (term where x =>  auxiliary_symbolic_var!*(j)*x_lcm);
262               unsubs := algebraic(auxiliary_symbolic_var!*(j) = x/x_lcm)
263                  . unsubs
264            >>;
265         >>
266      end;
267      return {term, 'list . unsubs}
268   end;
269
270symbolic procedure behandle ex;
271   begin scalar n, d;
272      %% FJW: Force (exp x)^n + (exp(-x))^n to simplify:
273      ex := algebraic(ex where pow2quot!*);
274      %% (Appears to have been unnecessary before REDUCE 3.7.)
275      ex := simp!* ex;
276      n := mk!*sq (numr ex ./ 1);
277      d := mk!*sq (denr ex ./ 1);
278      return algebraic((n where exp2trig1!*)/(d where exp2trig2!*))
279   end;
280
281
282%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283% General support routines %
284%%%%%%%%%%%%%%%%%%%%%%%%%%%%
285
286symbolic procedure get_trig_arguments(term, args);
287   %% Return a SET of all the arguments of the trig functions in the
288   %% expression.  (Note that trig functions are unary!)  The
289   %% arguments may themselves be general expressions -- they need not
290   %% be kernels!
291   if atom term then args else
292   begin scalar f, r;
293      f := car term;                    % function or operator
294      % Winfried Neun, 1 May 2008: you might in very special cases
295      % enter with equations which contain *SQs. These equations are
296      % not perfectly reval'ed to prefix form. This is special for
297      % equations and intentional, I think.  So...
298      if (f = '!*sq) then << term := reval term; f := car term >>;
299      r := cdr term;                    % arguments or operands
300      if f memq '(sin cos sinh cosh) then return
301         if not member(r := car r, args) then r . args else args;
302      for each j in r do args := get_trig_arguments(j, args);
303      return args
304   end;
305
306put('numberget, 'simpfn, 'numberget!-simpfn)$
307symbolic procedure numberget!-simpfn p;
308   %% Return the rational numeric content of a rational expression.
309   %% Algebraic-mode interface.
310   %% Cannot assume a numeric denominator!
311   numberget simp!* car p;
312
313symbolic procedure numberget p;
314   %% Return the rational numeric content of a rational expression.
315   %% Input and output in standard quotient form.
316   %% Assume integer domain mode.
317   %% Cannot assume a numeric denominator!
318   begin scalar n, d, g;
319      n := integer_content numr p;
320      d := integer_content denr p;
321      g := gcdn(n,d);
322      return (n/g) ./ (d/g)
323   end;
324
325% FJW: The following numeric content code is modelled on that in
326% poly/heugcd by James Davenport & Julian Padget.
327
328symbolic procedure integer_content p;
329   %% Extract INTEGER content of (multivariate) polynomial p in
330   %% standard form, assuming default (integer) domain mode!
331   if atom p then
332      p or 0
333   else if atom red p then
334      if red p then gcdn(integer_content lc p, red p)
335      else integer_content lc p
336   else integer_content1(red red p,
337      gcdn(integer_content lc p, integer_content lc red p));
338
339symbolic procedure integer_content1(p,a);
340   if a=1 then 1
341   else if atom p then
342      if p then gcdn(p,a) else a
343   else integer_content1(red p,
344      gcdn(remainder(integer_content lc p,a), a));
345
346
347%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
348% TrigGCD and TrigFactorize %
349%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
350
351symbolic operator degree$
352symbolic procedure degree(p,x);
353   %% cf. deg in poly/polyop
354   if numr(p := simp!* p) then
355      numrdeg(numr p, x) - numrdeg(denr p, x)
356   else 'inf;
357
358symbolic procedure balanced(p,x);
359   %% cf. deg in poly/polyop
360   <<
361      p := simp!* p;
362      numrdeg(numr p, x) = 2*numrdeg(denr p, x)
363   >>;
364
365symbolic procedure coordinated(p, x);
366   %% FJW: Returns t if p contains only even or only odd degree terms
367   %% wrt x; returns nil if p contains both even and odd degree terms.
368   begin scalar kord!*, evendeg, coord;
369      kord!* := {x := !*a2k x};
370      p := reorder numr simp!* p;
371      if domainp p or not(mvar p eq x) then return t; % degree = 0
372      evendeg := remainder(ldeg p, 2) = 0; % leading degree is even
373      p := red p;
374      coord := t;
375      while p and coord do
376         if domainp p or not(mvar p eq x) then <<
377            coord := evendeg eq t;         % degree = 0
378            p := nil
379         >> else <<
380            coord := evendeg eq (remainder(ldeg p, 2) = 0);
381            p := red p
382         >>;
383      return coord
384   end;
385
386flag ('(balanced coordinated), 'boolean)$
387
388algebraic procedure trig2ord(p,x,y);
389   if not balanced(p,x) or not balanced(p,y) then
390      rederr "trig2ord error: polynomial not balanced."
391   else if not coordinated(p,x) or not coordinated(p,y) then
392      rederr "trig2ord error: polynomial not coordinated."
393   else sub(x=sqrt(x), y=sqrt(y), x**degree(p,x)*y**degree(p,y)*p);
394
395algebraic procedure ord2trig(p,x,y);
396   x**(-degree(p,x))*y**(-degree(p,y))*sub(x=x**2, y=y**2, p);
397
398algebraic procedure subpoly2trig(p,x);
399   begin scalar r, d;
400      d := degree(den(p),x);
401      r := sub(x=cos(x)+i*sin(x), p*x**d);
402      return r*(cos(x)-i*sin(x))**d
403   end;
404
405algebraic procedure subpoly2hyp(p,x);
406   begin scalar r, d;
407      d := degree(den(p),x);
408      r := sub(x=cosh(x)+sinh(x), p*x**d);
409      return r*(cosh(x)-sinh(x))**d
410   end;
411
412algebraic procedure varget(p);
413   %% FJW: This procedure returns the variable `x' from an argument
414   %% `p' of the form `n*x', where `n' must be numeric and `x' must be
415   %% a kernel.
416   begin scalar var;
417      if not(var := mainvar num p) then rederr
418         "TrigGCD/Factorize error: no variable specified.";
419      if not numberp(p/var) then rederr
420         "TrigGCD/Factorize error: last arg must be [number*]variable.";
421      return var
422   end;
423
424symbolic procedure trigargcheck(p, var, nu);
425   %% Check validity of trig arguments.  Note that nu may be rational!
426   begin scalar df_arg_var;
427      for each arg in get_trig_arguments(p, nil) do algebraic
428         if (df_arg_var := df(arg,var)) and
429            not fixp(df_arg_var/nu) then
430            rederr "TrigGCD/Factorize error: basis not possible."
431   end;
432
433symbolic operator sub2poly$
434symbolic procedure sub2poly(p, var, nu, x, y);
435   <<
436      trigargcheck(p, var, nu);
437      p := trigsimp1(p, nil);
438      p := algebraic sub(
439         sin var = sin(var/nu),
440         cos var = cos(var/nu),
441         sinh var = sinh(var/nu),
442         cosh var = cosh(var/nu), p);
443      p := trigsimp1(p, nil);
444      algebraic sub(
445         sin var = (x-1/x)/(2i),
446         cos var = (x+1/x)/2,
447         sinh var = (y-1/y)/2,
448         cosh var = (y+1/y)/2, p)
449   >>;
450
451algebraic procedure triggcd(p, q, x);
452   begin scalar not_complex, var, nu, f;
453      symbolic if (not_complex := not !*complex) then on complex;
454      var := varget x;  nu := numberget x;
455      %% xx_x, yy_y should be gensyms?
456      p := sub2poly(p, var, nu, xx_x, yy_y);
457      q := sub2poly(q, var, nu, xx_x, yy_y);
458      if not and(balanced(p,xx_x), balanced(q,xx_x),
459         coordinated(p,xx_x), coordinated(q,xx_x),
460         balanced(p,yy_y), balanced(q,yy_y),
461         coordinated(p,yy_y), coordinated(q,yy_y))
462      then f := 1
463      else
464      begin scalar h, !*nopowers, !*ifactor;
465         symbolic(!*nopowers := t);
466         p := trig2ord(p, xx_x, yy_y);
467         q := trig2ord(q, xx_x, yy_y);
468         h := gcd(num p, num q);
469         h := ord2trig(h, xx_x, yy_y) / lcm(den p, den q);
470         h := subpoly2trig(h, xx_x);
471         h := subpoly2hyp(h, yy_y);
472         h := sub(xx_x=var*nu, yy_y=var*nu, h);
473         h := symbolic trigsimp1(h, nil);
474         %% What follows is an expensive way to extract the primitive
475         %% part!  Try using `integer_content' defined above or
476         %% `comfac' in alg/gcd?
477         h := factorize(num h);
478         if numberp first h then h := rest h;
479         f := for each r in h product r
480      end;
481      symbolic if not_complex then off complex;
482      return f
483   end;
484
485algebraic procedure trigfactorize(p, x);
486   begin scalar not_complex, var, nu, q, factors;
487      symbolic if (not_complex := not !*complex) then on complex;
488      var := varget x;  nu := numberget x;
489      %% xx_x, yy_y should be gensyms?
490      q := sub2poly(p, var, nu, xx_x, yy_y);
491      if not(balanced(q,xx_x) and coordinated(q,xx_x) and
492         balanced(q,yy_y) and coordinated(q,yy_y))
493      then factors := if symbolic !*nopowers then {p} else {{p,1}}
494      else
495      begin scalar pow, content;
496         %% Handle desired factorized form:
497         if symbolic(not !*nopowers) then pow := 1;
498         q := trig2ord(q, xx_x, yy_y);
499         content := 1/den q;
500         factors := {};
501         for each fac in factorize num q do <<
502            if pow then << pow := second fac;  fac := first fac >>;
503            fac := ord2trig(fac, xx_x, yy_y);
504            fac := subpoly2trig(fac, xx_x);
505            fac := subpoly2hyp(fac, yy_y);
506            fac := sub(xx_x=var*nu, yy_y=var*nu, fac);
507            fac := symbolic trigsimp1(fac, nil);
508            if fac freeof var then
509               content := content*(if pow > 1 then fac^pow else fac)
510            else
511            begin scalar !*nopowers;
512               for each u in factorize(fac) do
513                  if u freeof var then <<
514                     u := first u ^ second u;
515                     content := content*(if pow > 1 then u^pow else u);
516                     fac := fac/u
517                  >>;
518               factors := (if pow then {fac,pow} else fac) . factors
519            end
520         >>;
521         if content neq 1 then
522            factors := (if symbolic !*nopowers then content else {content,1})
523               . factors
524      end;
525      symbolic if not_complex then off complex;
526      return factors
527   end;
528
529endmodule;
530
531end;
532