Bug fix: gcfactor(x*%i) => lisp error
[maxima.git] / tests / rtest_trig.mac
blob5f85fa7a27205f6c4c5685c25e82c40c998ce858
1 (kill (all), 0);
2 0;
4 (trig : [cos(x), sin(x), tan(x), csc(x), sec(x), cot(x)],0);
5 0$
7 (htrig : [cosh(x), sinh(x), tanh(x), csch(x), sech(x), coth(x)],0);
8 0$
10 (complexfloatp(x) := block([xr,xi],
11   x : rectform(x),
12   xr : realpart(x),
13   xi : imagpart(x),
14   (?floatp(xr) and (?floatp(xi) or integerp(xi))) or (?floatp(xi) and integerp(xr))), 0);
16   
17 /*  Make sure that each trig function evaluates to a float for a float argument. */
19 (ok : true,0);
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);
29 ok;
30 true$
32 (complexbigfloatp(x) := block([xr,xi],
33   x : rectform(x),
34   xr : realpart(x),
35   xi : imagpart(x),
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. */
41 (ok : true,0);
44 (for f in append(trig,htrig) do (
45   ok : ok and complexbigfloatp(subst(x = 1.2b0,f))),0);
48 ok;
49 true$
51 /*  Test %piargs  */
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);
56 (for f in trig do (
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);
61 true$
63 /* Test %iargs   */
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  */
92 (trigsign : true,0);
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 intervall (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))))$
125 (forget(x>0,x<3),0);
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);
195 block ([buggy : []],
196   for f in invtrig do (
197     triginverses : true,
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.
204      */
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))),
211   buggy);
214 (invhtrig : [acosh(x), asinh(x), atanh(x), acsch(x), asech(x), acoth(x)],0);
217 block ([buggy : []],
218   for f in invhtrig do (
219     triginverses : true,
220     xtrig : subst(f,x,htrig),
221     triginverses : false,
222     xtrig : xtrig - subst(f,x,htrig),
223   
224     /* Paste expressions into a lambda to protect them from simplification and evaluation.
225      * Evaluate the lambda (no arguments) to get the test result.
226      */
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))),
233   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);
242 block ([buggy : []],
243   for f in alltrig do
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))))),
246       if e() > 1.0e-13
247         then buggy : cons (e, buggy)), buggy);
250 /* Like above, but for big floats.  We expect to be within 3 digits of the answer.
251    
252    Note that we might fail because the rectform and the function might
253    produce answers because of the branch cut!
254  */
255 block ([buggy : []],
256   for f in alltrig do
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),
276  0);
279 float_trig_test(cos(x)^2 + sin(x)^2 - 1, 1.0e-15);
280 true$
282 float_trig_test(cot(2*x) - (cot(x)^2-1)/(2*cot(x)), 1.0e-15);
283 true$
285 bfloat_trig_test(asin(x) + acos(x) - %pi/2, 1.0e-15);
286 true$
288 (declare(n,integer),0);
291 /* see SF bug 1754072 */
293 subst(n = 1, cos(%pi * n * 31 / 37));
294 cos(%pi * 31 / 37);
296 subst(n = 1, sin(%pi * n * 31 / 37));
297 sin(%pi * 31 / 37);
299 subst(n = 1, tan(%pi * n * 31 / 37));
300 tan(%pi * 31 / 37);
302 /* [ 580721 ] trigexpand bug */
303 tan(%pi*x+%pi/2),trigexpand;
304 -cot(%pi*x);
306 tan(n*2*%pi);
309 /* tan(%pi*integer) simplification - ID: 3058290 */
310 tan(n*%pi);
313 (kill(n,declare),0);
316 /* [ 1553866 ] %piargs inconsistent behavior */
317 cos(%pi*x + %pi);
318 -cos(%pi*x);
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) */
328 acot(0.0);
329 1.570796326794897;
331 /* [ 620246 ] carg(complex) */
332 (declare(u, complex),0);
335 rectform(log(u));
336 log(abs(u)) + %i*carg(u);
338 cabs(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 */
343 carg([1]);
344 [0];
346 /* Tests for half-angles simplification */
348 kill(all);
349 done;
351 (old_halfangles:halfangles,halfangles:true,old_%iargs:%iargs,%iargs:false,done);
352 done;
354 /* Sin function: The general result for real argument */
356 sin(x/2);
357 (-1)^floor(x/(2*%pi))*sqrt(1-cos(x))/sqrt(2);
359 /* The square simplifies correctly */
361 sin(x/2)^2;
362 (1-cos(x))/2;
364 /* Correct sign for negative and positive real argument */
366 (assume(x1>0,x1<2*%pi),done);
367 done;
369 sin(x1/2);
370 sqrt(1-cos(x1))/sqrt(2);
372 (assume(x2>-2*%pi,x2<0),done);
373 done;
375 sin(x2/2);
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);
382 done;
384 sin(%i*y1/2);
385 sqrt(1-cos(%i*y1))/sqrt(2);
387 sin(%i*y2/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);
398 done;
400 sin((x1+%i*yneg)/2);
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);
406 sin((x1+%i*ypos)/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 */
414 cos(x/2);
415 (-1)^floor((x+%pi)/(2*%pi))*sqrt(1+cos(x))/sqrt(2);
417 /* The square simplifies correctly */
419 cos(x/2)^2;
420 (1+cos(x))/2;
422 /* Correct sign for real argument in (-%pi,%pi) */
424 (forget(x1<2*%pi,x1>0),assume(x1>-%pi,x1<%pi),done);
425 done;
427 cos(x1/2);
428 sqrt(1+cos(x1))/sqrt(2);
430 /* Correct sign for pure imaginary argument
431    and complex argument with realpart a multiple of %pi */
433 cos(%i*y1/2);
434 sqrt(1+cos(%i*y1))/sqrt(2);
436 cos(%i*y2/2);
437 sqrt(1+cos(%i*y2))/sqrt(2);
439 cos((%pi+%i*y2)/2);
440 sqrt(1+cos(%pi+%i*y2))/sqrt(2);
442 cos((%pi+%i*y1)/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);
447 done;
449 cos((x1+%i*yneg)/2);
450 sqrt(1+cos(x1+%i*yneg))/sqrt(2);
452 cos((%pi+%i*yneg)/2);
453 sqrt(1+cos(%pi+%i*yneg))/sqrt(2);
455 cos((x1+%i*ypos)/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);
464 done;
466 sinh(x/2);
467 abs(x)/x*sqrt((cosh(x)-1)/2);
469 sinh(x/2)^2;
470 (cosh(x)-1)/2;
472 sinh(xpos/2);
473 sqrt((cosh(xpos)-1)/2);
475 sinh(xneg/2);
476 -sqrt((cosh(xneg)-1)/2);
478 /* Sinh function with Complex arguments */
480 /* x1 is in (-%pi,%pi) */
481 sinh(%i*x1/2);
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 */
492 cosh(x/2);
493 sqrt((cosh(x)+1)/2);
495 cosh(x/2)^2;
496 (cosh(x)+1)/2;
498 cosh(xpos/2);
499 sqrt((cosh(xpos)+1)/2);
501 cosh(xneg/2);
502 sqrt((cosh(xneg)+1)/2);
504 /* Cosh function with Complex arguments */
506 /* x1 is in (-%pi,%pi) */
507 cosh(%i*x1/2);
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))]);
535 [%pi/8, 3*%pi/8]$
537 block([%iargs : true], atan(%i*x));
538 %i*atanh(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]$
543 asin(sqrt(3)/2);
544 %pi/3$
546 asin(-sqrt(3)/2);
547 -%pi/3$
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))$
555 sin(641*%pi-2015);
556 sin(2015)$
558 asin(cos(2015));
559 2015-1283*%pi/2$
561 block([triginverses : true], asin(sin(1/5)));
562 1/5$
564 acos(sin(1066));
565 1066-677*%pi/2$
567 acot(cot(42));
568 42-13*%pi$
570 acsc(csc(191/8));
571 191/8-8*%pi$
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)];
584 [%pi/6,5*%pi/6]$
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);
596 done;
598 /* SF bug # 2570: "Make acos(cos(x)) simplify to x when on correct interval" */
600 (kill (foo),
601  acos (cos (foo)));
602 acos (cos (foo));
604 (ctxt : newcontext (),
605  assume (0 <= foo and foo <= %pi),
606  acos (cos (foo)));
607 foo;
609 killcontext (ctxt);
610 done;
612 /* mailing list 2015-08-25: "Strange symmetry of acoth(x), area cotangens hyperbolicus function (#552)" */
614 (kill (x), - acoth (- x));
615 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));
626 - atan2 (x, r);
628 atan2 (-n, r);
629 - atan2 (n, r);
631 atan2 (-p, r);
632 - atan2 (p, r);
634 forget (n < 0, p > 0);
635 [n < 0, p > 0];
637 /* mailing list 2017-11-27: "trigsimp fails with pderivop" */
639 (load (pdiff),
640  kill (f, y),
641  expr:pderivop(f,1)((y))^2);
642 pderivop(f,1)(y)^2;
644 trigsimp (expr);
645 pderivop(f,1)(y)^2;
647 (reset (use_pdiff), 0); /* disables pdiff */
650 (kill (foo, x),
651  trigsimp (foo(1)(x)^3));
652 foo(1)(x)^3;