reduce-historical/r35/lib/camal.red

%%off comp;
lisp;
%%module foursupport;

%% This section is to define macros and simple functions to handle the
%% data structures for harmonic forms.
%% The structure is a vector:
%% Coeff | FN | Angle | Next
%
%% This version only allows 8 angles.  Consider extending this later.

switch fourier;


%% A vector and counter to record link between angle names and index

global '(next!-angle!* fourier!-name!*);
next!-angle!* := 0;
if vectorp fourier!-name!* then <<
        for i :=0:7 do remprop(getv(fourier!-name!*, i), 'fourier!-angle)
>>;
fourier!-name!* := mkvect 7;

%% For non Cambridge LISP add 
smacro procedure putv!.unsafe(x,y,z); putv(x,y,z);
smacro procedure getv!.unsafe(x,y); getv(x,y);

%% Data abtraction says that we should define macros for access to 
%% the parts of the Fourier structure

smacro procedure fs!:set!-next(f,p); putv!.unsafe(f, 3, p);

smacro procedure fs!:next(f); getv!.unsafe(f,3);

smacro procedure fs!:set!-coeff(f,p); putv!.unsafe(f, 0, p);

smacro procedure fs!:coeff(f); getv!.unsafe(f, 0);

smacro procedure fs!:set!-fn(f,p); putv!.unsafe(f, 1, p);

smacro procedure fs!:fn(f); getv!.unsafe(f, 1);

smacro procedure fs!:set!-angle(f,p); putv!.unsafe(f, 2, p);

smacro procedure fs!:angle(f); getv!.unsafe(f, 2);

%% Some support functions for angle expressions

symbolic procedure fs!:make!-nullangle();
begin scalar ans;
    ans := mkvect 7;
    for i:=0:7 do putv!.unsafe(ans,i,0);
    return ans;
end;

symbolic procedure fs!:null!-angle!: u;
  fs!:null!-angle cdr u;

symbolic procedure fs!:null!-angle u;
begin scalar ans, i, x;
    x := fs!:angle u;
    ans := t;
    i := 0;
top:
    if not(getv!.unsafe(x,i)=0) then return nil;
    i := i+1;
    if (i<8) then go to top;
    return ans;
end;


%%module fourdom; % Domain definitions for angles and fourier series

% Authors: John Fitch 1991

global '(domainlist!*);

domainlist!*:=union('(!:fs!:),domainlist!*);

put('fourier,'tag,'!:fs!:);
put('!:fs!:,'dname,'fourier);
flag('(!:fs!:),'field); %% Should be ring really
put('!:fs!:,'i2d,'i2fourier);
put('!:fs!:,'minusp,'fs!:minusp!:);
put('!:fs!:,'plus,'fs!:plus!:);
put('!:fs!:,'times,'fs!:times!:);
put('!:fs!:, 'expt,'fs!:expt!:);
put('!:fs!:,'difference,'fs!:difference!:);
put('!:fs!:,'quotient,'fs!:quotient!:);
put('!:fs!:, 'divide, 'fs!:divide!:);
put('!:fs!:, 'gcd, 'fs!:gcd!:);
put('!:fs!:,'zerop,'fs!:zerop!:);
put('!:fs!:,'onep,'fs!:onep!:);
put('!:fs!:,'prepfn,'fs!:prepfn!:);
put('!:fs!:,'specprn,'fs!:prin!:);
put('!:fs!:,'prifn,'fs!:prin!:);
put('!:fs!:,'intequivfn,'fs!:intequiv!:);
flag('(!:fs!:),'ratmode);
% conversion functions

put('!:fs!:,'!:mod!:,mkdmoderr('!:fs!:,'!:mod!:));
% put('!:fs!:,'!:gi!:,mkdmoderr('!:fs!:,'!:gi!:));
% put('!:fs!:,'!:rn!:,mkdmoderr('!:fs!:,'!:rn!:));
put('!:rn!:,'!:fs!:,'!*d2fourier);
put('!:ft!:,'!:fs!:,'cdr);
put('!:gi!:,'!:fs!:,'!*d2fourier);
put('!:gf!:,'!:fs!:,'!*d2fourier);

put('expt, '!:fs!:, 'fs!:expt!:);

% Conversion functions

symbolic procedure i2fourier u; 
  if dmode!*='!:fs!: then !*d2fourier u else u;

symbolic procedure !*d2fourier u;
if null u then nil else
begin scalar fourier;
      fourier:=mkvect 3;
      fs!:set!-coeff(fourier,(u . 1)); 
      fs!:set!-fn(fourier,'cos);
      fs!:set!-angle(fourier,fs!:make!-nullangle()); 
      fs!:set!-next(fourier,nil); 
     return get('fourier,'tag) . fourier
end;

symbolic procedure !*sq2fourier u;
if null car u then nil else
begin scalar fourier;
      fourier:=mkvect 3;
      fs!:set!-coeff(fourier,u); 
      fs!:set!-fn(fourier,'cos);
      fs!:set!-angle(fourier,fs!:make!-nullangle()); 
      fs!:set!-next(fourier,nil); 
     return get('fourier,'tag) . fourier
end;

symbolic procedure fs!:minusp!:(x); fs!:minusp cdr x;

symbolic procedure fs!:minusp x; 
if null x then nil else 
   if null fs!:next x then minusf car fs!:coeff x
   else fs!:minusp fs!:next x;

%% Basic algebraic operations

symbolic procedure fs!:times!:(x,y);
% This function seems to be called with numeric values as well
   if null x then nil else if null y then nil
   else if numberp y then get('fourier,'tag) . fs!:timescoeff(y ./ 1, cdr x)
   else if numberp x then get('fourier,'tag) . fs!:timescoeff(x ./ 1, cdr y)
   else if not eqcar(x, get('fourier,'tag)) then
        get('fourier,'tag) . fs!:timescoeff(x,cdr y)
   else if not eqcar(y, get('fourier,'tag)) then
        get('fourier,'tag) . fs!:timescoeff(y,cdr x)
   else get('fourier,'tag) . fs!:times(cdr x, cdr y);

symbolic procedure fs!:timescoeff(x, y);
if null y then nil
   else begin scalar ans, coeff;
      coeff := multsq(x,fs!:coeff y); 
      if coeff = '(nil . 1) then <<
        print "zero in times";
        return fs!:timescoeff(x, fs!:next y) >>;
      ans := mkvect 3;
      fs!:set!-coeff(ans,coeff);
      fs!:set!-fn(ans,fs!:fn y);
      fs!:set!-angle(ans,fs!:angle y); 
      fs!:set!-next(ans, fs!:timescoeff(x, fs!:next y));
      return ans
   end;

symbolic procedure fs!:times(x,y);
if null x then nil else if null y then nil else
begin scalar ans;
        ans := fs!:timesterm(x, y);
        return fs!:plus(ans, fs!:times(fs!:next  x, y));
end;

symbolic procedure fs!:timesterm(x,y);
% Treat x as a term and y as a tree
if null y then nil else if null x then nil else
begin scalar ans;
        ans := fs!:timestermterm(x,y);
        return fs!:plus(ans, fs!:timesterm(x, fs!:next y));
end;

symbolic procedure fs!:timestermterm(x,y);
% x and y are terms.  Generate the two answer terms.
begin scalar sum, diff, ans, xv, yv, coeff;
        sum := mkvect 7;
        xv := fs!:angle x;
        yv := fs!:angle y;
        for i:=0:7 do putv!.unsafe(sum,i,
                                   getv!.unsafe(xv,i)+getv!.unsafe(yv,i));
        diff := mkvect 7;
        for i:=0:7 do putv!.unsafe(diff,i, 
                                   getv!.unsafe(xv,i)-getv!.unsafe(yv,i));
        coeff := multsq(fs!:coeff x, fs!:coeff y);
        coeff := multsq(coeff, '(1 . 2));
        if null car coeff then return nil;
        if fs!:fn x = 'sin then
            if fs!:fn y = 'sin then
                % sin x*sin y => [-cos(x+y)+cos(x-y)]/2
                return fs!:plus(make!-term('cos, sum, negsq coeff),
                                make!-term('cos,diff, coeff))
            else % fs!:fn y = 'cos
                % sin x * cos y => [sin(x+y)+sin(x-y)]/2
                return fs!:plus(make!-term('sin, sum, coeff),
                                make!-term('sin, diff,coeff))
        else % fs!:fn x='cos
            if fs!:fn y = 'sin then
                % cos x*sin y => [sin(x+y)-sin(x-y)]/2
                return fs!:plus(make!-term('sin, sum, coeff),
                                make!-term('sin,diff, negsq coeff))
            else % fs!:fn y = 'cos
                % cos x * cos y => [cos(x+y)+cos(x-y)]/2
                return fs!:plus(make!-term('cos, sum, coeff),
                                make!-term('cos, diff,coeff))
            
end;

symbolic procedure fs!:expt!:(x,n);
begin scalar ans, xx;
    ans := cdr !*d2fourier 1;
    x := cdr x;
    for i:=1:n do ans := fs!:times(ans,x);
    return get('fourier,'tag) . ans;
end;

symbolic procedure make!-term(fn, ang, coeff);
begin scalar fourier, sign, i;
      sign := 0;
      i:=0;
top:  if getv!.unsafe(ang,i)<0 then sign := -1
      else if getv!.unsafe(ang,i)>0 then sign := 1
      else if i=7 then <<
        if fn ='sin then return nil >>
      else << i := i #+ 1; goto top >>;
      fourier:=mkvect 3;
      if sign = 1 or fn = 'cos then fs!:set!-coeff(fourier,coeff)
      else fs!:set!-coeff(fourier, multsq('(-1 . 1), coeff));
      fs!:set!-fn(fourier,fn);
      if sign = -1 then << sign := mkvect 7;
        for i:=0:7 do putv!.unsafe(sign,i,-getv!.unsafe(ang,i));
        ang := sign
      >>;
      fs!:set!-angle(fourier,ang); 
      fs!:set!-next(fourier,nil); 
     return fourier
end;

symbolic procedure fs!:quotient!:(x,y);
if numberp y then fs!:times!:(x, !*sq2fourier (1 ./ y))
else rerror(fourier, 98, "Unimplemented");

symbolic procedure fs!:divide!:(x,y);
rerror(fourier, 98, "Unimplemented");

symbolic procedure fs!:gcd!:(x,y);
rerror(fourier, 98, "Unimplemented");

symbolic procedure fs!:difference!:(x,y);
   fs!:plus!:(x, fs!:negate!: y);

symbolic procedure fs!:negate!: x;
  get('fourier,'tag) . fs!:negate cdr x;

symbolic procedure fs!:negate x;
   if null x then nil
   else begin scalar ans;
      ans := mkvect 3;
      fs!:set!-coeff(ans,negsq fs!:coeff x); 
      fs!:set!-fn(ans,fs!:fn x);
      fs!:set!-angle(ans,fs!:angle x); 
      fs!:set!-next(ans, fs!:negate fs!:next x); 
      return ans
   end;

symbolic procedure fs!:zerop!:(u);
  null u or
  (not numberp u and
   null cdr u or
   (null fs!:next cdr u and 
   ((numberp v and zerop v) where v=fs!:coeff cdr u)));

symbolic procedure fs!:onep!:(u); fs!:onep cdr u;

symbolic procedure fs!:onep u;
  null fs!:next u and 
  onep fs!:coeff u and fs!:null!-angle u and fs!:fn(u) = 'cos;

symbolic procedure fs!:prepfn!:(x); x;

symbolic procedure simpfs u; u;

put('!:fs!:,'simpfn,'simpfs);

%% PRINTING FUNCTIONS

%% We have all the usual problems of unit coefficients, and zero angles

smacro procedure zeroterm x; fs!:coeff x = '(nil  . 1);

symbolic procedure fs!:prin!:(x);
  << prin2!* "["; fs!:prin cdr x; prin2!* "]" >>;

symbolic procedure fs!:prin x;
   if null x then prin2!* " 0 " else <<
   while x do <<
     fs!:prin1 x;
     x := fs!:next x;
     if x then prin2!* " + "
   >>
>>;

symbolic procedure fs!:prin1 x;
begin scalar first, u, v;
   first := t;
   if not(fs!:coeff x = '(1 . 1)) then <<
      prin2!* "("; sqprint fs!:coeff x;
      prin2!* ")" >>;
   if not(fs!:null!-angle x) then <<
     prin2!* fs!:fn x;
     prin2!* "[";
     u := fs!:angle x;
     for i:=0:7 do
         if not((v := getv!.unsafe(u,i)) = 0) then <<
            if v<0 then << first := t; prin2!* "-"; v := -v >>;
            if not first then prin2!* "+";
            if not(v=1) then prin2!* v;
            first := nil;
            prin2!* getv!.unsafe(fourier!-name!*, i)
     >>;
     prin2!* "]"
  >>
  else if fs!:coeff x = '(1 . 1) then prin2!* "1"
end;

symbolic procedure fs!:intequiv!:(u);
   null fs!:next x and 
   fs!:null!-angle x and
   fs!:fn(x) = 'cos and 
   fixp car fs!:coeff x and
   cdr fs!:coeff x = 1
        where x = cdr u;

%%module fourplus;

%% ARITHMETIC

%% Addition of Fourier expressionsis really a merge operation

symbolic procedure fs!:plus!:(x,y);
 %% Top level addition of two fourier series
    if fs!:zerop!: y then x
    else if fs!:zerop!: x then y
    else get('fourier,'tag) . fs!:plus(copy!-tree cdr x, copy!-tree cdr y);

% I cannot rely on the CAMAL selective copy, so I take the coward's way out
symbolic procedure copy!-tree x;
   if null x then nil
   else begin scalar ans;
      ans := mkvect 3;
      fs!:set!-coeff(ans,fs!:coeff x); 
      fs!:set!-fn(ans,fs!:fn x);
      fs!:set!-angle(ans,fs!:angle x); 
      fs!:set!-next(ans, copy!-tree fs!:next x); 
      return ans
   end;

symbolic procedure fs!:plus(x, y);
  %% The real addition.  x is a new tree to which y must be merged.
  if null y then x
  else if null x then y
  else if fs!:fn x = fs!:fn y and angles!-equal(fs!:angle x, fs!:angle y) then
        begin scalar coef;
            coef := addsq(fs!:coeff x, fs!:coeff y);
        % Really I should deal with the zero case here
            if null car coef then return fs!:plus(fs!:next x, fs!:next y);
            fs!:set!-coeff(x, coef);
            fs!:set!-next(x, fs!:plus(fs!:next x, fs!:next y));
            return x 
        end
    else if fs!:angle!-order(x, y) then <<
          fs!:set!-next(x, fs!:plus(fs!:next x, y));
          x >>
    else <<
          fs!:set!-next(y, fs!:plus(fs!:next y,x));
          y >>;

symbolic procedure angles!-equal(x, y);
% Are all angles the same?
begin scalar i;
    i := 0;
top:
    if not(getv!.unsafe(x,i)=getv!.unsafe(y,i)) then return nil;
    i := i+1;
    if (i<8) then go to top;
    return t;
end;

symbolic procedure fs!:angle!-order(x, y);
% Ordering function for angle expressions, also taking account of angle.
begin scalar ans, i, xx, yy;
    i := 0;
    xx := fs!:angle x;
    yy := fs!:angle y;
top:
    ans := (getv!.unsafe(xx,i)-getv!.unsafe(yy,i));
    if not(ans = 0) then return ans>0;
    i := i+1;
    if (i<8) then go to top;
    return
      if fs!:fn x = fs!:fn y then nil else if fs!:fn x = 'sin then nil else t;
end;

%%module makefour;

%% User interface; all rather iffy at present

symbolic procedure harmonicp u; get(u, 'fourier!-angle);

symbolic procedure harmonic u; 
<<
    for each x in u do if not(get(x, 'fourier!-angle)) then <<
        if (next!-angle!* > 7) then rerror(fourier, 3, "Too many angles");
        put(x, 'fourier!-angle, next!-angle!*);
        putv!.unsafe(fourier!-name!*, next!-angle!*, x);
        next!-angle!* := next!-angle!* #+ 1;
      >>
>>;

put('harmonic, 'stat, 'rlis);

symbolic procedure simpfourier u;
%% Handle the form fourier(...)  with treating sin and cos as special
begin 
  if not(length u = 1) then 
        rerror(fourier,1,"Argument should be single expression");
  return simpfourier1 prepsq simp!* car u;;
end;

symbolic procedure simpfourier1 u;
begin scalar ff;
  if atom u then <<
        if harmonicp u then rerror(fourier,2,"Secular angle not allowed");
        return (!*sq2fourier simp u) . 1;
  >>
  else if eqcar(u, '!:fs!:) then return u
  else if (ff := get(car u, 'simpfour)) then return apply1(ff, cdr u)
  else <<
    rerror(fourier,4,"Unknown function" . car u);
    return (!*sq2fourier u) . 1;
  >>
end;

put('fourier, 'simpfn, 'simpfourier);

symbolic procedure simpfouriersin u;
% Creation of a simple angle expression and function
begin scalar ans, vv;
  u := car u;
  if atom u then
        if harmonicp u then <<
              ans:=mkvect 3;
              fs!:set!-coeff(ans,(1 . 1)); 
              fs!:set!-fn(ans,'sin);
              vv := mkvect 7;
              for i:=0:7 do putv!.unsafe(vv,i,0);
              putv!.unsafe(vv, get(u, 'fourier!-angle), 1);
              fs!:set!-angle(ans,vv);
              fs!:set!-next(ans,nil); 
              return (get('fourier,'tag) . ans) . 1 >>
        else return !*sq2fourier(simp list('sin, u)) . 1;
  if angle!-expression!-p u then <<
      ans:=mkvect 3;
      fs!:set!-coeff(ans,(1 . 1)); 
      fs!:set!-fn(ans,'sin);
      vv := mkvect 7;
      for i:=0:7 do putv!.unsafe(vv,i,0);
      compile!-angle!-expression(u,vv);
      fs!:set!-angle(ans,vv);
      fs!:set!-next(ans,nil); 
      return (get('fourier,'tag) . ans) . 1 >>;
  rerror(fourier,99,"Not finished yet");
end;

put('sin, 'simpfour, 'simpfouriersin);

symbolic procedure simpfouriercos u;
% Creation of a simple angle expression and function
begin scalar ans, vv;
  u := car u;
  if atom u then
        if harmonicp u then <<
              ans:=mkvect 3;
              fs!:set!-coeff(ans,(1 . 1)); 
              fs!:set!-fn(ans,'cos);
              vv := mkvect 7;
              for i:=0:7 do putv!.unsafe(vv,i,0);
              putv!.unsafe(vv, get(u, 'fourier!-angle), 1);
              fs!:set!-angle(ans,vv);
              fs!:set!-next(ans,nil); 
              return (get('fourier,'tag) . ans) . 1 >>
        else return !*sq2fourier(simp list('cos, u)) . 1;
  if angle!-expression!-p u then <<
      ans:=mkvect 3;
      fs!:set!-coeff(ans,(1 . 1)); 
      fs!:set!-fn(ans,'cos);
      vv := mkvect 7;
      for i:=0:7 do putv!.unsafe(vv,i,0);
      compile!-angle!-expression(u,vv);
      fs!:set!-angle(ans,vv);
      fs!:set!-next(ans,nil); 
      return (get('fourier,'tag) . ans) . 1 >>;
  rerror(fourier,99,"Not finished yet");
end;

put('cos, 'simpfour, 'simpfouriercos);

%% Is the prefix expression u a sum of angles??

symbolic procedure angle!-expression!-p u;
if atom u and harmonicp u then t
else if eqcar(u,'plus) or eqcar(u,'difference) then 
    angle!-expression!-p cadr u and angle!-expression!-p caddr u
else if eqcar(u,'minus) then angle!-expression!-p cadr u
else if eqcar(u,'times) then
  if numberp cadr u then angle!-expression!-p caddr u
  else angle!-expression!-p cadr u and numberp caddr u
else nil;

%% We know that u is a sum of angles, so create the vector of coefficients;

symbolic procedure compile!-angle!-expression(u,v);
if atom u and harmonicp u then 
    putv!.unsafe(v, get(u, 'fourier!-angle),
                 1+getv!.unsafe(v, get(u, 'fourier!-angle)))
else if eqcar(u,'plus) then <<
        u := cdr u;
        while u do <<
            compile!-angle!-expression(car u,v);
            u := cdr u
        >>;
        v  >>
else if eqcar(u,'difference) then begin scalar vv;
    compile!-angle!-expression(cadr u,v);
    vv := mkvect 7;
    for i:=0:7 do putv!.unsafe(vv,i,0);
    compile!-angle!-expression(caddr u,vv);
    for i:=0:7 do putv!.unsafe(v,i,getv!.unsafe(v,i)-getv!.unsafe(vv,i));
    return v
  end
else if eqcar(u,'minus) then 
  begin scalar vv;
        vv := mkvect 7;
        for i:=0:7 do putv!.unsafe(vv,i,0);
        compile!-angle!-expression(cadr u,vv);
        for i:=0:7 do putv!.unsafe(v,i,getv!.unsafe(v,i)-getv!.unsafe(vv,i));
        return v;
  end
else if eqcar(u,'times) then
  if numberp cadr u then begin scalar vv;
      vv := mkvect 7;
      for i:=0:7 do putv!.unsafe(vv,i,0);
      compile!-angle!-expression(caddr u,vv);
      for i:=0:7 do putv!.unsafe(v, i, 
                        cadr u * getv!.unsafe(vv, i) + getv!.unsafe(v,i))
  end
  else begin scalar vv;
      vv := mkvect 7;
      for i:=0:7 do putv!.unsafe(vv,i,0);
      compile!-angle!-expression(cadr u,vv);
      for i:=0:7 do putv!.unsafe(v, i, 
                        caddr u * getv!.unsafe(vv, i) + getv!.unsafe(v,i))
  end
else nil;

symbolic procedure simpfouriertimes(u);
begin scalar z;
     z := car simpfourier1 car u;
     u := cdr u;
  a: if null u then return z ./ 1;
     z := fs!:times!:(car simpfourier1 car u,z);
     u := cdr u;
     go to a
   end;

put('times, 'simpfour, 'simpfouriertimes);

symbolic procedure simpfourierexpt(u);
  fs!:expt!:(car simpfourier1 car u, cadr u) . 1;

put('expt, 'simpfour, 'simpfourierexpt);

symbolic procedure simpfourierplus(u);
begin scalar z;
     z := car simpfourier1 car u;
     u := cdr u;
  a: if null u then return z ./ 1;
     z := fs!:plus!:(car simpfourier1 car u,z);
     u := cdr u;
     go to a
   end;

put('plus, 'simpfour, 'simpfourierplus);

symbolic procedure simpfourierdifference(u);
  fs!:difference!:(car simpfourier1 car u, car simpfourier1 cadr u) ./ 1;

put('difference, 'simpfour, 'simpfourierdifference);

symbolic procedure simpfourierminus(u);
  fs!:negate!:(car simpfourier1 car u) . 1;

put('minus, 'simpfour, 'simpfourierminus);

symbolic procedure simpfourierquot(u);
begin scalar v;
  v := simp!* cadr u;
  v := cdr v . car v;
  return fs!:times!:(car simpfourier1 car u, !*sq2fourier v) ./ 1
end;

put('quotient, 'simpfour, 'simpfourierquot);

symbolic procedure simphsin u;
begin 
  if not(length u = 1) then 
        rerror(fourier,5,"Argument should be single expression");
  return simpfouriersin list(u := prepsq simp!* car u)
end;

put('hsin, 'simpfn, 'simphsin);

symbolic procedure simphcos u;
begin 
  if not(length u = 1) then 
        rerror(fourier,6,"Argument should be single expression");
  return simpfouriercos list(u := prepsq simp!* car u)
end;

put('hcos, 'simpfn, 'simphcos);

%%endmodule;

%%module hsub

%% Harmonic substitution: the CAMAL HSUB operation, as well as other
%% substitutions.

fluid '(!*trharm);
switch trham;

symbolic procedure hsub(x,u,v,A,n);
%% Substitute v+A for u in x to order n
begin scalar ans, c, tmp, fs!:zero!-generated;
%%    fs!:zero!-generated := 0;
    ans := fs!:subang(x, u, v);
%    c := ensure!-fourier A;
    c := car A;
    if c then c := cdr c;
    A := c;
    if !*trham then << print "A"; if null A then print 0 else fs!:prin A >>;
    for i:=1:n do <<
        if !*trham then << print "i="; print i >>;
        x := hdiff(x, u);
        if !*trham then << prin2!* "df(x,u,i)="; fs!:prin x; terpri!* t;
                prin2!* "A^i ="; fs!:prin c; terpri!* t >>;
        c := fs!:times(cdr !*sq2fourier (1 ./ i), c);
        if !*trham then << prin2!* "A^i/fact(i) ="; fs!:prin c; terpri!* t>>;
        tmp := fs!:times(fs!:subang(x, u, v), c);
        if !*trham then <<
            prin2!* "f'(0)*A^i/fact i = "; fs!:prin tmp; terpri!* t>>;
        ans := fs!:plus(ans, tmp);
        if !*trham then << prin2!* "partial sum ="; fs!:prin ans; terpri!* t>>;
        if not(i=n) then c := fs!:times(c,A);
    >>;
    return ans
end;

symbolic procedure fs!:subang(x, u, v);
if null x then nil 
else begin scalar vv, n;
    vv := mkvect 7;
    n := getv!.unsafe(fs!:angle x, u);
    for i:=0:7 do if i = u then putv!.unsafe(vv, i, n*getv!.unsafe(v,i))
        else putv!.unsafe(vv, i, 
                getv!.unsafe(fs!:angle x,i) + n*getv!.unsafe(v,i));
    return fs!:plus(fs!:subang(fs!:next x, u, v),
                    make!-term(fs!:fn x, vv, fs!:coeff x));
end;

symbolic procedure fs!:sub(x,u);
if null x then nil else
begin scalar ans;
    ans := aeval prepsq fs!:coeff x;
    if not fixp ans then ans := subsq(cadr ans, u) 
    else ans := fs!:coeff x;
    if eqcar(numr ans, '!:fs!:) then ans := cdar ans
        else ans := cdr !*sq2fourier ans;
    ans := fs!:times(make!-term(fs!:fn x, fs!:angle x, 1 ./ 1), ans);
    return fs!:plus(fs!:sub(fs!:next x, u), ans);
end;

symbolic procedure simphsub uu;
begin scalar x, u, v, vv, A, n, dmode!*;
    dmode!* := '!:fs!:;
    if (length uu = 5) then <<
        x := car uu; uu := cdr uu;
        u := car uu; uu := cdr uu;
        v := car uu; uu := cdr uu;
        A := car uu; uu := cdr uu;
        n := car uu
    >>
    else if (length uu = 3) then <<
        x := car uu; uu := cdr uu;
        u := car uu; uu := cdr uu;
        v := car uu; uu := cdr uu;
        if not harmonicp u then <<
            A := ( ((get('fourier, 'tag) . 
                         fs!:sub(cdar simp x, list(u . v))) ./ 1)
            )  where wtl!*=delasc(u,wtl!*);
            return A;
        >>;
        A := 0;
        n := 0
    >>;
    if not harmonicp u then
        rerror(fourier, 7, "Not an angle in HSUB");
    x := cdar simp x;
    if not angle!-expression!-p v then
        rerror(fourier, 8, "Not an angle expression in HSUB");
    vv := mkvect 7;
    for i:=0:7 do putv!.unsafe(vv,i,0);
    compile!-angle!-expression(v, vv);
    A := simp!* A;
    n := simp!* n;
    if null car n then n := 0 ./ 1
    else if not(fixp car n and cdr n = 1)  then
        rerror(fourier, 9, "Non integer expansion in HSUB");
    n := car n;
    return (get('fourier, 'tag) . hsub(x,get(u,'fourier!-angle),vv,A,n)) ./ 1;
end;

put('hsub, 'simpfn, 'simphsub);

%%endmodule

%%module hdiff

%% Harmonic differentiation and Integration

symbolic procedure hdiff(x, u);
if null x then nil
    else fs!:plus(hdiff(fs!:next x,u), hdiffterm(x,u));

symbolic procedure hdiffterm(x, u);
begin scalar n;
   n := getv!.unsafe(fs!:angle x, u);
   if n = 0 then return nil;
   n := multsq( n . 1, fs!:coeff x);
   if fs!:fn x = 'cos then return make!-term('sin, fs!:angle x, negsq n)
   else                    return make!-term('cos, fs!:angle x, n)
end;

symbolic procedure hdiff1(x, u);
if null x then nil
else begin scalar ans, aaa;
        ans := diffsq(fs!:coeff x, u);
        if ans then <<
            aaa := mkvect 3;
            fs!:set!-coeff(aaa, ans); 
            fs!:set!-fn(aaa, fs!:fn x);
            fs!:set!-angle(aaa,fs!:angle x);
            fs!:set!-next(aaa, hdiff1(fs!:next x, u));
            return aaa >>
        else return hdiff1(fs!:next x, u)
end;

symbolic procedure simphdiff uu;
begin scalar x, u;
    if not (length uu = 2) then
        rerror(fourier, 10, "Improper number of arguments to HDIFF");
    x := car uu; uu := cdr uu;
    u := car uu;
    x := simp x;
    if not eqcar(car x, '!:fs!:) then x := !*sq2fourier x ./ 1;
    if not harmonicp u then 
        return (get('fourier, 'tag) . hdiff1(cdar x, u)) ./ 1;
    x := hdiff(cdar x,get(u,'fourier!-angle));
    if null x then return nil ./ 1;
    return (get('fourier, 'tag) . x) ./ 1
end;

put('hdiff, 'simpfn, 'simphdiff);

symbolic procedure hint(x, u);
if null x then nil
%% Bind fs!:zero!-generated ??
    else fs!:plus(hint(fs!:next x,u), hintterm(x,u));

symbolic procedure hintterm(x, u);
begin scalar n;
   n := getv!.unsafe(fs!:angle x, u);
   if n = 0 then return make!-term(fs!:fn x, fs!:angle x, fs!:coeff x);
   n := multsq( 1 ./ n, fs!:coeff x);
   if fs!:fn x = 'cos then return make!-term('sin, fs!:angle x, n)
   else                    return make!-term('cos, fs!:angle x, negsq n)
end;

symbolic procedure hint1(x , u);
if null x then nil
else begin scalar aaa;
        aaa := mkvect 3;
        fs!:set!-coeff(aaa, simpint list(prepsq fs!:coeff x, u)); 
        fs!:set!-fn(aaa, fs!:fn x);
        fs!:set!-angle(aaa,fs!:angle x);
        fs!:set!-next(aaa, hint1(fs!:next x, u));
        return aaa
end;

symbolic procedure simphint uu;
begin scalar x, u;
    if not (length uu = 2) then
        rerror(fourier, 11, "Improper number of arguments to HINT");
    x := car uu; uu := cdr uu;
    u := car uu;
    x := simp x;
    if not eqcar(car x, '!:fs!:) then x := !*sq2fourier x ./ 1;
    if not harmonicp u then 
        return (get('fourier, 'tag) . hint1(cdar x, u)) ./ 1;
    x := hint(cdar x,get(u,'fourier!-angle));
    if null x then return nil ./ 1;
    return (get('fourier, 'tag) . x) ./ 1
end;

put('hint, 'simpfn, 'simphint);

initdmode 'fourier;

algebraic;

end;


==John