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