FORM  4.2.1
diagrams.c
Go to the documentation of this file.
1 
5 /* #[ License : */
6 /*
7  * Copyright (C) 1984-2017 J.A.M. Vermaseren
8  * When using this file you are requested to refer to the publication
9  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
10  * This is considered a matter of courtesy as the development was paid
11  * for by FOM the Dutch physics granting agency and we would like to
12  * be able to track its scientific use to convince FOM of its value
13  * for the community.
14  *
15  * This file is part of FORM.
16  *
17  * FORM is free software: you can redistribute it and/or modify it under the
18  * terms of the GNU General Public License as published by the Free Software
19  * Foundation, either version 3 of the License, or (at your option) any later
20  * version.
21  *
22  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
23  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25  * details.
26  *
27  * You should have received a copy of the GNU General Public License along
28  * with FORM. If not, see <http://www.gnu.org/licenses/>.
29  */
30 /* #] License : */
31 /*
32  #[ Includes : diagrams.c
33 */
34 
35 #include "form3.h"
36 
37 static WORD one = 1;
38 
39 /*
40  #] Includes :
41  #[ CoCanonicalize :
42 
43  Syntax:
44  Canonicalize,mainoption,...
45  With the main options currently:
46  Canonicalize,topology,vertexfunction,edgefunction,OutDollar,extraoptions;
47  Canonicalize,polynomial,InDollar,set_or_setname_or_dollar,OutDollar,extraoptions;
48  The vertex function needs to have the format (assume it is called vx):
49  vx(p1,p2,-p3) or vx(-p1,p2,p3,-p4) etc.
50  The external lines have a vertex with only one line.
51  All momenta that form connections should be unique.
52  In principle the - signs are not relevant for the topology, but they
53  may exist already in the remaining part of the diagram. They are also
54  part of the canonical form.
55  The edge function can be used in different ways, depending on the options.
56  The extraoption(s) should be nonnegative integers or $variables that evaluate
57  into nonnegative integers (integers less than 2^31).
58  The Indollar variable contains the polynomial to be canonicalized.
59  The OutDollar variable should be the name of a $-variable (as in $out) which
60  will be filled with a replace_ function. The canonicalization can then
61  be executed in the whole term with the Multiply $out; command.
62 */
63 
64 int CoCanonicalize(UBYTE *s)
65 {
66  WORD args[10], *a, num;
67  UBYTE *t, c;
68  args[0] = TYPECANONICALIZE;
69  a = args+2;
70  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
71  t = s; while ( FG.cTable[*s] == 0 ) s++;
72  c = *s; *s = 0;
73  if ( StrICmp(t,(UBYTE *)("topology")) == 0 ) {
74  *s = c; *a++ = 0;
75  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
76  s = GetFunction(s,a);
77  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
78  s = GetFunction(s,a+1);
79  if ( *a == 0 || a[1] == 0 ) return(1);
80  a += 2;
81  }
82  else if ( StrICmp(t,(UBYTE *)("polynomial")) == 0 ) {
83  *s = c; *a++ = 1;
84  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
85 /*
86  Now get the name of the input $-variable
87 */
88  if ( *s != '$' ) {
89  MesPrint("&Canonicalize statement needs a $-variable for its input.");
90  return(1);
91  }
92  s++; t = s; while ( FG.cTable[*s] < 2 ) s++;
93  c = *s; *s = 0;
94  if ( GetName(AC.dollarnames,t,&num,NOAUTO) == CDOLLAR ) *a++ = num;
95  else { *a++ = AddDollar(t,DOLINDEX,&one,1); }
96  *s = c;
97 /*
98  Now the set
99 */
100  if ( *s == '{' ) {
101  t = s+1; SKIPBRA2(s)
102  c = *s; *s = 0;
103  *a++ = DoTempSet(t,s);
104  *s++ = c;
105  }
106  else if ( FG.cTable[*s] == 0 || *s == '[' ) {
107  t = s;
108  if ( ( s = SkipAName(s) ) == 0 ) {
109  MesPrint("&Illegal name for set in Canonicalize statement: %s",t);
110  return(1);
111  }
112  c = *s; *s = 0;
113  if ( GetName(AC.varnames,t,a,WITHAUTO) == CSET ) {
114  if ( Sets[*a].type != CSYMBOL ) {
115  MesPrint("&In Canonicalize: %s is not a set of symbols.",t);
116  return(1);
117  }
118  }
119  else {
120  MesPrint("&In Canonicalize: %s is not a set.",t);
121  return(1);
122  }
123  *s = c; a++;
124  }
125  else if ( *s == '$' ) {
126  s++; t = s; while ( FG.cTable[*s] < 2 ) s++;
127  c = *s; *s = 0;
128  if ( GetName(AC.dollarnames,t,&num,NOAUTO) == CDOLLAR ) *a++ = -num-2;
129  else {
130  MesPrint("&In Canonicalize: %s is undefined.",t-1);
131  return(1);
132  }
133  *s = c;
134  }
135  else {
136  MesPrint("&In Canonicalize: Illegal third(=set) argument.");
137  return(1);
138  }
139  }
140  else {
141  MesPrint("&Unrecognized option in Canonicalize statement: %s",t);
142  return(1);
143  }
144 /*
145  Now get the name of the output $-variable
146 */
147  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
148  if ( *s != '$' ) {
149  MesPrint("&Canonicalize statement needs a $-variable for its output.");
150  return(1);
151  }
152  s++; t = s; while ( FG.cTable[*s] < 2 ) s++;
153  c = *s; *s = 0;
154  if ( GetName(AC.dollarnames,t,&num,NOAUTO) == CDOLLAR ) *a++ = num;
155  else { *a++ = AddDollar(t,DOLINDEX,&one,1); }
156  *s = c;
157 /*
158  Now the options. At the moment we just do one of them.
159  (the first extra option is relevant to determine the use of the edge function)
160  In the future we may have to be more flexible.
161 */
162  a[0] = 0; /* default value */
163  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
164  if ( *s != 0 ) {
165  s = GetNumber(s,a);
166  if ( *a == -1 ) return(1);
167  while ( *s == ',' || *s == '\t' || *s == ' ' ) s++;
168  a++;
169  }
170 /*
171  Now complete the args string and put it in the compiler buffer
172 */
173  args[1] = a-args;
174  AddNtoL(args[1],args);
175  return(0);
176 }
177 
178 /*
179  #] CoCanonicalize :
180  #[ DoCanonicalize :
181 
182  Does the canonicalization. The output term overwrites the input term.
183 */
184 
185 int DoCanonicalize(PHEAD WORD *term, WORD *params)
186 {
187  WORD args[10];
188  int i;
189 /*
190  First check whether we need to expand dollars;
191 */
192  for ( i = 0; i < params[1]; i++ ) args[i] = params[i];
193  if ( args[2] == 0 ) { /* topology */
194  for ( i = 3; i < 5; i++ ) {
195  if ( args[i] < 0 ) { /* This is a dollar */
196  args[i] = DolToFunction(BHEAD -args[i]-2);
197  if ( args[i] == 0 ) {
198  MLOCK(ErrorMessageLock);
199  MesPrint("Value of $-variable in Canonicalize statement should be a function.");
200  MUNLOCK(ErrorMessageLock);
201  Terminate(-1);
202  }
203  }
204  }
205  for ( i = 6; i < args[1]; i++ ) { /* Extra options */
206  if ( args[i] < 0 ) { /* This is a dollar */
207  args[i] = DolToNumber(BHEAD -args[i]-2);
208  if ( args[i] < 0 ) {
209  MLOCK(ErrorMessageLock);
210  MesPrint("Value of $-variable in Canonicalize statement should be a nonnegative number < %l.",(LONG)MAXPOSITIVE);
211  MUNLOCK(ErrorMessageLock);
212  Terminate(-1);
213  }
214  }
215  }
216  switch ( args[6] ) {
217  case 1: {/* pass the vertex and edge functions. */
218  WORD *tstop, *t, *tedge, *te;
219  tstop = term + *term; tstop -= ABS(tstop[-1]);
220  t = term+1;
221  tedge = AT.WorkPointer; te = tedge+1;
222  while ( t < tstop ) {
223  if ( *t != args[3] && *t != args[4] ) { t += t[1]; continue; }
224  for ( i = 0; i < t[1]; i++ ) te[i] = t[i];
225  te += t[1]; t += t[1];
226  }
227  *te++ = 1; *te++ = 1; *te++ = 3;
228  tedge[0] = te-tedge;
229  AT.WorkPointer = te;
230 /*
231  DoVertexCanonicalize(BHEAD term,tedge,args[3],args[4],args[5],args[6]);
232 */
233  AT.WorkPointer = tedge;
234  } break;
235  case 2: {/* pass the edge functions only */
236  WORD *tstop, *t, *tedge, *te;
237  tstop = term + *term; tstop -= ABS(tstop[-1]);
238  t = term+1;
239  tedge = AT.WorkPointer; te = tedge+1;
240  while ( t < tstop ) {
241  if ( *t != args[4] ) { t += t[1]; continue; }
242  for ( i = 0; i < t[1]; i++ ) te[i] = t[i];
243  te += t[1]; t += t[1];
244  }
245  *te++ = 1; *te++ = 1; *te++ = 3;
246  tedge[0] = te-tedge;
247  AT.WorkPointer = te;
248 /*
249  DoEdgeCanonicalize(BHEAD term,tedge,args[5]);
250 */
251  AT.WorkPointer = tedge;
252  } break;
253  default: {
254  DoTopologyCanonicalize(BHEAD term,args[3],args[4],args+5);
255  } break;
256  }
257 /*
258  Call here the topology canonicalization
259  We will have the arguments:
260  args[3]: The function used as vertex function
261  args[4]: The function used as edge function
262  args[5]: The number of the dollar to be used for the output
263  args[6]: Potentially other options (like saying how to use args[4]).
264  term: The term in which the topology resides.
265 */
266 
267  }
268  else if ( args[2] == 1 ) { /* polynomial */
269  WORD *symlist, nsymlist;
270  for ( i = 6; i < args[1]; i++ ) { /* Extra options */
271  if ( args[i] < 0 ) { /* This is a dollar */
272  args[i] = DolToNumber(BHEAD -args[i]-2);
273  if ( args[i] < 0 ) {
274  MLOCK(ErrorMessageLock);
275  MesPrint("Value of $-variable in Canonicalize statement should be a nonnegative number < %l.",(LONG)MAXPOSITIVE);
276  MUNLOCK(ErrorMessageLock);
277  Terminate(-1);
278  }
279  }
280  }
281 /*
282  Now we sort out the set. We create a pointer to the list of set
283  elements, and we determine the number of elements in the set.
284 */
285  symlist = AT.WorkPointer;
286  if ( args[4] < -1 ) { /* Dollar that should expand into a list of symbols */
287  DOLLARS d = Dollars - args[4] - 2;
288  WORD *ds, *insym;
289  if ( d->type != DOLWILDARGS ) {
290 NoWildArg:
291  MLOCK(ErrorMessageLock);
292  MesPrint("Value of $-variable in Canonicalize statement should be a argument wildcard of symbol arguments.");
293  MUNLOCK(ErrorMessageLock);
294  Terminate(-1);
295  }
296  insym = symlist; ds = d->where+1;
297  while ( *ds ) {
298  if ( *ds != -SYMBOL ) goto NoWildArg;
299  *insym++ = ds[1];
300  ds += 2;
301  }
302  nsymlist = insym-symlist;
303  }
304  else { /*if ( args[4] >= 0 ) */
305  WORD *ss, *sy, n;
306  ss = (WORD *)(AC.SetElementList.lijst)+Sets[args[4]].first;
307  nsymlist = n = Sets[args[4]].last-Sets[args[4]].first;
308  sy = symlist = AT.WorkPointer;
309  NCOPY(sy,ss,n);
310  }
311  AT.WorkPointer = symlist+nsymlist;
312 /*
313  Call here the polynomial canonicalization
314  We will have the arguments:
315  args[3]: The number of the dollar to be used for the input
316  symlist: an array of symbols
317  nsymlist: the number of symbols in symlist
318  args[5]: The number of the dollar to be used for the output
319  args[6]: Potentially other options.
320 */
321 /*
322  DoPolynomialCanonicalize(BHEAD args[3],symlist,nsymlist,args[5],args[6]);
323 */
324  AT.WorkPointer = symlist;
325  }
326  return(0);
327 }
328 
329 /*
330  #] DoCanonicalize :
331  #[ GenTopologies :
332 
333  This function has the syntax
334  topologies_(nloops, Number of loops
335  nlegs, Number of legs
336  setvertexsizes, A set which tells which vertices are allowed like {3,4}.
337  set_extmomenta, The name of a set with the external momenta
338  set_intmomenta The name of a set with the internal momenta
339  [,options])
340  The output will be using the built in functions vertex_ and edge_.
341 
342  The test for whether this function can be evaluated is in TestSub (inside
343  file proces.c) (search for the string TOPOLOGIES).
344  This passes the code -15 in AN.TeInFun to Generator, which then calls
345  the GenTopologies routine.
346 */
347 
348 WORD GenTopologies(PHEAD WORD *term,WORD level)
349 {
350  WORD *t1, *tt1, *tstop, *t, *tt;
351  WORD *oldworkpointer = AT.WorkPointer;
352  WORD option1 = 0, option2 = 0, setoption = -1;
353  WORD retval;
354 /*
355 
356  We have to go through the testing procedure again, because there could
357  be more than one topologies_ function and not all have to be expandable.
358 */
359  tstop = term+*term; tstop -= ABS(tstop[-1]);
360  tt = term+1;
361  while ( tt < tstop ) {
362  t = tt; tt = t+t[1];
363  if ( *t != TOPOLOGIES ) continue;
364  tt = t + t[1]; t1 = t + FUNHEAD;
365  if ( t1+10 > tt || *t1 != -SNUMBER || t1[1] < 0 || /* loops */
366  t1[2] != -SNUMBER || ( t1[3] < 0 && t1[3] != -2 ) ||/* legs */
367  t1[4] != -SETSET || Sets[t1[5]].type != CNUMBER || /* set vertices */
368  t1[6] != -SETSET || Sets[t1[7]].type != CVECTOR || /* outvectors */
369  t1[8] != -SETSET || Sets[t1[9]].type != CVECTOR ) continue;
370  tt1 = t1 + 10;
371  if ( tt1+2 <= tt && tt1[0] == -SETSET ) {
372  if ( Sets[t1[5]].last-Sets[t1[5]].first !=
373  Sets[tt1[1]].last-Sets[tt1[1]].first ) continue;
374  setoption = tt1[1]; tt1 += 2;
375  }
376  if ( tt1+2 <= tt && tt1[0] == -SNUMBER ) { option1 = tt1[1]; tt1 += 2; }
377  if ( tt1+2 <= tt && tt1[0] == -SNUMBER ) { option2 = tt1[1]; tt1 += 2; }
378  AT.setinterntopo = t1[9];
379  AT.setexterntopo = t1[7];
380  AT.TopologiesTerm = term;
381  AT.TopologiesStart = t;
382  AT.TopologiesLevel = level;
383  AT.TopologiesOptions[0] = option1;
384  AT.TopologiesOptions[1] = option2;
385  retval = GenerateTopologies(BHEAD t1[1],t1[3],t1[5],setoption);
386  AT.WorkPointer = oldworkpointer;
387  return(retval);
388  }
389  MLOCK(ErrorMessageLock);
390  MesPrint("Internal error: topologies_ function not encountered.");
391  MUNLOCK(ErrorMessageLock);
392  return(-1);
393 
394 }
395 
396 /*
397  #] GenTopologies :
398  #[ GenDiagrams :
399 */
400 
401 WORD GenDiagrams(PHEAD WORD *term,WORD level)
402 {
403 #ifdef WITHPTHREADS
404  DUMMYUSE(B)
405 #endif
406  DUMMYUSE(term)
407  DUMMYUSE(level)
408  return(0);
409 }
410 
411 /*
412  #] GenDiagrams :
413  #[ DoTopologyCanonicalize :
414 
415  term: The term
416  vert: the vertex function
417  edge: the edge function
418  args[0]: the number of the output dollar
419  args[1]: options
420  return value should be zero if all is correct.
421 
422  The external lines connect to an 'external vertex' which has only one line.
423  The vertices are of a type: vertex_(p1,p2,-p3)*vertex_(p3,p4,p5) etc.
424  The edges indicate noninteger powers of the lines:
425  edge_(p1,2) means 1/p1.p1^(2*ep)
426 */
427 
428 int DoTopologyCanonicalize(PHEAD WORD *term,WORD vert,WORD edge,WORD *args)
429 {
430  int nvert = 0, nvert2, i, ii, jj, flipnames = 0, nparts, level, num;
431  WORD *tstop, *t, *tt, *tend, *td;
432  WORD *oldworkpointer = AT.WorkPointer;
433  WORD *termcopy = TermMalloc("TopologyCanonize1");
434  WORD *vet= TermMalloc("TopologyCanonize2");
435  WORD *partition, *environ, *connect, *pparts, *p;
436 /*
437  WORD *pparts;
438 */
439  WORD momenta[150],flips[50],nmomenta = 0, nflips = 0;
440 /*
441  Step one: the vertices should get a number. We copy the term for this.
442  We need a high number for the vertex function to make sure
443  that it comes after the edge function in the sorting.
444 */
445  if ( args[0] < args[1] ) { flipnames = 1; }
446  tend = term + *term; tend -= ABS(tend[-1]); t = term+1; tt = termcopy+1;
447  while ( t < tend ) {
448  if ( *t == vert ) {
449  for ( i = FUNHEAD; i < t[1]; i += 2 ) {
450  if ( t[i] == -VECTOR || ( t[i] == -INDEX && t[i+1] < 0 ) ) {
451  momenta[nmomenta++] = -VECTOR;
452  momenta[nmomenta++] = t[i+1];
453  }
454  else if ( t[i] == -MINVECTOR ) {
455  momenta[nmomenta++] = -MINVECTOR;
456  momenta[nmomenta++] = t[i+1];
457  }
458  else goto notgoodvertex;
459  momenta[nmomenta++] = nvert;
460  }
461  ii = FUNHEAD; i = t[1]-FUNHEAD;
462  NCOPY(tt,t,ii)
463  if ( flipnames ) tt[-FUNHEAD] = edge;
464  tt[-FUNHEAD+1] += 2;
465  *tt++ = -CNUMBER; *tt++ = nvert++;
466  }
467  else if ( *t == edge && flipnames ) {
468  i = t[1] - 1; *tt++ = vert; t++;
469  }
470  else {
471 notgoodvertex:
472  i = t[1];
473  }
474  NCOPY(tt,t,i)
475  }
476  while ( t < tend ) *tt++ = *t++;
477  termcopy[0] = tt - termcopy;
478  if ( flipnames ) EXCH(edge,vert)
479  nvert2 = nvert*nvert;
480 /*
481  Sort the momenta. Keep the sign order.
482 */
483  for ( i = 0; i < nmomenta-3; i+=3 ) {
484  jj = i;
485  while ( jj >= 0 && momenta[jj+4] > momenta[jj+1] ) {
486  EXCH(momenta[jj+5],momenta[jj+2])
487  EXCH(momenta[jj+4],momenta[jj+1])
488  EXCH(momenta[jj+3],momenta[jj])
489  jj -= 3;
490  }
491  }
492 /*
493  Step two: make now the edge functions in the proper notation.
494 */
495  t = vet+1;
496  for ( i = 0; i < nmomenta; i += 6 ) {
497  if ( momenta[i] == -VECTOR && momenta[i+3] == -MINVECTOR
498  && momenta[i+1] == momenta[i+4] ) {
499  }
500  else if ( momenta[i] == -MINVECTOR && momenta[i+3] == -VECTOR
501  && momenta[i+1] == momenta[i+4] ) {
502  flips[nflips++] = momenta[i+1];
503  DUMMYUSE(flips[nflips-1]);
504  }
505  else { /* something wrong with the momenta */
506  MLOCK(ErrorMessageLock);
507  MesPrint("No momentum conservation or wrong momenta in Canonicalize statement");
508  MUNLOCK(ErrorMessageLock);
509  Terminate(-1);
510  }
511  *t++ = EDGE; *t++ = FUNHEAD+10; FILLFUN(t)
512  *t++ = -SNUMBER; *t++ = momenta[i+2];
513  *t++ = -SNUMBER; *t++ = momenta[i+5];
514  *t++ = -VECTOR; *t++ = momenta[i+1];
515  *t++ = -SNUMBER; *t++ = 0; /* provisional power/color, multiple of ep */
516  *t++ = -SNUMBER; *t++ = 0; /* provisional power/color, integer */
517  }
518  tend = t;
519  *t++ = 1; *t++ = 1; *t++ = 3; vet[0] = t-vet; *t = 0;
520 /*
521  Now the powers of the denominators
522 */
523  tstop = termcopy+*termcopy; tstop -= ABS(tstop[-1]); td = termcopy+1;
524  while ( td < tstop ) {
525  if ( *td == edge && td[1] == FUNHEAD+4 ) {
526  if ( td[FUNHEAD+2] == -SNUMBER && ( td[FUNHEAD] == -VECTOR || td[FUNHEAD] == -INDEX
527  || td[FUNHEAD] == -MINVECTOR ) ) {}
528  else {
529  MLOCK(ErrorMessageLock);
530  MesPrint("Illegal argument in edge function in Canonicalize statement");
531  MUNLOCK(ErrorMessageLock);
532  Terminate(-1);
533  }
534  tt = vet+1;
535  while ( tt < tend ) {
536  if ( tt[FUNHEAD+5] == td[FUNHEAD+1] ) { tt[FUNHEAD+7] = td[FUNHEAD+3]; break; }
537  tt += tt[1];
538  }
539  }
540  else if ( *td == DOTPRODUCT ) break;
541  td += td[1];
542  }
543  if ( td < tstop ) {
544  tt = vet+1;
545  while ( tt < tend ) {
546 /*
547  tt[FUNHEAD+5] is a vector. We look for dotproducts with twice tt[FUNHEAD+5]
548 */
549  for ( i = 2; i < td[1]; i += 3 ) {
550  if ( td[i] == tt[FUNHEAD+5] && td[i+1] == tt[FUNHEAD+5] ) {
551  tt[FUNHEAD+9] = td[i+2];
552  break;
553  }
554  }
555  tt += tt[1];
556  }
557  }
558  Normalize(BHEAD vet);
559 /*
560  Now we have a term `vet' in the proper notation and we can start.
561  To keep track of the shattering we use an array of 2*nvert.
562  each entry is Number,marker
563  When the marker is zero, the vertices are in the same partition.
564  For the environment we need a matrix that is nvert x nvert
565  At the same time we keep the connectivity matrix, because that will
566  save much time later.
567  The partitions are stored in a matrix as well. This allows them to be
568  treated as a stack. The entries are separated by 0 if they belong to
569  the same part, and by a 1 when they belong to different parts.
570 */
571  partition = AT.WorkPointer; AT.WorkPointer += 2*nvert2;
572  for ( i = 0; i < nvert; i++ ) { partition[2*i] = i; partition[2*i+1] = 0; }
573  partition[2*i-1] = -1; /* end of the first part which is currently all vertices */
574  nparts = 1;
575 
576  connect = AT.WorkPointer; AT.WorkPointer += nvert2;
577  for ( i = 0; i < nvert2; i++ ) connect[i] = 0;
578  tstop = vet+*vet; tstop -= ABS(tstop[-1]); t = vet+1;
579  while ( t < tstop ) {
580  if ( *t == EDGE ) {
581  connect[t[FUNHEAD+1]*nvert+t[FUNHEAD+3]]++;
582  connect[t[FUNHEAD+3]*nvert+t[FUNHEAD+1]]++;
583  }
584  t += t[1];
585  }
586 for ( i = 0; i < nvert; i++ ) {
587 MesPrint("connectivity: %d -- %a",i,nvert,connect+i*nvert);
588 }
589 /*
590  Create the environment matrix and sort it.
591 */
592  environ = AT.WorkPointer; AT.WorkPointer += nvert2;
593 /*
594  And now the refinement process starts.
595 */
596  WantAddPointers(nvert+1);
597  for ( i = 0; i < nvert2; i++ ) environ[i] = 0;
598  level = 0;
599  pparts = partition;
600  while ( nparts < nvert ) {
601  nparts = DoShattering(BHEAD connect,environ,pparts,nvert);
602  if ( nparts < nvert ) { /* raise level and make a copy and split a part */
603  p = pparts + 2*nvert;
604  level++;
605  for ( i = 0; i < 2*nvert; i++ ) p[i] = pparts[i];
606  for ( ii = 0; ii < 2*nvert; ii += 2 ) {
607  if ( p[ii+1] == 0 ) { /* found a part with more than one */
608  num = 2; i = ii+2;
609  while ( p[i+1] == 0 ) { num++; i += 2; }
610 
611 
612  p[ii+1] = -1; pparts = p;
613  break;
614  }
615  }
616 
617 
618  }
619 MesPrint("partition: %d -- %a",nparts,2*nvert,pparts);
620 
621  }
622 /*
623  Just for now
624 */
625  PutTermInDollar(vet,args[0]);
626 
627 
628  TermFree(vet,"TopologyCanonize2");
629  TermFree(termcopy,"TopologyCanonize1");
630  AT.WorkPointer = oldworkpointer;
631  return(0);
632 }
633 
634 /*
635  #] DoTopologyCanonicalize :
636  #[ DoShattering :
637 */
638 
639 int DoShattering(PHEAD WORD *connect, WORD *environ, WORD *partitions, WORD nvert)
640 {
641  int nparts, i, j, ii, jj, iii, jjj, newmarker;
642  WORD **p = AT.pWorkSpace + AT.pWorkPointer, *part, *endpart;
643  WORD *poin1, *poin2;
644 #ifdef SHATBUG
645 MesPrint("Entering DoShattering. partitions = %a",2*nvert,partitions);
646 #endif
647 restart:
648 /*
649  Determine the number of parts
650  p will be an array with pointers to the parts.
651  We made space for this array in the calling routine and because this
652  routine is not calling any other routines we do not need to raise
653  the pointer in this stack (AT.pWorkPointer).
654 */
655  nparts = 0; newmarker = 0;
656  part = partitions; endpart = part + 2*nvert;
657  p[0] = part;
658  while ( part < endpart ) {
659  if ( part[1] != 0 ) { p[++nparts] = part+2; }
660  part += 2;
661  }
662  for ( i = 0; i < nparts; i++ )
663  AT.WorkPointer[i] = (p[i+1]-p[i])/2;
664 #ifdef SHATBUG
665 MesPrint("DoShattering: calculated the pointers");
666 MesPrint("DoShattering: sizes: %a",nparts,AT.WorkPointer);
667 MesPrint("DoShattering: p[0]: %a, p[1]: %a",6,p[0],6,p[1]);
668 #endif
669  for ( i = 0; i < nparts; i++ ) {
670  if ( AT.WorkPointer[i] > 1 ) {
671  for ( j = 0; j < nparts; j++ ) {
672 /*
673  Shatter part i wrt part j.
674  if there is action, go to restart.
675 */
676  for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
677  for ( jj = 0; jj < AT.WorkPointer[j]; jj++ ) {
678  environ[ii*AT.WorkPointer[j]+jj] += connect[p[i][2*ii]*nvert+p[j][2*jj]];
679  }
680  }
681 #ifdef SHATBUG
682 for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
683 MesPrint("Environ(%d,%d): %a",i,j,AT.WorkPointer[j],environ+ii*AT.WorkPointer[j]);
684 }
685 #endif
686 /*
687  Sort the rows internally, then sort the rows wrt each other
688  and finally place new markers. If a new marker, we restart.
689  Don't forget to clean up the environ array.
690 */
691  for ( ii = 0; ii < nvert; ii++ ) {
692  poin1 = environ+ii*AT.WorkPointer[j];
693  for ( jj = 0; jj < AT.WorkPointer[j]-1; jj++ ) {
694  jjj = jj;
695  while ( jjj >= 0 && poin1[jjj+1] > poin1[jjj] ) {
696  EXCH(poin1[jjj+1],poin1[jjj])
697  jjj--;
698  }
699  }
700  }
701 #ifdef SHATBUG
702 for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
703 MesPrint("environ(%d,%d): %a",i,j,AT.WorkPointer[j],environ+ii*AT.WorkPointer[j]);
704 }
705 #endif
706  for ( ii = 0; ii < AT.WorkPointer[i]-1; ii++ ) {
707  poin2 = environ+ii*AT.WorkPointer[j]; poin1 = poin2+AT.WorkPointer[j];
708  iii = ii;
709  while ( iii >= 0 && ( CmpArray(poin1,poin2,AT.WorkPointer[j]) < 0 ) ) {
710  EXCHN(poin2,poin1,AT.WorkPointer[j])
711  EXCH(p[i][2*iii+2],p[i][2*iii])
712  iii--; poin1 = poin2; poin2 = poin1-AT.WorkPointer[j];
713  }
714  }
715 #ifdef SHATBUG
716 for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
717 MesPrint("environ(%d,%d): %a",i,j,AT.WorkPointer[j],environ+ii*AT.WorkPointer[j]);
718 }
719 MesPrint("partitions = %a",2*nvert,partitions);
720 #endif
721  for ( ii = 0; ii < AT.WorkPointer[i]-1; ii++ ) {
722  poin2 = environ+ii*AT.WorkPointer[j]; poin1 = poin2+AT.WorkPointer[j];
723  if ( CmpArray(poin1,poin2,AT.WorkPointer[j]) == 0 ) continue;
724  p[i][2*ii+1] = -1; nparts++; newmarker++;
725  }
726 #ifdef SHATBUG
727 MesPrint("partitions = %a",2*nvert,partitions);
728 #endif
729 /*
730  Clear environ. This is probably faster than just clearing the whole array.
731  Maybe in the future a test could be done on nvert to decide how to clear.
732 */
733  for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
734  for ( jj = 0; jj < AT.WorkPointer[j]; jj++ ) {
735  environ[ii*AT.WorkPointer[j]+jj] = 0;
736  }
737  }
738  if ( newmarker ) { goto restart; }
739  }
740  }
741  }
742  return(nparts);
743 }
744 
745 /*
746  #] DoShattering :
747 */
int AddNtoL(int n, WORD *array)
Definition: comtool.c:288