FORM  4.1
optim.c
Go to the documentation of this file.
1 
5 /* #[ License : */
6 /*
7  * Copyright (C) 1984-2013 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:
33 
34  Routines for the optimization of FORTRAN or C output.
35  We have to do a number of things.
36  1: map all objects to a special type of scalars. In terms of these
37  scalars the terms are a bit simpler:
38  size numsca sca1,pow1,sca2,pow2,...,coeff
39  If size < 0 we have to skip -size words
40  A zero ends a subexpression, two zeroes the whole expression
41  2: Apply the routines from the codegenerator of spider.
42  The problem is where to leave the intermediate results.
43  We put them in a buffer and try to be careful.
44  Also the temporary assignments are put in a buffer.
45  We need a system of pointers for this.
46  The name of the vector with intermediate results should be selected
47  as a variable with as few letters as possible. Try which 1 character
48  name is not in use. Next try two character names etc. Assume z.
49  After that all intermediate variables have the name z(i).
50 */
51 
52 #include "form3.h"
53 
54 static char *useprop =
55 "Only Symbols, DotProducts, User-defined Functions and vectors with indices\n\
56 are allowed in optimized FORTRAN/C output";
57 
58 typedef struct ScAlAr {
59  WORD *buffer;
60  WORD *pointer;
61  WORD *top;
62  LONG bufsize;
63  LONG numterms;
64 } SCALAR;
65 
66 static SCALAR *scabuffer = 0;
67 static LONG scasize = 0, scanumber = -1;
68 
69 static LONG *objoffset = 0, objsize = 0;
70 static WORD *objbuffer = 0, *objpointer = 0, *objtop = 0;
71 static int numobjects = 0, numobjoffsets = 0;
72 static int iniobjects = 0;
73 
74 static UBYTE scratchname[10];
75 static WORD **helppointers = 0, **sortpointers = 0;
76 static LONG sizepointers = 0;
77 
78 static LONG *multiplicities = 0, multiplicitysize = 0;
79 
80 /*
81  #] Includes:
82  #[ Optimize:
83 */
84 
85 int Optimize(WORD numexpr)
86 {
87  int retval;
88 /* WORD powmax; */
89  LONG j;
90  if ( FindScratchName() ) return(-1);
91  if ( LoadOpti(numexpr) < 0 ) return(-1);
92  iniobjects = numobjects;
93 /*
94  Now we have the terms in the buffer.
95  We pass the various optimization steps.
96 */
97  HuntNumFactor(0L,0,0);
98 /*
99  powmax = 0;
100  {
101  WORD i;
102  for ( j = 0; j < scanumber; j++ ) {
103  SortOpti(j);
104  i = MaxPowerOpti(j);
105  if ( i > powmax ) powmax = i;
106  }
107  for ( i = powmax; i > 0; i-- ) {
108  for ( j = 0; j < scanumber; j++ ) {
109  HuntPairs(j,i);
110  }
111  }
112  }
113 */
114  for ( j = 0; j < scanumber; j++ ) HuntBrackets(j);
115  CombiOpti();
116 /*
117  Finally we print the results.
118 */
119  retval = PrintOptima(numexpr);
120  CleanOptiBuffer();
121  return(retval);
122 }
123 
124 /*
125  #] Optimize:
126  #[ LoadOpti:
127 */
128 
129 int LoadOpti(WORD numexpr)
130 {
131  GETIDENTITY
132  WORD *term, *t, *tstop, *m, *mstop, *oldwork = AT.WorkPointer, *tbuf, *w, nc;
133  WORD pow = 0;
134  LONG objnum = 0;
135  int i;
136  tbuf = oldwork;
137  AT.WorkPointer = (WORD *)(((UBYTE *)(tbuf)) + 2*AM.MaxTer);
138  term = AT.WorkPointer;
139  while ( GetTerm(BHEAD term) > 0 ) {
140  AT.WorkPointer = term + *term;
141  Normalize(BHEAD term);
142  AT.WorkPointer = term + *term;
143 /*
144  First hunt objects. Some we don't accept!
145  We do accept symbols, dotproducts, vectors, user defined functions.
146  We do not accept indices, system functions and the rest.
147 */
148  w = tbuf + 1;
149  GETSTOP(term,tstop);
150  t = term + 1; mstop = term + *term - 1;
151  nc = *mstop;
152  *w++ = LNUMBER; *w++ = ABS(nc)+2; *w++ = nc; m = tstop;
153  while ( m < mstop ) *w++ = *m++;
154  while ( t < tstop ) {
155  mstop = t + t[1];
156  if ( *t == SYMBOL ) {
157  m = t + 2;
158  while ( m < mstop ) {
159  objnum = PutObject(m,SYMBOL);
160  m += 2; pow = m[-1];
161  *w++ = objnum >> BITSINWORD;
162  *w++ = objnum & WORDMASK;
163  *w++ = pow;
164  }
165  }
166  else if ( *t == DOTPRODUCT ) {
167  m = t + 2;
168  while ( m < mstop ) {
169  objnum = PutObject(m,DOTPRODUCT);
170  m += 3; pow = m[-1];
171  *w++ = objnum >> BITSINWORD;
172  *w++ = objnum & WORDMASK;
173  *w++ = pow;
174  }
175  }
176  else if ( *t == VECTOR ) {
177  m = t + 2;
178  while ( m < mstop ) {
179  objnum = PutObject(m,VECTOR);
180  m += 2;
181  pow = 1;
182  while ( m < mstop && m[0] == m[-2] && m[1] == m[-1] ) {
183  pow++; m += 2;
184  }
185  *w++ = objnum >> BITSINWORD;
186  *w++ = objnum & WORDMASK;
187  *w++ = pow;
188  }
189  }
190  else if ( *t > MAXBUILTINFUNCTION ) {
191  objnum = PutObject(t,FUNCTION);
192  pow = 1;
193  while ( mstop < tstop ) {
194  m = mstop;
195  for ( i = 0; i < t[1]; i++ ) {
196  if ( m[i] != t[i] ) break;
197  }
198  if ( i < t[1] ) break;
199  pow++;
200  t = mstop; mstop = t + t[1];
201  }
202  *w++ = objnum >> BITSINWORD;
203  *w++ = objnum & WORDMASK;
204  *w++ = pow;
205  }
206  else {
207  MLOCK(ErrorMessageLock);
208  MesPrint("Problems with expression %s:",EXPRNAME(numexpr));
209  MesPrint(useprop);
210  MUNLOCK(ErrorMessageLock);
211  AT.WorkPointer = oldwork;
212  return(-1);
213  }
214  t = mstop;
215  }
216  tbuf[0] = w-tbuf;
217  AddToOpti(tbuf,0);
218  }
219  AT.WorkPointer = oldwork;
220 /*
221  We sort the terms internally. This makes life easier later.
222 */
223  t = scabuffer->buffer;
224  while ( t < scabuffer->pointer ) {
225  NormOpti(t);
226  t += *t;
227  }
228  return(0);
229 }
230 
231 /*
232  #] LoadOpti:
233  #[ CleanOptiBuffer:
234 */
235 
236 void CleanOptiBuffer()
237 {
238  LONG i;
239  if ( scanumber > 0 ) {
240  for ( i = 0; i < scanumber; i++ ) {
241  if ( scabuffer[i].buffer ) M_free(scabuffer[i].buffer,"scabuf2");
242  scabuffer[i].buffer = scabuffer[i].pointer = scabuffer[i].top = 0;
243  scabuffer[i].bufsize = scabuffer[i].numterms = 0;
244  }
245  if ( scabuffer ) M_free(scabuffer,"scabuffer");
246  scabuffer = 0;
247  scasize = 0; scanumber = -1;
248  }
249  if ( objbuffer ) M_free(objbuffer,"objbuffer");
250  objbuffer = objtop = 0;
251  objsize = 0;
252  if ( objoffset ) M_free(objoffset,"objoffset");
253  objoffset = 0;
254  numobjoffsets = 0;
255  if ( sortpointers ) M_free(sortpointers,"optisort");
256  sortpointers = 0; helppointers = 0; sizepointers = 0;
257  if ( multiplicities ) M_free(multiplicities,"multiplicities");
258  multiplicities = 0;
259  multiplicitysize = 0;
260 }
261 
262 /*
263  #] CleanOptiBuffer:
264  #[ PutObject:
265 */
266 
267 int PutObject(WORD *object, int type)
268 {
269  int i, k;
270  WORD *obj, *o, *oo, *t, *newobjbuffer;
271  LONG size = 2, j, newobjsize, *newoffsets;
272  switch ( type ) {
273  case SYMBOL: size = 2; break;
274  case DOTPRODUCT: size = 3; break;
275  case VECTOR: size = 3; break;
276  case SUBEXPRESSION:
277  t = object; while ( *t ) t += *t; size = t-object+2; break;
278  case FUNCTION: size = object[1]+1; break;
279  }
280  for ( i = 1; i <= numobjects; i++ ) {
281  obj = objbuffer + objoffset[i];
282  if ( *obj != type ) continue;
283  switch ( type ) {
284  case SYMBOL: if ( obj[1] == *object ) return(i); break;
285  case DOTPRODUCT:
286  case VECTOR:
287  if ( obj[1] == *object && obj[2] == object[1] ) return(i);
288  break;
289  case FUNCTION:
290  k = object[1]; o = obj+1; oo = object;
291  while ( --k >= 0 ) {
292  if ( *o != *oo ) break;
293  o++; oo++;
294  }
295  if ( k < 0 ) return(i);
296  break;
297  case SUBEXPRESSION:
298  if ( size == objoffset[i+1]-objoffset[i] ) {
299  o = obj + 1; oo = object; j = size-1;
300  while ( --j >= 0 ) {
301  if ( *o != *oo ) break;
302  o++; oo++;
303  }
304  if ( j < 0 ) return(i);
305  }
306  break;
307  }
308  }
309 /*
310  Here we should add an object. We have to check two buffers.
311 */
312  while ( objpointer + size >= objtop ) {
313  if ( objsize == 0 ) newobjsize = 2*size;
314  else newobjsize = 2*objsize;
315  if ( newobjsize < 200 ) newobjsize = 200;
316  newobjbuffer = (WORD *)Malloc1(newobjsize*sizeof(WORD),"objbuffer");
317  if ( objsize ) {
318  for ( j = 0; j < objsize; j++ ) newobjbuffer[j] = objbuffer[j];
319  objpointer = (objpointer-objbuffer) + newobjbuffer;
320  if ( objbuffer ) M_free(objbuffer,"objbuffer");
321  objbuffer = newobjbuffer;
322  }
323  else {
324  objpointer = objbuffer = newobjbuffer;
325  }
326  objsize = newobjsize;
327  objtop = objbuffer + objsize;
328  }
329  if ( numobjects + 3 >= numobjoffsets ) {
330  if ( numobjoffsets == 0 ) newobjsize = 50;
331  else newobjsize = 2*numobjoffsets;
332  newoffsets = (LONG *)Malloc1(newobjsize*sizeof(LONG),"newoffsets");
333  if ( numobjects > 0 ) {
334  for ( j = numobjects+1; j >= 0; j-- ) newoffsets[j] = objoffset[j];
335  if ( objoffset ) M_free(objoffset,"objoffset");
336  }
337  else {
338  newoffsets[0] = newoffsets[1] = newoffsets[2] = 0;
339  }
340  numobjoffsets = newobjsize; objoffset = newoffsets;
341  }
342  numobjects++;
343  o = objbuffer + objoffset[numobjects]; *o++ = type;
344  switch ( type ) {
345  case SYMBOL: *o++ = *object; break;
346  case VECTOR:
347  case DOTPRODUCT: *o++ = *object++; *o++ = *object; break;
348  case FUNCTION: k = object[1]; oo = object;
349  while ( --k >= 0 ) *o++ = *oo++;
350  break;
351  case SUBEXPRESSION:
352  oo = object; j = size;
353  while ( --j >= 0 ) *o++ = *oo++;
354  break;
355  }
356  objoffset[numobjects+1] = o - objbuffer;
357  return(numobjects);
358 }
359 
360 /*
361  #] PutObject:
362  #[ AddToOpti:
363 */
364 
365 int AddToOpti(WORD *term, int num)
366 {
367  LONG newnumber, i;
368  WORD *w;
369  int j;
370  SCALAR *newsca;
371  if ( num >= scasize ) { /* increase scalar buffer */
372  if ( scasize <= 0 ) newnumber = 100;
373  else newnumber = scasize * 2;
374  newsca = (SCALAR *)Malloc1(newnumber * sizeof(SCALAR),"scabuffer");
375  if ( scanumber > 0 ) {
376  for ( i = 0; i < scanumber; i++ ) newsca[i] = scabuffer[i];
377  for ( ; i < newnumber; i++ ) {
378  newsca[i].buffer = newsca[i].pointer = newsca[i].top = 0;
379  newsca[i].numterms = newsca[i].bufsize = 0;
380  }
381  if ( scabuffer ) M_free(scabuffer,"scabuffer");
382  }
383  scabuffer = newsca;
384  scasize = newnumber;
385  }
386  if ( num >= scanumber ) { /* create new entry */
387  for ( i = scanumber+1; i <= num; i++ ) {
388  newsca = scabuffer + i;
389  if ( newsca->bufsize == 0 ) {
390  newsca->bufsize = 40;
391  newsca->buffer = (WORD *)Malloc1(newsca->bufsize*sizeof(WORD),"scabuf2");
392  newsca->pointer = newsca->buffer;
393  newsca->top = newsca->buffer + newsca->bufsize;
394  newsca->numterms = 0;
395  }
396  }
397  scanumber = num+1;
398  }
399  newsca = scabuffer + num;
400  while ( newsca->pointer + term[0]+1 >= newsca->top ) {
401  newnumber = newsca->bufsize * 2;
402  w = (WORD *)Malloc1(newnumber*sizeof(WORD),"newscabuffer");
403  i = newsca->pointer - newsca->buffer;
404  while ( --i >= 0 ) w[i] = newsca->buffer[i];
405  newsca->pointer = ( newsca->pointer - newsca->buffer ) + w;
406  newsca->bufsize = newnumber;
407  newsca->top = w + newsca->bufsize;
408  if ( newsca->buffer ) M_free(newsca->buffer,"newscabuffer");
409  newsca->buffer = w;
410  }
411  w = term;
412  j = *term;
413  while ( --j >= 0 ) *newsca->pointer++ = *w++;
414  *(newsca->pointer) = 0;
415  newsca->numterms++;
416  return(0);
417 }
418 
419 /*
420  #] AddToOpti:
421  #[ FindScratchName:
422 */
423 
424 int FindScratchName()
425 {
426  int i, j, k, l, m;
427  WORD number;
428  UBYTE sname[10];
429 /*
430  First attempt. Single characters. Be carefull with the case.
431 */
432  scratchname[0] = 'z'; scratchname[1] = 0;
433  sname[0] = 'Z'; sname[1] = 0;
434  for ( i = 25; i >= 0; i--, scratchname[0]--, sname[0]-- ) {
435  if ( ( GetName(AC.varnames,scratchname,&number,NOAUTO) == NAMENOTFOUND )
436  && ( GetName(AC.exprnames,scratchname,&number,NOAUTO) == NAMENOTFOUND )
437  && ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
438  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
439  return(0);
440  }
441  scratchname[0] = 'z'; scratchname[2] = 0; sname[2] = 0;
442  for ( i = 25; i >= 0; i--, scratchname[0]-- ) {
443  scratchname[1] = 'z';
444  for ( j = 25; j >= 0; j--, scratchname[1]-- ) {
445  sname[0] = scratchname[0]; sname[1] = scratchname[1];
446  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
447  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
448  sname[0] = scratchname[0]; sname[1] = (UBYTE)(scratchname[1]-'a'+'A');
449  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
450  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
451  sname[0] = (UBYTE)(scratchname[0]-'a'+'A'); sname[1] = (UBYTE)(scratchname[1]-'a'+'A');
452  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
453  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
454  sname[0] = (UBYTE)(scratchname[0]-'a'+'A'); sname[1] = scratchname[1];
455  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
456  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
457  return(0);
458  }
459  }
460  }
461  }
462  scratchname[1] = '9';
463  for ( j = 9; j >= 0; j--, scratchname[1]-- ) {
464  sname[0] = scratchname[0]; sname[1] = scratchname[1];
465  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
466  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
467  sname[0] = (UBYTE)(scratchname[0]-'a'+'A');
468  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
469  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
470  return(0);
471  }
472  }
473  }
474  scratchname[0] = 'z'; scratchname[3] = 0; sname[3] = 0;
475  for ( i = 25; i >= 0; i--, scratchname[0]-- ) {
476  scratchname[1] = '9';
477  for ( j = 9; j >= 0; j--, scratchname[1]-- ) {
478  scratchname[2] = '9';
479  for ( k = 9; k >= 0; k--, scratchname[2]-- ) {
480  sname[0] = scratchname[0]; sname[1] = scratchname[1];
481  sname[2] = scratchname[2];
482  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
483  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
484  sname[0] = (UBYTE)(scratchname[0]-'a'+'A');
485  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
486  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
487  return(0);
488  }
489  }
490  }
491  }
492  scratchname[0] = 'z'; scratchname[4] = 0; sname[4] = 0;
493  for ( i = 25; i >= 0; i--, scratchname[0]-- ) {
494  scratchname[1] = '9';
495  for ( j = 9; j >= 0; j--, scratchname[1]-- ) {
496  scratchname[2] = '9';
497  for ( k = 9; k >= 0; k--, scratchname[2]-- ) {
498  scratchname[3] = '9';
499  for ( l = 9; l >= 0; l--, scratchname[3]-- ) {
500  sname[0] = scratchname[0]; sname[1] = scratchname[1];
501  sname[2] = scratchname[2]; sname[3] = scratchname[3];
502  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
503  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
504  sname[0] = (UBYTE)(scratchname[0]-'a'+'A');
505  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
506  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
507  return(0);
508  }
509  }
510  }
511  }
512  }
513  scratchname[0] = 'z'; scratchname[5] = 0; sname[5] = 0;
514  for ( i = 25; i >= 0; i--, scratchname[0]-- ) {
515  scratchname[1] = '9';
516  for ( j = 9; j >= 0; j--, scratchname[1]-- ) {
517  scratchname[2] = '9';
518  for ( k = 9; k >= 0; k--, scratchname[2]-- ) {
519  scratchname[3] = '9';
520  for ( l = 9; l >= 0; l--, scratchname[3]-- ) {
521  scratchname[4] = '9';
522  for ( m = 9; m >= 0; m--, scratchname[4]-- ) {
523  sname[0] = scratchname[0]; sname[1] = scratchname[1];
524  sname[2] = scratchname[2]; sname[3] = scratchname[3];
525  sname[4] = scratchname[4];
526  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
527  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) ) {
528  sname[0] = (UBYTE)(scratchname[0]-'a'+'A');
529  if ( ( GetName(AC.varnames,sname,&number,NOAUTO) == NAMENOTFOUND )
530  && ( GetName(AC.exprnames,sname,&number,NOAUTO) == NAMENOTFOUND ) )
531  return(0);
532  }
533  }
534  }
535  }
536  }
537  }
538  MLOCK(ErrorMessageLock);
539  MesPrint("Could not find a decent name for the scratch variable in Optimize");
540  MUNLOCK(ErrorMessageLock);
541  return(-1);
542 }
543 
544 /*
545  #] FindScratchName:
546  #[ PrintOptima:
547 */
548 
549 int PrintOptima(WORD numexpr)
550 {
551  UBYTE obuffer[80], lbuf[24];
552  WORD *obj, stermbuf[10], *t, *m, n, oldskip = AO.OutSkip;
553  int i, first, fsym, *used;
554  SCALAR *sca;
555  LONG num, totnum = numobjects + scanumber, j;
556  AC.OutputMode = FORTRANMODE;
557  used = (int *)Malloc1((totnum+1)*sizeof(int),"PrintOptima");
558  for ( j = 0; j < totnum; j++ ) used[j] = 0;
559 /*
560  First we print the objects for as far as they should be printed
561  (we do not convert the symbols).
562  Because these routines are called from the writing routines we can
563  safely assume that the writing environment is active.
564 */
565  FiniLine();
566  NumToStr(lbuf,numobjects+scanumber-1);
567  sprintf((char *)obuffer,"DOUBLE PRECISION %s(%s)",(char *)scratchname,lbuf);
568  TokenToLine(obuffer);
569  FiniLine();
570  for ( i = 1; i <= iniobjects; i++ ) {
571  obj = objbuffer + objoffset[i];
572  switch ( *obj ) {
573  case SYMBOL:
574  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i);
575  TokenToLine(obuffer);
576  stermbuf[0] = SYMBOL;
577  stermbuf[1] = 4;
578  stermbuf[2] = obj[1];
579  stermbuf[3] = 1;
580  WriteSubTerm(stermbuf,1);
581  break;
582  case DOTPRODUCT:
583  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i);
584  TokenToLine(obuffer);
585  stermbuf[0] = DOTPRODUCT;
586  stermbuf[1] = 5;
587  stermbuf[2] = obj[1];
588  stermbuf[3] = obj[2];
589  stermbuf[4] = 1;
590  WriteSubTerm(stermbuf,1);
591  break;
592  case VECTOR:
593  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i);
594  TokenToLine(obuffer);
595  stermbuf[0] = VECTOR;
596  stermbuf[1] = 4;
597  stermbuf[2] = obj[1];
598  stermbuf[3] = obj[2];
599  WriteSubTerm(stermbuf,1);
600  break;
601 /*
602  case LNUMBER:
603  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i);
604  TokenToLine(obuffer);
605  n = REDLENG(obj[2]);
606  if ( n < 0 ) { n = -n; TokenToLine((UBYTE *)" - "); }
607  RatToLine((UWORD *)(obj+3),n);
608  break;
609 */
610  case SUBEXPRESSION:
611  break;
612  case FUNCTION:
613  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i);
614  TokenToLine(obuffer);
615  AO.OutSkip = 7;
616  WriteSubTerm(obj+1,1);
617  AO.OutSkip = oldskip;
618  break;
619  }
620  FiniLine();
621  used[i] = 1;
622  }
623  for (;;) {
624  for ( i = scanumber-1; i >= 0; i-- ) {
625  if ( used[i+numobjects] && i != 0 ) continue;
626  sca = scabuffer+i;
627  t = sca->buffer;
628  while ( t < sca->pointer ) {
629  m = t + t[2] + 1; t += *t;
630  while ( m < t ) {
631  num = ( (LONG)(m[0]) << BITSINWORD ) + m[1];
632  if ( used[num] == 0 ) goto nexti;
633  m += 3;
634  }
635  }
636 
637  if ( i == 0 ) {
638  sprintf((char *)obuffer,"%s=",(char *)(EXPRNAME(numexpr)));
639  TokenToLine(obuffer);
640  }
641  else {
642  sprintf((char *)obuffer,"%s(%d)=",(char *)scratchname,i+numobjects);
643  TokenToLine(obuffer);
644  used[i+numobjects] = 1;
645  }
646  sca = scabuffer+i;
647  t = sca->buffer; first = 1;
648  AO.OutSkip = 7;
649  while ( t < sca->pointer ) {
650  m = t + 1; t += *t;
651  n = REDLENG(m[2]);
652  fsym = 1;
653  if ( n < 0 ) { n = -n; TokenToLine((UBYTE *)"-"); }
654  else if ( !first ) TokenToLine((UBYTE *)"+");
655  first = 0;
656  if ( n != 1 || m[3] != 1 || m[4] != 1 || t == m + m[1] ) {
657  fsym = 0; RatToLine((UWORD *)(m+3),n); }
658  m += m[1];
659  while ( m < t ) {
660  if ( fsym == 0 ) {
661  if ( m[2] < 0 ) TokenToLine((UBYTE *)"/");
662  else TokenToLine((UBYTE *)"*");
663  }
664  else if ( m[2] < 0 ) TokenToLine((UBYTE *)"1/");
665  fsym = 0;
666  num = ( (LONG)(m[0]) << BITSINWORD ) + m[1];
667  NumToStr(lbuf,num);
668  sprintf((char *)obuffer,"%s(%s)",(char *)scratchname,lbuf);
669  TokenToLine(obuffer);
670  if ( m[2] > 1 ) {
671  sprintf((char *)obuffer,"**%d",m[2]);
672  TokenToLine(obuffer);
673  }
674  else if ( m[2] < -1 ) {
675  sprintf((char *)obuffer,"**%d",-m[2]);
676  TokenToLine(obuffer);
677  }
678  m += 3;
679  }
680  }
681  AO.OutSkip = oldskip;
682  FiniLine();
683 nexti:;
684  }
685  for ( j = 1; j < totnum; j++ ) {
686  if ( used[j] == 0 ) break;
687  }
688  if ( j >= totnum ) break;
689  }
690  AC.OutputMode = VORTRANMODE;
691  AO.OutSkip = oldskip;
692  M_free(used,"PrintOptima");
693  return(0);
694 }
695 
696 /*
697  #] PrintOptima:
698  #[ MaxPowerOpti:
699 */
700 
701 WORD MaxPowerOpti(LONG number)
702 {
703  SCALAR *sca = scabuffer + number;
704  WORD pow = 0, *t, *m;
705  t = sca->buffer;
706  while ( t < sca->pointer ) {
707  m = t + t[2] + 1;
708  t += *t;
709  while ( m < t ) {
710  if ( m[2] > pow ) pow = m[2];
711  else if ( -m[2] > pow ) pow = -m[2];
712  m += 3;
713  }
714  }
715  return(pow);
716 }
717 
718 /*
719  #] MaxPowerOpti:
720  #[ HuntNumFactor:
721 */
722 
723 WORD HuntNumFactor(LONG number, WORD *coef, int par)
724 {
725  GETIDENTITY
726  SCALAR *sca = scabuffer + number, *ss;
727  WORD *t, *tt, *ttt, *m, *mm, *coef2, ncoef, n, nn, n1, n2, nt;
728  int i;
729  LONG numnewsca, a;
730  if ( par != 1 ) { coef = AT.WorkPointer + 4; }
731  t = sca->buffer;
732  m = t + 4; i = t[2]; mm = coef; ncoef = t[3];
733  NCOPY(mm,m,i);
734  if ( ncoef < 0 ) ncoef = -ncoef;
735  t += *t;
736  while ( t < sca->pointer ) {
737  if ( AccumGCD(BHEAD (UWORD *)coef,&ncoef,(UWORD *)(t+4),t[3]) ) goto ExitHunt;
738  if ( ncoef == 3 && coef[0] == 1 && coef[1] == 1 ) return(0);
739  t += *t;
740  }
741 /*
742  We have now the gcd in coef. First test for triviality:
743 */
744  if ( ncoef < 0 ) ncoef = -ncoef;
745  if ( ncoef == 3 && coef[0] == 1 && coef[1] == 1 ) return(0);
746 /*
747  Next we divide out the gcd. First flip it so we can use MulRat.
748 */
749  n = (ncoef-1)/2; t = coef; m = coef+n;
750  for ( i = 0; i < n; i++ ) { nn = *t; *t++ = *m; *m++ = nn; }
751  t = m = sca->buffer;
752  coef2 = coef + ncoef;
753  while ( t < sca->pointer ) {
754  nt = *t; ttt = t + t[2] + 1;
755  n1 = REDLENG(t[3]);
756  mm = m+1; *mm++ = LNUMBER;
757  if ( MulRat(BHEAD (UWORD *)(t+4),n1,(UWORD *)coef,n,(UWORD *)coef2,&n2) ) goto ExitHunt;
758  i = INCLENG(n2);
759  *mm++ = ABS(i)+2;
760  *mm++ = i; if ( i < 0 ) i = -i; i--;
761  tt = coef2; while ( --i >= 0 ) *mm++ = *tt++;
762  t += nt;
763  while ( ttt < t ) *mm++ = *ttt++;
764  *m = mm-m; m = mm;
765  }
766  *m = 0;
767  sca->pointer = m;
768 /*
769  Now we make a flip. We put the buffer of this scalar in a new scalar.
770  In the buffer of the new scalar we put the factor and the pointer to
771  the now reduced term.
772 */
773  n = (ncoef-1)/2; t = coef; m = coef+n;
774  for ( i = 0; i < n; i++ ) { nn = *t; *t++ = *m; *m++ = nn; }
775  if ( par == 1 ) return(ncoef);
776 
777  numnewsca = scanumber;
778  AT.WorkPointer[1] = LNUMBER;
779  AT.WorkPointer[2] = ncoef+2;
780  AT.WorkPointer[3] = ncoef;
781  m = AT.WorkPointer + ncoef + 3;
782  *m++ = (numnewsca+numobjects) >> BITSINWORD;
783  *m++ = (numnewsca+numobjects) & WORDMASK;
784  *m++ = 1;
785  AT.WorkPointer[0] = m - AT.WorkPointer;
786  AddToOpti(AT.WorkPointer,numnewsca);
787 /*
788  Note: AddToOpti can move the scabuffer. Hence we need to redefine sca
789 */
790  ss = scabuffer + numnewsca;
791  sca = scabuffer + number;
792  m = sca->buffer; sca->buffer = ss->buffer; ss->buffer = m;
793  m = sca->pointer; sca->pointer = ss->pointer; ss->pointer = m;
794  m = sca->top; sca->top = ss->top; ss->top = m;
795  a = sca->numterms; sca->numterms = ss->numterms; ss->numterms = a;
796  a = sca->bufsize; sca->bufsize = ss->bufsize; ss->bufsize = a;
797 
798  return(ncoef);
799 ExitHunt:
800  MLOCK(ErrorMessageLock);
801  MesCall("HuntNumFactor");
802  MUNLOCK(ErrorMessageLock);
803  Terminate(-1);
804  return(0);
805 }
806 
807 /*
808  #] HuntNumFactor:
809  #[ HuntFactor:
810 
811  Hunts for factors in scalar 'number'
812  The action depends of the value of 'par'
813  par = 0 Symbolic and numeric factors. Extra sca is made.
814  par = 1 Symbolic and numeric factors. Factor is divided out.
815  par = 2 Symbolic and numeric factors. Factor is only determined.
816  par = -1 Symbolic factors only. Factor is divided out.
817  par = -2 Symbolic factors only. Factor is only determined.
818  The factor is returned as a regular term in 'factor'
819  In the case of par = 0, one can put 'factor' to zero.
820  If an extra sca is made, the original sca and the new one are exchanged.
821  The return value is the size of the factor, or zero if the factor is 1.
822 */
823 
824 WORD HuntFactor(LONG number, WORD *factor, int par)
825 {
826  GETIDENTITY
827  SCALAR *sca, *scb;
828  WORD *t, *m, *ft, *fm, *fr, *frr, *fact, *coef, ncoef, retval;
829  int i, size;
830  LONG newnum, a;
831  if ( factor == 0 && par != 0 ) {
832  MLOCK(ErrorMessageLock);
833  MesPrint("Internal error: wrong value (%d) of par in HuntFactor",
834  (WORD)par);
835  MUNLOCK(ErrorMessageLock);
836  Terminate(-1);
837  }
838  if ( par >= 0 ) { /* We need the coefficient also */
839  if ( factor == 0 ) factor = AT.WorkPointer;
840  coef = factor + 4;
841  ncoef = HuntNumFactor(number,coef,1);
842  if ( ncoef == 0 || ( ncoef == 3 && coef[0] == 1
843  && coef[1] == 1 ) ) {
844  retval = 0; ncoef = 3; coef[0] = coef[1] = 1; }
845  else retval = ABS(ncoef)+3;
846  }
847  else {
848  coef = factor + 4; ncoef = 3; coef[0] = coef[1] = 1;
849  retval = 6;
850  }
851  factor[1] = LNUMBER; factor[2] = ABS(ncoef)+2; factor[3] = ncoef;
852  factor[0] = factor[2]+1;
853  fact = factor + factor[0];
854  sca = scabuffer+number;
855  t = sca->buffer;
856  m = t + t[2] + 1; t += *t;
857  fm = fact;
858  while ( m < t ) *fm++ = *m++;
859  ft = fm; size = ft - fact;
860  if ( size == 0 ) { /* no extra factor */
861  return(retval);
862  }
863  while ( t < sca->pointer ) {
864  m = t + t[2] + 1; t += *t;
865  fm = fact;
866  while ( fm < ft && m < t ) {
867  if ( fm[0] == m[0] && fm[1] == m[1] ) {
868  if ( fm[2] == m[2] ) { m += 3; fm += 3; }
869  else if ( ( fm[2] > 0 && m[2] < 0 )
870  || ( fm[2] < 0 && m[2] > 0 ) ) goto loosethis;
871  else if ( ( m[2] > 0 && fm[2] > m[2] )
872  || ( m[2] < 0 && fm[2] < m[2] ) ) {
873  fm[2] = m[2]; m += 3; fm += 3;
874  }
875  else { m += 3; fm += 3; }
876  }
877  else if ( fm[0] < m[0] || ( fm[0] == m[0] && fm[1] < m[1] ) ) {
878 loosethis: fr = fm + 3; frr = fm;
879  while ( fr < ft ) *frr++ = *fr++;
880  ft = frr; size -= 3;
881  if ( size == 0 ) { /* no symbolic factor left */
882  return(retval);
883  }
884  m += 3;
885  }
886  else m += 3;
887  }
888  if ( fm < ft ) {
889  ft = fm; size = ft - fact;
890  if ( size == 0 ) { /* no symbolic factor left */
891  return(retval);
892  }
893  }
894  }
895 /*
896  We have a factor now. We can divide it out if the parameter asks for it.
897 */
898  factor[0] += size;
899  size = factor[0];
900  if ( par > 1 || par < -1 ) return(size);
901  fr = t = sca->buffer;
902  while ( t < sca->pointer ) {
903  frr = fr; i = t[2]+1; m = t; t += *t;
904  while ( --i >= 0 ) *fr++ = *m++; /* coefficient */
905  fm = fact;
906  while ( fm < ft && m < t ) {
907  if ( fm[0] == m[0] && fm[1] == m[1] ) {
908  if ( fm[2] == m[2] ) { m += 3; fm += 3; }
909  else {
910  *fr++ = *m++; *fr++ = *m++; *fr++ = *m++ - fm[2];
911  fm += 3;
912  }
913  }
914  else { *fr++ = *m++; *fr++ = *m++; *fr++ = *m++; }
915  }
916  while ( m < t ) *fr++ = *m++;
917  *frr = fr - frr;
918  NormOpti(frr);
919  }
920  sca->pointer = fr;
921  if ( par == 1 || par == -1 ) return(size);
922 /*
923  Now we create a term with the factor and put it as a scalar.
924  We have to exchange the order of the scalars!
925 */
926  t = AT.WorkPointer;
927  *t++ = size + 9; *t++ = LNUMBER; *t++ = 5; *t++ = 3; *t++ = 1; *t++ = 1;
928  t += size;
929  *t++ = ( scanumber + numobjects ) >> BITSINWORD;
930  *t++ = ( scanumber + numobjects ) & WORDMASK;
931  *t++ = 1;
932  newnum = scanumber;
933  AddToOpti(AT.WorkPointer,newnum);
934  sca = scabuffer + number;
935  scb = scabuffer + newnum;
936  t = sca->buffer; sca->buffer = scb->buffer; scb->buffer = t;
937  t = sca->pointer; sca->pointer = scb->pointer; scb->pointer = t;
938  t = sca->top; sca->top = scb->top; scb->top = t;
939  a = sca->numterms; sca->numterms = scb->numterms; scb->numterms = a;
940  a = sca->bufsize; sca->bufsize = scb->bufsize; scb->bufsize = a;
941  return(size);
942 }
943 
944 /*
945  #] HuntFactor:
946  #[ HuntPairs:
947 
948  Routine looks for an object (x) to a power of at least 'power'.
949  It looks then for terms with are similar except for x^power.
950  It will also try to look for (x+y)^power.
951  This is more complicated than it looks because y can involve the
952  coefficient of the term.
953 */
954 
955 void HuntPairs(LONG number, WORD power)
956 {
957  GETIDENTITY
958  SCALAR *sca = scabuffer + number;
959  WORD *t, *tt, *m, *mm, *w, *w1, *pattern, *p, *pp, *patstop,
960  *coef, ncoef, nf, nc2, nc, *newter;
961  int patsize, i, first, pushback = 0, nons, action = 0;
962  LONG numnewsca = 0, ns;
963  pattern = AT.WorkPointer;
964  t = sca->buffer;
965  while ( t < sca->pointer ) {
966  if ( *t < 0 ) { t -= *t; continue; }
967  w1 = t; m = t + 1 + t[2]; t += *t;
968  while ( m < t ) {
969  if ( m[2] >= power || -m[2] >= power ) {
970  w = w1 + 1 + w1[2];
971  p = pattern;
972  while ( w < t ) {
973  if ( w == m ) {
974  if ( power == m[2] ) { w += 3; }
975  else if ( power == -m[2] ) { w += 3; }
976  else { *p++ = *w++; *p++ = *w++;
977  if ( m[2] < 0 ) *p++ = *w++ + power;
978  else *p++ = *w++ - power;
979  }
980  }
981  else { *p++ = *w++; *p++ = *w++; *p++ = *w++; }
982  }
983  patsize = p - pattern;
984  if ( patsize == 0 ) { pushback++; goto nextm; }
985  AT.WorkPointer = patstop = p;
986  first = 1;
987 /*
988  Now collect all terms which contain this pattern and put
989  them into a new subexpression. Do this however only
990  once there are at least two terms!
991 */
992  tt = sca->buffer;
993  while ( tt < sca->pointer ) {
994  if ( *tt < 0 ) { tt -= *tt; continue; }
995  w = tt; mm = tt + tt[2] + 1; tt += *tt;
996  if ( w == w1 ) continue;
997  p = pattern; pp = AT.WorkPointer + w[2] + 1;
998  while ( p < patstop && mm < tt ) {
999  if ( mm[0] == p[0] && mm[1] == p[1] ) {
1000  if ( ( p[2] > 0 && mm[2] >= p[2] )
1001  || ( p[2] < 0 && mm[2] <= p[2] ) ) {
1002  if ( mm[2] == p[2] ) { mm += 3; p += 3; }
1003  else {
1004  *pp++ = *mm++; *pp++ = *mm++;
1005  *pp++ = *mm++ - p[2];
1006  p += 3;
1007  }
1008  }
1009  else goto nexttt;
1010  }
1011  else { *pp++ = *mm++; *pp++ = *mm++; *pp++ = *mm++; }
1012  }
1013  if ( p >= patstop ) { /* match! */
1014  while ( mm < tt ) *pp++ = *mm++;
1015  p = AT.WorkPointer; mm = w+1; i = mm[1];
1016  *p++ = pp - AT.WorkPointer;
1017  while ( --i >= 0 ) *p++ = *mm++;
1018  if ( first ) {
1019  p = pp+1; mm = w1+1; i = mm[1];
1020  while ( --i >= 0 ) *p++ = *mm++;
1021  *p++ = m[0]; *p++ = m[1];
1022  if ( m[2] < 0 ) *p++ = -power;
1023  else *p++ = power;
1024  *pp = p - pp;
1025  NormOpti(pp);
1026  numnewsca = scanumber;
1027  AddToOpti(pp,numnewsca);
1028  sca = scabuffer + number;
1029  first = 0;
1030  *w1 = -*w1;
1031  }
1032  NormOpti(AT.WorkPointer);
1033  AddToOpti(AT.WorkPointer,numnewsca);
1034  *w = -*w;
1035  sca->numterms--;
1036  }
1037 nexttt:;
1038  }
1039  if ( *w1 > 0 ) { /* no second term! */
1040  goto nextm;
1041  }
1042  action++;
1043 /*
1044  First simplify the expression
1045 */
1046  nf = HuntFactor(numnewsca,AT.WorkPointer,1);
1047  SortOpti(numnewsca);
1048 /*
1049  Check whether we have this scalar already.
1050  The factor could still be a problem.
1051 */
1052  coef = AT.WorkPointer + 4 + nf;
1053  ns = TestNewSca(numnewsca,coef,&ncoef);
1054  if ( ns != numnewsca ) {
1055  nons = 1;
1056  scanumber--;
1057  }
1058  else nons = 0;
1059 /*
1060  See whether we have to combine things
1061  After this we only look for nf != 0
1062 */
1063  if ( nons ) {
1064  if ( nf ) {
1065  newter = coef + 2*ABS(ncoef) + 2;
1066  ncoef = REDLENG(ncoef);
1067  nc2 = REDLENG(AT.WorkPointer[3]);
1068  if ( MulRat(BHEAD (UWORD *)coef,ncoef,(UWORD *)(AT.WorkPointer+4)
1069  ,nc2,(UWORD *)(newter+4),&nc) ) {
1070  MLOCK(ErrorMessageLock);
1071  MesCall("HuntPairs");
1072  MUNLOCK(ErrorMessageLock);
1073  Terminate(-1);
1074  }
1075  nc = INCLENG(nc);
1076  p = AT.WorkPointer + AT.WorkPointer[2] + 1;
1077  t = AT.WorkPointer + nf;
1078  newter[2] = ABS(nc) + 2;
1079  newter[3] = nc;
1080  newter[1] = LNUMBER;
1081  m = newter + newter[2] + 1;
1082  while ( p < t ) *m++ = *p++;
1083  nf = i = newter[0] = m-newter;
1084  m = newter; t = AT.WorkPointer;
1085  while ( --i >= 0 ) *t++ = *m++;
1086  }
1087  else {
1088  nf = ABS(ncoef);
1089  AT.WorkPointer[0] = nf + 3;
1090  AT.WorkPointer[1] = LNUMBER;
1091  AT.WorkPointer[2] = nf + 2;
1092  AT.WorkPointer[3] = ncoef;
1093  nf += 3;
1094  }
1095  }
1096 /*
1097  Compress the buffer and continue from the next term
1098 */
1099  m = t = sca->buffer; p = AT.WorkPointer + 1;
1100  while ( t < sca->pointer ) {
1101  if ( *t < 0 ) {
1102  if ( t == w1 ) {
1103  w = m;
1104  if ( nf == 0 ) {
1105  *p++ = LNUMBER; *p++ = 5; *p++ = 3; *p++ = 1;
1106  *p++ = 1;
1107  }
1108  else {
1109  p = AT.WorkPointer + nf;
1110  }
1111  mm = pattern;
1112  while ( mm < patstop ) *p++ = *mm++;
1113  *p++ = (ns+numobjects) >> BITSINWORD;
1114  *p++ = (ns+numobjects) & WORDMASK;
1115  *p++ = 1;
1116  *(AT.WorkPointer) = p - AT.WorkPointer;
1117  NormOpti(AT.WorkPointer);
1118  }
1119  t -= *t;
1120  }
1121  else if ( t > m ) {
1122  i = *t; while ( --i >= 0 ) *m++ = *t++;
1123  }
1124  else { t += *t; m += *m; }
1125  }
1126  i = *(AT.WorkPointer); p = AT.WorkPointer;
1127  while ( --i >= 0 ) *m++ = *p++;
1128  sca->pointer = m; *m = 0;
1129  t = w;
1130  goto nextterm;
1131  }
1132 nextm:
1133  m += 3;
1134  }
1135 nextterm:;
1136  }
1137  if ( action ) SortOpti(number);
1138 /*
1139  Next we have to look in the new scalars for things of the type
1140  y^power + power*a*y^(power-1) + ... + a^power
1141  We trigger on three terms.
1142  We do this term by term and try to accumulate powers ie:
1143  (a+b+c)^3 = a^3+b^3+c^3+3*a^2*b+3*a^2*c+3*a*b^2+3*a*c^2+3*b^2*c+3*b*c^2
1144  +6*a*b*c -> (a^3+3*a^2*b +b^3)
1145  -> ab^3+c^3+3*a^2*c+3*a*c^2+3*b^2*c+3*b*c^2+6*a*b*c;
1146  Next try to find more ab combinations (in squares etc) ->
1147  ab^3+c^3+3*ab^2*c+3*ab*c^2
1148  and then we start from the beginning which gives directly
1149  abc^3
1150  If an object simplifies completely, we adjust the power where it is
1151  used so that more tricks will be possible later.
1152 */
1153  AT.WorkPointer = pattern;
1154 }
1155 
1156 /*
1157  #] HuntPairs:
1158  #[ HuntBrackets:
1159 
1160  Routine goes through a sca and tries to find the most popular object.
1161  Then breaks up the expression in the style a*x+b.
1162  It keeps doing this till there are no more objects like this.
1163 */
1164 
1165 void HuntBrackets(LONG number)
1166 {
1167  GETIDENTITY
1168  SCALAR *sca = scabuffer + number;
1169  WORD *t, *m, *tt, *mm, *left = 0, mostpopular[2], *coef, nf, *newter;
1170  WORD ncoef, nc2, nc;
1171  LONG mostmultiple, n, i, num, newscanum, ns;
1172  int j, neg, nons;
1173  for (;;) {
1174  if ( sca->numterms <= 2 ) return;
1175  n = scanumber + numobjects;
1176  while ( 2*n >= multiplicitysize ) {
1177  if ( multiplicitysize == 0 ) multiplicitysize = 500;
1178  else multiplicitysize *= 2;
1179  if ( multiplicities ) M_free(multiplicities,"multiplicities");
1180  multiplicities = (LONG *)Malloc1(multiplicitysize*sizeof(LONG),"multiplicities");
1181  }
1182  for ( i = 0; i <= n; i++ ) {
1183  multiplicities[i] = 0; multiplicities[i+n] = 0;
1184  }
1185  t = sca->buffer;
1186  while ( t < sca->pointer ) {
1187  m = t + t[2] + 1; t += *t;
1188  while ( m < t ) {
1189  num = ( (LONG)(m[0]) << BITSINWORD ) + m[1];
1190  if ( m[2] > 0 ) multiplicities[num]++;
1191  else multiplicities[n+num]++;
1192  m += 3;
1193  }
1194  }
1195  for ( i = 0, num = -1, mostmultiple = 0; i <= n; i++ ) {
1196  if ( multiplicities[i] > mostmultiple ) {
1197  mostmultiple = multiplicities[i]; num = i;
1198  }
1199  if ( multiplicities[i+n] > mostmultiple ) {
1200  mostmultiple = multiplicities[i+n]; num = i+n;
1201  }
1202  }
1203  if ( mostmultiple <= 1 ) return;
1204  if ( num > n ) { neg = 1; num -= n; }
1205  else neg = 0;
1206  mostpopular[0] = num >> BITSINWORD;
1207  mostpopular[1] = num & WORDMASK;
1208 /*
1209  Now all terms containing num will be taken out ( with one power of num less )
1210  and sent off to a new expression.
1211  First find the first to start the system.
1212 */
1213  t = sca->buffer; newscanum = scanumber;
1214  while ( t < sca->pointer ) {
1215  m = t + t[2] + 1; tt = t; t += *t;
1216  while ( m < t ) {
1217  if ( m[0] == mostpopular[0] && m[1] == mostpopular[1] ) {
1218  if ( ( neg && m[2] > 0 ) || ( !neg && m[2] < 0 ) ) break;
1219  mm = AT.WorkPointer;
1220  left = tt;
1221  while ( tt < m ) *mm++ = *tt++;
1222  if ( neg ) {
1223  if ( m[2] < -1 ) {
1224  *mm++ = *tt++; *mm++ = *tt++; *mm++ = 1 + *tt++;
1225  }
1226  else tt += 3;
1227  }
1228  else {
1229  if ( m[2] > 1 ) {
1230  *mm++ = *tt++; *mm++ = *tt++; *mm++ = *tt++ - 1;
1231  }
1232  else tt += 3;
1233  }
1234  while ( tt < t ) *mm++ = *tt++;
1235  AT.WorkPointer[0] = mm - AT.WorkPointer;
1236  sca->numterms--;
1237  AddToOpti(AT.WorkPointer,newscanum);
1238  goto dorest;
1239  }
1240  m += 3;
1241  }
1242  }
1243 dorest:
1244  sca = scabuffer+number;
1245  while ( t < sca->pointer ) {
1246  m = t + t[2] + 1; tt = t; t += *t;
1247  while ( m < t ) {
1248  if ( m[0] == mostpopular[0] && m[1] == mostpopular[1] ) {
1249  if ( ( neg && m[2] > 0 ) || ( !neg && m[2] < 0 ) ) break;
1250  mm = AT.WorkPointer;
1251  while ( tt < m ) *mm++ = *tt++;
1252  if ( neg ) {
1253  if ( m[2] < -1 ) {
1254  *mm++ = *tt++; *mm++ = *tt++; *mm++ = 1 + *tt++;
1255  }
1256  else tt += 3;
1257  }
1258  else {
1259  if ( m[2] > 1 ) {
1260  *mm++ = *tt++; *mm++ = *tt++; *mm++ = *tt++ - 1;
1261  }
1262  else tt += 3;
1263  }
1264  while ( tt < t ) *mm++ = *tt++;
1265  AT.WorkPointer[0] = mm - AT.WorkPointer;
1266  sca->numterms--;
1267  AddToOpti(AT.WorkPointer,newscanum);
1268  goto nextt;
1269  }
1270  m += 3;
1271  }
1272  while ( tt < t ) *left++ = *tt++; /* leftovers */
1273 nextt:;
1274  }
1275 /*
1276  Now the remaining term contains the popular factor, the new sca
1277  and possibly factors from the new sca that we can take out.
1278  First find the factors.
1279 */
1280  nf = HuntFactor(newscanum,AT.WorkPointer,1);
1281  SortOpti(newscanum);
1282 /*
1283  Check whether we have this scalar already.
1284  The factor could still be a problem.
1285 */
1286  coef = AT.WorkPointer + 4 + nf;
1287  ns = TestNewSca(newscanum,coef,&ncoef);
1288  if ( ns != newscanum ) {
1289  nons = 1;
1290  scanumber--;
1291  }
1292  else nons = 0;
1293 /*
1294  See whether we have to combine things
1295  After this we only look for nf != 0
1296 */
1297  if ( nons ) {
1298  if ( nf ) {
1299  newter = coef + 2*ABS(ncoef) + 2;
1300  ncoef = REDLENG(ncoef);
1301  nc2 = REDLENG(AT.WorkPointer[3]);
1302  if ( MulRat(BHEAD (UWORD *)coef,ncoef,(UWORD *)(AT.WorkPointer+4),nc2
1303  ,(UWORD *)(newter+4),&nc) ) {
1304  MLOCK(ErrorMessageLock);
1305  MesCall("HuntPairs");
1306  MUNLOCK(ErrorMessageLock);
1307  Terminate(-1);
1308  }
1309  nc = INCLENG(nc);
1310  tt = AT.WorkPointer + AT.WorkPointer[2] + 1;
1311  t = AT.WorkPointer + nf;
1312  newter[1] = LNUMBER;
1313  newter[2] = ABS(nc) + 2;
1314  newter[3] = nc;
1315  m = newter + newter[2] + 1;
1316  while ( tt < t ) *m++ = *tt++;
1317  nf = j = newter[0] = m-newter;
1318  m = newter; t = AT.WorkPointer;
1319  while ( --j >= 0 ) *t++ = *m++;
1320  }
1321  else {
1322  nf = ABS(ncoef);
1323  AT.WorkPointer[0] = nf + 3;
1324  AT.WorkPointer[1] = LNUMBER;
1325  AT.WorkPointer[2] = nf + 2;
1326  AT.WorkPointer[3] = ncoef;
1327  nf += 3;
1328  }
1329  }
1330  if ( nf == 0 ) {
1331  m = AT.WorkPointer + 1;
1332  *m++ = LNUMBER; *m++ = 5; *m++ = 3; *m++ = 1; *m++ = 1;
1333  *m++ = mostpopular[0];
1334  *m++ = mostpopular[1];
1335  if ( neg ) *m++ = -1;
1336  else *m++ = 1;
1337  }
1338  else {
1339  t = AT.WorkPointer + nf;
1340  m = AT.WorkPointer + AT.WorkPointer[2] + 1;
1341  while ( m < t ) {
1342  if ( m[0] == mostpopular[0] && m[1] == mostpopular[1] ) {
1343  if ( neg ) m[2]--;
1344  else m[2]++;
1345  break;
1346  }
1347  m += 3;
1348  }
1349  if ( m >= t ) {
1350  *m++ = mostpopular[0];
1351  *m++ = mostpopular[1];
1352  if ( neg ) *m++ = -1;
1353  else *m++ = 1;
1354  }
1355  else m = t;
1356  }
1357  *m++ = (ns+numobjects) >> BITSINWORD;
1358  *m++ = (ns+numobjects) & WORDMASK;
1359  *m++ = 1;
1360  *(AT.WorkPointer) = m - AT.WorkPointer;
1361  NormOpti(AT.WorkPointer);
1362  j = *(AT.WorkPointer); m = AT.WorkPointer;
1363  while ( --j >= 0 ) *left++ = *m++;
1364  sca->pointer = left; *left = 0;
1365  sca->numterms++;
1366  SortOpti(number);
1367  }
1368 }
1369 
1370 /*
1371  #] HuntBrackets:
1372  #[ HuntNumBrackets:
1373 
1374  Hunts for terms with an identical coefficient
1375  Makes a bracket of them. in the style of 2*(a+b)+c.
1376 */
1377 
1378 void HuntNumBrackets(LONG number)
1379 {
1380  DUMMYUSE(number);
1381 }
1382 
1383 /*
1384  #] HuntNumBrackets:
1385  #[ HuntPowers:
1386 
1387  Tries to look for something of the type
1388  a^n + n*a^(n-1)*b + b^n
1389  This pattern will trigger the substitution
1390  a^n + n*a^(n-1)*b + b^n = (a+b)^n - rest
1391  a+b will become a new sca
1392 */
1393 
1394 void HuntPowers(LONG number, WORD power)
1395 {
1396  GETIDENTITY
1397  SCALAR *sca = scabuffer + number;
1398  WORD *t1, *m1, *r1, *t2, *m2, *r2, *t3, *m3, *r3, *quotient, *extra
1399  , *q, n1, n2, n3, nq;
1400  int i;
1401  quotient = AT.WorkPointer;
1402  t1 = sca->buffer;
1403  while ( t1 < sca->pointer ) {
1404  m1 = t1 + t1[2] + 1; r1 = t1; t1 += *t1;
1405  while ( m1 < t1 ) {
1406  if ( m1[2] >= power ) {
1407 /*
1408  Now we have to find a second term with the property
1409  that term1/term2 is a clean power (like a^n/b^n)
1410 */
1411  n1 = REDLENG(r1[3]); t2 = t1; m1 = r1 + r1[2] + 1;
1412  while ( t2 < sca->pointer ) {
1413  m2 = t2 + t2[2] + 1; r2 = t2; t2 += *t2;
1414  n2 = REDLENG(r2[3]);
1415  if ( DivRat(BHEAD (UWORD *)(r2+4),n2,(UWORD *)(r1+4),n1,
1416  (UWORD *)(quotient+4),&nq) ) goto callHP;
1417  if ( TakeRatRoot((UWORD *)(quotient+4),&nq,power) ) continue;
1418  n2 = INCLENG(nq);
1419  quotient[2] = ABS(n2)+2; quotient[3] = n2;
1420  quotient[1] = LNUMBER; q = quotient + quotient[2] + 1;
1421  while ( m2 < t2 && m1 < t1 ) {
1422  if ( *m1 == *m2 && m1[1] == m2[1] ) {
1423  if ( ((m2[2]-m1[2])%power) != 0 ) goto nextt2;
1424  *q++ = *m1++; *q++ = *m1++;
1425  *q++ = (m2[2]-*m1++)/power;
1426  m2 += 3;
1427  }
1428  else if ( *m1 < *m2 || ( *m1 == *m2 && m1[1] < m2[1] ) ) {
1429  if ( ( m1[2] % power ) != 0 ) goto nextt2;
1430  *q++ = *m1++; *q++ = *m1++; *q++ = -*m1++ / power;
1431  }
1432  else {
1433  if ( ( m2[2] % power ) != 0 ) goto nextt2;
1434  *q++ = *m2++; *q++ = *m2++; *q++ = *m2++ / power;
1435  }
1436  }
1437  while ( m1 < t1 ) {
1438  if ( ( m1[2] % power ) != 0 ) goto nextt2;
1439  *q++ = *m1++; *q++ = *m1++; *q++ = -*m1++ / power;
1440  }
1441  while ( m2 < t2 ) {
1442  if ( ( m2[2] % power ) != 0 ) goto nextt2;
1443  *q++ = *m2++; *q++ = *m2++; *q++ = *m2++ / power;
1444  }
1445  quotient[0] = q - quotient;
1446 /*
1447  Now we have the a^power and the b^power
1448  +/- a/b is in quotient
1449  if the coefficient is negative it has to be negative.
1450  if it is positive and power is odd it has to be positive.
1451  otherwise it can be either positive or negative.
1452  Construct now the third term.
1453  It is number 1 divided by quotient times the power.
1454 */
1455  extra = quotient + quotient[0];
1456  if ( DivRat(BHEAD (UWORD *)(r1+4),n1,(UWORD *)(quotient+4),nq
1457  ,(UWORD *)(extra+4),&n3) ) goto callHP;
1458  if ( Mully(BHEAD (UWORD *)(extra+4),&n3,(UWORD *)(&power),1) ) goto callHP;
1459  n2 = INCLENG(n3);
1460  extra[1] = LNUMBER; extra[2] = ABS(n2)+2; extra[3] = n2;
1461  m3 = extra + extra[2] + 1;
1462  m1 = r1 + r1[2] + 1;
1463  m2 = quotient + quotient[2] + 1;
1464  while ( m1 < t1 && m2 < q ) {
1465  if ( m1[0] == m2[0] && m1[1] == m2[1] ) {
1466  if ( m1[2] == m2[2] ) { m1+=3; }
1467  else {
1468  *m3++ = *m1++; *m3++ = *m1++; *m3++ = *m1++ - m2[2];
1469  }
1470  m2 += 3;
1471  }
1472  else if ( m1[0] < m2[0] || ( m1[0] == m2[0]
1473  && m1[1] < m2[1] ) ) {
1474  *m3++ = *m1++; *m3++ = *m1++; *m3++ = *m1++;
1475  }
1476  else {
1477  *m3++ = *m2++; *m3++ = *m2++; *m3++ = *m2++;
1478  }
1479  }
1480  while ( m1 < t1 ) {
1481  *m3++ = *m1++; *m3++ = *m1++; *m3++ = *m1++;
1482  }
1483  while ( m2 < q ) {
1484  *m3++ = *m2++; *m3++ = *m2++; *m3++ = *m2++;
1485  }
1486  extra[0] = m3 - extra;
1487 /*
1488  Now look for +/- extra
1489 */
1490  t3 = sca->buffer;
1491  while ( t3 < sca->pointer ) {
1492  r3 = t3; t3 += *t3;
1493  if ( *r3 != *extra ) continue;
1494  if ( r3[2] != extra[2] ) continue;
1495  for ( i = 4; i < *r3; i++ ) {
1496  if ( extra[i] != r3[i] ) break;
1497  }
1498  if ( i < *r3 ) continue;
1499 /*
1500  Now we can finally clear up the sign.
1501 */
1502  if ( r3[3] < 0 ) {
1503  if ( ( power & 1 ) != 0 && quotient[3] > 0 ) goto nextt2;
1504  if ( quotient[3] > 0 ) quotient[3] = -quotient[3];
1505  }
1506  else if ( quotient[3] < 0 ) goto nextt2;
1507 /*
1508  We have everything now
1509 */
1510 
1511  }
1512 nextt2:;
1513  }
1514  }
1515  m1 += 3;
1516  }
1517  }
1518  return;
1519 callHP:
1520  MLOCK(ErrorMessageLock);
1521  MesCall("HuntPowers");
1522  MUNLOCK(ErrorMessageLock);
1523  Terminate(-1);
1524  return;
1525 }
1526 
1527 /*
1528  #] HuntPowers:
1529  #[ CombiOpti:
1530 
1531  We try to express the bigger sca's into smaller sca's
1532 */
1533 
1534 void CombiOpti()
1535 {
1536  LONG i1, i2;
1537  SCALAR *sca1, *sca2;
1538  WORD *t1, *m1, *tt1, *t2, *m2, *tt2, *fill;
1539  int sign = 1, j;
1540  for ( i1 = 0; i1 < scanumber; i1++ ) {
1541  sca1 = scabuffer + i1;
1542  for ( i2 = 0; i2 < scanumber; i2++ ) {
1543  sca2 = scabuffer + i2;
1544  if ( sca2->numterms <= sca1->numterms ) continue;
1545  t1 = sca1->buffer; t2 = sca2->buffer;
1546  tt1 = t1; t1 += *t1;
1547 /*
1548  First find the first term and determine the factor
1549 */
1550  while ( t2 < sca2->pointer ) {
1551  if ( *tt1 != *t2 ) { t2 += *t2; continue; }
1552  m1 = tt1 + tt1[2] + 1;
1553  m2 = t2 + t2[2] + 1; tt2 = t2; t2 += *t2;
1554  while ( m1 < t1 && m2 < t2 ) {
1555  if ( *m1 != *m2 || m1[1] != m2[1] || m1[2] != m2[2] ) break;
1556  m1 += 3; m2 += 3;
1557  }
1558  if ( m1 >= t1 && m2 >= t2 ) { /* we do have a match */
1559 /*
1560  For now we tolerate just a sign
1561 */
1562  if ( tt1[2] != tt2[2] ) goto nexti2;
1563  if ( tt1[3] > 0 ) {
1564  if ( tt2[3] > 0 ) sign = 1;
1565  else sign = -1;
1566  }
1567  else {
1568  if ( tt2[3] < 0 ) sign = 1;
1569  else sign = -1;
1570  }
1571  m1 = tt1 + tt1[2] + 1;
1572  tt2 += 4; tt1 += 4;
1573  while ( tt1 < m1 ) {
1574  if ( *tt1 != *tt2 ) goto nexti2;
1575  tt1++; tt2++;
1576  }
1577  break;
1578  }
1579  }
1580  if ( t2 >= sca2->pointer ) goto nexti2;
1581  while ( t1 < sca1->pointer && t2 < sca2->pointer ) {
1582  if ( *t1 != *t2 ) { t2 += *t2; continue; }
1583  tt1 = t1 + *t1; m1 = t1 + t1[2] + 1;
1584  tt2 = t2 + *t2; m2 = t2 + t2[2] + 1;
1585  while ( m1 < tt1 && m2 < tt2 ) {
1586  if ( *m1 != *m2 || m1[1] != m2[1] || m1[2] != m2[2] ) break;
1587  m1 += 3; m2 += 3;
1588  }
1589  if ( m1 >= tt1 && m2 >= tt2 ) {
1590  if ( t1[2] != t2[2] ) goto nexti2;
1591  if ( sign*t1[3] != t2[3] ) goto nexti2;
1592  m1 = t1 + t1[2] + 1; t1 += 4; t2 += 4;
1593  while ( t1 < m1 ) {
1594  if ( *t1 != *t2 ) goto nexti2;
1595  t1++; t2++;
1596  }
1597  t1 = tt1; t2 = tt2;
1598  }
1599  else t2 = tt2;
1600  }
1601  if ( t1 < sca1->pointer ) goto nexti2;
1602 /*
1603  Now we made it. We have to remove the terms and add the new term.
1604 */
1605  t1 = sca1->buffer; fill = t2 = sca2->buffer;
1606  sca2->numterms = 0;
1607  while ( t1 < sca1->pointer && t2 < sca2->pointer ) {
1608  if ( *t1 != *t2 ) { t2 += *t2; continue; }
1609  m1 = t1 + t1[2] + 1; tt1 = t1 + *t1;
1610  m2 = t2 + t2[2] + 1; tt2 = t2 + *t2;
1611  while ( m1 < tt1 && m2 < tt2 ) {
1612  if ( *m1 != *m2 || m1[1] != m2[1] || m1[2] != m2[2] ) break;
1613  m1 += 3; m2 += 3;
1614  }
1615  if ( m1 >= tt1 && m2 >= tt2 ) {
1616  t1 = tt1; t2 = tt2;
1617  }
1618  else {
1619  while ( t2 < tt2 ) *fill++ = *t2++;
1620  sca2->numterms++;
1621  }
1622  }
1623  while ( t2 < sca2->pointer ) {
1624  j = *t2; while ( --j >= 0 ) *fill++ = *t2++;
1625  sca2->numterms++;
1626  }
1627  *fill++ = 9; *fill++ = LNUMBER; *fill++ = 5;
1628  *fill++ = sign*3; *fill++ = 1; *fill++ = 1;
1629  *fill++ = ( i1+numobjects ) >> BITSINWORD;
1630  *fill++ = ( i1+numobjects ) & WORDMASK; *fill++ = 1;
1631  *fill = 0;
1632  sca2->numterms++;
1633  sca2->pointer = fill;
1634  SortOpti(i2);
1635 nexti2:;
1636  }
1637  }
1638 }
1639 
1640 /*
1641  #] CombiOpti:
1642  #[ TestNewSca:
1643 */
1644 
1645 LONG TestNewSca(LONG number, WORD *coef, WORD *ncoef)
1646 {
1647  GETIDENTITY
1648  SCALAR *sca = scabuffer + number, *s;
1649  WORD *t1, *t2, *m1, *m2, *coef2, ncoef2, n1, n2;
1650  int no;
1651  LONG i, ii;
1652  for ( i = 0, s = scabuffer; i < scanumber; i++, s++ ) {
1653  if ( i == number ) continue;
1654  if ( sca->numterms != s->numterms ) continue;
1655 /*
1656  Test 1: equal up to coefficient
1657 */
1658  t1 = sca->buffer;
1659  t2 = s->buffer;
1660  no = 0;
1661  while ( t1 < sca->pointer && t2 < s->pointer ) {
1662  m1 = t1 + t1[2] + 1; t1 += *t1;
1663  m2 = t2 + t2[2] + 1; t2 += *t2;
1664  while ( m1 < t1 && m2 < t2 ) {
1665  if ( *m1 != *m2 || m1[1] != m2[1] || m1[2] != m2[2] ) {
1666  no = 1; break; }
1667  m1 += 3; m2 += 3;
1668  }
1669  if ( m1 < t1 || m2 < t2 ) { no = 1; break; }
1670  }
1671  if ( !no ) {
1672 /*
1673  Now we have to be more precise with the coefficients
1674 */
1675  t1 = sca->buffer;
1676  t2 = s->buffer;
1677  n1 = REDLENG(t1[3]); n2 = REDLENG(t2[3]);
1678  DivRat(BHEAD (UWORD *)(t1+4),n1,(UWORD *)(t2+4),n2,(UWORD *)coef,ncoef);
1679  t1 += *t1; t2 += *t2;
1680  coef2 = coef + 2*ABS(*ncoef)+2;
1681  while ( t1 < sca->pointer && t2 < s->pointer ) {
1682  n1 = REDLENG(t1[3]); n2 = REDLENG(t2[3]);
1683  DivRat(BHEAD (UWORD *)(t1+4),n1,(UWORD *)(t2+4),n2,(UWORD *)coef2,&ncoef2);
1684  if ( *ncoef != ncoef2 ) break;
1685  ii = 2*ABS(*ncoef);
1686  while ( --ii >= 0 ) {
1687  if ( coef[ii] != coef2[ii] ) break;
1688  }
1689  if ( ii >= 0 ) break;
1690  t1 += *t1; t2 += *t2;
1691  }
1692  if ( t1 >= sca->pointer && t2 >= s->pointer ) {
1693 /*
1694  There is a match! And we can keep the factor.
1695 */
1696  sca->pointer = sca->buffer;
1697  sca->numterms = 0;
1698  *ncoef = INCLENG(*ncoef);
1699  return(i);
1700  }
1701  }
1702  }
1703  return(number);
1704 }
1705 
1706 /*
1707  #] TestNewSca:
1708  #[ NormOpti:
1709 */
1710 
1711 void NormOpti(WORD *term)
1712 {
1713  WORD *t, *m, *w, *tt, a;
1714  tt = m = term + term[2] + 4; t = term + *term;
1715  while ( m < t ) {
1716  w = m;
1717  while ( w >= tt && ( w[0] < w[-3]
1718  || ( w[0] == w[-3] && w[1] < w[-2] ) ) ) {
1719  a = w[0]; w[0] = w[-3]; w[-3] = a;
1720  a = w[1]; w[1] = w[-2]; w[-2] = a;
1721  a = w[2]; w[2] = w[-1]; w[-1] = a;
1722  w -= 3;
1723  }
1724  m += 3;
1725  }
1726 }
1727 
1728 /*
1729  #] NormOpti:
1730  #[ SortOpti:
1731 */
1732 
1733 void SortOpti(LONG number)
1734 {
1735  SCALAR *sca = scabuffer + number;
1736  WORD *newbuffer, *t, *m, **p, j;
1737  LONG i, newsize, num = 0;
1738  if ( sca->numterms <= 1 ) return;
1739  if ( (sca->numterms+sca->numterms/2) > sizepointers ) {
1740  if ( sortpointers ) M_free(sortpointers,"optisort");
1741  sizepointers = sca->numterms + sca->numterms/2 + 10;
1742  sortpointers = (WORD **)Malloc1(sizepointers*sizeof(WORD *),"optisort");
1743  }
1744  p = sortpointers;
1745  t = sca->buffer;
1746  while ( t < sca->pointer ) { *p++ = t; t += *t; num++; }
1747  if ( num != sca->numterms ) {
1748  MLOCK(ErrorMessageLock);
1749  MesPrint("Help! sca->numterms = %l, actual number = %l",
1750  sca->numterms, num);
1751  MUNLOCK(ErrorMessageLock);
1752  }
1753  helppointers = p;
1754  SplitOpti(sortpointers,num);
1755  newsize = ( sca->pointer - sca->buffer ) + 2;
1756  newbuffer = (WORD *)Malloc1(newsize*sizeof(WORD),"optisort2");
1757  m = newbuffer;
1758  for ( i = 0; i < num; i++ ) {
1759  t = sortpointers[i]; j = *t;
1760  while ( --j >= 0 ) *m++ = *t++;
1761  }
1762  sca->numterms = num;
1763  *m = 0;
1764  M_free(sca->buffer,"optisort2");
1765  sca->buffer = newbuffer; sca->pointer = m; sca->top = sca->buffer + newsize;
1766  sca->bufsize = newsize;
1767 }
1768 
1769 /*
1770  #] SortOpti:
1771  #[ SplitOpti:
1772 
1773  SplitMerge for the SortOpti routine.
1774 */
1775 
1776 void SplitOpti(WORD **pointers, LONG number)
1777 {
1778  WORD *t1, *t2, *m1, *m2, **p;
1779  int n;
1780  LONG left, right, i, j;
1781  if ( number <= 1 ) return;
1782  if ( number == 2 ) {
1783  t1 = pointers[0]; t2 = pointers[1];
1784  m1 = t1 + t1[2] + 1; m2 = t2 + t2[2] + 1;
1785  t1 += *t1; t2 += *t2;
1786  n = 0;
1787  while ( m1 < t1 && m2 < t2 ) {
1788  if ( m1[0] > m2[0] ) { n = 1; break; }
1789  else if ( m1[0] < m2[0] ) return;
1790  if ( m1[1] > m2[1] ) { n = 1; break; }
1791  else if ( m1[1] < m2[1] ) return;
1792  if ( m1[2] < m2[2] ) { n = 1; break; }
1793  else if ( m1[2] > m2[2] ) return;
1794  m1 += 3; m2 += 3;
1795  }
1796  if ( n > 0 || m1 < t1 ) {
1797  t1 = pointers[0]; pointers[0] = pointers[1]; pointers[1] = t1;
1798  }
1799  return;
1800  }
1801  left = number / 2;
1802  right = number - left;
1803  SplitOpti(pointers,left);
1804  SplitOpti(pointers+left,right);
1805  for ( i = 0; i < left; i++ ) helppointers[i] = pointers[i];
1806  i = 0; j = left; p = pointers;
1807  while ( i < left && j < number ) {
1808  t1 = helppointers[i]; t2 = pointers[j];
1809  m1 = t1 + t1[2] + 1; m2 = t2 + t2[2] + 1;
1810  t1 += *t1; t2 += *t2;
1811  n = 0;
1812  while ( m1 < t1 && m2 < t2 ) {
1813  if ( m1[0] > m2[0] ) { n = 1; break; }
1814  else if ( m1[0] < m2[0] ) { n = -1; break; }
1815  if ( m1[1] > m2[1] ) { n = 1; break; }
1816  else if ( m1[1] < m2[1] ) { n = -1; break; }
1817  if ( m1[2] < m2[2] ) { n = 1; break; }
1818  else if ( m1[2] > m2[2] ) { n = -1; break; }
1819  m1 += 3; m2 += 3;
1820  }
1821  if ( n > 0 || ( n == 0 && m1 < t1 ) ) { *p++ = pointers[j++]; }
1822  else *p++ = helppointers[i++];
1823  }
1824  while ( i < left ) *p++ = helppointers[i++];
1825 }
1826 
1827 /*
1828  #] SplitOpti:
1829 */
1830 
Definition: optim.c:58