FORM 4.3
diagrams.c
Go to the documentation of this file.
1
5/* #[ License : */
6/*
7 * Copyright (C) 1984-2022 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
37static 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
64int 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
185int 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 ) {
290NoWildArg:
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
348WORD 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
401WORD 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
428int 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 {
471notgoodvertex:
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 }
586for ( i = 0; i < nvert; i++ ) {
587MesPrint("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 }
619MesPrint("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
639int 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
645MesPrint("Entering DoShattering. partitions = %a",2*nvert,partitions);
646#endif
647restart:
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
665MesPrint("DoShattering: calculated the pointers");
666MesPrint("DoShattering: sizes: %a",nparts,AT.WorkPointer);
667MesPrint("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
682for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
683MesPrint("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
702for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
703MesPrint("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
716for ( ii = 0; ii < AT.WorkPointer[i]; ii++ ) {
717MesPrint("environ(%d,%d): %a",i,j,AT.WorkPointer[j],environ+ii*AT.WorkPointer[j]);
718}
719MesPrint("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
727MesPrint("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