1module elem; % Simplification rules for elementary functions. 2 3% Author: Anthony C. Hearn. 4% Modifications by: Herbert Melenk, Rainer Schoepf and others. 5 6% Copyright (c) 1993 The RAND Corporation. All rights reserved. 7 8% Redistribution and use in source and binary forms, with or without 9% modification, are permitted provided that the following conditions are met: 10% 11% * Redistributions of source code must retain the relevant copyright 12% notice, this list of conditions and the following disclaimer. 13% * Redistributions in binary form must reproduce the above copyright 14% notice, this list of conditions and the following disclaimer in the 15% documentation and/or other materials provided with the distribution. 16% 17% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 19% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 20% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 21% CONTRIBUTORS 22% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 23% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 24% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 25% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 26% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 27% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 28% POSSIBILITY OF SUCH DAMAGE. 29% 30 31 32fluid '(!*!*sqrt !*complex !*keepsqrts !*precise !*precise_complex 33 !*rounded dmode!* !*elem!-inherit); 34 35% No references to RPLAC-based functions in this module. 36 37% For a proper bootstrapping the following order of operator 38% declarations is essential: 39 40% sqrt 41% sign with reference to sqrt 42% trigonometrical functions using abs which uses sign 43 44algebraic; 45 46% Square roots. 47 48deflist('((sqrt outer!-simpsqrt)),'simpfn); 49 50% for all x let sqrt x**2=x; 51 52% !*!*sqrt: used to indicate that SQRTs have been used. 53 54% !*keepsqrts: causes SQRT rather than EXPT to be used. 55 56symbolic procedure mksqrt u; 57 if not !*keepsqrts then list('expt,u,list('quotient,1,2)) 58 else <<if null !*!*sqrt then <<!*!*sqrt := t; 59 algebraic for all x let sqrt x**2=x>>; 60 list('sqrt,u)>>; 61 62for all x let df(sqrt x,x)=sqrt x/(2*x); 63 64 65% SIGN operator. 66 67symbolic procedure sign!-of u; 68 % Returns -1,0 or 1 if the sign of u is known. Otherwise nil. 69 (numberp s and s) where s = numr simp!-sign{u}; 70 71symbolic procedure simp!-sign1 u; 72 begin scalar s,n; 73 s:=if atom u then 74 << if !*modular then simpiden {'sign,u} 75 else if numberp u then 76 (if u > 0 then 1 else if u < 0 then -1 else nil) ./ 1 77 else simp!-sign2 u >> 78 else if car u eq '!:rd!: then 79 (if rd!:zerop u then nil else if rd!:minusp u then -1 else 1) ./ 1 80% This may be wrong in general, eg. sign(abs(y)) would always return 1, even if y is later set to 0. 81% else if car u eq 'abs then 1 ./ 1 82 else if car u eq 'minus then negsq simp!-sign1 cadr u 83 else if car u eq 'times then simp!-sign!-times u 84 else if car u eq 'quotient then simp!-sign!-quot u 85 else if car u eq 'plus then simp!-sign!-plus u 86 else if car u eq 'expt then simp!-sign!-expt u 87 else if car u eq 'sqrt then simp!-sign!-sqrt u 88 else simp!-sign2 u; 89 n:=numr s; 90 if not numberp n or n=1 or n=-1 or n=0 then return s; 91 typerr(n,"sign value"); 92 end; 93 94symbolic procedure simp!-sign2 u; 95 % fallback to standard rules: 1. try numeric evaluation, 2. call simpiden 96 (if x then x ./ 1 else simpiden{'sign,u}) where x := rd!-sign u; 97 98symbolic procedure simp!-sign u; 99 simp!-sign1 reval car u; 100 101symbolic inline procedure sq!-is!-sign u; 102 % Returns t is s.q. u is either 1, -1, or 0 103 denr u = 1 and (nu=1 or nu=-1 or nu=0) where nu=numr u; 104 105symbolic procedure simp!-sign!-times w; 106 % Factor all known signs out of the product. 107 begin scalar n,s,x; 108 n:=1; 109 for each f in cdr w do 110 << x:=simp!-sign1 f; 111 % Make sure that only values 1,-1, and 0 are used here, everything else is left in place 112 if sq!-is!-sign x then n:=n * numr x else s:=f.s>>; 113 s:=if null s then '(1 . 1) 114 else simp!-sign2(if cdr s then 'times.reversip s else car s); 115 return multsq (n ./ 1,s) 116 end; 117 118symbolic procedure simp!-sign!-quot w; 119 % Factor known signs out of the quotient. 120 begin scalar x,y,flg,z; 121 if eqcar(cadr w,'minus) then << flg:=t; w := {car w,cadr cadr w, caddr w} >>; 122 x := simp!-sign1 cadr w; % numerator 123 y := simp!-sign1 caddr w; % denominator 124 if sq!-is!-sign x or sq!-is!-sign y then z := quotsq (x,y) 125 else z := simp!-sign2 w; 126 return if flg then negsq z else z 127 end; 128 129symbolic procedure simp!-sign!-plus w; 130 % Stop sign evaluation as soon as two different signs 131 % or one unknown sign were found. 132 begin scalar n,m,x,q; 133 for each f in cdr w do if null q then 134 <<x:=simp!-sign1 f; 135 m:=if sq!-is!-sign x then numr x; 136 if null m or n and m neq n then q:=t; 137 n:=m>>; 138 return if null q then n ./ 1 else simp!-sign2 w; 139 end; 140 141symbolic procedure simp!-sign!-expt w; 142 (if fixp ex and evenp ex and not(!*complex or !*precise_complex) then (1 ./ 1) 143 else ( 144 if fixp ex and sq!-is!-sign sb 145 then (if not evenp ex and numr sb < 0 then -1 else 1) ./ 1 146 else if fixp ex and not(!*complex or !*precise_complex) then sb 147 else if ex = '(quotient 1 2) and sb = (1 ./ 1) then 1 ./ 1 148 else if sb = (1 ./ 1) and realvaluedp ex then (1 ./ 1) 149 else simp!-sign2 w 150 ) where sb := simp!-sign1 cadr w 151 ) where ex := caddr w; 152 153symbolic procedure simp!-sign!-sqrt w; 154 (if u = (1 ./ 1) then u else simp!-sign2 w) 155 where u := simp!-sign!-expt {'expt,cadr w,'(quotient 1 2)}; 156 157fluid '(rd!-sign!*); 158 159symbolic procedure rd!-sign u; 160 % if U is constant evaluable return sign of u. 161 % the value is set aside. 162 if pairp rd!-sign!* and u=car rd!-sign!* then cdr rd!-sign!* 163 else 164 if !*complex or !*rounded or not constant_exprp u then nil 165 else 166 (begin scalar x,y,dmode!*; 167 setdmode('rounded,t); 168 x := aeval u; 169 if evalnumberp x and 0=reval {'impart,x} 170 then y := if evalgreaterp(x,0) then 1 else 171 if evalequal(x,0) then 0 else -1; 172 setdmode('rounded,nil); 173 rd!-sign!*:=(u.y); 174 return y 175 end) where alglist!*=alglist!*; 176 177symbolic operator rd!-sign; 178 179operator sign; 180 181put('sign,'simpfn,'simp!-sign); 182 183% The rules for products and sums are covered by the routines 184% below in order to avoid a combinatoric explosion. Abs (sign ~x) 185% cannot be defined by a rule because the evaluation of abs needs 186% sign. 187 188sign_rules := 189 { 190%% sign ~x => (if x>0 then 1 else if x<0 then -1 else 0) 191%% when numberp x and impart x=0, 192%% sign(e) => 1, 193%% sign(pi) => 1, 194%% sign(-~x) => -sign(x), 195%% sign( ~x * ~y) => sign x * sign y 196%% when numberp sign x or numberp sign y, 197%% sign( ~x / ~y) => sign x * sign y 198%% when y neq 1 and (numberp sign x or numberp sign y), 199%% sign( ~x + ~y) => sign x when sign x = sign y, 200%% sign( ~x ^ ~n) => 1 when fixp (n/2) and 201%% lisp(not (!*complex or !*precise_complex)), 202%% sign( ~x ^ ~n) => sign x^n when fixp n and numberp sign x, 203%% sign( ~x ^ ~n) => sign x when fixp n and 204%% lisp(not (!*complex or !*precise_complex)), 205%% sign(sqrt ~a) => 1 when sign a=1, 206 sign(sinh ~x) => sign(x) when numberp sign(x) or realvaluedp x, 207 sign(cosh ~x) => 1 when realvaluedp x, 208 sign(tanh ~x) => sign(x) when numberp sign(x) or realvaluedp x, 209%% sign( ~a ^ ~x) => 1 when sign a=1 and realvaluedp x, 210%% sign(abs ~a) => 1, 211%% sign ~a => rd!-sign a when rd!-sign a, 212 % Next rule here for convenience. 213 abs(~x)^2 => x^2 when symbolic not !*precise or realvaluedp x 214 }$ 215 % $ above needed for bootstrap. 216 217let sign_rules; 218 219% Rule for I**2. 220 221remflag('(i),'reserved); 222 223let i**2= -1; 224 225flag('(e i nil pi),'reserved); % Leave out T for now. 226 227% Some more reserved ids (for realroot) 228 229flag('(positive negative infinity),'reserved); 230 231% Logarithms. 232 233%%let log(e)= 1, 234%% log(1)= 0, 235%% log10(10) = 1, 236%% log10(1) = 0; 237%% 238%%for all x let logb(1,x) = 0, 239%% logb(x,x) = 1, 240%% logb(x,e) = log(x), 241%% logb(x,10) = log10(x); 242%% 243%%for all x let log(e**x)=x; % e**log x=x now done by simpexpt. 244%% 245%%for all x let logb(a**x,a)=x; 246%% 247%%for all x let log10(10**x)=x; 248 249for all x let 10^log10(x)=x; 250 251%% Remove these rules as they return wrong results in complex mode 252%% and are superceded by the code in poly/compopr.red 253%let impart(log(~a)) => 0 when numberp a and a>0; 254%let repart(log(~a)) => log(a) when numberp a and a>0; 255 256%% The following rule interferes with HE vector simplification, 257%% until the problem is resolved use a more complicated rule 258%for all a,x let a^logb(x,a)=x; 259for all a,b,x such that a=b let a^logb(x,b)=x; 260 261 262% The next rule is implemented via combine/expand logs. 263 264% for all x,y let log(x*y) = log x + log y, log(x/y) = log x - log y; 265 266let df(log(~x),~x) => 1/x; 267 268let df(log(~x/~y),~z) => df(log x,z) - df(log y,z); 269 270let df(log10(~x),~x) => 1/(x*log(10)); 271 272let df(logb(~x,~a),~x) => 1/(x*log(a)) - logb(x,a)/(a*log(a))*df(a,x); 273 274% Trigonometrical functions. 275 276deflist('((acos simpiden) (asin simpiden) (atan simpiden) 277 (acosh simpiden) (asinh simpiden) (atanh simpiden) 278 (acot simpiden) (cos simpiden) (sin simpiden) (tan simpiden) 279 (sec simpiden) (sech simpiden) (csc simpiden) (csch simpiden) 280 (cot simpiden)(acot simpiden)(coth simpiden)(acoth simpiden) 281 (cosh simpiden) (sinh simpiden) (tanh simpiden) 282 (asec simpiden) (acsc simpiden) 283 (asech simpiden) (acsch simpiden) 284 ),'simpfn); 285 286% The following declaration causes the simplifier to pass the full 287% expression (including the function) to simpiden. 288 289flag ('(acos asin atan acosh acot asinh atanh cos sin tan cosh sinh tanh 290 csc csch sec sech cot acot coth acoth asec acsc asech acsch), 291 'full); 292 293% flag ('(atan),'oddreal); 294 295flag('(acoth acsc acsch asin asinh atan atanh sin tan csc csch sinh 296 tanh cot coth), 297 'odd); 298 299flag('(cos sec sech cosh),'even); 300 301flag('(cot coth csc csch acoth),'nonzero); 302 303% added by A Barnes 2021 304deflist('((sec 1) (sech 1) (csc 1) (csch 1) (cot 1) (coth 1)), 305 'number!-of!-args); 306 307% In the following rules, it is not necessary to let f(0)=0, when f 308% is odd, since simpiden already does this. 309 310% Some value have been commented out since these can be computed from 311% other functions. 312 313let cos(0)= 1, 314 % sec(0)= 1, 315 % cos(pi/12)=sqrt(2)/4*(sqrt 3+1), 316 sin(pi/12)=sqrt(2)/4*(sqrt 3-1), 317 sin(5pi/12)=sqrt(2)/4*(sqrt 3+1), 318 % cos(pi/6)=sqrt 3/2, 319 sin(pi/6)= 1/2, 320 % cos(pi/4)=sqrt 2/2, 321 sin(pi/4)=sqrt 2/2, 322 % cos(pi/3) = 1/2, 323 sin(pi/3) = sqrt(3)/2, 324 cos(pi/2)= 0, 325 sin(pi/2)= 1, 326 sin(pi)= 0, 327 cos(pi)=-1, 328 cosh 0=1, 329 sech(0) =1, 330 sinh(i) => i*sin(1), 331 cosh(i) => cos(1), 332 acosh(0) => i*pi/2, 333 acosh(1) => 0, 334 acosh(-1) => i*pi, 335 acoth(0) => i*pi/2 336 % acos(0)= pi/2, 337 % acos(1)=0, 338 % acos(1/2)=pi/3, 339 % acos(sqrt 3/2) = pi/6, 340 % acos(sqrt 2/2) = pi/4, 341 % acos(1/sqrt 2) = pi/4 342 % asin(1/2)=pi/6, 343 % asin(-1/2)=-pi/6, 344 % asin(1)=pi/2, 345 % asin(-1)=-pi/2 346 ; 347 348for all x let cos acos x=x, sin asin x=x, tan atan x=x, 349 cosh acosh x=x, sinh asinh x=x, tanh atanh x=x, 350 cot acot x=x, coth acoth x=x, sec asec x=x, 351 csc acsc x=x, sech asech x=x, csch acsch x=x; 352 353for all x let acos(-x)=pi-acos(x), 354 asec(-x)=pi-asec(x), 355 acot(-x)=pi-acot(x); 356 357% Fold the elementary trigonometric functions down to the origin. 358 359let 360 361 sin( (~~w + ~~k*pi)/~~d) 362 => (if evenp fix repart(k/d) then 1 else -1) 363 * sin(w/d + ((k/d)-fix repart(k/d))*pi) 364 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)), 365 366 sin( ~~k*pi/~~d) => sin((1-k/d)*pi) 367 when ratnump(k/d) and k/d > 1/2, 368 369 cos( (~~w + ~~k*pi)/~~d) 370 => (if evenp fix repart(k/d) then 1 else -1) 371 * cos(w/d + ((k/d)-fix repart(k/d))*pi) 372 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)), 373 374 cos( ~~k*pi/~~d) => -cos((1-k/d)*pi) 375 when ratnump(k/d) and k/d > 1/2, 376 377% next two rules needed with current definition of knowledge_about 378 csc( (~~w + ~~k*pi)/~~d) 379 => (if evenp fix repart(k/d) then 1 else -1) 380 * csc(w/d + ((k/d)-fix repart(k/d))*pi) 381 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)), 382 383 csc( ~~k*pi/~~d) => csc((1-k/d)*pi) 384 when ratnump(k/d) and k/d > 1/2, 385 386 sec( (~~w + ~~k*pi)/~~d) 387 => (if evenp fix repart(k/d) then 1 else -1) 388 * sec(w/d + ((k/d)-fix repart(k/d))*pi) 389 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)), 390 391 sec( ~~k*pi/~~d) => -sec((1-k/d)*pi) 392 when ratnump(k/d) and k/d > 1/2, 393 394 tan( (~~w + ~~k*pi)/~~d) 395 => tan(w/d + ((k/d)-fix repart(k/d))*pi) 396 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)), 397 398 cot( (~~w + ~~k*pi)/~~d) 399 => cot(w/d + ((k/d)-fix repart(k/d))*pi) 400 when ((ratnump(rp) and abs(rp) >= 1) where rp => repart(k/d)); 401 402% The following rules follow the pattern 403% sin(~x + pi/2)=> cos(x) when x freeof pi 404% however allowing x to be a quotient and a negative pi/2 shift. 405% We need to handleonly pi/2 shifts here because 406% the bigger shifts are already covered by the rules above. 407% 408% Note the use of ~~d instead of ~d in the denominator for rational k. 409 410let sin((~x + ~~k*pi)/~~d) => sign repart(k/d)*cos(x/d + i*pi*impart(k/d)) 411 when abs repart(k/d) = 1/2, 412 413 cos((~x + ~~k*pi)/~~d) => -sign repart(k/d)*sin(x/d + i*pi*impart(k/d)) 414 when abs repart(k/d) = 1/2, 415 416 csc((~x + ~~k*pi)/~~d) => sign repart(k/d)*sec(x/d + i*pi*impart(k/d)) 417 when abs repart(k/d) = 1/2, 418 419 sec((~x + ~~k*pi)/~~d) => -sign repart(k/d)*csc(x/d + i*pi*impart(k/d)) 420 when abs repart(k/d) = 1/2, 421 422 tan((~x + ~~k*pi)/~~d) => -cot(x/d + i*pi*impart(k/d)) 423 when abs repart(k/d) = 1/2, 424 425 cot((~x + ~~k*pi)/~~d) => -tan(x/d + i*pi*impart(k/d)) 426 when x freeof pi and abs repart(k/d) = 1/2; 427 428% Inherit function values. 429 430symbolic (!*elem!-inherit := t); 431 432symbolic procedure knowledge_about(op,arg,top); 433 % True if the form '(op arg) can be formally simplified. 434 % Avoiding recursion from rules for the target operator top by 435 % a local remove of the property opmtch. 436 % The internal switch !*elem!-inherit!* allows us to turn the 437 % inheritage temporarily off. 438 if dmode!* eq '!:rd!: or dmode!* eq '!:cr!: 439 or null !*elem!-inherit then nil else 440 (begin scalar r,old; 441 old:=get(top,'opmtch); put(top,'opmtch,nil); 442 unwind!-protect(r:= errorset!*({'aeval,mkquote{op,arg}},nil), 443 put(top,'opmtch,old)); 444 return not errorp r and not smemq(op,car r) 445 and not smemq(top,car r); 446 end) where varstack!*=nil; 447 448symbolic operator knowledge_about; 449 450symbolic procedure trigquot(n,d); 451 % Form a quotient n/d, replacing sin and cos by tan/cot 452 % whenver possible. 453 begin scalar m,u,v,w; 454 u:=if eqcar(n,'minus) then <<m:=t; cadr n>> else n; 455 v:=if eqcar(d,'minus) then <<m:=not m; cadr d>> else d; 456 if pairp u and pairp v then 457 if car u eq 'sin and car v eq 'cos and cadr u=cadr v 458 then w:='tan else 459 if car u eq 'cos and car v eq 'sin and cadr u=cadr v 460 then w:='cot; 461 if null w then return{'quotient,n,d}; 462 w:={w,cadr u}; 463 return if m then {'minus,w} else w; 464 end; 465 466symbolic operator trigquot; 467 468% cos, tan, cot, sec, csc inherit from sin. 469 470 471let cos(~x)=>sin(x+pi/2) 472 when not(x=0) and (x+pi/2)/pi freeof pi and knowledge_about(sin,x+pi/2,cos), 473 cos(~x)=>-sin(x-pi/2) 474 when not(x=0) and (x-pi/2)/pi freeof pi and knowledge_about(sin,x-pi/2,cos), 475 tan(~x)=>trigquot(sin(x),cos(x)) when knowledge_about(sin,x,tan), 476 cot(~x)=>trigquot(cos(x),sin(x)) when knowledge_about(sin,x,cot), 477 sec(~x)=>1/cos(x) when knowledge_about(cos,x,sec), 478 csc(~x)=>1/sin(x) when knowledge_about(sin,x,csc); 479 480% area functions 481 482let asin(~x)=>pi/2 - acos(x) when knowledge_about(acos,x,asin), 483 acot(~x)=>pi/2 - atan(x) when knowledge_about(atan,x,acot), 484 acsc(~x) => asin(1/x) when knowledge_about(asin,1/x,acsc), 485 asec(~x) => acos(1/x) when knowledge_about(acos,1/x,asec), 486 acsch(~x) => acsc(-i*x)/i when knowledge_about(acsc,-i*x,acsch), 487 asech(~x) => asec(x)/i when knowledge_about(asec,x,asech); 488 489% hyperbolic functions 490 491let sinh(i*~x)=>i*sin(x) when knowledge_about(sin,x,sinh), 492 sinh(i*~x/~n)=>i*sin(x/n) when knowledge_about(sin,x/n,sinh), 493 cosh(i*~x)=>cos(x) when knowledge_about(cos,x,cosh), 494 cosh(i*~x/~n)=>cos(x/n) when knowledge_about(cos,x/n,cosh), 495 cosh(~x)=>-i*sinh(x+i*pi/2) 496 when not(x=0) and (x+i*pi/2)/pi freeof pi 497 and knowledge_about(sinh,x+i*pi/2,cosh), 498 cosh(~x)=>i*sinh(x-i*pi/2) 499 when not(x=0) and (x-i*pi/2)/pi freeof pi 500 and knowledge_about(sinh,x-i*pi/2,cosh), 501 tanh(~x)=>sinh(x)/cosh(x) when knowledge_about(sinh,x,tanh), 502 coth(~x)=>cosh(x)/sinh(x) when knowledge_about(sinh,x,coth), 503 sech(~x)=>1/cosh(x) when knowledge_about(cosh,x,sech), 504 csch(~x)=>1/sinh(x) when knowledge_about(sinh,x,csch); 505 506let acsch(~x) => asinh(1/x) when knowledge_about(asinh,1/x,acsch), 507 asech(~x) => acosh(1/x) when knowledge_about(acosh,1/x,asech), 508 asinh(~x) => -i*asin(i*x) when i*x freeof i 509 and knowledge_about(asin,i*x,asinh); 510 511 512% hyperbolic functions 513 514let 515 516 sinh( (~~w + ~~k*pi)/~~d) 517 => (if evenp fix impart(k/d) then 1 else -1) 518 * sinh(w/d + (k/d-i*fix impart(k/d))*pi) 519 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)), 520 521 sinh( ~~k*pi/~~d) => sinh((i-k/d)*pi) 522 when ratnump(i*k/d) and abs(i*k/d) > 1/2, 523 524 cosh( (~~w + ~~k*pi)/~~d) 525 => (if evenp fix impart(k/d) then 1 else -1) 526 * cosh(w/d + (k/d-i*fix impart(k/d))*pi) 527 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)), 528 529 cosh( ~~k*pi/~~d) => -cosh((i-k/d)*pi) 530 when ratnump(i*k/d) and abs(i*k/d) > 1/2, 531 532% next two rules needed with current definition of knowledge_about 533 csch( (~~w + ~~k*pi)/~~d) 534 => (if evenp fix impart(k/d) then 1 else -1) 535 * csch(w/d + (k/d-i*fix impart(k/d))*pi) 536 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)), 537 538 csch( ~~k*pi/~~d) => csch((i-k/d)*pi) 539 when ratnump(i*k/d) and abs(i*k/d) > 1/2, 540 541 sech( (~~w + ~~k*pi)/~~d) 542 => (if evenp fix impart(k/d) then 1 else -1) 543 * sech(w/d + (k/d-i*fix impart(k/d))*pi) 544 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)), 545 546 sech( ~~k*pi/~~d) => -sech((i-k/d)*pi) 547 when ratnump(i*k/d) and abs(i*k/d) > 1/2, 548 549 tanh( (~~w + ~~k*pi)/~~d) 550 => tanh(w/d + (k/d-i*fix impart(k/d))*pi) 551 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)), 552 553 coth( (~~w + ~~k*pi)/~~d) 554 => coth(w/d + (k/d-i*fix impart(k/d))*pi) 555 when ((ratnump(ip) and abs(ip) >= 1) where ip => impart(k/d)); 556 557% The following rules follow the pattern 558% sinh(~x + i*pi/2)=> cosh(x) when x freeof pi 559% however allowing x to be a quotient and a negative i*pi/2 shift. 560% We need to handle only pi/2 shifts here because 561% the bigger shifts are already covered by the rules above. 562 563let 564 sinh((~~x + ~~k*pi)/~~d) => i*sign impart(k/d)*cosh(x/d + pi*repart(k/d)) 565 when abs impart(k/d) = 1/2, 566 567 cosh((~~x + ~~k*pi)/~~d) => i*sign impart(k/d)*sinh(x/d+pi*repart(k/d)) 568 when abs impart(k/d) = 1/2, 569 570% next 2 rules could be omitted and use inheritance 571 csch((~~x + ~~k*pi)/~~d) => -i*sign impart(k/d)*sech(x/d + pi*repart(k/d)) 572 when abs impart(k/d) = 1/2, 573 574 sech((~~x + ~~k*pi)/~~d) => -i*sign impart(k/d)*csch(x/d+pi*repart(k/d)) 575 when abs impart(k/d) = 1/2, 576 577 tanh((~~x + ~~k*pi)/~~d) => coth(x/d+pi*repart(k/d)) 578 when abs impart(k/d) = 1/2, 579 580 coth((~~x + ~~k*pi)/~~d) => tanh(x/d+pi*repart(k/d)) 581 when abs impart(k/d) = 1/2; 582 583 584% Transfer inverse function values from cos to acos and tan to atan. 585% Negative values not needed. 586 587acos_rules := 588 symbolic( 589 'list . for j:=0:12 join 590 (if eqcar(q,'acos) and cadr q=w then {{'replaceby,q,u}}) 591 where q=reval{'acos,w} 592 where w=reval{'cos,u} 593 where u=reval{'quotient,{'times,'pi,j},12})$ 594 595let acos_rules; 596 597clear acos_rules; 598 599atan_rules := 600 symbolic( 601 'list . for j:=0:5 join 602 (if eqcar(q,'atan) and cadr q=w then {{'replaceby,q,u}}) 603 where q= reval{'atan,w} 604 where w= reval{'tan,u} 605 where u= reval{'quotient,{'times,'pi,j},12})$ 606 607let atan_rules; 608 609clear atan_rules; 610 611 612repart(pi) := pi$ % $ used for bootstrapping purposes. 613impart(pi) := 0$ 614 615% ***** Differentiation rules *****. 616 617for all x let df(acos(x),x)= -sqrt(1-x**2)/(1-x**2), 618 df(asin(x),x)= sqrt(1-x**2)/(1-x**2), 619 df(atan(x),x)= 1/(1+x**2), 620 df(acosh(x),x)= sqrt(x-1)*sqrt(x+1)/(x**2-1), 621 df(acot(x),x)= -1/(1+x**2), 622 df(acoth(x),x)= -1/(1-x**2), 623 df(asinh(x),x)= sqrt(x**2+1)/(x**2+1), 624 df(atanh(x),x)= 1/(1-x**2), 625 df(acoth(x),x)= 1/(1-x**2), 626 df(cos x,x)= -sin(x), 627 df(sin x,x)= cos(x), 628 df(sec x,x) = sec(x)*tan(x), 629 df(csc x,x) = -csc(x)*cot(x), 630 df(tan x,x)=1 + tan x**2, 631 df(sinh x,x)=cosh x, 632 df(cosh x,x)=sinh x, 633 df(sech x,x) = -sech(x)*tanh(x), 634% df(tanh x,x)=sech x**2, 635 % J.P. Fitch prefers this one for integration purposes 636 df(tanh x,x)=1-tanh(x)**2, 637 df(csch x,x)= -csch x*coth x, 638 df(cot x,x)=-1-cot x**2, 639 df(coth x,x)=1-coth x**2; 640 641% Beware we cannot take an x factor inside the sqrt in the 4 formulae below 642% as this changes the cuts and introduces sign errors in part of the domain. 643let df(acsc(~x),x) => -1/(x^2*sqrt(1-1/x^2)), 644 df(asec(~x),x) => 1/(x^2*sqrt(1-1/x^2)), 645 df(acsch(~x),x)=> -1/(x^2*sqrt(1+ 1/x^2)), 646 df(asech(~x),x)=> -1/(x^2*sqrt(1/x-1)*sqrt(1/x+1)); 647 648% rules for atan2 649% This rule must appear before simp!-atan2 during bootstrap. 650% Otherwise, the simplication of the right hand side calls ezgcdf which 651% isn't defined yet. 652% Temporarily set up simiden as simpfn (overwritten below) 653put('atan2,'simpfn,'simpiden); %temporary, see below 654let df(atan2(~y,~x),~z) => (x*df(y, z)-y*df(x, z))/(x^2+y^2); 655 656% This procedure now works for complex arguments and gives results compatible 657% with the numerical results returned by cratan2!* (and rdatan2!*) 658% It may currently be somewhat inefficient as it frequently inter-converts 659% between prefix and standard quotient forms 660symbolic procedure simp!-atan2 u; 661<< if length u neq 2 then 662 rerror(alg,17,list("Wrong number of arguments to", 'atan2)); 663 (if val then val % where val=valuechk('atan2, u) 664 else % NB some simplifications seem to fail if !*complex=t 665 begin scalar x,y,z,v,w, !*complex; 666 y := reval car u; 667 x := reval cadr u; 668 669 % deal with usual case of real arguments separately as more 670 % simp!-sign1 succeeds more often with simpler arguments 671 if null numr simpimpart list x and null numr simpimpart list y then 672 return simpatan2r(y, x); 673 674 % save simplified original arguments for default return value 675 u := {y, x}; 676 677 if x=0 then << 678 z:= simp!-sign1 prepsq simprepart list y; 679 if null numr z then z:= simp!-sign1 reval {'quotient, y, 'i}; 680 if denr z=1 and fixp numr z then 681 return multsq(z, simp {'quotient,'pi, 2})>> 682 else if y=0 then << 683 z:= simp!-sign1 prepsq simprepart list x; 684 if null numr z then z:= simp!-sign1 reval {'quotient, x, 'i}; 685 if denr z=1 and fixp numr z then 686 if numr z=1 then return nil ./ 1 687 else return simp 'pi>> 688 else << 689 z := simp!* {'plus, {'expt, x, 2}, {'expt, y, 2}}; 690 if null numr z then 691 rerror(alg, 212, "Essential singularity encountered in atan"); 692 x := {'quotient, x, {'sqrt, z:=prepsq z}}; 693 y := {'quotient, y, {'sqrt, z}}; 694 z:= simp!-sign1 prepsq simprepart list x; 695 if null numr z then z := simp!-sign1 reval {'quotient, x, 'i}; 696 v := simp!-sign1 prepsq simprepart list y; 697 if null numr v then v := simp!-sign1 reval {'quotient, y, 'i}; 698 if denr z=1 and fixp numr z and denr v=1 and fixp numr v then << 699 w := simp {'atan, prepsq rationalizesq simp {'quotient, y, x}}; 700 if numr z=1 then return w 701 else if numr v=1 then return addsq(simp 'pi, w) 702 else return subtrsq(w, simp 'pi) >> 703 >>; 704 return simpiden {'atan2, car u, cadr u} 705 end) where val = valuechk('atan2, u) 706>>; 707 708put('atan2,'simpfn,'simp!-atan2); 709 710symbolic procedure simpatan2r(y, x); 711begin scalar z,v,w; 712% Arguments are real and simplified 713 if x=0 then << 714 z := simp!-sign1 y; 715 if null numr z then rerror(alg, 211, "atan2(0, 0) formed") 716 else return quotsq(simp {'quotient, 'pi, 2}, z)>> 717 else if y=0 then << 718 z := simp!-sign1 x; 719 return multsq(addsq(1 ./ 1, quotsq((-1) ./ 1, z)), 720 simp {'quotient, 'pi, 2})>> 721 else << 722 z := simp!-sign1 x; v := simp!-sign1 y; 723 if denr z=1 and fixp numr z and denr v=1 and fixp numr v then << 724 w := simp {'atan, {'quotient, y, x}}; 725 if numr z=1 then return w 726 else if numr v=1 then return addsq(simp 'pi, w) 727 else return subtrsq(w, simp 'pi) >> 728 >>; 729 return simpiden {'atan2, y, x}; 730end; 731 732 733%for all x let e**log x=x; % Requires every power to be checked. 734 735for all x,y let df(x**y,x)= y*x**(y-1), 736 df(x**y,y)= log x*x**y; 737 738 739% erf, erfc, erfi, exp and dilog. 740 741operator dilog,erf,erfi,erfc; 742 743let { 744 dilog(0) => pi**2/6, 745 dilog(1) => 0, 746 dilog(2) => -pi^2/12, 747 dilog(-1) => pi^2/4-i*pi*log(2) 748}; 749 750let df(dilog(~x),(~x)) => -log(x)/(x-1); 751 752 753 754let erf 0=0; 755 756for all x let erf(-x)=-erf x; 757 758for all x let df(erf x,x)=2*sqrt(pi)*e**(-x**2)/pi; 759 760let erf (~x) => compute!:int!:functions(x,erf) 761 when numberp x and abs(x)<5 and lisp !*rounded; 762 763let erfc(~x) => 1 - erf(x); 764let erfi(~z) => -i * erf(i*z); 765 766for all x let exp(x)=e**x; 767 768% Supply missing argument and simplify 1/4 roots of unity. 769 770let e**(i*pi/2) = i, 771 e**(i*pi) = -1; 772% e**(3*i*pi/2)=-i; 773 774 775 776% Rule for derivative of absolute value. 777 778for all x let df(abs x,x)=abs x/x; 779 780% More trigonometrical rules. 781 782invtrigrules := { 783 sin(atan ~u) => u/sqrt(1+u^2), 784 cos(atan ~u) => 1/sqrt(1+u^2), 785 sin(2*atan ~u) => 2*u/(1+u^2), 786 cos(2*atan ~u) => (1-u^2)/(1+u^2), 787 sin(~n*atan ~u) => sin((n-2)*atan u) * (1-u^2)/(1+u^2) + 788 cos((n-2)*atan u) * 2*u/(1+u^2) 789 when fixp n and n>2, 790 cos(~n*atan ~u) => cos((n-2)*atan u) * (1-u^2)/(1+u^2) - 791 sin((n-2)*atan u) * 2*u/(1+u^2) 792 when fixp n and n>2, 793 sin(acos ~u) => sqrt(1-u^2), 794 cos(asin ~u) => sqrt(1-u^2), 795 sin(2*acos ~u) => 2 * u * sqrt(1-u^2), 796 cos(2*acos ~u) => 2*u^2 - 1, 797 sin(2*asin ~u) => 2 * u * sqrt(1-u^2), 798 cos(2*asin ~u) => 1 - 2*u^2, 799 sin(~n*acos ~u) => sin((n-2)*acos u) * (2*u^2 - 1) + 800 cos((n-2)*acos u) * 2 * u * sqrt(1-u^2) 801 when fixp n and n>2, 802 cos(~n*acos ~u) => cos((n-2)*acos u) * (2*u^2 - 1) - 803 sin((n-2)*acos u) * 2 * u * sqrt(1-u^2) 804 when fixp n and n>2, 805 sin(~n*asin ~u) => sin((n-2)*asin u) * (1 - 2*u^2) + 806 cos((n-2)*asin u) * 2 * u * sqrt(1-u^2) 807 when fixp n and n>2, 808 cos(~n*asin ~u) => cos((n-2)*asin u) * (1 - 2*u^2) - 809 sin((n-2)*asin u) * 2 * u * sqrt(1-u^2) 810 when fixp n and n>2 811% Next rule causes a simplification loop in solve(atan y=y). 812% atan(~x) => acos((1-x^2)/(1+x^2)) * sign (x) / 2 813% when symbolic(not !*complex) and x^2 neq -1 814% and acos((1-x^2)/(1+x^2)) freeof acos 815}$ 816 817% The following are useful for simplifying sqrts of complex values 818% Added by A Barnes, Aug 2015 819invtrigrules2 := { 820 sin(atan(~x)/2) => sin(atan((sqrt(1+x^2)-1)/x)), 821 cos(atan(~x)/2) => cos(atan((sqrt(1+x^2)-1)/x)) 822% tan(atan(~x)/2) => (sqrt(1+x^2)-1)/x 823}$ 824 825let invtrigrules2; 826 827invhyprules := { 828 sinh(atanh ~u) => u/sqrt(1-u^2), 829 cosh(atanh ~u) => 1/sqrt(1-u^2), 830 sinh(2*atanh ~u) => 2*u/(1-u^2), 831 cosh(2*atanh ~u) => (1+u^2)/(1-u^2), 832 sinh(~n*atanh ~u) => sinh((n-2)*atanh u) * (1+u^2)/(1-u^2) + 833 cosh((n-2)*atanh u) * 2*u/(1-u^2) 834 when fixp n and n>2, 835 cosh(~n*atanh ~u) => cosh((n-2)*atanh u) * (1+u^2)/(1-u^2) + 836 sinh((n-2)*atanh u) * 2*u/(1-u^2) 837 when fixp n and n>2, 838 sinh(acosh ~u) => sqrt(u-1)*sqrt(u+1), 839 cosh(asinh ~u) => sqrt(1+u^2), 840 sinh(2*acosh ~u) => 2 * u * sqrt(u-1)*sqrt(u+1), 841 cosh(2*acosh ~u) => 2*u^2 - 1, 842 sinh(2*asinh ~u) => 2 * u * sqrt(1+u^2), 843 cosh(2*asinh ~u) => 1 + 2*u^2, 844 sinh(~n*acosh ~u) => sinh((n-2)*acosh u) * (2*u^2 - 1) + 845 cosh((n-2)*acosh u) * 2 * u * sqrt(u-1)*sqrt(u+1) 846 when fixp n and n>2, 847 cosh(~n*acosh ~u) => cosh((n-2)*acosh u) * (2*u^2 - 1) + 848 sinh((n-2)*acosh u) * 2 * u * sqrt(u-1)*sqrt(u+1) 849 when fixp n and n>2, 850 sinh(~n*asinh ~u) => sinh((n-2)*asinh u) * (1 + 2*u^2) + 851 cosh((n-2)*asinh u) * 2 * u * sqrt(1+u^2) 852 when fixp n and n>2, 853 cosh(~n*asinh ~u) => cosh((n-2)*asinh u) * (1 + 2*u^2) + 854 sinh((n-2)*asinh u) * 2 * u * sqrt(1+u^2) 855 when fixp n and n>2, 856 atanh(~x) => acosh((1+x^2)/(1-x^2)) * sign (x) / 2 857 when symbolic(not !*complex) 858 and x^2 neq 1 859 and acosh((1+x^2)/(1-x^2)) freeof acosh 860}$ 861 862let invtrigrules,invhyprules; 863 864trig_imag_rules := { 865 sin(i * ~~x / ~~y) => i * sinh(x/y) when impart(y)=0, 866 cos(i * ~~x / ~~y) => cosh(x/y) when impart(y)=0, 867 sinh(i * ~~x / ~~y) => i * sin(x/y) when impart(y)=0, 868 cosh(i * ~~x / ~~y) => cos(x/y) when impart(y)=0, 869 asin(i * ~~x / ~~y) => i * asinh(x/y) when impart(y)=0, 870 atan(i * ~~x / ~~y) => i * atanh(x/y) when impart(y)=0 871 and not((x=1 or x=-1) and y=1), 872 asinh(i * ~~x / ~~y) => i * asin(x/y) when impart(y)=0, 873 atanh(i * ~~x / ~~y) => i * atan(x/y) when impart(y)=0 874}$ 875 876let trig_imag_rules; 877 878% Generalized periodicity rules for trigonometric functions. 879% FJW, 16 October 1996. 880% exp rule corrected and others generalised to work for negative n 881% by AB March 2015 (negative n used to give error) 882 883operator arbint; 884 885let { 886 cos(~n*pi*arbint(~i) + ~~x) => cos((if evenp n then 0 else 1)*pi*arbint(i) + x) 887 when fixp n, 888 sin(~n*pi*arbint(~i) + ~~x) => sin((if evenp n then 0 else 1)*pi*arbint(i) + x) 889 when fixp n, 890 tan(~n*pi*arbint(~i) + ~~x) => tan(x) when fixp n, 891 sec(~n*pi*arbint(~i) + ~~x) => sec((if evenp n then 0 else 1)*pi*arbint(i) + x) 892 when fixp n, 893 csc(~n*pi*arbint(~i) + ~~x) => csc((if evenp n then 0 else 1)*pi*arbint(i) + x) 894 when fixp n, 895 cot(~n*pi*arbint(~i) + ~~x) => cot(x) when fixp n, 896 exp(~n*i*pi*arbint(~k) + ~~x) => 897 exp((if evenp n then 0 else 1)*i*pi*arbint(k) + x) when fixp n 898}; 899 900endmodule; 901 902end; 903 904