FORM  4.1
factor.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 : factor.c
33 */
34 
35 #include "form3.h"
36 
37 /*
38  #] Includes :
39  #[ FactorIn :
40 
41  This routine tests for a factor in a dollar expression.
42 
43  Note that unlike with regular active or hidden expressions we cannot
44  add memory as easily as dollars are rather volatile.
45 */
46 int FactorIn(PHEAD WORD *term, WORD level)
47 {
48  GETBIDENTITY
49  WORD *t, *tstop, *m, *mm, *oldwork, *mstop, *n1, *n2, *n3, *n4, *n1stop, *n2stop;
50  WORD *r1, *r2, *r3, *r4, j, k, kGCD, kGCD2, kLCM, jGCD, kkLCM, jLCM, size;
51  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMb, *LCMc;
52  int fromwhere = 0, i;
53  DOLLARS d;
54  t = term; GETSTOP(t,tstop); t++;
55  while ( ( t < tstop ) && ( *t != FACTORIN || ( ( *t == FACTORIN )
56  && ( t[FUNHEAD] != -DOLLAREXPRESSION || t[1] != FUNHEAD+2 ) ) ) ) t += t[1];
57  if ( t >= tstop ) {
58  MLOCK(ErrorMessageLock);
59  MesPrint("Internal error. Could not find proper factorin_ function.");
60  MUNLOCK(ErrorMessageLock);
61  return(-1);
62  }
63  oldwork = AT.WorkPointer;
64  d = Dollars + t[FUNHEAD+1];
65 #ifdef WITHPTHREADS
66  {
67  int nummodopt, dtype = -1;
68  if ( AS.MultiThreaded ) {
69  for ( nummodopt = 0; nummodopt < NumModOptdollars; nummodopt++ ) {
70  if ( t[FUNHEAD+1] == ModOptdollars[nummodopt].number ) break;
71  }
72  if ( nummodopt < NumModOptdollars ) {
73  dtype = ModOptdollars[nummodopt].type;
74  if ( dtype == MODLOCAL ) {
75  d = ModOptdollars[nummodopt].dstruct+AT.identity;
76  }
77  }
78  }
79  }
80 #endif
81 
82  if ( d->type == DOLTERMS ) {
83  fromwhere = 1;
84  }
85  else if ( ( d = DolToTerms(BHEAD t[FUNHEAD+1]) ) == 0 ) {
86 /*
87  The variable cannot convert to an expression
88  We replace the function by 1.
89 */
90  m = oldwork; n1 = term;
91  while ( n1 < t ) *m++ = *n1++;
92  n1 = t + t[1]; tstop = term + *term;
93  while ( n1 < tstop ) *m++ = *n1++;
94  *oldwork = m - oldwork;
95  AT.WorkPointer = m;
96  if ( Generator(BHEAD oldwork,level) ) return(-1);
97  AT.WorkPointer = oldwork;
98  return(0);
99  }
100  if ( d->where[0] == 0 ) {
101  if ( fromwhere == 0 ) {
102  if ( d->factors ) M_free(d->factors,"Dollar factors");
103  M_free(d,"Dollar in FactorIn_");
104  }
105  return(0);
106  }
107 /*
108  Now we have an expression in d->where. Find the symbolic factor that
109  divides the expression and the numerical factor that makes all
110  coefficients integer.
111 
112  For the symbolic factor we make a copy of the first term, and then
113  go through all terms, scratching in the copy the objects that do not
114  occur in the terms.
115 */
116  m = oldwork;
117  mm = d->where;
118  k = *mm - ABS((mm[*mm-1]));
119  for ( j = 0; j < k; j++ ) *m++ = *mm++;
120  mstop = m;
121  *oldwork = k;
122 /*
123  The copy is in place. Now search through the terms. Start at the second term
124 */
125  mm = d->where + d->where[0];
126  while ( *mm ) {
127  m = oldwork+1;
128  r2 = mm+*mm;
129  r2 -= ABS(r2[-1]);
130  r1 = mm+1;
131  while ( m < mstop ) {
132  while ( r1 < r2 ) {
133  if ( *r1 != *m ) {
134  r1 += r1[1]; continue;
135  }
136 /*
137  Now the various cases
138  #[ SYMBOL :
139 */
140  if ( *m == SYMBOL ) {
141  n1 = m+2; n1stop = m+m[1];
142  n2stop = r1+r1[1];
143  while ( n1 < n1stop ) {
144  n2 = r1+2;
145  while ( n2 < n2stop ) {
146  if ( *n1 != *n2 ) { n2 += 2; continue; }
147  if ( n1[1] > 0 ) {
148  if ( n2[1] < 0 ) { n2 += 2; continue; }
149  if ( n2[1] < n1[1] ) n1[1] = n2[1];
150  }
151  else {
152  if ( n2[1] > 0 ) { n2 += 2; continue; }
153  if ( n2[1] > n1[1] ) n1[1] = n2[1];
154  }
155  break;
156  }
157  if ( n2 >= n2stop ) { /* scratch symbol */
158  if ( m[1] == 4 ) goto scratch;
159  m[1] -= 2;
160  n3 = n1; n4 = n1+2;
161  while ( n4 < mstop ) *n3++ = *n4++;
162  *oldwork = n3 - oldwork;
163  mstop -= 2; n1stop -= 2;
164  continue;
165  }
166  n1 += 2;
167  }
168  break;
169  }
170 /*
171  #] SYMBOL :
172  #[ DOTPRODUCT :
173 */
174  else if ( *m == DOTPRODUCT ) {
175  n1 = m+2; n1stop = m+m[1];
176  n2stop = r1+r1[1];
177  while ( n1 < n1stop ) {
178  n2 = r1+2;
179  while ( n2 < n2stop ) {
180  if ( *n1 != *n2 || n1[1] != n2[1] ) { n2 += 3; continue; }
181  if ( n1[2] > 0 ) {
182  if ( n2[2] < 0 ) { n2 += 3; continue; }
183  if ( n2[2] < n1[2] ) n1[2] = n2[2];
184  }
185  else {
186  if ( n2[2] > 0 ) { n2 += 3; continue; }
187  if ( n2[2] > n1[2] ) n1[2] = n2[2];
188  }
189  break;
190  }
191  if ( n2 >= n2stop ) { /* scratch symbol */
192  if ( m[1] == 5 ) goto scratch;
193  m[1] -= 3;
194  n3 = n1; n4 = n1+3;
195  while ( n4 < mstop ) *n3++ = *n4++;
196  *oldwork = n3 - oldwork;
197  mstop -= 3; n1stop -= 3;
198  continue;
199  }
200  n1 += 3;
201  }
202  break;
203  }
204 /*
205  #] DOTPRODUCT :
206  #[ VECTOR :
207 */
208  else if ( *m == VECTOR ) {
209 /*
210  Here we have to be careful if there is more than
211  one of the same
212 */
213  n1 = m+2; n1stop = m+m[1];
214  n2 = r1+2;n2stop = r1+r1[1];
215  while ( n1 < n1stop ) {
216  while ( n2 < n2stop ) {
217  if ( *n1 == *n2 && n1[1] == n2[1] ) {
218  n2 += 2; goto nextn1;
219  }
220  n2 += 2;
221  }
222  if ( n2 >= n2stop ) { /* scratch symbol */
223  if ( m[1] == 4 ) goto scratch;
224  m[1] -= 2;
225  n3 = n1; n4 = n1+2;
226  while ( n4 < mstop ) *n3++ = *n4++;
227  *oldwork = n3 - oldwork;
228  mstop -= 2; n1stop -= 2;
229  continue;
230  }
231  n2 = r1+2;
232 nextn1: n1 += 2;
233  }
234  break;
235  }
236 /*
237  #] VECTOR :
238  #[ REMAINDER :
239 */
240  else {
241 /*
242  Again: watch for multiple occurrences of the same object
243 */
244  if ( m[1] != r1[1] ) { r1 += r1[1]; continue; }
245  for ( j = 2; j < m[1]; j++ ) {
246  if ( m[j] != r1[j] ) break;
247  }
248  if ( j < m[1] ) { r1 += r1[1]; continue; }
249  r1 += r1[1]; /* to restart at the next potential match */
250  goto nextm; /* match */
251  }
252 /*
253  #] REMAINDER :
254 */
255  }
256  if ( r1 >= r2 ) { /* no factor! */
257 scratch:;
258  r3 = m + m[1]; r4 = m;
259  while ( r3 < mstop ) *r4++ = *r3++;
260  *oldwork = r4 - oldwork;
261  if ( *oldwork == 1 ) goto nofactor;
262  mstop = r4;
263  r1 = mm + 1;
264  continue;
265  }
266  r1 = mm + 1;
267 nextm: m += m[1];
268  }
269  mm = mm + *mm;
270  }
271 
272 nofactor:;
273 /*
274  For the coefficient we have to determine the LCM of the denominators
275  and the GCD of the numerators.
276 */
277  GCDbuffer = NumberMalloc("FactorIn"); GCDbuffer2 = NumberMalloc("FactorIn");
278  LCMbuffer = NumberMalloc("FactorIn"); LCMb = NumberMalloc("FactorIn"); LCMc = NumberMalloc("FactorIn");
279  r1 = d->where;
280 /*
281  First take the first term to load up the LCM and the GCD
282 */
283  r2 = r1 + *r1;
284  j = r2[-1];
285  r3 = r2 - ABS(j);
286  k = REDLENG(j);
287  if ( k < 0 ) k = -k;
288  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
289  for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
290  k = REDLENG(j);
291  if ( k < 0 ) k = -k;
292  r3 += k;
293  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
294  for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
295  r1 = r2;
296 /*
297  Now go through the rest of the terms in this dollar buffer.
298 */
299  while ( *r1 ) {
300  r2 = r1 + *r1;
301  j = r2[-1];
302  r3 = r2 - ABS(j);
303  k = REDLENG(j);
304  if ( k < 0 ) k = -k;
305  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
306  if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
307 /*
308  GCD is already 1
309 */
310  }
311  else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
312  if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
313  goto onerror;
314  }
315  kGCD = kGCD2;
316  for ( i = 0; i < kGCD; i++ ) GCDbuffer[i] = GCDbuffer2[i];
317  }
318  else {
319  kGCD = 1; GCDbuffer[0] = 1;
320  }
321  k = REDLENG(j);
322  if ( k < 0 ) k = -k;
323  r3 += k;
324  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
325  if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
326  for ( kLCM = 0; kLCM < k; kLCM++ )
327  LCMbuffer[kLCM] = r3[kLCM];
328  }
329  else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
330  if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
331  goto onerror;
332  }
333  DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
334  MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
335  for ( kLCM = 0; kLCM < jLCM; kLCM++ )
336  LCMbuffer[kLCM] = LCMc[kLCM];
337  }
338  else {} /* LCM doesn't change */
339  r1 = r2;
340  }
341 /*
342  Now put the factor together: GCD/LCM
343 */
344  r3 = (WORD *)(GCDbuffer);
345  if ( kGCD == kLCM ) {
346  for ( jGCD = 0; jGCD < kGCD; jGCD++ )
347  r3[jGCD+kGCD] = LCMbuffer[jGCD];
348  k = kGCD;
349  }
350  else if ( kGCD > kLCM ) {
351  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
352  r3[jGCD+kGCD] = LCMbuffer[jGCD];
353  for ( jGCD = kLCM; jGCD < kGCD; jGCD++ )
354  r3[jGCD+kGCD] = 0;
355  k = kGCD;
356  }
357  else {
358  for ( jGCD = kGCD; jGCD < kLCM; jGCD++ )
359  r3[jGCD] = 0;
360  for ( jGCD = 0; jGCD < kLCM; jGCD++ )
361  r3[jGCD+kLCM] = LCMbuffer[jGCD];
362  k = kLCM;
363  }
364  j = 2*k+1;
365  mm = m = oldwork + oldwork[0];
366 /*
367  Now compose the new term
368 */
369  n1 = term;
370  while ( n1 < t ) *m++ = *n1++;
371  n1 += n1[1];
372  n2 = oldwork+1;
373  while ( n2 < mm ) *m++ = *n2++;
374  while ( n1 < tstop ) *m++ = *n1++;
375 /*
376  And the coefficient
377 */
378  size = term[*term-1];
379  size = REDLENG(size);
380  if ( MulRat(BHEAD (UWORD *)tstop,size,(UWORD *)r3,k,
381  (UWORD *)m,&size) ) goto onerror;
382  size = INCLENG(size);
383  k = size < 0 ? -size: size;
384  m[k-1] = size;
385  m += k;
386  *mm = (WORD)(m - mm);
387  AT.WorkPointer = m;
388  if ( Generator(BHEAD mm,level) ) goto onerror;
389  AT.WorkPointer = oldwork;
390  if ( fromwhere == 0 ) {
391  if ( d->factors ) M_free(d->factors,"Dollar factors");
392  M_free(d,"Dollar in FactorIn");
393  }
394  NumberFree(GCDbuffer,"FactorIn"); NumberFree(GCDbuffer2,"FactorIn");
395  NumberFree(LCMbuffer,"FactorIn"); NumberFree(LCMb,"FactorIn"); NumberFree(LCMc,"FactorIn");
396  return(0);
397 onerror:
398  AT.WorkPointer = oldwork;
399  MLOCK(ErrorMessageLock);
400  MesCall("FactorIn");
401  MUNLOCK(ErrorMessageLock);
402  NumberFree(GCDbuffer,"FactorIn"); NumberFree(GCDbuffer2,"FactorIn");
403  NumberFree(LCMbuffer,"FactorIn"); NumberFree(LCMb,"FactorIn"); NumberFree(LCMc,"FactorIn");
404  return(-1);
405 }
406 
407 /*
408  #] FactorIn :
409  #[ FactorInExpr :
410 
411  This routine tests for a factor in an active or hidden expression.
412 
413  The factor from the last call is stored in a cache.
414  Main problem here is whether the cache is global or private to each thread.
415  A global cache gives most likely the same traffic jam for the computation
416  as a local cache. The local cache may stay valid longer.
417  In the future we may make it such that threads can look at the cache
418  of the others, or even whether a certain factor is under construction.
419 */
420 
421 int FactorInExpr(PHEAD WORD *term, WORD level)
422 {
423  GETBIDENTITY
424  WORD *t, *tstop, *m, *oldwork, *mstop, *n1, *n2, *n3, *n4, *n1stop, *n2stop;
425  WORD *r1, *r2, *r3, *r4, j, k, kGCD, kGCD2, kLCM, jGCD, kkLCM, jLCM, size, sign;
426  WORD *newterm, expr = 0;
427  WORD olddeferflag = AR.DeferFlag, oldgetfile = AR.GetFile, oldhold = AR.KeptInHold;
428  WORD newgetfile, newhold;
429  int i;
430  EXPRESSIONS e;
431  FILEHANDLE *file = 0;
432  POSITION position, oldposition, startposition;
433  WORD *oldcpointer = AR.CompressPointer;
434  UWORD *GCDbuffer, *GCDbuffer2, *LCMbuffer, *LCMb, *LCMc;
435  GCDbuffer = NumberMalloc("FactorInExpr"); GCDbuffer2 = NumberMalloc("FactorInExpr");
436  LCMbuffer = NumberMalloc("FactorInExpr"); LCMb = NumberMalloc("FactorInExpr"); LCMc = NumberMalloc("FactorInExpr");
437  t = term; GETSTOP(t,tstop); t++;
438  while ( t < tstop ) {
439  if ( *t == FACTORIN && t[1] == FUNHEAD+2 && t[FUNHEAD] == -EXPRESSION ) {
440  expr = t[FUNHEAD+1];
441  break;
442  }
443  t += t[1];
444  }
445  if ( t >= tstop ) {
446  MLOCK(ErrorMessageLock);
447  MesPrint("Internal error. Could not find proper factorin_ function.");
448  MUNLOCK(ErrorMessageLock);
449  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
450  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
451  return(-1);
452  }
453  oldwork = AT.WorkPointer;
454  if ( AT.previousEfactor && ( expr == AT.previousEfactor[0] ) ) {
455 /*
456  We have a hit in the cache. Construct the new term.
457  At the moment AT.previousEfactor[1] is reserved for future flags
458 */
459  goto PutTheFactor;
460  }
461 /*
462  No hit. We have to do the work. We start with constructing the factor
463  in the WorkSpace. Later we will move it to the cache.
464  Finally we will jump to PutTheFactor.
465 */
466  e = Expressions + expr;
467  switch ( e->status ) {
468  case LOCALEXPRESSION:
469  case SKIPLEXPRESSION:
470  case DROPLEXPRESSION:
471  case GLOBALEXPRESSION:
472  case SKIPGEXPRESSION:
473  case DROPGEXPRESSION:
474 /*
475  Expression is to be found in the input Scratch file.
476  Set the file handle and the position.
477  The rest is done by GetTerm.
478 */
479  newhold = 0;
480  newgetfile = 1;
481  file = AR.infile;
482  break;
483  case HIDDENLEXPRESSION:
484  case HIDDENGEXPRESSION:
485  case HIDELEXPRESSION:
486  case HIDEGEXPRESSION:
487  case DROPHLEXPRESSION:
488  case DROPHGEXPRESSION:
489  case UNHIDELEXPRESSION:
490  case UNHIDEGEXPRESSION:
491 /*
492  Expression is to be found in the hide Scratch file.
493  Set the file handle and the position.
494  The rest is done by GetTerm.
495 */
496  newhold = 0;
497  newgetfile = 2;
498  file = AR.hidefile;
499  break;
500  case STOREDEXPRESSION:
501 /*
502  This is an 'illegal' case
503 */
504  MLOCK(ErrorMessageLock);
505  MesPrint("Error: factorin_ cannot determine factors in stored expressions.");
506  MUNLOCK(ErrorMessageLock);
507  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
508  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
509  return(-1);
510  case DROPPEDEXPRESSION:
511 /*
512  We replace the function by 1.
513 */
514  m = oldwork; n1 = term;
515  while ( n1 < t ) *m++ = *n1++;
516  n1 = t + t[1]; tstop = term + *term;
517  while ( n1 < tstop ) *m++ = *n1++;
518  *oldwork = m - oldwork;
519  AT.WorkPointer = m;
520  if ( Generator(BHEAD oldwork,level) ) {
521  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
522  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
523  return(-1);
524  }
525  AT.WorkPointer = oldwork;
526  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
527  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
528  return(0);
529  default:
530  MLOCK(ErrorMessageLock);
531  MesPrint("Error: Illegal expression in factorinexpr.");
532  MUNLOCK(ErrorMessageLock);
533  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
534  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
535  return(-1);
536  }
537 /*
538  Before we start with the file we set the buffers for the coefficient
539  For the coefficient we have to determine the LCM of the denominators
540  and the GCD of the numerators.
541 */
542  position = AS.OldOnFile[expr];
543  AR.DeferFlag = 0; AR.KeptInHold = newhold; AR.GetFile = newgetfile;
544  SeekScratch(file,&oldposition);
545  SetScratch(file,&position);
546  if ( GetTerm(BHEAD oldwork) <= 0 ) {
547  MLOCK(ErrorMessageLock);
548  MesPrint("(5) Expression %d has problems in scratchfile",expr);
549  MUNLOCK(ErrorMessageLock);
550  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
551  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
552  return(-1);
553  }
554  SeekScratch(file,&startposition);
555  SeekScratch(file,&position);
556 /*
557  Load the first term in the workspace
558 */
559  if ( GetTerm(BHEAD oldwork) == 0 ) {
560  SetScratch(file,&oldposition); /* We still need this untill Processor is clean */
561  AR.DeferFlag = olddeferflag;
562  oldwork[0] = 4; oldwork[1] = 1; oldwork[2] = 1; oldwork[3] = 3;
563  goto Complete;
564  }
565  SeekScratch(file,&position);
566  AR.DeferFlag = olddeferflag; AR.KeptInHold = oldhold; AR.GetFile = oldgetfile;
567 
568  r2 = m = oldwork + *oldwork;
569  j = m[-1];
570  m -= ABS(j);
571  *oldwork = (WORD)(m-oldwork);
572  AT.WorkPointer = newterm = mstop = m;
573 /*
574  Now take the coefficient of the first term to load up the LCM and the GCD
575 */
576  r3 = m;
577  k = REDLENG(j);
578  if ( k < 0 ) { k = -k; sign = -1; }
579  else { sign = 1; }
580  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
581  for ( kGCD = 0; kGCD < k; kGCD++ ) GCDbuffer[kGCD] = r3[kGCD];
582  k = REDLENG(j);
583  if ( k < 0 ) k = -k;
584  r3 += k;
585  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
586  for ( kLCM = 0; kLCM < k; kLCM++ ) LCMbuffer[kLCM] = r3[kLCM];
587 /*
588  The copy and the coefficient are in place. Now search through the terms.
589 */
590  for (;;) {
591  AR.DeferFlag = 0; AR.KeptInHold = newhold; AR.GetFile = newgetfile;
592  SetScratch(file,&position);
593  size = GetTerm(BHEAD newterm);
594  SeekScratch(file,&position);
595  AR.DeferFlag = olddeferflag; AR.KeptInHold = oldhold; AR.GetFile = oldgetfile;
596  if ( size == 0 ) break;
597  m = oldwork+1;
598  r2 = newterm + *newterm;
599  r2 -= ABS(r2[-1]);
600  r1 = newterm+1;
601  while ( m < mstop ) {
602  while ( r1 < r2 ) {
603  if ( *r1 != *m ) {
604  r1 += r1[1]; continue;
605  }
606 /*
607  Now the various cases
608  #[ SYMBOL :
609 */
610  if ( *m == SYMBOL ) {
611  n1 = m+2; n1stop = m+m[1];
612  n2stop = r1+r1[1];
613  while ( n1 < n1stop ) {
614  n2 = r1+2;
615  while ( n2 < n2stop ) {
616  if ( *n1 != *n2 ) { n2 += 2; continue; }
617  if ( n1[1] > 0 ) {
618  if ( n2[1] < 0 ) { n2 += 2; continue; }
619  if ( n2[1] < n1[1] ) n1[1] = n2[1];
620  }
621  else {
622  if ( n2[1] > 0 ) { n2 += 2; continue; }
623  if ( n2[1] > n1[1] ) n1[1] = n2[1];
624  }
625  break;
626  }
627  if ( n2 >= n2stop ) { /* scratch symbol */
628  if ( m[1] == 4 ) goto scratch;
629  m[1] -= 2;
630  n3 = n1; n4 = n1+2;
631  while ( n4 < mstop ) *n3++ = *n4++;
632  *oldwork = n3 - oldwork;
633  mstop -= 2; n1stop -= 2;
634  continue;
635  }
636  n1 += 2;
637  }
638  break;
639  }
640 /*
641  #] SYMBOL :
642  #[ DOTPRODUCT :
643 */
644  else if ( *m == DOTPRODUCT ) {
645  n1 = m+2; n1stop = m+m[1];
646  n2stop = r1+r1[1];
647  while ( n1 < n1stop ) {
648  n2 = r1+2;
649  while ( n2 < n2stop ) {
650  if ( *n1 != *n2 || n1[1] != n2[1] ) { n2 += 3; continue; }
651  if ( n1[2] > 0 ) {
652  if ( n2[2] < 0 ) { n2 += 3; continue; }
653  if ( n2[2] < n1[2] ) n1[2] = n2[2];
654  }
655  else {
656  if ( n2[2] > 0 ) { n2 += 3; continue; }
657  if ( n2[2] > n1[2] ) n1[2] = n2[2];
658  }
659  break;
660  }
661  if ( n2 >= n2stop ) { /* scratch dotproduct */
662  if ( m[1] == 5 ) goto scratch;
663  m[1] -= 3;
664  n3 = n1; n4 = n1+3;
665  while ( n4 < mstop ) *n3++ = *n4++;
666  *oldwork = n3 - oldwork;
667  mstop -= 3; n1stop -= 3;
668  continue;
669  }
670  n1 += 3;
671  }
672  break;
673  }
674 /*
675  #] DOTPRODUCT :
676  #[ VECTOR :
677 */
678  else if ( *m == VECTOR ) {
679 /*
680  Here we have to be careful if there is more than
681  one of the same
682 */
683  n1 = m+2; n1stop = m+m[1];
684  n2 = r1+2;n2stop = r1+r1[1];
685  while ( n1 < n1stop ) {
686  while ( n2 < n2stop ) {
687  if ( *n1 == *n2 && n1[1] == n2[1] ) {
688  n2 += 2; goto nextn1;
689  }
690  n2 += 2;
691  }
692  if ( n2 >= n2stop ) { /* scratch vector */
693  if ( m[1] == 4 ) goto scratch;
694  m[1] -= 2;
695  n3 = n1; n4 = n1+2;
696  while ( n4 < mstop ) *n3++ = *n4++;
697  *oldwork = n3 - oldwork;
698  mstop -= 2; n1stop -= 2;
699  continue;
700  }
701  n2 = r1+2;
702 nextn1: n1 += 2;
703  }
704  break;
705  }
706 /*
707  #] VECTOR :
708  #[ REMAINDER :
709 */
710  else {
711 /*
712  Again: watch for multiple occurrences of the same object
713 */
714  if ( m[1] != r1[1] ) { r1 += r1[1]; continue; }
715  for ( j = 2; j < m[1]; j++ ) {
716  if ( m[j] != r1[j] ) break;
717  }
718  if ( j < m[1] ) { r1 += r1[1]; continue; }
719  r1 += r1[1]; /* to restart at the next potential match */
720  goto nextm; /* match */
721  }
722 /*
723  #] REMAINDER :
724 */
725  }
726  if ( r1 >= r2 ) { /* no factor! */
727 scratch:;
728  r3 = m + m[1]; r4 = m;
729  while ( r3 < mstop ) *r4++ = *r3++;
730  *oldwork = r4 - oldwork;
731  if ( *oldwork == 1 ) goto nofactor;
732  mstop = r4;
733  r1 = newterm + 1;
734  continue;
735  }
736  r1 = newterm + 1;
737 nextm: m += m[1];
738  }
739 nofactor:;
740 /*
741  Now the coefficient
742 */
743  r2 = newterm + *newterm;
744  j = r2[-1];
745  r3 = r2 - ABS(j);
746  k = REDLENG(j);
747  if ( k < 0 ) k = -k;
748  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
749  if ( ( ( GCDbuffer[0] == 1 ) && ( kGCD == 1 ) ) ) {
750 /*
751  GCD is already 1
752 */
753  }
754  else if ( ( ( k != 1 ) || ( r3[0] != 1 ) ) ) {
755  if ( GcdLong(BHEAD GCDbuffer,kGCD,(UWORD *)r3,k,GCDbuffer2,&kGCD2) ) {
756  goto onerror;
757  }
758  kGCD = kGCD2;
759  for ( i = 0; i < kGCD; i++ ) GCDbuffer[i] = GCDbuffer2[i];
760  }
761  else {
762  kGCD = 1; GCDbuffer[0] = 1;
763  }
764  k = REDLENG(j);
765  if ( k < 0 ) k = -k;
766  r3 += k;
767  while ( ( k > 1 ) && ( r3[k-1] == 0 ) ) k--;
768  if ( ( ( LCMbuffer[0] == 1 ) && ( kLCM == 1 ) ) ) {
769  for ( kLCM = 0; kLCM < k; kLCM++ )
770  LCMbuffer[kLCM] = r3[kLCM];
771  }
772  else if ( ( k != 1 ) || ( r3[0] != 1 ) ) {
773  if ( GcdLong(BHEAD LCMbuffer,kLCM,(UWORD *)r3,k,LCMb,&kkLCM) ) {
774  goto onerror;
775  }
776  DivLong((UWORD *)r3,k,LCMb,kkLCM,LCMb,&kkLCM,LCMc,&jLCM);
777  MulLong(LCMbuffer,kLCM,LCMb,kkLCM,LCMc,&jLCM);
778  for ( kLCM = 0; kLCM < jLCM; kLCM++ )
779  LCMbuffer[kLCM] = LCMc[kLCM];
780  }
781  else {} /* LCM doesn't change */
782  }
783  SetScratch(file,&oldposition); /* Needed until Processor is thread safe */
784  AR.DeferFlag = olddeferflag;
785 /*
786  Now put the term together in oldwork: GCD/LCM
787  We have already the algebraic contents there.
788 */
789  r3 = (WORD *)(GCDbuffer);
790  r4 = (WORD *)(LCMbuffer);
791  r1 = oldwork + *oldwork;
792  if ( kGCD == kLCM ) {
793  for ( jGCD = 0; jGCD < kGCD; jGCD++ ) *r1++ = *r3++;
794  for ( jGCD = 0; jGCD < kGCD; jGCD++ ) *r1++ = *r4++;
795  k = 2*kGCD+1;
796  }
797  else if ( kGCD > kLCM ) {
798  for ( jGCD = 0; jGCD < kGCD; jGCD++ ) *r1++ = *r3++;
799  for ( jGCD = 0; jGCD < kLCM; jGCD++ ) *r1++ = *r4++;
800  for ( jGCD = kLCM; jGCD < kGCD; jGCD++ ) *r1++ = 0;
801  k = 2*kGCD+1;
802  }
803  else {
804  for ( jGCD = 0; jGCD < kGCD; jGCD++ ) *r1++ = *r3++;
805  for ( jGCD = kGCD; jGCD < kLCM; jGCD++ ) *r1++ = 0;
806  for ( jGCD = 0; jGCD < kLCM; jGCD++ ) *r1++ = *r4++;
807  k = 2*kLCM+1;
808  }
809  if ( sign < 0 ) *r1++ = -k;
810  else *r1++ = k;
811  *oldwork = (WORD)(r1-oldwork);
812 /*
813  Now put the new term in the cache
814 */
815 Complete:;
816  if ( AT.previousEfactor ) M_free(AT.previousEfactor,"Efactor cache");
817  AT.previousEfactor = (WORD *)Malloc1((*oldwork+2)*sizeof(WORD),"Efactor cache");
818  AT.previousEfactor[0] = expr;
819  r1 = oldwork; r2 = AT.previousEfactor + 2; k = *oldwork;
820  NCOPY(r2,r1,k)
821  AT.previousEfactor[1] = 0;
822 /*
823  Now we construct the new term in the workspace.
824 */
825 PutTheFactor:;
826  if ( AT.WorkPointer + AT.previousEfactor[2] >= AT.WorkTop ) {
827  MLOCK(ErrorMessageLock);
828  MesWork();
829  MesPrint("Called from factorin_");
830  MUNLOCK(ErrorMessageLock);
831  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
832  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
833  return(-1);
834  }
835  n1 = oldwork; n2 = term; while ( n2 < t ) *n1++ = *n2++;
836  n2 = AT.previousEfactor+2; GETSTOP(n2,n2stop); n3 = n2 + *n2;
837  n2++; while ( n2 < n2stop ) *n1++ = *n2++;
838  n2 = t + t[1]; while ( n2 < tstop ) *n1++ = *n2++;
839  size = term[*term-1];
840  size = REDLENG(size);
841  k = n3[-1]; k = REDLENG(k);
842  if ( MulRat(BHEAD (UWORD *)tstop,size,(UWORD *)n2stop,k,
843  (UWORD *)n1,&size) ) goto onerror;
844  size = INCLENG(size);
845  k = size < 0 ? -size: size;
846  n1 += k; n1[-1] = size;
847  *oldwork = n1 - oldwork;
848  AT.WorkPointer = n1;
849  if ( Generator(BHEAD oldwork,level) ) {
850  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
851  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
852  return(-1);
853  }
854  AT.WorkPointer = oldwork;
855  AR.CompressPointer = oldcpointer;
856  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
857  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
858  return(0);
859 onerror:
860  AT.WorkPointer = oldwork;
861  AR.CompressPointer = oldcpointer;
862  MLOCK(ErrorMessageLock);
863  MesCall("FactorInExpr");
864  MUNLOCK(ErrorMessageLock);
865  NumberFree(GCDbuffer,"FactorInExpr"); NumberFree(GCDbuffer2,"FactorInExpr");
866  NumberFree(LCMbuffer,"FactorInExpr"); NumberFree(LCMb,"FactorInExpr"); NumberFree(LCMc,"FactorInExpr");
867  return(-1);
868 }
869 
870 /*
871  #] FactorInExpr :
872 */
873 
Definition: structs.h:618
#define PHEAD
Definition: ftypes.h:56
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:2865