declare([a,b],special)$
define_variable(SINNPIFLAG,TRUE,boolean)$
define_variable(COSNPIFLAG,TRUE,boolean)$

REMFUN1(FUN,EXP):=SCANMAP(LAMBDA([Q],DELFUN1(FUN,Q)),EXP)$

DELFUN1(FUN,EXP):=IF NOT ATOM(EXP) AND INPART(EXP,0) = FUN
	 THEN FIRST(ARGS(EXP)) ELSE EXP$

REMFUNN1(FUN,EXP):=SCANMAP(LAMBDA([Q],DELFUNN1(FUN,Q)),EXP)$

DELFUNN1(FUN,EXP):=IF NOT ATOM(EXP) AND INPART(EXP,0) = FUN
	  THEN -FIRST(ARGS(EXP)) ELSE EXP$

REMFUN2(FUN,EXP,VAR):=SCANMAP(LAMBDA([Q],DELFUN2(FUN,Q,VAR)),EXP)$

DELFUN2(FUN,EXP,VAR):=
     IF	NOT ATOM(EXP) AND INPART(EXP,0) = FUN AND MEMBER(VAR,LISTOFVARS(EXP))
	 THEN FIRST(ARGS(EXP)) ELSE EXP$

REMFUNN2(FUN,EXP,VAR):=SCANMAP(LAMBDA([Q],DELFUNN2(FUN,Q,VAR)),EXP)$

DELFUNN2(FUN,EXP,VAR):=
      IF NOT ATOM(EXP) AND INPART(EXP,0) = FUN
		       AND MEMBER(VAR,LISTOFVARS(EXP))
	  THEN -FIRST(ARGS(EXP)) ELSE EXP$

REMFUN(FUN,EXP,[VAR]):=IF VAR = [] THEN REMFUN1(FUN,EXP)
        ELSE (IF LENGTH(VAR) = 1 THEN REMFUN2(FUN,EXP,FIRST(VAR))
		  ELSE ERROR("TOO MANY ARGUMENTS TO REMFUN"))$

REMFUNN(FUN,EXP,[VAR]):=IF VAR = [] THEN REMFUNN1(FUN,EXP)
	 ELSE (IF LENGTH(VAR) = 1 THEN REMFUNN2(FUN,EXP,FIRST(VAR))
		   ELSE ERROR("TOO MANY ARGUMENTS TO REMFUNN"))$

FUNP1(FUN,EXP):=BLOCK([INFLAG],INFLAG:TRUE,
      IF MAPATOM(EXP) THEN FALSE
	  ELSE (IF INPART(EXP,0) = FUN THEN TRUE
		    ELSE member(true,MAPLIST(LAMBDA([Q],FUNP1(FUN,Q)),EXP))))$
 /* commented out of DOE MACSYMA
; APPLY("OR",MAPLIST(LAMBDA([Q],FUNP1(FUN,Q)),EXP))))$ */


FUNP2(FUN,EXP,VAR):=BLOCK([INFLAG],INFLAG:TRUE,
      IF MAPATOM(EXP) THEN FALSE
	  ELSE (IF INPART(EXP,0) = FUN AND MEMBER(VAR,LISTOFVARS(EXP))
		    THEN TRUE
		    ELSE member(true,
                                MAPLIST(LAMBDA([Q],FUNP2(FUN,Q,VAR)),EXP))))$
/* commented out of DOE MACSYMA
APPLY("OR",MAPLIST(LAMBDA([Q],FUNP2(FUN,Q,VAR)),EXP))))$ */

FUNP(FUN,EXP,[VAR]):=IF VAR = [] THEN FUNP1(FUN,EXP)
      ELSE (IF LENGTH(VAR) = 1 THEN FUNP2(FUN,EXP,FIRST(VAR))
	        ELSE ERROR("TOO MANY ARGUMENTS TO FUNP"))$

EQUALP(X,Y):=BLOCK([PREDERROR],PREDERROR:FALSE,
       IF IS(EQUAL(X,Y)) = TRUE THEN TRUE)$

PARITY(F,X):=IF EVENFUNP(F,X) THEN 'EVEN
        ELSE (IF ODDFUNP(F,X) THEN 'ODD ELSE 'NEITHER)$

