fourdom.red 8.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282
  1. module fourdom; % Domain definitions for angles and fourier series
  2. % Author: John Fitch 1991.
  3. global '(domainlist!*);
  4. domainlist!*:=union('(!:fs!:),domainlist!*);
  5. put('fourier,'tag,'!:fs!:);
  6. put('!:fs!:,'dname,'fourier);
  7. flag('(!:fs!:),'field); %% Should be ring really
  8. put('!:fs!:,'i2d,'i2fourier);
  9. put('!:fs!:,'minusp,'fs!:minusp!:);
  10. put('!:fs!:,'plus,'fs!:plus!:);
  11. put('!:fs!:,'times,'fs!:times!:);
  12. put('!:fs!:, 'expt,'fs!:expt!:);
  13. put('!:fs!:,'difference,'fs!:difference!:);
  14. put('!:fs!:,'quotient,'fs!:quotient!:);
  15. put('!:fs!:, 'divide, 'fs!:divide!:);
  16. put('!:fs!:, 'gcd, 'fs!:gcd!:);
  17. put('!:fs!:,'zerop,'fs!:zerop!:);
  18. put('!:fs!:,'onep,'fs!:onep!:);
  19. put('!:fs!:,'prepfn,'fs!:prepfn!:);
  20. put('!:fs!:,'specprn,'fs!:prin!:);
  21. put('!:fs!:,'prifn,'fs!:prin!:);
  22. put('!:fs!:,'intequivfn,'fs!:intequiv!:);
  23. flag('(!:fs!:),'ratmode);
  24. % conversion functions
  25. put('!:fs!:,'!:mod!:,mkdmoderr('!:fs!:,'!:mod!:));
  26. % put('!:fs!:,'!:gi!:,mkdmoderr('!:fs!:,'!:gi!:));
  27. % put('!:fs!:,'!:rn!:,mkdmoderr('!:fs!:,'!:rn!:));
  28. put('!:rn!:,'!:fs!:,'!*d2fourier);
  29. put('!:ft!:,'!:fs!:,'cdr);
  30. put('!:gi!:,'!:fs!:,'!*d2fourier);
  31. put('!:gf!:,'!:fs!:,'!*d2fourier);
  32. put('expt, '!:fs!:, 'fs!:expt!:);
  33. % Conversion functions
  34. symbolic procedure i2fourier u;
  35. if dmode!*='!:fs!: then !*d2fourier u else u;
  36. symbolic procedure !*d2fourier u;
  37. if null u then nil else
  38. begin scalar fourier;
  39. fourier:=mkvect 3;
  40. fs!:set!-coeff(fourier,(u . 1));
  41. fs!:set!-fn(fourier,'cos);
  42. fs!:set!-angle(fourier,fs!:make!-nullangle());
  43. fs!:set!-next(fourier,nil);
  44. return get('fourier,'tag) . fourier
  45. end;
  46. symbolic procedure !*sq2fourier u;
  47. if null car u then nil else
  48. begin scalar fourier;
  49. fourier:=mkvect 3;
  50. fs!:set!-coeff(fourier,u);
  51. fs!:set!-fn(fourier,'cos);
  52. fs!:set!-angle(fourier,fs!:make!-nullangle());
  53. fs!:set!-next(fourier,nil);
  54. return get('fourier,'tag) . fourier
  55. end;
  56. symbolic procedure fs!:minusp!:(x); fs!:minusp cdr x;
  57. symbolic procedure fs!:minusp x;
  58. if null x then nil else
  59. if null fs!:next x then minusf car fs!:coeff x
  60. else fs!:minusp fs!:next x;
  61. %% Basic algebraic operations
  62. symbolic procedure fs!:times!:(x,y);
  63. % This function seems to be called with numeric values as well
  64. if null x then nil else if null y then nil
  65. else if numberp y
  66. then get('fourier,'tag) . fs!:timescoeff(y ./ 1, cdr x)
  67. else if numberp x
  68. then get('fourier,'tag) . fs!:timescoeff(x ./ 1, cdr y)
  69. else if not eqcar(x, get('fourier,'tag)) then
  70. get('fourier,'tag) . fs!:timescoeff(x,cdr y)
  71. else if not eqcar(y, get('fourier,'tag)) then
  72. get('fourier,'tag) . fs!:timescoeff(y,cdr x)
  73. else get('fourier,'tag) . fs!:times(cdr x, cdr y);
  74. symbolic procedure fs!:timescoeff(x, y);
  75. if null y then nil
  76. else begin scalar ans, coeff;
  77. coeff := multsq(x,fs!:coeff y);
  78. if coeff = '(nil . 1) then <<
  79. print "zero in times";
  80. return fs!:timescoeff(x, fs!:next y) >>;
  81. ans := mkvect 3;
  82. fs!:set!-coeff(ans,coeff);
  83. fs!:set!-fn(ans,fs!:fn y);
  84. fs!:set!-angle(ans,fs!:angle y);
  85. fs!:set!-next(ans, fs!:timescoeff(x, fs!:next y));
  86. return ans
  87. end;
  88. symbolic procedure fs!:times(x,y);
  89. if null x then nil else if null y then nil else
  90. begin scalar ans;
  91. ans := fs!:timesterm(x, y);
  92. return fs!:plus(ans, fs!:times(fs!:next x, y));
  93. end;
  94. symbolic procedure fs!:timesterm(x,y);
  95. % Treat x as a term and y as a tree
  96. if null y then nil else if null x then nil else
  97. begin scalar ans;
  98. ans := fs!:timestermterm(x,y);
  99. return fs!:plus(ans, fs!:timesterm(x, fs!:next y));
  100. end;
  101. symbolic procedure fs!:timestermterm(x,y);
  102. % x and y are terms. Generate the two answer terms.
  103. begin scalar sum, diff, ans, xv, yv, coeff;
  104. sum := mkvect 7;
  105. xv := fs!:angle x;
  106. yv := fs!:angle y;
  107. for i:=0:7 do putv!.unsafe(sum,i,
  108. getv!.unsafe(xv,i)+getv!.unsafe(yv,i));
  109. diff := mkvect 7;
  110. for i:=0:7 do putv!.unsafe(diff,i,
  111. getv!.unsafe(xv,i)-getv!.unsafe(yv,i));
  112. coeff := multsq(fs!:coeff x, fs!:coeff y);
  113. coeff := multsq(coeff, '(1 . 2));
  114. if null car coeff then return nil;
  115. if fs!:fn x = 'sin then
  116. if fs!:fn y = 'sin then
  117. % sin x*sin y => [-cos(x+y)+cos(x-y)]/2
  118. return fs!:plus(make!-term('cos, sum, negsq coeff),
  119. make!-term('cos,diff, coeff))
  120. else % fs!:fn y = 'cos
  121. % sin x * cos y => [sin(x+y)+sin(x-y)]/2
  122. return fs!:plus(make!-term('sin, sum, coeff),
  123. make!-term('sin, diff,coeff))
  124. else % fs!:fn x='cos
  125. if fs!:fn y = 'sin then
  126. % cos x*sin y => [sin(x+y)-sin(x-y)]/2
  127. return fs!:plus(make!-term('sin, sum, coeff),
  128. make!-term('sin,diff, negsq coeff))
  129. else % fs!:fn y = 'cos
  130. % cos x * cos y => [cos(x+y)+cos(x-y)]/2
  131. return fs!:plus(make!-term('cos, sum, coeff),
  132. make!-term('cos, diff,coeff))
  133. end;
  134. symbolic procedure fs!:expt!:(x,n);
  135. begin scalar ans, xx;
  136. ans := cdr !*d2fourier 1;
  137. x := cdr x;
  138. for i:=1:n do ans := fs!:times(ans,x);
  139. return get('fourier,'tag) . ans;
  140. end;
  141. symbolic procedure make!-term(fn, ang, coeff);
  142. begin scalar fourier, sign, i;
  143. sign := 0;
  144. i:=0;
  145. top: if getv!.unsafe(ang,i)<0 then sign := -1
  146. else if getv!.unsafe(ang,i)>0 then sign := 1
  147. else if i=7 then <<
  148. if fn ='sin then return nil >>
  149. else << i := i #+ 1; goto top >>;
  150. fourier:=mkvect 3;
  151. if sign = 1 or fn = 'cos then fs!:set!-coeff(fourier,coeff)
  152. else fs!:set!-coeff(fourier, multsq('(-1 . 1), coeff));
  153. fs!:set!-fn(fourier,fn);
  154. if sign = -1 then << sign := mkvect 7;
  155. for i:=0:7 do putv!.unsafe(sign,i,-getv!.unsafe(ang,i));
  156. ang := sign
  157. >>;
  158. fs!:set!-angle(fourier,ang);
  159. fs!:set!-next(fourier,nil);
  160. return fourier
  161. end;
  162. symbolic procedure fs!:quotient!:(x,y);
  163. if numberp y then fs!:times!:(x, !*sq2fourier (1 ./ y))
  164. else rerror(fourier, 98, "Unimplemented");
  165. symbolic procedure fs!:divide!:(x,y);
  166. rerror(fourier, 98, "Unimplemented");
  167. symbolic procedure fs!:gcd!:(x,y);
  168. rerror(fourier, 98, "Unimplemented");
  169. symbolic procedure fs!:difference!:(x,y);
  170. fs!:plus!:(x, fs!:negate!: y);
  171. symbolic procedure fs!:negate!: x;
  172. get('fourier,'tag) . fs!:negate cdr x;
  173. symbolic procedure fs!:negate x;
  174. if null x then nil
  175. else begin scalar ans;
  176. ans := mkvect 3;
  177. fs!:set!-coeff(ans,negsq fs!:coeff x);
  178. fs!:set!-fn(ans,fs!:fn x);
  179. fs!:set!-angle(ans,fs!:angle x);
  180. fs!:set!-next(ans, fs!:negate fs!:next x);
  181. return ans
  182. end;
  183. symbolic procedure fs!:zerop!:(u);
  184. null u or
  185. (not numberp u and
  186. null cdr u or
  187. (null fs!:next cdr u and
  188. ((numberp v and zerop v) where v=fs!:coeff cdr u)));
  189. symbolic procedure fs!:onep!:(u); fs!:onep cdr u;
  190. symbolic procedure fs!:onep u;
  191. null fs!:next u and
  192. onep fs!:coeff u and fs!:null!-angle u and fs!:fn(u) = 'cos;
  193. symbolic procedure fs!:prepfn!:(x); x;
  194. symbolic procedure simpfs u; u;
  195. put('!:fs!:,'simpfn,'simpfs);
  196. %% PRINTING FUNCTIONS
  197. %% We have all the usual problems of unit coefficients, and zero angles
  198. smacro procedure zeroterm x; fs!:coeff x = '(nil . 1);
  199. symbolic procedure fs!:prin!:(x);
  200. << prin2!* "["; fs!:prin cdr x; prin2!* "]" >>;
  201. symbolic procedure fs!:prin x;
  202. if null x then prin2!* " 0 " else <<
  203. while x do <<
  204. fs!:prin1 x;
  205. x := fs!:next x;
  206. if x then prin2!* " + "
  207. >>
  208. >>;
  209. symbolic procedure fs!:prin1 x;
  210. begin scalar first, u, v;
  211. first := t;
  212. if not(fs!:coeff x = '(1 . 1)) then <<
  213. prin2!* "("; sqprint fs!:coeff x;
  214. prin2!* ")" >>;
  215. if not(fs!:null!-angle x) then <<
  216. prin2!* fs!:fn x;
  217. prin2!* "[";
  218. u := fs!:angle x;
  219. for i:=0:7 do
  220. if not((v := getv!.unsafe(u,i)) = 0) then <<
  221. if v<0 then << first := t; prin2!* "-"; v := -v >>;
  222. if not first then prin2!* "+";
  223. if not(v=1) then prin2!* v;
  224. first := nil;
  225. prin2!* getv!.unsafe(fourier!-name!*, i)
  226. >>;
  227. prin2!* "]"
  228. >>
  229. else if fs!:coeff x = '(1 . 1) then prin2!* "1"
  230. end;
  231. symbolic procedure fs!:intequiv!:(u);
  232. null fs!:next x and
  233. fs!:null!-angle x and
  234. fs!:fn(x) = 'cos and
  235. fixp car fs!:coeff x and
  236. cdr fs!:coeff x = 1
  237. where x = cdr u;
  238. endmodule;
  239. end;