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