EVENFUNP(F,X):=IF EQUALP(F,SUBST(-X,X,F)) THEN TRUE$

ODDFUNP(F,X):=IF EQUALP(F,-SUBST(-X,X,F)) THEN TRUE$

PARINT(F,X,A,B):=
    IF NOT (EQUALP(-A,B) OR A = 'MINF AND B = 'INF OR A = 'INF AND B = 'MINF)
        THEN F
        ELSE (IF ATOM(F) OR INPART(F,0) # "+" THEN PARINT1(F,X)
		  ELSE MAP(LAMBDA([Q],PARINT1(Q,X)),F))$

PARINT1(F,X):=IF ODDFUNP(F,X) THEN 0 ELSE F$

ADEFINT(F,X,A,B):=BLOCK([ASIGN,BSIGN],
	IF EQUALP(A,B) THEN 0
	    ELSE (IF NOT (FREEOF(%I,F) AND FREEOF(%I,A) AND FREEOF(%I,B))
		      THEN LDEFINT(F,X,A,B)
		      ELSE (F:PARINT(F,X,A,B),
			    IF EQUALP(F,0) THEN 0
			        ELSE (IF NOT FUNP2('ABS,F,X)
					  THEN LDEFINT(F,X,A,B)
					  ELSE (ASIGN:ASKSIGN(A),
						BSIGN:ASKSIGN(B),
						IF 
						 (ASIGN = 'NEG
						   OR ASIGN = 'ZERO)
						   AND (BSIGN = 'NEG
						    OR BSIGN = 'ZERO)
						  THEN LDEFINT(
						  REMFUNN2('ABS,F,X),X,A,B)
						  ELSE (IF 
						  (ASIGN = 'ZERO
						    OR ASIGN = 'POS)
						    AND (BSIGN = 'ZERO
						     OR BSIGN = 'POS)
						   THEN LDEFINT(
						   REMFUN2('ABS,F,X),X,A,B)
						   ELSE (IF ASIGN = 'NEG
							     THEN RATSIMP(
							     LDEFINT(
							      REMFUNN2(
							       'ABS,F,X),X,A,
							      0)
							      +LDEFINT(
							       REMFUN2(
								'ABS,F,X),X,
							       0,B))
							     ELSE RATSIMP(
							     LDEFINT(
							      REMFUN2(
							       'ABS,F,X),X,A,
							      0)
							      +LDEFINT(
							       REMFUNN2(
								'ABS,F,X),X,
							       0,B)))))))))$

INDEFINT(F,X,HALFPLANE):=IF HALFPLANE = 'POS
	  THEN INTEGRATE(REMFUN2('ABS,F,X),X)
	  ELSE (IF HALFPLANE = 'NEG THEN INTEGRATE(REMFUNN2('ABS,F,X),X)
		    ELSE (IF HALFPLANE = 'BOTH
			      THEN APPEND(
			      LDISP(INTEGRATE(REMFUNN2('ABS,F,X),X)),
			      LDISP(INTEGRATE(REMFUN2('ABS,F,X),X)))
			      ELSE ERROR("INVALID HALFPLANE:",HALFPLANE)))$

ABSINT(F,X,[RANGE]):=IF RANGE = [] THEN INDEFINT(F,X,'POS)
        ELSE (IF LENGTH(RANGE) = 1 THEN INDEFINT(F,X,FIRST(RANGE))
		  ELSE (IF LENGTH(RANGE) = 2
			    THEN ADEFINT(F,X,RANGE[1],RANGE[2])
			    ELSE ERROR("TOO MANY ARGUMENTS TO ABSINT")))$

FOURINT(F,X):=block([z],
IF EVENFUNP(F,X) THEN APPEND(FOURINTCOS(F,X),LDISP(B[Z] = 0))
	 ELSE (IF ODDFUNP(F,X) THEN APPEND(LDISP(A[Z] = 0),FOURINTSIN(F,X))
		   ELSE FOURINTCOEFF(F,X)))$

FOURINTCOEFF(F,X):=BLOCK([AZ,BZ,Z],ASSUME(Z > 0),
	     AZ:ADEFINT(F*COS(Z*X),X,'MINF,'INF)/%PI,
	     BZ:ADEFINT(F*SIN(Z*X),X,'MINF,'INF)/%PI,
	     APPEND(LDISP(A[Z] = AZ),LDISP(B[Z] = BZ)))$

FOURINTCOS(F,X):=BLOCK([AZ,Z],ASSUME(Z > 0),
	   AZ:2*ADEFINT(F*COS(Z*X),X,0,'INF)/%PI,LDISP(A[Z] = AZ))$

FOURINTSIN(F,X):=BLOCK([BZ,Z],ASSUME(Z > 0),
	   BZ:2*ADEFINT(F*SIN(Z*X),X,0,'INF)/%PI,LDISP(B[Z] = BZ))$

FOURIER(F,X,P):=block([n],
IF EVENFUNP(F,X) THEN APPEND(FOURCOS(F,X,P),LDISP(B[N] = 0))
	 ELSE (IF ODDFUNP(F,X)
		   THEN APPEND(LDISP(A[0] = 0),LDISP(A[N] = 0),
			       FOURSIN(F,X,P)) ELSE FOURCOEFF(F,X,P)))$

FOURCOEFF(F,X,P):=BLOCK([A0,AN,BN,N],ASSUME(N > 0),
	  A0:ADEFINT(F,X,-P,P)/(2*P),AN:ADEFINT(F*COS(N*%PI*X/P),X,-P,P)/P,
	  BN:ADEFINT(F*SIN(N*%PI*X/P),X,-P,P)/P,
	  APPEND(LDISP(A[0] = A0),LDISP(A[N] = AN),LDISP(B[N] = BN)))$

FOURCOS(F,X,P):=BLOCK([A0,AN,N],ASSUME(N > 0),A0:ADEFINT(F,X,0,P)/P,
	AN:2*ADEFINT(F*COS(N*%PI*X/P),X,0,P)/P,
	APPEND(LDISP(A[0] = A0),LDISP(A[N] = AN)))$

FOURSIN(F,X,P):=BLOCK([BN,N],ASSUME(N > 0),
	BN:2*ADEFINT(F*SIN(N*%PI*X/P),X,0,P)/P,LDISP(B[N] = BN))$

FOURSIMP(EXP):=IF LISTP(EXP)
	  THEN MAP(LAMBDA([Q],FIRST(LDISP(FOURSIMP(apply('EV,[Q]))))),EXP)
	  ELSE (IF NOT FREEOF("=",EXP) THEN LHS(EXP) = FOURSIMPLE(RHS(EXP))
		    ELSE FOURSIMPLE(EXP))$

FOURSIMPLE(EXP):=BLOCK([n],
	   IF FUNP1('INTEGRATE,EXP) THEN EXP
	       ELSE (IF SINNPIFLAG THEN EXP:RATSUBST(0,SIN(N*%PI),EXP),
		     IF COSNPIFLAG THEN EXP:RATSUBST((-1)^N,COS(N*%PI),EXP),
		     FACTOR(EXP)))$

FOUREXPAND(L,X,P,NN):=BLOCK([SIMPSUM,SERIES,L1,LHSL1,n],
	   IF NOT LISTP(L) THEN ERROR("FIRST ARGUMENT NOT A LIST")
	       ELSE (IF L = [] THEN ERROR("ARGUMENT LIST IS EMPTY")
			 ELSE (L:apply('EV,[L]),SIMPSUM:TRUE,SERIES:0,
			       UNLESS L = [] DO
				      (L1:FIRST(L),L:REST(L),LHSL1:LHS(L1),
				       IF LHSL1 = A[0]
					   THEN SERIES:SERIES+RHS(L1)
					   ELSE (IF LHSL1 = A[N]
						     THEN SERIES
						     :SERIES
						      +APPLY('SUM,
							     [RHS(L1)
							       *COS(
								N*%PI*X/P),N,
							      1,NN])
						     ELSE (IF LHSL1 = B[N]
							       THEN SERIES
							       :SERIES
								+APPLY(
								 'SUM,
								 [
								  RHS(L1)
								   *SIN(
								    N*%PI*X
								     /P),N,1,
								  NN])
							       ELSE ERROR(
							       "INVALID EQUATI
ON IN ARGUMENT LIST:",
							       L1)))),
			       SERIES)))$

TOTALFOURIER(F,X,P):=FOUREXPAND(FOURSIMP(FOURIER(F,X,P)),X,P,'INF)$