4 (trig : [cos(x), sin(x), tan(x), csc(x), sec(x), cot(x)],0);
7 (htrig : [cosh(x), sinh(x), tanh(x), csch(x), sech(x), coth(x)],0);
10 (complexfloatp(x) := block([xr,xi],
14 (?floatp(xr) and (?floatp(xi) or integerp(xi))) or (?floatp(xi) and integerp(xr))), 0);
17 /* Make sure that each trig function evaluates to a float for a float argument. */
22 (for f in append(trig,htrig) do (
23 ok : ok and complexfloatp(subst(x = 1.2 + %i,f)),
24 ok : ok and complexfloatp(subst(x = 1 + %i/7.0,f)),
25 ok : ok and complexfloatp(subst(x = -2.3 + %i/7.0,f)),
26 ok : ok and complexfloatp(subst(x = -2.3 - %i,f))),0);
32 (complexbigfloatp(x) := block([xr,xi],
36 (bfloatp(xr) and (bfloatp(xi) or integerp(xi))) or (bfloatp(xi) and integerp(xr))), 0);
39 /* Make sure that each trig function evaluates to a big float for a big float argument. */
44 (for f in append(trig,htrig) do (
45 ok : ok and complexbigfloatp(subst(x = 1.2b0,f))),0);
53 (maxerror : 0.0, %piargs : true, save_flags_names : '[%piargs, %iargs, trigsign, listarith, exponentialize, halfangles, trigexpand, trigexpandplus, trigexpandtimes, triginverses], save_flags_values : ev (save_flags_names), 0);
57 for i : -24 thru 24 do (
58 z : errcatch(subst(x = %pi * i / 24, f)),
59 if z # [ ] then maxerror : max(maxerror, abs(float(first(z) - subst(x = %pi * i / 24.0, f)))))),
60 if is(maxerror < 1.0e-13) then true else maxerror);
65 (%iargs : true, trigsign : true, 0);
68 subst(-%i*x,x, subst(%i*x,x,append(trig,htrig))) - append(trig,htrig);
69 ''(makelist(0,i,1,length(append(trig,htrig))));
71 (itrig : subst(%i*x,x,append(trig,htrig)),0);
74 (%iargs : false, listarith : true,0);
77 (itrig : itrig - subst(%i*x,x,append(trig,htrig)),0);
80 ratsimp(taylor(itrig,x,0,5));
81 ''(makelist(0,i,1,length(append(trig,htrig))))$
83 (exponentialize : false,0);
86 taylor(exponentialize(append(trig,htrig)) - append(trig,htrig),x,0,5);
87 ''(makelist(taylor(0,x,0,5), i,1,length(append(trig,htrig))))$
90 /* Test reflection rules */
95 subst(-x,x, subst(-x,x, append(trig,htrig)));
96 ''(append(trig,htrig))$
98 (mtrig : subst(x=5.6, subst(-x,x,append(trig,htrig))),0);
101 (trigsign : false,0);
104 mtrig - subst(x=-5.6, append(trig,htrig));
105 ''(makelist(0.0,i,1,length(append(trig,htrig))))$
107 /* Test half angles */
109 /* Because we have generalized the half-angle-transformation, we have to
110 restrict the arguments to a positive interval (0,%pi). We take not %pi
111 but the number 3 as upper limit, because of problems with assume. */
113 (assume(x>0,x<3),halfangles : true, 0);
116 (xtrig : subst(x/2,x, append(trig,htrig)),0);
119 (halfangles : false, 0);
122 xtrig : taylor(xtrig - subst(x/2,x, append(trig, htrig)),x,0,5);
123 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
128 /* Test trig expand */
130 (trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
133 (xtrig : subst(x=x+y, append(trig,htrig)),0);
136 (trigexpand : false,0);
139 xtrig : taylor(xtrig - subst(x=x+y, append(trig, htrig)),[x,y],0,5);
140 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
142 (trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
145 (xtrig : subst(x=x-y, append(trig,htrig)),0);
148 (trigexpand : false,0);
151 xtrig : taylor(xtrig - subst(x=x-y, append(trig, htrig)),[x,y],0,5);
152 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
154 (trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
157 (xtrig : subst(x = 2*x, append(trig,htrig)),0);
160 (trigexpand : false,0);
163 xtrig : taylor(xtrig - subst(x=2*x, append(trig, htrig)),[x,y],0,5);
164 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
166 (trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
169 (xtrig : subst(x = 3*x, append(trig,htrig)),0);
172 (trigexpand : false,0);
175 xtrig : taylor(xtrig - subst(x=3*x, append(trig, htrig)),[x,y],0,5);
176 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
178 (trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
181 (xtrig : subst(x = 7*x, append(trig,htrig)),0);
184 (trigexpand : false,0);
187 xtrig : taylor(xtrig - subst(x=7*x, append(trig, htrig)),[x,y],0,5);
188 ''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$
190 /* Test triginverses */
192 (invtrig : [acos(x), asin(x), atan(x), acsc(x), asec(x), acot(x)],0);
196 for f in invtrig do (
198 xtrig : subst(f,x,trig),
199 triginverses : false,
200 xtrig : xtrig - subst(f,x,trig),
202 /* Paste expressions into a lambda to protect them from simplification and evaluation.
203 * Evaluate the lambda (no arguments) to get the test result.
205 xtrig : buildq ([L : map (nounify (cabs), xtrig)], lambda ([], L)),
206 for i : 1 thru 9 do block ([xi : i/10.0, xtrigi],
207 xtrigi : apply (buildq, [[funmake (":", ['x, xi])], xtrig]),
208 if apply ('max, ev (xtrigi(), nouns)) >= 1.0e-12
209 then buggy : cons (xtrigi, buggy))),
214 (invhtrig : [acosh(x), asinh(x), atanh(x), acsch(x), asech(x), acoth(x)],0);
218 for f in invhtrig do (
220 xtrig : subst(f,x,htrig),
221 triginverses : false,
222 xtrig : xtrig - subst(f,x,htrig),
224 /* Paste expressions into a lambda to protect them from simplification and evaluation.
225 * Evaluate the lambda (no arguments) to get the test result.
227 xtrig : buildq ([L : map (nounify (cabs), xtrig)], lambda ([], L)),
228 for i : 1 thru 9 do block ([xi : i/10.0, xtrigi],
229 xtrigi : apply (buildq, [[funmake (":", ['x, xi])], xtrig]),
230 if apply ('max, ev (xtrigi(), nouns)) >= 1.0e-12
231 then buggy : cons (xtrigi, buggy))),
236 (alltrig : append(trig, htrig, invtrig, invhtrig),0);
239 (pts : [2,2+%i, 2-%i,-2,-2+%i,-2-%i,1/2,1/2 +%i/2,1/2-%i/2,-1/2,-1/2+%i/2,-1/2-%i/2],0);
244 for p in pts do block ([e, fop : op(f)],
245 e : buildq ([p, fop], lambda ([], cabs (float (rectform (fop (p))) - fop (float (p))))),
247 then buggy : cons (e, buggy)), buggy);
250 /* Like above, but for big floats. We expect to be within 3 digits of the answer.
252 Note that we might fail because the rectform and the function might
253 produce answers because of the branch cut!
257 for p in pts do block ([e, fop : op(f)],
258 e : buildq ([p, fop], lambda ([], cabs (bfloat (rectform (fop (p))) - fop (bfloat (p))))),
259 if e() > 10b0^(3-fpprec)
260 then buggy : cons (e, buggy)), buggy);
264 /* Numerically check some identities */
266 (float_trig_test(id, threshold) :=
267 block([maxerror : 0],
268 for p in pts do block([numer : true],
269 maxerror : max(maxerror, cabs(expand(rectform(subst(x = p, id)))))),
270 if is(maxerror < threshold) then true else maxerror),
271 bfloat_trig_test(id, threshold) :=
272 block([maxerror : 0],
273 for p in pts do block([numer : true],
274 maxerror : max(maxerror, cabs(expand(subst(x = bfloat(p), id))))),
275 if is(maxerror < threshold) then true else maxerror),
279 float_trig_test(cos(x)^2 + sin(x)^2 - 1, 1.0e-15);
282 float_trig_test(cot(2*x) - (cot(x)^2-1)/(2*cot(x)), 1.0e-15);
285 bfloat_trig_test(asin(x) + acos(x) - %pi/2, 1.0e-15);
288 (declare(n,integer),0);
291 /* see SF bug 1754072 */
293 subst(n = 1, cos(%pi * n * 31 / 37));
296 subst(n = 1, sin(%pi * n * 31 / 37));
299 subst(n = 1, tan(%pi * n * 31 / 37));
302 /* [ 580721 ] trigexpand bug */
303 tan(%pi*x+%pi/2),trigexpand;
309 /* tan(%pi*integer) simplification - ID: 3058290 */
316 /* [ 1553866 ] %piargs inconsistent behavior */
320 (remfunction(complexfloatp, complexbigfloatp),0);
323 /* Probably we should beef up reset slightly to handle this, but for now this is OK. */
324 (map (lambda ([a, b], a :: b), save_flags_names, save_flags_values), 0);
327 /* [ 1644575 ] acot(0.0) vs acot(0) */
331 /* [ 620246 ] carg(complex) */
332 (declare(u, complex),0);
336 log(abs(u)) + %i*carg(u);
339 abs(u); /* for a symbol we get a noun form with abs */
340 /*sqrt(?%realpart(u)^2+?%imagpart(u)^2);*/
342 /* [ 617699 ] carg([1]) is bogus */
346 /* Tests for half-angles simplification */
351 (old_halfangles:halfangles,halfangles:true,old_%iargs:%iargs,%iargs:false,done);
354 /* Sin function: The general result for real argument */
357 (-1)^floor(x/(2*%pi))*sqrt(1-cos(x))/sqrt(2);
359 /* The square simplifies correctly */
364 /* Correct sign for negative and positive real argument */
366 (assume(x1>0,x1<2*%pi),done);
370 sqrt(1-cos(x1))/sqrt(2);
372 (assume(x2>-2*%pi,x2<0),done);
376 -sqrt(1-cos(x2))/sqrt(2);
378 /* Correct sign for pure imaginary argument
379 and complex argument with realpart a multiple of 2*%pi */
381 (assume(y1>0,y2<0),done);
385 sqrt(1-cos(%i*y1))/sqrt(2);
388 -sqrt(1-cos(%i*y2))/sqrt(2);
390 sin((2*%pi+%i*y2)/2);
391 sqrt(1-cos(2*%pi+%i*y2))/sqrt(2);
393 sin((2*%pi+%i*y1)/2);
394 -sqrt(1-cos(2*%pi+%i*y1))/sqrt(2);
396 /* Simplification for negative or positive imaginary part */
397 (assume(yneg<0,ypos>0),done);
401 sqrt(1-cos(x1+%i*yneg))/sqrt(2);
403 sin((2*%pi+%i*yneg)/2);
404 sqrt(1-cos(2*%pi+%i*yneg))/sqrt(2);
407 sqrt(1-cos(x1+%i*ypos))/sqrt(2);
409 sin((2*%pi+%i*ypos)/2);
410 -sqrt(1-cos(2*%pi+%i*ypos))/sqrt(2); /* for this case an -1 */
412 /* Cos function: The general result for real argument */
415 (-1)^floor((x+%pi)/(2*%pi))*sqrt(1+cos(x))/sqrt(2);
417 /* The square simplifies correctly */
422 /* Correct sign for real argument in (-%pi,%pi) */
424 (forget(x1<2*%pi,x1>0),assume(x1>-%pi,x1<%pi),done);
428 sqrt(1+cos(x1))/sqrt(2);
430 /* Correct sign for pure imaginary argument
431 and complex argument with realpart a multiple of %pi */
434 sqrt(1+cos(%i*y1))/sqrt(2);
437 sqrt(1+cos(%i*y2))/sqrt(2);
440 sqrt(1+cos(%pi+%i*y2))/sqrt(2);
443 -sqrt(1+cos(%pi+%i*y1))/sqrt(2);
445 /* Simplification for negative or positive imaginary part */
446 (assume(yneg<0,ypos>0),done);
450 sqrt(1+cos(x1+%i*yneg))/sqrt(2);
452 cos((%pi+%i*yneg)/2);
453 sqrt(1+cos(%pi+%i*yneg))/sqrt(2);
456 sqrt(1+cos(x1+%i*ypos))/sqrt(2);
458 cos((%pi+%i*ypos)/2);
459 -sqrt(1+cos(%pi+%i*ypos))/sqrt(2); /* for this case an -1 */
461 /* Sinh function with Real arguments */
463 (assume(xpos>0,xneg<0),done);
467 abs(x)/x*sqrt((cosh(x)-1)/2);
473 sqrt((cosh(xpos)-1)/2);
476 -sqrt((cosh(xneg)-1)/2);
478 /* Sinh function with Complex arguments */
480 /* x1 is in (-%pi,%pi) */
482 abs(x1)/x1*sqrt((cosh(%i*x1)-1)/2);
484 sinh((xpos+%i*x1)/2);
485 sqrt((xpos+%i*x1)^2)/(xpos+%i*x1)*sqrt((cosh(xpos+%i*x1)-1)/2);
487 sinh((xneg+%i*x1)/2);
488 sqrt((xneg+%i*x1)^2)/(xneg+%i*x1)*sqrt((cosh(xneg+%i*x1)-1)/2);
490 /* Cosh function with Real arguments */
499 sqrt((cosh(xpos)+1)/2);
502 sqrt((cosh(xneg)+1)/2);
504 /* Cosh function with Complex arguments */
506 /* x1 is in (-%pi,%pi) */
508 sqrt((cosh(%i*x1)+1)/2);
510 cosh((xpos+%i*x1)/2);
511 sqrt((cosh(xpos+%i*x1)+1)/2);
513 cosh((xneg+%i*x1)/2);
514 sqrt((cosh(xneg+%i*x1)+1)/2);
516 /* two tests for atan2(sqrt(1-u)*(u-1),1) - ID: 3521596 */
517 atan2(sqrt(1-u)*(u-1),1);
518 atan2(sqrt(1-u)*(u-1),1)$
520 atan2((1-u)^(1/2^10)*(u-1),1);
521 atan2((1-u)^(1/2^10)*(u-1),1)$
523 /* adding tests for some previously uncovered code in trigi.lisp and trigo.lisp */
525 block([triginverses : true], [sin(atan2(y,x)), cos(atan2(y,x)),tan(atan2(y,x)), cot(atan2(y,x)), csc(atan2(y,x)), sec(atan2(y,x))]);
526 [y/sqrt(y^2+x^2), x/sqrt(y^2+x^2), y/x, x/y, sqrt(y^2+x^2)/y, sqrt(y^2+x^2)/x]$
528 block([exponentialize : true], [sec(x),csch(x), sech(x)]);
529 [2/(%e^(%i*x)+%e^-(%i*x)), 2/(%e^x-%e^-x), 2/(%e^x+%e^-x)]$
531 [atan(inf), atan(-minf), atan(minf), atan(-inf)];
532 [%pi/2, %pi/2, -%pi/2, -%pi/2]$
534 block([%piargs : true], [atan(-1 + sqrt(2)), atan(1+sqrt(2))]);
537 block([%iargs : true], atan(%i*x));
540 [atan(tan(2015)), atan(tan(-2014)), atan(cot(-2014)), atan(cot(2015))];
541 [2015-641*%pi,641*%pi-2014,2014-1283*%pi/2,1283*%pi/2-2015]$
549 map(lambda([x], asin(sin(x)) - (-min(2*%pi*floor(x/(2*%pi))-x+%pi,-2*%pi*floor(x/(2*%pi))+x-2*%pi)-max(x-2*%pi*floor(x/(2*%pi)),2*%pi*floor(x/(2*%pi))-x+%pi)-2*%pi*floor(x/(2*%pi))+x-%pi)), makelist(i/3,i,0,43));
550 ''(makelist(0,i,0,43))$
552 map(lambda([x], asec(sec(x)) - (%pi-abs(2*%pi*floor(x/(2*%pi))-x+%pi))), makelist(i/19,i,0,25));
553 ''(makelist(0,i,0,25))$
561 block([triginverses : true], asin(sin(1/5)));
573 block([triginverses : all],
574 [acsc(csc(x)), acot(cot(x)), acos(cos(x)), asinh(sinh(x)), acosh(cosh(x)), atanh(tanh(x)),acoth(coth(x)), acsch(csch(x)), asech(sech(x))]);
575 [x,x, x, x, x, x, x, x, x]$
577 [acot(0), acot(1),acot(-1),acot(1/sqrt(3)), acot(sqrt(3))];
578 [%pi/2,%pi/4,-%pi/4,%pi/3,%pi/6]$
580 [atan(inf), atan(1),atan(-1),atan(sqrt(3)), atan(1/sqrt(3))];
581 [%pi/2,%pi/4,-%pi/4,%pi/3,%pi/6]$
583 [acos(sqrt(3)/2), acos(-sqrt(3)/2)];
586 [acsc(1), acsc(-1), acsc(sqrt(2)), acsc(2/sqrt(3))];
587 [%pi/2,-%pi/2,%pi/4,%pi/3]$
589 [asec(-1),asec(sqrt(2)), 2/sqrt(3)];
590 [%pi,%pi/4,2/sqrt(3)]$
592 /* end of tests for some previously uncovered code in trigi.lisp and trigo.lisp */
594 /* Restore global flags */
595 (halfangles:old_halfangles,%iargs:old_%iargs,done);
598 /* SF bug # 2570: "Make acos(cos(x)) simplify to x when on correct interval" */
604 (ctxt : newcontext (),
605 assume (0 <= foo and foo <= %pi),
612 /* mailing list 2015-08-25: "Strange symmetry of acoth(x), area cotangens hyperbolicus function (#552)" */
614 (kill (x), - acoth (- x));
617 /* SF bug #2620: "atan2(y,x)+atan2(-y,x) doesn't always return 0 " */
619 (kill (x, y), atan2(y,x)+atan2(-y,x));
622 (assume (y > 0), atan2(y,x)+atan2(-y,x));
625 (kill (n, p, r), forget (y > 0), assume(n<0,p>0), atan2 (- x, r));
634 forget (n < 0, p > 0);
637 /* mailing list 2017-11-27: "trigsimp fails with pderivop" */
641 expr:pderivop(f,1)((y))^2);
647 (reset (use_pdiff), 0); /* disables pdiff */
651 trigsimp (foo(1)(x)^3));