FORM 4.3
factor.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 : 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*/
46int 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;
232nextn1: 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! */
257scratch:;
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;
267nextm: m += m[1];
268 }
269 mm = mm + *mm;
270 }
271
272nofactor:;
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);
397onerror:
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
421int 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;
702nextn1: 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! */
727scratch:;
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;
737nextm: m += m[1];
738 }
739nofactor:;
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*/
815Complete:;
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*/
825PutTheFactor:;
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);
859onerror:
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*/
WORD Generator(PHEAD WORD *, WORD)
Definition: proces.c:3101
Definition: structs.h:633