1module scripts;
2
3% Redistribution and use in source and binary forms, with or without
4% modification, are permitted provided that the following conditions are met:
5%
6%    * Redistributions of source code must retain the relevant copyright
7%      notice, this list of conditions and the following disclaimer.
8%    * Redistributions in binary form must reproduce the above copyright
9%      notice, this list of conditions and the following disclaimer in the
10%      documentation and/or other materials provided with the distribution.
11%
12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
16% CONTRIBUTORS
17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
23% POSSIBILITY OF SUCH DAMAGE.
24%
25
26
27COMMENT
28
29               ######################
30               ##                  ##
31               ##     ADVANCED     ##
32               ##   APPLICATIONS   ##
33               ##                  ##
34               ######################
35
36This module contains several additional advanced applications of
37standard basis computations, inspired partly by the scripts
38distributed with the commutative algebra package MACAULAY
39(Bayer/Stillman/Eisenbud).
40
41The following topics are currently covered :
42        - [BGK]'s heuristic variable optimization
43        - certain stuff on maps (preimage, ratpreimage)
44        - ideals of points (in affine and proj. spaces)
45        - ideals of (affine and proj.) monomial curves
46        - General Rees rings, associated graded rings, and related
47                topics (analytic spread, symmetric algebra)
48        - several short scripts (minimal generators, symbolic powers
49                of primes, singular locus)
50
51
52END COMMENT;
53
54%---------- [BGK]'s heuristic variable optimization ----------
55
56symbolic operator varopt;
57symbolic procedure varopt m;
58  if !*mode='algebraic then makelist varopt!* dpmat_from_a m
59  else varopt!* m;
60
61symbolic procedure varopt!* m;
62% Find a heuristically optimal variable order.
63  begin scalar c; c:=mo_zero();
64  for each x in dpmat_list m do
65    for each y in bas_dpoly x do c:=mo_lcm(c,car y);
66  return
67    for each x in sort(mo_2list c,function greaterpcdr) collect car x;
68  end;
69
70% ----- Certain stuff on maps -------------
71
72% A ring map is represented as a list
73%   {preimage_ring, image_ring, subst_list},
74% where subst_list is a substitution list {v1=ex1,v2=ex2,...} in
75% algebraic prefix form, i.e. looks like (list (equal var image) ...)
76
77symbolic operator preimage;
78symbolic procedure preimage(m,map);
79% Compute the preimage of an ideal m under a (polynomial) ring map.
80  if !*mode='algebraic then
81  begin map:=cdr reval map;
82     return preimage!*(reval m,
83        {ring_from_a first map, ring_from_a second map, third map});
84  end
85  else preimage!*(m,map);
86
87symbolic procedure preimage!*(m,map);
88% m and the result are given and returned in algebraic prefix form.
89  if not !*noetherian then
90        rederr"PREIMAGE only for noetherian term orders"
91  else begin scalar u,oldring,newring,oldnames;
92  if not eqcar(m,'list) then rederr"PREIMAGE only for ideals";
93  oldring:=first map; newring:=second map;
94  oldnames:=ring_names oldring;
95  setring!* ring_sum(newring,oldring);
96  u:=bas_renumber for each x in cdr third map collect
97  << if not member(second x,oldnames) then
98            typerr(second x,"var. name");
99     bas_make(0,dp_diff(dp_from_a second x,dp_from_a third x))
100  >>;
101  m:=matsum!* {dpmat_from_a m,dpmat_make(length u,0,u,nil,nil)};
102  m:=dpmat_2a eliminate!*(m,ring_names newring);
103  setring!* oldring;
104  return m;
105  end;
106
107symbolic operator ratpreimage;
108symbolic procedure ratpreimage(m,map);
109% Compute the preimage of an ideal m under a rational ring map.
110  if !*mode='algebraic then
111  begin map:=cdr reval map;
112  return ratpreimage!*(reval m,
113        {ring_from_a first map, ring_from_a second map, third map});
114  end
115  else ratpreimage!*(m,map);
116
117symbolic procedure ratpreimage!*(m,map);
118% m and the result are given and returned in algebraic prefix form.
119  if not !*noetherian then
120        rederr"RATPREIMAGE only for noetherian term orders"
121  else begin scalar u,oldring,newnames,oldnames,f,g,v,g0;
122  if not eqcar(m,'list) then rederr"RATPREIMAGE only for ideals";
123  oldring:=first map; v:=make_cali_varname();
124  newnames:=v . ring_names second map;
125  oldnames:=ring_names oldring; u:=append(oldnames,newnames);
126  setring!* ring_define(u,nil,'lex,for each x in u collect 1);
127  g0:=dp_fi 1;
128  u:=bas_renumber for each x in cdr third map collect
129  << if not member(second x,oldnames) then
130            typerr(second x,"var. name");
131     f:=simp third x; g:=dp_from_a prepf denr f;
132     f:=dp_from_a prepf numr f; g0:=dp_prod(g,g0);
133     bas_make(0,dp_diff(dp_prod(g,dp_from_a second x),f))
134  >>;
135  u:=bas_make(0,dp_diff(dp_prod(g0,dp_from_a v),dp_fi 1)) . u;
136  m:=matsum!* {dpmat_from_a m,dpmat_make(length u,0,u,nil,nil)};
137  m:=dpmat_2a eliminate!*(m,newnames);
138  setring!* oldring;
139  return m;
140  end;
141
142% ---- The ideals of affine resp. proj. points. The old stuff, but the
143% ---- algebraic interface now uses the linear algebra approach.
144
145symbolic procedure affine_points1!* m;
146  begin scalar names;
147  if length(names:=ring_names cali!=basering) neq length cadr m then
148        typerr(m,"coordinate matrix");
149  m:=for each x in cdr m collect
150         'list . for each y in pair(names,x) collect
151                {'plus,car y,{'minus,reval cdr y}};
152  m:=for each x in m collect dpmat_from_a x;
153  m:=matintersect!* m;
154  return m;
155  end;
156
157symbolic procedure scripts!=ideal u;
158  'list . for each x in cali_choose(u,2) collect
159        {'plus,{'times, car first x,cdr second x},
160        {'minus,{'times, car second x,cdr first x}}};
161
162symbolic procedure proj_points1!* m;
163  begin scalar names;
164  if length(names:=ring_names cali!=basering) neq length cadr m then
165        typerr(m,"coordinate matrix");
166  m:=for each x in cdr m collect scripts!=ideal pair(names,x);
167  m:=for each x in m collect interreduce!* dpmat_from_a x;
168  m:=matintersect!* m;
169  return m;
170  end;
171
172% ----- Affine and proj. monomial curves ------------
173
174symbolic operator affine_monomial_curve;
175symbolic procedure affine_monomial_curve(l,r);
176% l is a list of integers, R contains length l ring var. names.
177% Returns the generators of the monomial curve (t^i : i\in l) in R.
178  if !*mode='algebraic then
179        dpmat_2a affine_monomial_curve!*(cdr reval l,cdr reval r)
180  else affine_monomial_curve!*(l,r);
181
182symbolic procedure affine_monomial_curve!*(l,r);
183  if not numberlistp l then typerr(l,"number list")
184  else if length l neq length r then
185        rederr"number of variables doesn't match"
186  else begin scalar u,t0,v;
187    v:=list make_cali_varname();
188    r:=ring_define(r,{l},'revlex,l);
189    setring!* ring_sum(r,ring_define(v,degreeorder!* v,'lex,'(1)));
190    t0:=dp_from_a car v;
191    u:=bas_renumber for each x in pair(l,ring_names r) collect
192        bas_make(0,dp_diff(dp_from_a cdr x,dp_power(t0,car x)));
193    u:=dpmat_make(length u,0,u,nil,nil);
194    u:=(eliminate!*(u,v) where cali!=monset=ring_names cali!=basering);
195    setring!* r;
196    return dpmat_neworder(u,dpmat_gbtag u);
197    end;
198
199symbolic operator proj_monomial_curve;
200symbolic procedure proj_monomial_curve(l,r);
201% l is a list of integers, R contains length l ring var. names.
202% Returns the generators of the monomial curve
203% (s^(d-i)*t^i : i\in l) in R where d = max { x : x \in l}
204  if !*mode='algebraic then
205        dpmat_2a proj_monomial_curve!*(cdr reval l,cdr reval r)
206  else proj_monomial_curve!*(l,r);
207
208symbolic procedure proj_monomial_curve!*(l,r);
209  if not numberlistp l then typerr(l,"number list")
210  else if length l neq length r then
211        rederr"number of variables doesn't match"
212  else begin scalar u,t0,t1,v,d;
213    t0:=make_cali_varname(); t1:=make_cali_varname(); v:={t0,t1};
214    d:=listexpand(function max2,l);
215    r:=ring_define(r,degreeorder!* r,'revlex,for each x in r collect 1);
216    setring!* ring_sum(r,ring_define(v,degreeorder!* v,'lex,'(1 1)));
217    t0:=dp_from_a t0; t1:=dp_from_a t1;
218    u:=bas_renumber for each x in pair(l,ring_names r) collect
219        bas_make(0,dp_diff(dp_from_a cdr x,
220                dp_prod(dp_power(t0,car x),dp_power(t1,d-car x))));
221    u:=dpmat_make(length u,0,u,nil,nil);
222    u:=(eliminate!*(u,v) where cali!=monset=ring_names cali!=basering);
223    setring!* r;
224    return dpmat_neworder(u,dpmat_gbtag u);
225    end;
226
227% -- General Rees rings, associated graded rings, and related topics --
228
229symbolic operator blowup;
230symbolic procedure blowup(m,n,vars);
231% vars is a list of var. names for the ring R
232%       of the same length as dpmat_list n.
233% Returns an ideal J such that (S+R)/J == S/M [ N.t ]
234%       ( with S = the current ring )
235% is the blow up ring of the ideal N over S/M.
236% (S+R) is the new current ring.
237  if !*mode='algebraic then
238        dpmat_2a blowup!*(dpmat_from_a reval m,dpmat_from_a reval n,
239                cdr reval vars)
240  else blowup!*(m,n,vars);
241
242symbolic procedure blowup!*(m,n,vars);
243  if (dpmat_cols m > 0)or(dpmat_cols n > 0) then
244        rederr"BLOWUP defined only for ideals"
245  else if not !*noetherian then
246        rederr"BLOWUP only for noetherian term orders"
247  else begin scalar u,s,t0,v,r1;
248    if length vars neq dpmat_rows n then
249        rederr {"ring must have",dpmat_rows n,"variables"};
250    u:=for each x in dpmat_rowdegrees n collect mo_ecart cdr x;
251    r1:=ring_define(vars,list u,'revlex,u);
252    s:=ring_sum(cali!=basering,r1); v:=list(make_cali_varname());
253    setring!* ring_sum(s,ring_define(v,degreeorder!* v,'lex,'(1)));
254    t0:=dp_from_a car v;
255    n:=for each x in
256            pair(vars,for each y in dpmat_list n collect bas_dpoly y)
257            collect dp_diff(dp_from_a car x,
258                            dp_prod(dp_neworder cdr x,t0));
259    m:=bas_renumber append(bas_neworder dpmat_list m,
260            for each x in n collect bas_make(0,x));
261    m:=(eliminate!*(interreduce!* dpmat_make(length m,0,m,nil,nil),v)
262        where cali!=monset=nil);
263    setring!* s;
264    return dpmat_neworder(m,dpmat_gbtag m);
265    end;
266
267symbolic operator assgrad;
268symbolic procedure assgrad(m,n,vars);
269% vars is a list of var. names for the ring T
270%       of the same length as dpmat_list n.
271% Returns an ideal J such that (S+T)/J == (R/N + N/N^2 + ... )
272%       ( with R=S/M and S the current ring )
273% is the associated graded ring of the ideal N over R.
274% (S+T) is the new current ring.
275  if !*mode='algebraic then
276        dpmat_2a assgrad!*(dpmat_from_a reval m,dpmat_from_a reval n,
277                cdr reval vars)
278  else assgrad!*(m,n,vars);
279
280symbolic procedure assgrad!*(m,n,vars);
281  if (dpmat_cols m > 0)or(dpmat_cols n > 0) then
282        rederr"ASSGRAD defined only for ideals"
283  else begin scalar u;
284    u:=blowup!*(m,n,vars);
285    return matsum!* {u,dpmat_neworder(n,nil)};
286    end;
287
288symbolic operator analytic_spread;
289symbolic procedure analytic_spread m;
290% Returns the analytic spread of the ideal m.
291  if !*mode='algebraic then analytic_spread!* dpmat_from_a reval m
292  else analytic_spread!* m;
293
294symbolic procedure analytic_spread!* m;
295   if (dpmat_cols m>0) then rederr"ANALYTIC SPREAD only for ideals"
296   else (begin scalar r,m1,vars;
297   r:=ring_names cali!=basering;
298   vars:=for each x in dpmat_list m collect make_cali_varname();
299   m1:=blowup!*(dpmat_from_dpoly nil,m,vars);
300   return dim!* gbasis!* matsum!*{m1,dpmat_from_a('list . r)};
301   end) where cali!=basering=cali!=basering;
302
303symbolic operator sym;
304symbolic procedure sym(m,vars);
305% vars is a list of var. names for the ring R
306%       of the same length as dpmat_list M.
307% Returns an ideal J such that (S+R)/J == Sym(M)
308%       ( with S = the current ring )
309% is the symmetric algebra of M over S.
310% (S+R) is the new current ring.
311  if !*mode='algebraic then
312        dpmat_2a sym!*(dpmat_from_a m,cdr reval vars)
313  else sym!*(m,vars);
314
315symbolic procedure sym!*(m,vars);
316% The symmetric algebra of the dpmat m.
317   if not !*noetherian then
318        rederr"SYM only for noetherian term orders"
319   else begin scalar n,u,r1;
320    if length vars neq dpmat_rows m then
321        rederr {"ring must have",dpmat_rows m,"variables"};
322    cali!=degrees:=dpmat_coldegs m;
323    u:=for each x in dpmat_rowdegrees m collect mo_ecart cdr x;
324    r1:=ring_define(vars,list u,'revlex,u); n:=syzygies!* m;
325    setring!* ring_sum(cali!=basering,r1);
326    return mat2list!* interreduce!*
327                dpmat_mult(dpmat_neworder(n,nil),
328                        ideal2mat!* dpmat_from_a('list . vars));
329    end;
330
331% ----- Several short scripts ----------
332
333% ------ Minimal generators of an ideal or module.
334symbolic operator minimal_generators;
335symbolic procedure minimal_generators m;
336  if !*mode='algebraic then
337        dpmat_2a minimal_generators!* dpmat_from_a reval m
338  else minimal_generators!* m;
339
340symbolic procedure minimal_generators!* m;
341  car groeb_minimize(m,syzygies!* m);
342
343% ------- Symbolic powers of prime (or unmixed) ideals
344symbolic operator symbolic_power;
345symbolic procedure symbolic_power(m,d);
346  if !*mode='algebraic then
347        dpmat_2a symbolic_power!*(dpmat_from_a m,reval d)
348  else symbolic_power!*(m,d);
349
350symbolic procedure symbolic_power!*(m,d);
351  eqhull!* idealpower!*(m,d);
352
353% ---- non zero divisor property -----------
354
355put('nzdp,'psopfn,'scripts!=nzdp);
356symbolic procedure scripts!=nzdp m;
357  if length m neq 2 then rederr"Syntax : nzdp(dpoly,dpmat)"
358  else begin scalar f,b;
359    f:=reval car m; intf_get second m;
360    if null(b:=get(second m,'gbasis)) then
361        put(second m,'gbasis,b:=gbasis!* get(second m,'basis));
362    return if nzdp!*(dp_from_a f,b) then 'yes else 'no;
363    end;
364
365symbolic procedure nzdp!*(f,m);
366% Test dpoly f for a non zero divisor on coker m. m must be a gbasis.
367  submodulep!*(matqquot!*(m,f),m);
368
369endmodule; % scripts
370
371end;
372