xaux.red 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263
  1. module xaux;
  2. % Auxiliary functions for XIDEAL
  3. % Author: David Hartley
  4. Comment. The routines in EXCALC sometimes use a new type, here called
  5. wedgepf, internally. It has the same structure as a pf, but the powers
  6. are lists of factors in an implicit wedge product. The WEDGE tag may
  7. or may not be present. A pf, typically a 0- or 1-form, can be
  8. converted to this type using mkunarywedge. More general routines for
  9. converting pf <-> wedgepf are provided here.
  10. It is not necessary for the WEDGE kernels passed to the EXCALC product
  11. routines to be unique (and the output is not), hence two conversions
  12. lpow wedgepf -> lpow pf are given below: mkuwedge constructs a unique
  13. kernel, while mknwedge may be non-unique. The results of the product
  14. routine wedgepf defined here are unique.
  15. endcomment;
  16. symbolic procedure !*wedgepf2pf f;
  17. % f:wedgepf -> !*wedgepf2pf:pf
  18. if null f then nil
  19. else mkuwedge lpow f .* lc f .+ !*wedgepf2pf red f;
  20. symbolic procedure !*pf2wedgepf f;
  21. % f:wedgepf -> !*pf2wedgepf:pf
  22. if null f then nil
  23. else wedgefax lpow f .* lc f .+ !*pf2wedgepf red f;
  24. symbolic procedure mkuwedge u;
  25. % u:list of kernel -> mkuwedge:lpow pf
  26. % result is a unique kernel
  27. if cdr u then car fkern('wedge . u) else car u;
  28. symbolic procedure mknwedge u;
  29. % u:list of kernel -> mknwedge:lpow pf
  30. % result is a non-unique kernel
  31. if cdr u then 'wedge . u else car u;
  32. symbolic procedure wedgefax u;
  33. % u:lpow pf -> wedgefax:list of kernel
  34. if eqcar(u,'wedge) then cdr u else {u};
  35. symbolic procedure wedgepf(u,v);
  36. % u,v:pf -> wedgepf:pf
  37. !*wedgepf2pf wedgepf2(u,!*pf2wedgepf v);
  38. Comment. The list xvars!* is used to decide which 0-form kernels are
  39. counted as parameters and which as variables ("xvars") in partitioned
  40. pf's. The xvars statement allows this list to be set.
  41. endcomment;
  42. fluid '(xvars!*);
  43. rlistat '(xvars);
  44. symbolic procedure xvars u;
  45. % u:list of prefix -> xvars:nil
  46. begin
  47. xvars!* := if u = {nil} then t else xvarlist u;
  48. end;
  49. symbolic procedure xvarlist u;
  50. % u:list of prefix -> xvarlist:list of kernel
  51. % recursively evaluate and expand lists
  52. for each x in u join
  53. if eqcar(x := reval x,'list) then xvarlist cdr x
  54. else {!*a2k x};
  55. symbolic procedure xpartitsq u;
  56. % u:sq -> xpartitsq:pf
  57. % Leaves unexpanded structure if possible
  58. (if null x then nil
  59. else if domainp x then 1 .* u .+ nil
  60. else addpf(if sfp mvar x then
  61. wedgepf(xexptpf(xpartitsq(mvar x ./ 1),ldeg x),
  62. xpartitsq cancel(lc x ./ y))
  63. else if xvarp mvar x then
  64. wedgepf(xexptpf(xpartitk mvar x,ldeg x),
  65. xpartitsq cancel(lc x ./ y))
  66. else
  67. multpfsq(xpartitsq cancel(lc x ./ y),
  68. !*p2q lpow x),
  69. xpartitsq(red x ./ y)))
  70. where x = numr u, y = denr u;
  71. symbolic procedure xpartitk k;
  72. % k:kernel -> xpartitk:pf
  73. % k is an xvar. If k is not a variable (eg a wedge product)
  74. % then its arguments may need reordering if they've been through subf1.
  75. if memqcar(k,'(wedge partdf)) then
  76. (if j=k then !*k2pf k else xpartitop j) where j=reval k
  77. else !*k2pf k;
  78. symbolic procedure xpartitop u;
  79. xpartitsq simp!* u;
  80. symbolic procedure xexptpf(u,n);
  81. % u:pf,n:posint -> xexptpf:pf
  82. if n = 1 then u
  83. else wedgepf(u,xexptpf(u,n-1));
  84. symbolic procedure xvarp u;
  85. % u:kernel -> xvarp:bool
  86. % Test for exterior variables: p-forms (incl. p=0) and vectors
  87. % xvars!* controls whether 0-forms are included: if t, then all
  88. % 0-forms are included, otherwise only those in xvars!*. Forms of
  89. % degree other than 0 are always included. If xvars!* contains x,
  90. % then sin(x) is not an xvar (unless explicitly listed) since it is
  91. % algebraically independent.
  92. % Should the last line be exformp u?
  93. if xvars!* neq t then
  94. xdegree u neq 0 or u memq xvars!*
  95. else if atom u then
  96. get(u,'fdegree)
  97. else if flagp(car u,'indexvar) then
  98. assoc(length cdr u,get(car u,'ifdegree))
  99. else
  100. car u memq '(wedge d partdf hodge innerprod liedf);
  101. symbolic operator excoeffs;
  102. symbolic procedure excoeffs u;
  103. begin scalar x;
  104. u := 1 .+ xpartitop u;
  105. while (u := red u) do
  106. x := mk!*sq lc u . x;
  107. return makelist reverse x;
  108. end;
  109. symbolic operator exvars;
  110. symbolic procedure exvars u;
  111. begin scalar x;
  112. u := 1 .+ xpartitop u;
  113. while (u := red u) do
  114. x := lpow u . x;
  115. return makelist reverse x;
  116. end;
  117. % Various auxilliary functions
  118. symbolic procedure xdegree f;
  119. % f:prefix -> xdegree:int
  120. % This procedure gives the degree of a homogeneous form (deg!*form in
  121. % excalc returns nil for 0-forms). Behaves erratically with
  122. % inhomogeneous forms.
  123. (if null x then 0 else x) where x = deg!*form f;
  124. symbolic procedure xhomogeneous f;
  125. % f:pf -> xhomogeneous:int or nil
  126. % Result is degree of f if homogeneous, otherwise nil.
  127. if null f then 0
  128. else if null red f then xdegree lpow f
  129. else (if d = xhomogeneous red f then d) where d = xdegree lpow f;
  130. symbolic procedure xmaxdegree f;
  131. % f:pf -> xmaxdegree:int
  132. % Returns the maximum degree among the terms of f
  133. if null f then 0
  134. else max(xdegree lpow f,xmaxdegree red f);
  135. symbolic procedure xnormalise f;
  136. % f:pf -> xnormalise:pf
  137. % rescale f so that the leading coefficient is 1
  138. if null f then nil
  139. else if lc f = (1 ./ 1) then f
  140. else multpfsq(f,invsq lc f);
  141. symbolic procedure subs2pf f;
  142. % f:pf -> subs2pf:pf
  143. % Power check for pf. Only leading term is guaranteed correct.
  144. if f then
  145. (if numr c then lpow f .* c .+ red f else subs2pf red f)
  146. where c = subs2 resimp lc f;
  147. symbolic procedure subs2pf!* f;
  148. % f:pf -> subs2pf!*:pf
  149. % Power check for pf. All terms guaranteed correct.
  150. if f then
  151. (if numr c then lpow f .* c .+ subs2pf!* red f else subs2pf!* red f)
  152. where c = subs2 resimp lc f;
  153. % Partitioned form printing
  154. symbolic procedure !*pf2a f;
  155. % f:pf -> !*pf2a:!*sq prefix
  156. % Returns 0-form ^ 0-form to 0-form * 0-form.
  157. mk!*sq !*pf2sq repartit f;
  158. symbolic procedure !*pf2a1(f,v);
  159. % f:pf, v:bool -> !*pf2a1:prefix
  160. % !*sq prefix if v null, else true prefix.
  161. % Returns 0-form ^ 0-form to 0-form * 0-form.
  162. !*q2a1(!*pf2sq repartit f,v);
  163. symbolic procedure preppf f;
  164. % f:pf -> preppf:prefix
  165. % produce a partitioned prefix form
  166. if null(f := preppf0 f) then 0
  167. else if length f = 1 then car f
  168. else 'plus . f;
  169. symbolic procedure preppf0 f;
  170. % f:pf -> preppf0:list of prefix
  171. % produce a list of prefix terms
  172. % prepsq!* takes out over minus signs
  173. if null f then nil
  174. else preppf1(lpow f,prepsq!* lc f) . preppf0 red f;
  175. symbolic procedure preppf1(k,c);
  176. % k:lpow pf, c:prefix -> preppf1:prefix
  177. % extract an overall minus sign, and expand an overall product
  178. if k = 1 then c
  179. else if c = 1 then k
  180. else if eqcar(c,'minus) then {'minus,preppf1(k,cadr c)}
  181. else if eqcar(c,'times) then append(c,{k})
  182. else if eqcar(c,'quotient) and eqcar(cadr c,'minus) then
  183. preppf1(k,{'minus,{'quotient,cadr cadr c,caddr c}})
  184. else {'times,c,k};
  185. symbolic procedure printpf f;
  186. % f:pf -> printpf:nil
  187. % A simple printing routine for use in tracing
  188. mathprint preppf f;
  189. endmodule;
  190. end;