(kill (all), 0);
0;

(trig : [cos(x), sin(x), tan(x), csc(x), sec(x), cot(x)],0);
0$

(htrig : [cosh(x), sinh(x), tanh(x), csch(x), sech(x), coth(x)],0);
0$

(complexfloatp(x) := block([xr,xi],
  x : rectform(x),
  xr : realpart(x),
  xi : imagpart(x),
  (?floatp(xr) and (?floatp(xi) or integerp(xi))) or (?floatp(xi) and integerp(xr))), 0);
0$
  
/*  Make sure that each trig function evaluates to a float for a float argument. */

(ok : true,0);
0$

(for f in append(trig,htrig) do (
  ok : ok and complexfloatp(subst(x = 1.2 + %i,f)),
  ok : ok and complexfloatp(subst(x = 1 + %i/7.0,f)),
  ok : ok and complexfloatp(subst(x = -2.3 + %i/7.0,f)),
  ok : ok and complexfloatp(subst(x = -2.3 - %i,f))),0);
0$

ok;
true$

(complexbigfloatp(x) := block([xr,xi],
  x : rectform(x),
  xr : realpart(x),
  xi : imagpart(x),
  (bfloatp(xr) and (bfloatp(xi) or integerp(xi))) or (bfloatp(xi) and integerp(xr))), 0);
0$

/*  Make sure that each trig function evaluates to a big float for a big float argument. */

(ok : true,0);
0$

(for f in append(trig,htrig) do (
  ok : ok and complexbigfloatp(subst(x = 1.2b0,f))),0);
0$

ok;
true$

/*  Test %piargs  */

(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);
0;

(for f in trig do (
   for i : -24 thru 24 do (
      z : errcatch(subst(x = %pi * i / 24, f)),
      if z # [ ] then maxerror : max(maxerror, abs(float(first(z) - subst(x = %pi * i / 24.0, f)))))),
 if is(maxerror < 1.0e-13) then true else maxerror);
true$

/* Test %iargs   */

(%iargs : true, trigsign : true, 0);
0$

subst(-%i*x,x, subst(%i*x,x,append(trig,htrig))) - append(trig,htrig);
''(makelist(0,i,1,length(append(trig,htrig))));

(itrig : subst(%i*x,x,append(trig,htrig)),0);
0$

(%iargs : false, listarith : true,0);
0$

(itrig : itrig - subst(%i*x,x,append(trig,htrig)),0);
0$

ratsimp(taylor(itrig,x,0,5));
''(makelist(0,i,1,length(append(trig,htrig))))$

(exponentialize : false,0);
0$

taylor(exponentialize(append(trig,htrig)) - append(trig,htrig),x,0,5);
''(makelist(taylor(0,x,0,5), i,1,length(append(trig,htrig))))$


/*  Test reflection rules  */

(trigsign : true,0);
0$

subst(-x,x, subst(-x,x, append(trig,htrig)));
''(append(trig,htrig))$

(mtrig : subst(x=5.6, subst(-x,x,append(trig,htrig))),0);
0$

(trigsign : false,0);
0$

mtrig - subst(x=-5.6, append(trig,htrig));
''(makelist(0.0,i,1,length(append(trig,htrig))))$

/* Test half angles */

/* Because we have generalized the half-angle-transformation, we have to
   restrict the arguments to a positive interval (0,%pi). We take not %pi
   but the number 3 as upper limit, because of problems with assume. */

(assume(x>0,x<3),halfangles : true, 0);
0$

(xtrig : subst(x/2,x, append(trig,htrig)),0);
0$

(halfangles : false, 0);
0$

xtrig : taylor(xtrig - subst(x/2,x, append(trig, htrig)),x,0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

(forget(x>0,x<3),0);
0;

/* Test trig expand */

(trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
0$

(xtrig : subst(x=x+y, append(trig,htrig)),0);
0$

(trigexpand : false,0);
0$

xtrig : taylor(xtrig - subst(x=x+y, append(trig, htrig)),[x,y],0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

(trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
0$

(xtrig : subst(x=x-y, append(trig,htrig)),0);
0$

(trigexpand : false,0);
0$

xtrig : taylor(xtrig - subst(x=x-y, append(trig, htrig)),[x,y],0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

(trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
0$

(xtrig : subst(x = 2*x, append(trig,htrig)),0);
0$

(trigexpand : false,0);
0$

xtrig : taylor(xtrig - subst(x=2*x, append(trig, htrig)),[x,y],0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

(trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
0$

(xtrig : subst(x = 3*x, append(trig,htrig)),0);
0$

(trigexpand : false,0);
0$

xtrig : taylor(xtrig - subst(x=3*x, append(trig, htrig)),[x,y],0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

(trigexpand : true, trigexpandplus : true, trigexpandtimes : true, 0);
0$

(xtrig : subst(x = 7*x, append(trig,htrig)),0);
0$

(trigexpand : false,0);
0$

xtrig : taylor(xtrig - subst(x=7*x, append(trig, htrig)),[x,y],0,5);
''(makelist(taylor(0,x,0,5),i,1,length(append(trig,htrig))))$

/* Test triginverses  */

(invtrig : [acos(x), asin(x), atan(x), acsc(x), asec(x), acot(x)],0);
0$

block ([buggy : []],
  for f in invtrig do (
    triginverses : true,
    xtrig : subst(f,x,trig),
    triginverses : false,
    xtrig : xtrig - subst(f,x,trig),

    /* Paste expressions into a lambda to protect them from simplification and evaluation.
     * Evaluate the lambda (no arguments) to get the test result.
     */
    xtrig : buildq ([L : map (nounify (cabs), xtrig)], lambda ([], L)),
    for i : 1 thru 9 do block ([xi : i/10.0, xtrigi],
      xtrigi : apply (buildq, [[funmake (":", ['x, xi])], xtrig]),
      if apply ('max, ev (xtrigi(), nouns)) >= 1.0e-12
        then buggy : cons (xtrigi, buggy))),

  buggy);
[]$

(invhtrig : [acosh(x), asinh(x), atanh(x), acsch(x), asech(x), acoth(x)],0);
0$

block ([buggy : []],
  for f in invhtrig do (
    triginverses : true,
    xtrig : subst(f,x,htrig),
    triginverses : false,
    xtrig : xtrig - subst(f,x,htrig),
  
    /* Paste expressions into a lambda to protect them from simplification and evaluation.
     * Evaluate the lambda (no arguments) to get the test result.
     */
    xtrig : buildq ([L : map (nounify (cabs), xtrig)], lambda ([], L)),
    for i : 1 thru 9 do block ([xi : i/10.0, xtrigi],
      xtrigi : apply (buildq, [[funmake (":", ['x, xi])], xtrig]),
      if apply ('max, ev (xtrigi(), nouns)) >= 1.0e-12
        then buggy : cons (xtrigi, buggy))),

  buggy);
[]$

(alltrig : append(trig, htrig, invtrig, invhtrig),0);
0$

(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);
0$

block ([buggy : []],
  for f in alltrig do
    for p in pts do block ([e, fop : op(f)],
      e : buildq ([p, fop], lambda ([], cabs (float (rectform (fop (p))) - fop (float (p))))),
      if e() > 1.0e-13
        then buggy : cons (e, buggy)), buggy);
[]$

/* Like above, but for big floats.  We expect to be within 3 digits of the answer.
   
   Note that we might fail because the rectform and the function might
   produce answers because of the branch cut!
 */
block ([buggy : []],
  for f in alltrig do
    for p in pts do block ([e, fop : op(f)],
      e : buildq ([p, fop], lambda ([], cabs (bfloat (rectform (fop (p))) - fop (bfloat (p))))),
      if e() > 10b0^(3-fpprec)
        then buggy : cons (e, buggy)), buggy);
[]$


/*  Numerically check some identities  */

(float_trig_test(id, threshold) :=
  block([maxerror : 0],
    for p in pts do block([numer : true],
      maxerror : max(maxerror, cabs(expand(rectform(subst(x = p, id)))))),
    if is(maxerror < threshold) then true else maxerror),
 bfloat_trig_test(id, threshold) :=
  block([maxerror : 0],
    for p in pts do block([numer : true],
      maxerror : max(maxerror, cabs(expand(subst(x = bfloat(p), id))))),
    if is(maxerror < threshold) then true else maxerror),
 0);
0;

float_trig_test(cos(x)^2 + sin(x)^2 - 1, 1.0e-15);
true$

float_trig_test(cot(2*x) - (cot(x)^2-1)/(2*cot(x)), 1.0e-15);
true$

bfloat_trig_test(asin(x) + acos(x) - %pi/2, 1.0e-15);
true$

(declare(n,integer),0);
0$

/* see SF bug 1754072 */

subst(n = 1, cos(%pi * n * 31 / 37));
cos(%pi * 31 / 37);

subst(n = 1, sin(%pi * n * 31 / 37));
sin(%pi * 31 / 37);

subst(n = 1, tan(%pi * n * 31 / 37));
tan(%pi * 31 / 37);

/* [ 580721 ] trigexpand bug */
tan(%pi*x+%pi/2),trigexpand;
-cot(%pi*x);

tan(n*2*%pi);
0;

/* tan(%pi*integer) simplification - ID: 3058290 */
tan(n*%pi);
0;

(kill(n,declare),0);
0$

/* [ 1553866 ] %piargs inconsistent behavior */
cos(%pi*x + %pi);
-cos(%pi*x);

(remfunction(complexfloatp, complexbigfloatp),0);
0$

/* Probably we should beef up reset slightly to handle this, but for now this is OK. */
(map (lambda ([a, b], a :: b), save_flags_names, save_flags_values), 0);
0;

/* [ 1644575 ] acot(0.0) vs acot(0) */
acot(0.0);
1.570796326794897;

/* [ 620246 ] carg(complex) */
(declare(u, complex),0);
0;

rectform(log(u));
log(abs(u)) + %i*carg(u);

cabs(u);
abs(u); /* for a symbol we get a noun form with abs */
/*sqrt(?%realpart(u)^2+?%imagpart(u)^2);*/

/* [ 617699 ] carg([1]) is bogus */
carg([1]);
[0];

/* Tests for half-angles simplification */

kill(all);
done;

(old_halfangles:halfangles,halfangles:true,old_%iargs:%iargs,%iargs:false,done);
done;

/* Sin function: The general result for real argument */

sin(x/2);
(-1)^floor(x/(2*%pi))*sqrt(1-cos(x))/sqrt(2);

/* The square simplifies correctly */

sin(x/2)^2;
(1-cos(x))/2;

/* Correct sign for negative and positive real argument */

(assume(x1>0,x1<2*%pi),done);
done;

sin(x1/2);
sqrt(1-cos(x1))/sqrt(2);

(assume(x2>-2*%pi,x2<0),done);
done;

sin(x2/2);
-sqrt(1-cos(x2))/sqrt(2);

/* Correct sign for pure imaginary argument
   and complex argument with realpart a multiple of 2*%pi */

(assume(y1>0,y2<0),done);
done;

sin(%i*y1/2);
sqrt(1-cos(%i*y1))/sqrt(2);

sin(%i*y2/2);
-sqrt(1-cos(%i*y2))/sqrt(2);

sin((2*%pi+%i*y2)/2);
sqrt(1-cos(2*%pi+%i*y2))/sqrt(2);

sin((2*%pi+%i*y1)/2);
-sqrt(1-cos(2*%pi+%i*y1))/sqrt(2);

/* Simplification for negative or positive imaginary part */
(assume(yneg<0,ypos>0),done);
done;

sin((x1+%i*yneg)/2);
sqrt(1-cos(x1+%i*yneg))/sqrt(2);

sin((2*%pi+%i*yneg)/2);
sqrt(1-cos(2*%pi+%i*yneg))/sqrt(2);

sin((x1+%i*ypos)/2);
sqrt(1-cos(x1+%i*ypos))/sqrt(2);

sin((2*%pi+%i*ypos)/2);
-sqrt(1-cos(2*%pi+%i*ypos))/sqrt(2); /* for this case an -1 */

/* Cos function: The general result for real argument */

cos(x/2);
(-1)^floor((x+%pi)/(2*%pi))*sqrt(1+cos(x))/sqrt(2);

/* The square simplifies correctly */

cos(x/2)^2;
(1+cos(x))/2;

/* Correct sign for real argument in (-%pi,%pi) */

(forget(x1<2*%pi,x1>0),assume(x1>-%pi,x1<%pi),done);
done;

cos(x1/2);
sqrt(1+cos(x1))/sqrt(2);

/* Correct sign for pure imaginary argument
   and complex argument with realpart a multiple of %pi */

cos(%i*y1/2);
sqrt(1+cos(%i*y1))/sqrt(2);

cos(%i*y2/2);
sqrt(1+cos(%i*y2))/sqrt(2);

cos((%pi+%i*y2)/2);
sqrt(1+cos(%pi+%i*y2))/sqrt(2);

cos((%pi+%i*y1)/2);
-sqrt(1+cos(%pi+%i*y1))/sqrt(2);

/* Simplification for negative or positive imaginary part */
(assume(yneg<0,ypos>0),done);
done;

cos((x1+%i*yneg)/2);
sqrt(1+cos(x1+%i*yneg))/sqrt(2);

cos((%pi+%i*yneg)/2);
sqrt(1+cos(%pi+%i*yneg))/sqrt(2);

cos((x1+%i*ypos)/2);
sqrt(1+cos(x1+%i*ypos))/sqrt(2);

cos((%pi+%i*ypos)/2);
-sqrt(1+cos(%pi+%i*ypos))/sqrt(2); /* for this case an -1 */

/* Sinh function with Real arguments */

(assume(xpos>0,xneg<0),done);
done;

sinh(x/2);
abs(x)/x*sqrt((cosh(x)-1)/2);

sinh(x/2)^2;
(cosh(x)-1)/2;

sinh(xpos/2);
sqrt((cosh(xpos)-1)/2);

sinh(xneg/2);
-sqrt((cosh(xneg)-1)/2);

/* Sinh function with Complex arguments */

/* x1 is in (-%pi,%pi) */
sinh(%i*x1/2);
abs(x1)/x1*sqrt((cosh(%i*x1)-1)/2);

sinh((xpos+%i*x1)/2);
sqrt((xpos+%i*x1)^2)/(xpos+%i*x1)*sqrt((cosh(xpos+%i*x1)-1)/2);

sinh((xneg+%i*x1)/2);
sqrt((xneg+%i*x1)^2)/(xneg+%i*x1)*sqrt((cosh(xneg+%i*x1)-1)/2);

/* Cosh function with Real arguments */

cosh(x/2);
sqrt((cosh(x)+1)/2);

cosh(x/2)^2;
(cosh(x)+1)/2;

cosh(xpos/2);
sqrt((cosh(xpos)+1)/2);

cosh(xneg/2);
sqrt((cosh(xneg)+1)/2);

/* Cosh function with Complex arguments */

/* x1 is in (-%pi,%pi) */
cosh(%i*x1/2);
sqrt((cosh(%i*x1)+1)/2);

cosh((xpos+%i*x1)/2);
sqrt((cosh(xpos+%i*x1)+1)/2);

cosh((xneg+%i*x1)/2);
sqrt((cosh(xneg+%i*x1)+1)/2);

/* two tests for atan2(sqrt(1-u)*(u-1),1) - ID: 3521596 */
atan2(sqrt(1-u)*(u-1),1);
atan2(sqrt(1-u)*(u-1),1)$

atan2((1-u)^(1/2^10)*(u-1),1);
atan2((1-u)^(1/2^10)*(u-1),1)$

/* adding tests for some previously uncovered code in trigi.lisp and trigo.lisp */

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))]);
[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]$

block([exponentialize : true], [sec(x),csch(x), sech(x)]);
[2/(%e^(%i*x)+%e^-(%i*x)), 2/(%e^x-%e^-x), 2/(%e^x+%e^-x)]$

[atan(inf), atan(-minf), atan(minf), atan(-inf)];
[%pi/2, %pi/2, -%pi/2, -%pi/2]$

block([%piargs : true], [atan(-1 + sqrt(2)), atan(1+sqrt(2))]);
[%pi/8, 3*%pi/8]$

block([%iargs : true], atan(%i*x));
%i*atanh(x)$

[atan(tan(2015)), atan(tan(-2014)),  atan(cot(-2014)), atan(cot(2015))];
[2015-641*%pi,641*%pi-2014,2014-1283*%pi/2,1283*%pi/2-2015]$

asin(sqrt(3)/2);
%pi/3$

asin(-sqrt(3)/2);
-%pi/3$

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));
''(makelist(0,i,0,43))$

map(lambda([x], asec(sec(x)) - (%pi-abs(2*%pi*floor(x/(2*%pi))-x+%pi))), makelist(i/19,i,0,25));
''(makelist(0,i,0,25))$

sin(641*%pi-2015);
sin(2015)$

asin(cos(2015));
2015-1283*%pi/2$

block([triginverses : true], asin(sin(1/5)));
1/5$

acos(sin(1066));
1066-677*%pi/2$

acot(cot(42));
42-13*%pi$

acsc(csc(191/8));
191/8-8*%pi$

block([triginverses : all],
   [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))]);
[x,x, x, x, x, x, x, x, x]$

[acot(0), acot(1),acot(-1),acot(1/sqrt(3)), acot(sqrt(3))];
[%pi/2,%pi/4,-%pi/4,%pi/3,%pi/6]$

[atan(inf), atan(1),atan(-1),atan(sqrt(3)), atan(1/sqrt(3))];
[%pi/2,%pi/4,-%pi/4,%pi/3,%pi/6]$

[acos(sqrt(3)/2), acos(-sqrt(3)/2)];
[%pi/6,5*%pi/6]$

[acsc(1), acsc(-1), acsc(sqrt(2)), acsc(2/sqrt(3))];
 [%pi/2,-%pi/2,%pi/4,%pi/3]$

[asec(-1),asec(sqrt(2)), 2/sqrt(3)];
[%pi,%pi/4,2/sqrt(3)]$

/* end of tests for some previously uncovered code in trigi.lisp and trigo.lisp */

/* Restore global flags */
(halfangles:old_halfangles,%iargs:old_%iargs,done);
done;

/* SF bug # 2570: "Make acos(cos(x)) simplify to x when on correct interval" */

(kill (foo),
 acos (cos (foo)));
acos (cos (foo));

(ctxt : newcontext (),
 assume (0 <= foo and foo <= %pi),
 acos (cos (foo)));
foo;

killcontext (ctxt);
done;

/* mailing list 2015-08-25: "Strange symmetry of acoth(x), area cotangens hyperbolicus function (#552)" */

(kill (x), - acoth (- x));
acoth (x);

/* SF bug #2620: "atan2(y,x)+atan2(-y,x) doesn't always return 0 " */
block([ans,trigsign : true],
  kill(x,y),
  assume(notequal(y,0)),
  ans : atan2(y,x)+atan2(-y,x),
  forget(notequal(y,0)),
  ans);
0$


(assume (y > 0), block([trigsign : true], atan2(y,x)+atan2(-y,x)));
0;

(kill (n, p, r), forget (y > 0), assume(n<0,p>0, notequal(x,0)), atan2 (- x, r));
- atan2 (x, r);

atan2 (-n, r);
- atan2 (n, r);

atan2 (-p, r);
- atan2 (p, r);

forget (n < 0, p > 0,notequal(x,0));
[n < 0, p > 0,notequal(x,0)];

/* mailing list 2017-11-27: "trigsimp fails with pderivop" */

(load (pdiff),
 kill (f, y),
 expr:pderivop(f,1)((y))^2);
pderivop(f,1)(y)^2;

trigsimp (expr);
pderivop(f,1)(y)^2;

(reset (use_pdiff), 0); /* disables pdiff */
0;

(kill (foo, x),
 trigsimp (foo(1)(x)^3));
foo(1)(x)^3;

/* maxima-discuss 2023-02-17: "acos(1) works; acos(1.0) doesn't !"
 * bug appears to be specific to GCL versions < 2.14.
 * more extensive tests for special cases would be terrific.
 */

[acos(1), acos(1.0), acos(1b0)];
[0, 0.0, 0b0];

/* Bug 4346:  csc(0b0) produces a funny message
 * Test that the error message has the expected value
 */
block([errormsg: false], errcatch(csc(0b0)), error);
["~A: argument ~:M isn't in the domain of ~A.", csc, 0.0b0, csc];

block([errormsg: false], errcatch(cot(0b0)), error);
["~A: argument ~:M isn't in the domain of ~A.", cot, 0.0b0, cot];

/* Bug 4461: funny error message from sec(%pi/2)
 *
 * Test that the error messages are reasonable and correct now.
 */
block([errormsg: false], errcatch(sec(%pi/2)), error);
["~A: argument ~:M isn't in the domain of ~A.",sec,%pi/2,sec];

block([errormsg: false], errcatch(cot(%pi)), error);
["~A: argument ~:M isn't in the domain of ~A.",cot,%pi,cot];