FORM v5.0.0-35-g6318119
flintinterface.cc
Go to the documentation of this file.
1
6/* #[ License : */
7/*
8 * Copyright (C) 1984-2026 J.A.M. Vermaseren
9 * When using this file you are requested to refer to the publication
10 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
11 * This is considered a matter of courtesy as the development was paid
12 * for by FOM the Dutch physics granting agency and we would like to
13 * be able to track its scientific use to convince FOM of its value
14 * for the community.
15 *
16 * This file is part of FORM.
17 *
18 * FORM is free software: you can redistribute it and/or modify it under the
19 * terms of the GNU General Public License as published by the Free Software
20 * Foundation, either version 3 of the License, or (at your option) any later
21 * version.
22 *
23 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
24 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
25 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
26 * details.
27 *
28 * You should have received a copy of the GNU General Public License along
29 * with FORM. If not, see <http://www.gnu.org/licenses/>.
30 */
31/* #] License : */
32
33#ifdef HAVE_CONFIG_H
34#ifndef CONFIG_H_INCLUDED
35#define CONFIG_H_INCLUDED
36#include <config.h>
37#endif
38#endif
39
40extern "C" {
41#include "form3.h"
42}
43
44#include "flintinterface.h"
45
46/*
47 #[ Types
48 * Use fixed-size integer types (int32_t, uint32_t, int64_t, uint64_t) except when we refer
49 * specifically to FORM term data, when we use WORD. This code has only been tested on 64bit
50 * systems where WORDs are int32_t, and some functions (fmpz_get_form, fmpz_set_form) require this
51 * to be the case. Enforce this:
52 */
53static_assert(sizeof(WORD) == sizeof(int32_t), "flint interface expects 32bit WORD");
54static_assert(sizeof(UWORD) == sizeof(uint32_t), "flint interface expects 32bit WORD");
55static_assert(BITSINWORD == 32, "flint interface expects 32bit WORD");
56static_assert(sizeof(ulong) == sizeof(uint64_t), "flint interface expects ulong is uint64_t");
57static_assert(sizeof(slong) == sizeof(int64_t), "flint interface expects slong is int64_t");
58
59/*
60 * Flint functions take arguments or return values which may be "slong" or "ulong" in its
61 * documentation, and these are int64_t and uint64_t respectively.
62 *
63 #] Types
64*/
65
66/*
67 * FLINT's univariate poly has a dense representation. For sufficiently sparse polynomials it is
68 * faster to use mpoly instead, which is sparse. For a density <= this threshold we switch, where
69 * the density defined as is "number of terms" / "maximum degree".
70 */
71#define UNIVARIATE_DENSITY_THR 0.02f
72
73/*
74 #[ flint::cleanup :
75*/
76void flint::cleanup(void) {
77 flint_cleanup();
78}
79/*
80 #] flint::cleanup :
81 #[ flint::cleanup_master :
82*/
83void flint::cleanup_master(void) {
84 flint_cleanup_master();
85}
86/*
87 #] flint::cleanup_master :
88
89 #[ flint::divmod_mpoly :
90*/
91WORD* flint::divmod_mpoly(PHEAD const WORD *a, const WORD *b, const bool return_rem,
92 const WORD must_fit_term, const var_map_t &var_map) {
93
94 flint::mpoly_ctx ctx(var_map.size());
95 flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
96
97 flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
98 flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
99
100 // The input won't have any symbols with negative powers, but there may be rational
101 // coefficients. Verify this:
102 if ( fmpz_mpoly_is_fmpz(denpa.d, ctx.d) != 1 ) {
103 MLOCK(ErrorMessageLock);
104 MesPrint("flint::divmod_mpoly: error: denpa is non-constant");
105 MUNLOCK(ErrorMessageLock);
106 Terminate(-1);
107 }
108 if ( fmpz_mpoly_is_fmpz(denpb.d, ctx.d) != 1 ) {
109 MLOCK(ErrorMessageLock);
110 MesPrint("flint::divmod_mpoly: error: denpb is non-constant");
111 MUNLOCK(ErrorMessageLock);
112 Terminate(-1);
113 }
114
115
116 flint::fmpz scale;
117 flint::mpoly div(ctx.d), rem(ctx.d);
118
119 fmpz_mpoly_quasidivrem(scale.d, div.d, rem.d, pa.d, pb.d, ctx.d);
120
121 // The quotient must be multiplied by the denominator of the divisor
122 fmpz_mpoly_mul(div.d, div.d, denpb.d, ctx.d);
123
124 // The overall denominator of both div and rem is given by scale*denpa. This we will pass to
125 // to_argument_mpoly's "denscale" argument which performs the division in the result. We have
126 // already checked denpa is just an fmpz.
127 fmpz_mul(scale.d, scale.d, fmpz_mpoly_term_coeff_ref(denpa.d, 0, ctx.d));
128
129
130 // Determine the size of the output by passing write = false.
131 const bool with_arghead = false;
132 bool write = false;
133 const uint64_t prev_size = 0;
134 const uint64_t out_size = return_rem ?
135 (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
136 rem.d, var_map, ctx.d, scale.d)
137 :
138 (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
139 div.d, var_map, ctx.d, scale.d)
140 ;
141 WORD* res = (WORD *)Malloc1(sizeof(WORD)*out_size, "flint::divrem_mpoly");
142
143
144 // Write out the result
145 write = true;
146 if ( return_rem ) {
147 (uint64_t)flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
148 rem.d, var_map, ctx.d, scale.d);
149 }
150 else {
151 (uint64_t)flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
152 div.d, var_map, ctx.d, scale.d);
153 }
154
155 return res;
156}
157/*
158 #] flint::divmod_mpoly :
159 #[ flint::divmod_poly :
160*/
161WORD* flint::divmod_poly(PHEAD const WORD *a, const WORD *b, const bool return_rem,
162 const WORD must_fit_term, const var_map_t &var_map) {
163
164 flint::poly pa, pb, denpa, denpb;
165
166 flint::from_argument_poly(pa.d, denpa.d, a, false);
167 flint::from_argument_poly(pb.d, denpb.d, b, false);
168
169 // The input won't have any symbols with negative powers, but there may be rational
170 // coefficients. Verify this:
171 if ( fmpz_poly_length(denpa.d) != 1 ) {
172 MLOCK(ErrorMessageLock);
173 MesPrint("flint::divmod_poly: error: denpa is non-constant");
174 MUNLOCK(ErrorMessageLock);
175 Terminate(-1);
176 }
177 if ( fmpz_poly_length(denpb.d) != 1 ) {
178 MLOCK(ErrorMessageLock);
179 MesPrint("flint::divmod_poly: error: denpb is non-constant");
180 MUNLOCK(ErrorMessageLock);
181 Terminate(-1);
182 }
183
184 flint::fmpz scale;
185 uint64_t scale_power = 0;
186 flint::poly div, rem;
187
188 // Get the leading coefficient of pb. fmpz_poly_pseudo_divrem returns only the scaling power.
189 fmpz_poly_get_coeff_fmpz(scale.d, pb.d, fmpz_poly_degree(pb.d));
190
191 fmpz_poly_pseudo_divrem(div.d, rem.d, (ulong*)&scale_power, pa.d, pb.d);
192
193 // The quotient must be multiplied by the denominator of the divisor
194 fmpz_poly_mul(div.d, div.d, denpb.d);
195
196 // The overall denominator of both div and rem is given by scale^scale_power*denpa. This we will
197 // pass to to_argument_poly's "denscale" argument which performs the division in the result. We
198 // have already checked denpa is just an fmpz.
199 fmpz_pow_ui(scale.d, scale.d, scale_power);
200 fmpz_mul(scale.d, scale.d, fmpz_poly_get_coeff_ptr(denpa.d, 0));
201
202
203 // Determine the size of the output by passing write = false.
204 const bool with_arghead = false;
205 bool write = false;
206 const uint64_t prev_size = 0;
207 const uint64_t out_size = return_rem ?
208 (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
209 rem.d, var_map, scale.d)
210 :
211 (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
212 div.d, var_map, scale.d)
213 ;
214 WORD* res = (WORD *)Malloc1(sizeof(WORD)*out_size, "flint::divrem_poly");
215
216
217 // Write out the result
218 write = true;
219 if ( return_rem ) {
220 (uint64_t)flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
221 rem.d, var_map, scale.d);
222 }
223 else {
224 (uint64_t)flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
225 div.d, var_map, scale.d);
226 }
227
228 return res;
229}
230/*
231 #] flint::divmod_poly :
232
233 #[ flint::factorize_mpoly :
234*/
235WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool with_arghead,
236 const bool is_fun_arg, const var_map_t &var_map) {
237
238 flint::mpoly_ctx ctx(var_map.size());
239 flint::mpoly arg(ctx.d), den(ctx.d), base(ctx.d);
240
241 flint::from_argument_mpoly(arg.d, den.d, argin, with_arghead, var_map, ctx.d);
242 // The denominator must be 1:
243 if ( fmpz_mpoly_is_one(den.d, ctx.d) != 1 ) {
244 MLOCK(ErrorMessageLock);
245 MesPrint("flint::factorize_mpoly error: den != 1");
246 MUNLOCK(ErrorMessageLock);
247 Terminate(-1);
248 }
249
250
251 // Now we can factor the mpoly:
252 flint::mpoly_factor arg_fac(ctx.d);
253 fmpz_mpoly_factor(arg_fac.d, arg.d, ctx.d);
254 const int64_t num_factors = fmpz_mpoly_factor_length(arg_fac.d, ctx.d);
255
256 flint::fmpz overall_constant;
257 fmpz_mpoly_factor_get_constant_fmpz(overall_constant.d, arg_fac.d, ctx.d);
258 // FORM should always have taken the overall constant out in the content. Thus this overall
259 // constant factor should be +- 1 here. Verify this:
260 if ( ! ( fmpz_equal_si(overall_constant.d, 1) || fmpz_equal_si(overall_constant.d, -1) ) ) {
261 MLOCK(ErrorMessageLock);
262 MesPrint("flint::factorize_mpoly error: overall constant factor != +-1");
263 MUNLOCK(ErrorMessageLock);
264 Terminate(-1);
265 }
266
267 // Construct the output. If argout is not NULL, we write the result there.
268 // Otherwise, allocate memory.
269 // The output is zero-terminated list of factors. If with_arghead, each has an arghead which
270 // contains its size. Otherwise, the factors are zero separated.
271
272 // We only need to determine the size of the output if we are allocating memory, but we need to
273 // loop through the factors to fix their signs anyway. Do both together in one loop:
274
275 // Initially 1, for the final trailing 0.
276 uint64_t output_size = 1;
277
278 // For finding the highest symbol, in FORM's lexicographic ordering
279 var_map_t var_map_inv;
280 for ( auto x: var_map ) {
281 var_map_inv[x.second] = x.first;
282 }
283
284 // Store whether we should flip the factor sign in the ouput:
285 vector<int32_t> base_sign(num_factors, 0);
286
287 for ( int64_t i = 0; i < num_factors; i++ ) {
288 fmpz_mpoly_factor_get_base(base.d, arg_fac.d, i, ctx.d);
289 const int64_t exponent = fmpz_mpoly_factor_get_exp_si(arg_fac.d, i, ctx.d);
290
291 // poly_factorize makes sure the highest power of the "highest symbol" (in FORM's
292 // lexicographic ordering) has a positive coefficient. Check this, update overall_constant
293 // of the factorization if necessary.
294 // Store the sign per factor, so that we can flip the signs in the output without re-checking
295 // the individual terms again.
296 uint32_t max_var = 0; // FORM symbols start at 20, 0 is a good initial value.
297 int32_t max_pow = -1;
298 vector<int64_t> base_term_exponents(var_map.size(), 0);
299
300 for ( int64_t j = 0; j < fmpz_mpoly_length(base.d, ctx.d); j++ ) {
301 fmpz_mpoly_get_term_exp_si((slong*)base_term_exponents.data(), base.d, (slong)j, ctx.d);
302
303 for ( size_t k = 0; k < var_map.size(); k++ ) {
304 if ( base_term_exponents[k] > 0 && ( var_map_inv.at(k) > max_var ||
305 ( var_map_inv.at(k) == max_var && base_term_exponents[k] > max_pow ) ) ) {
306
307 max_var = var_map_inv.at(k);
308 max_pow = base_term_exponents[k];
309 base_sign[i] = fmpz_sgn(fmpz_mpoly_term_coeff_ref(base.d, j, ctx.d));
310 }
311 }
312 }
313 // If this base's sign will be flipped an odd number of times, there is a contribution to
314 // the overall sign of the whole factorization:
315 if ( ( base_sign[i] == -1 ) && ( exponent % 2 == 1 ) ) {
316 fmpz_neg(overall_constant.d, overall_constant.d);
317 }
318
319 // Now determine the output size of the factor, if we are allocating the memory
320 if ( argout == NULL ) {
321 const bool write = false;
322 for ( int64_t j = 0; j < exponent; j++ ) {
323 output_size += (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead,
324 is_fun_arg, write, 0, base.d, var_map, ctx.d);
325 }
326 }
327 }
328 if ( fmpz_sgn(overall_constant.d) == -1 && argout == NULL ) {
329 // Add space for a fast-notation number or a normal-notation number and zero separator
330 output_size += with_arghead ? 2 : 4+1;
331 }
332
333 // Now make the allocation if necessary:
334 if ( argout == NULL ) {
335 argout = (WORD*)Malloc1(sizeof(WORD)*output_size, "flint::factorize_mpoly");
336 }
337
338
339 // And now comes the actual output:
340 WORD* old_argout = argout;
341
342 // If the overall sign is negative, first write a full-notation -1. It will be absorbed into the
343 // overall factor in the content by the caller.
344 if ( fmpz_sgn(overall_constant.d) == -1 ) {
345 if ( with_arghead ) {
346 // poly writes in fast notation in this case. Fast notation is expected by the caller, to
347 // properly merge it with the overall factor of the content.
348 *argout++ = -SNUMBER;
349 *argout++ = -1;
350 }
351 else {
352 *argout++ = 4; // term size
353 *argout++ = 1; // numerator
354 *argout++ = 1; // denominator
355 *argout++ = -3; // coeff size, negative number
356 *argout++ = 0; // factor separator
357 }
358 }
359
360 for ( int64_t i = 0; i < num_factors; i++ ) {
361 fmpz_mpoly_factor_get_base(base.d, arg_fac.d, i, ctx.d);
362 const int64_t exponent = fmpz_mpoly_factor_get_exp_si(arg_fac.d, i, ctx.d);
363
364 if ( base_sign[i] == -1 ) {
365 fmpz_mpoly_neg(base.d, base.d, ctx.d);
366 }
367
368 const bool write = true;
369 for ( int64_t j = 0; j < exponent; j++ ) {
370 argout += flint::to_argument_mpoly(BHEAD argout, with_arghead, is_fun_arg, write,
371 argout-old_argout, base.d, var_map, ctx.d);
372 }
373 }
374 // Final trailing zero to denote the end of the factors.
375 *argout++ = 0;
376
377 return old_argout;
378}
379/*
380 #] flint::factorize_mpoly :
381 #[ flint::factorize_poly :
382*/
383WORD* flint::factorize_poly(PHEAD const WORD *argin, WORD *argout, const bool with_arghead,
384 const bool is_fun_arg, const var_map_t &var_map) {
385
386 flint::poly arg, den;
387
388 flint::from_argument_poly(arg.d, den.d, argin, with_arghead);
389 // The denominator must be 1:
390 if ( fmpz_poly_is_one(den.d) != 1 ) {
391 MLOCK(ErrorMessageLock);
392 MesPrint("flint::factorize_poly error: den != 1");
393 MUNLOCK(ErrorMessageLock);
394 Terminate(-1);
395 }
396
397
398 // Now we can factor the poly:
399 flint::poly_factor arg_fac;
400 fmpz_poly_factor(arg_fac.d, arg.d);
401 // fmpz_poly_factor_t lacks some convenience functions which fmpz_mpoly_factor_t has.
402 // I have worked out how to get the factors by looking at how fmpz_poly_factor_print works.
403 const long num_factors = (arg_fac.d)->num;
404
405
406 // Construct the output. If argout is not NULL, we write the result there.
407 // Otherwise, allocate memory.
408 // The output is zero-terminated list of factors. If with_arghead, each has an arghead which
409 // contains its size. Otherwise, the factors are zero separated.
410 if ( argout == NULL ) {
411 // First we need to determine the size of the output. This is the same procedure as the
412 // loop below, but we don't write anything in to_argument_poly (write arg: false).
413 // Initially 1, for the final trailing 0.
414 uint64_t output_size = 1;
415
416 for ( long i = 0; i < num_factors; i++ ) {
417 fmpz_poly_struct* base = (arg_fac.d)->p + i;
418
419 const long exponent = (arg_fac.d)->exp[i];
420
421 const bool write = false;
422 for ( long j = 0; j < exponent; j++ ) {
423 output_size += (uint64_t)flint::to_argument_poly(BHEAD NULL, with_arghead,
424 is_fun_arg, write, 0, base, var_map);
425 }
426 }
427
428 argout = (WORD*)Malloc1(sizeof(WORD)*output_size, "flint::factorize_poly");
429 }
430
431 WORD* old_argout = argout;
432
433 for ( long i = 0; i < num_factors; i++ ) {
434 fmpz_poly_struct* base = (arg_fac.d)->p + i;
435
436 const long exponent = (arg_fac.d)->exp[i];
437
438 const bool write = true;
439 for ( long j = 0; j < exponent; j++ ) {
440 argout += flint::to_argument_poly(BHEAD argout, with_arghead, is_fun_arg, write,
441 argout-old_argout, base, var_map);
442 }
443 }
444 *argout = 0;
445
446 return old_argout;
447}
448/*
449 #] flint::factorize_poly :
450
451 #[ flint::form_sort :
452*/
453// Sort terms using form's sorting routines. Uses a custom (faster) compare routine, since here
454// only symbols can appear.
455// This is a modified poly_sort from polywrap.cc.
456void flint::form_sort(PHEAD WORD *terms) {
457
458 if ( terms[0] < 0 ) {
459 // Fast notation, there is nothing to do
460 return;
461 }
462
463 const WORD oldsorttype = AR.SortType;
464 AR.SortType = SORTHIGHFIRST;
465
466 const WORD in_size = terms[0] - ARGHEAD;
467 WORD out_size;
468
469 if ( NewSort(BHEAD0) ) {
470 Terminate(-1);
471 }
472 AR.CompareRoutine = (COMPAREDUMMY)(&CompareSymbols);
473
474 // Make sure the symbols are in the right order within the terms
475 for ( WORD i = ARGHEAD; i < terms[0]; i += terms[i] ) {
476 if ( SymbolNormalize(terms+i) < 0 || StoreTerm(BHEAD terms+i) ) {
477 AR.SortType = oldsorttype;
478 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
479 LowerSortLevel();
480 Terminate(-1);
481 }
482 }
483
484 if ( ( out_size = EndSort(BHEAD terms+ARGHEAD, 1) ) < 0 ) {
485 AR.SortType = oldsorttype;
486 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
487 Terminate(-1);
488 }
489
490 // Check the final size
491 if ( in_size != out_size ) {
492 MLOCK(ErrorMessageLock);
493 MesPrint("flint::form_sort: error: unexpected sorted arg length change %d->%d", in_size,
494 out_size);
495 MUNLOCK(ErrorMessageLock);
496 Terminate(-1);
497 }
498
499 AR.SortType = oldsorttype;
500 AR.CompareRoutine = (COMPAREDUMMY)(&Compare1);
501 terms[1] = 0; // set dirty flag to zero
502}
503/*
504 #] flint::form_sort :
505
506 #[ flint::from_argument_mpoly :
507*/
508// Convert a FORM argument (or 0-terminated list of terms: with_arghead == false) to a
509// (multi-variate) fmpz_mpoly_t poly. The "denominator" is return in denpoly, and contains the
510// overall negative-power numeric and symbolic factor.
511uint64_t flint::from_argument_mpoly(fmpz_mpoly_t poly, fmpz_mpoly_t denpoly, const WORD *args,
512 const bool with_arghead, const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx) {
513
514 // Some callers re-use their poly, denpoly to avoid calling init/clear unnecessarily.
515 // Make sure they are 0 to begin.
516 fmpz_mpoly_set_si(poly, 0, ctx);
517 fmpz_mpoly_set_si(denpoly, 0, ctx);
518
519 // First check for "fast notation" arguments:
520 if ( *args == -SNUMBER ) {
521 fmpz_mpoly_set_si(poly, *(args+1), ctx);
522 fmpz_mpoly_set_si(denpoly, 1, ctx);
523 return 2;
524 }
525
526 if ( *args == -SYMBOL ) {
527 // A "fast notation" SYMBOL has a power and coefficient of 1:
528 vector<uint64_t> exponents(var_map.size(), 0);
529 exponents[var_map.at(*(args+1))] = 1;
530
531 fmpz_mpoly_set_coeff_ui_ui(poly, (ulong)1, (ulong*)exponents.data(), ctx);
532 fmpz_mpoly_set_ui(denpoly, (ulong)1, ctx);
533 return 2;
534 }
535
536
537 // Now we can iterate through the terms of the argument. If we have
538 // an ARGHEAD, we already know where to terminate. Otherwise we'll have
539 // to loop until the terminating 0.
540 const WORD* arg_stop = with_arghead ? args+args[0] : (WORD*)UINT64_MAX;
541 uint64_t arg_size = 0;
542 if ( with_arghead ) {
543 arg_size = args[0];
544 args += ARGHEAD;
545 }
546
547
548 // Search for numerical or symbol denominators to create "denpoly".
549 flint::fmpz den_coeff, tmp;
550 fmpz_set_si(den_coeff.d, 1);
551 vector<uint64_t> neg_exponents(var_map.size(), 0);
552
553 for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
554 const WORD* term_stop = term+term[0];
555 const WORD coeff_size = (term_stop)[-1];
556 const WORD* symbol_stop = term_stop - ABS(coeff_size);
557 const WORD* t = term;
558
559 t++;
560 if ( t == symbol_stop ) {
561 // Just a number, no symbols
562 }
563 else {
564 t++; // this entry is SYMBOL
565 t++; // this entry just has the size of the symbol array, but we can use symbol_stop
566
567 for ( const WORD* s = t; s < symbol_stop; s += 2 ) {
568 if ( *(s+1) < 0 ) {
569 neg_exponents[var_map.at(*s)] =
570 MaX(neg_exponents[var_map.at(*s)], (uint64_t)(-(*(s+1))) );
571 }
572 }
573 }
574
575 // Now check for a denominator in the coefficient:
576 if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
577 flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
578 // Record the LCM of the coefficient denominators:
579 fmpz_lcm(den_coeff.d, den_coeff.d, tmp.d);
580 }
581
582 if ( !with_arghead && *term_stop == 0 ) {
583 // + 1 for the terminating 0
584 arg_size = term_stop - args + 1;
585 break;
586 }
587
588 }
589 // Assemble denpoly.
590 fmpz_mpoly_set_coeff_fmpz_ui(denpoly, den_coeff.d, (ulong*)neg_exponents.data(), ctx);
591
592
593 // For the term coefficients
594 flint::fmpz coeff;
595
596 for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
597 const WORD* term_stop = term+term[0];
598 const WORD coeff_size = (term_stop)[-1];
599 const WORD* symbol_stop = term_stop - ABS(coeff_size);
600 const WORD* t = term;
601
602 vector<uint64_t> exponents(var_map.size(), 0);
603
604 t++; // skip over the total size entry
605 if ( t == symbol_stop ) {
606 // Just a number, no symbols
607 }
608 else {
609 t++; // this entry is SYMBOL
610 t++; // this entry just has the size of the symbol array, but we can use symbol_stop
611 for ( const WORD* s = t; s < symbol_stop; s += 2 ) {
612 exponents[var_map.at(*s)] = *(s+1);
613 }
614 }
615 // Now read the coefficient
616 flint::fmpz_set_form(coeff.d, (UWORD*)symbol_stop, coeff_size/2);
617
618 // Multiply by denominator LCM
619 fmpz_mul(coeff.d, coeff.d, den_coeff.d);
620
621 // Shift by neg_exponents
622 for ( size_t i = 0; i < var_map.size(); i++ ) {
623 exponents[i] += neg_exponents[i];
624 }
625
626 // Read the denominator if there is one, and divide it out of the coefficient
627 if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
628 flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
629 // By construction, this is an exact division
630 fmpz_divexact(coeff.d, coeff.d, tmp.d);
631 }
632
633 // Push the term to the mpoly, remember to sort when finished! This is much faster than using
634 // fmpz_mpoly_set_coeff_fmpz_ui when the terms arrive in the "wrong order".
635 fmpz_mpoly_push_term_fmpz_ui(poly, coeff.d, (ulong*)exponents.data(), ctx);
636
637 if ( !with_arghead && *term_stop == 0 ) {
638 break;
639 }
640
641 }
642 // And now sort the mpoly
643 fmpz_mpoly_sort_terms(poly, ctx);
644
645 return arg_size;
646}
647/*
648 #] flint::from_argument_mpoly :
649 #[ flint::from_argument_poly :
650*/
651// Convert a FORM argument (or 0-terminated list of terms: with_arghead == false) to a
652// (uni-variate) fmpz_poly_t poly. The "denominator" is return in denpoly, and contains the
653// overall negative-power numeric and symbolic factor.
654uint64_t flint::from_argument_poly(fmpz_poly_t poly, fmpz_poly_t denpoly, const WORD *args,
655 const bool with_arghead) {
656
657 // Some callers re-use their poly, denpoly to avoid calling init/clear unnecessarily.
658 // Make sure they are 0 to begin.
659 fmpz_poly_set_si(poly, 0);
660 fmpz_poly_set_si(denpoly, 0);
661
662 // First check for "fast notation" arguments:
663 if ( *args == -SNUMBER ) {
664 fmpz_poly_set_si(poly, *(args+1));
665 fmpz_poly_set_si(denpoly, 1);
666 return 2;
667 }
668
669 if ( *args == -SYMBOL ) {
670 // A "fast notation" SYMBOL has a power and coefficient of 1:
671 fmpz_poly_set_coeff_si(poly, 1, 1);
672 fmpz_poly_set_si(denpoly, 1);
673 return 2;
674 }
675
676 // Now we can iterate through the terms of the argument. If we have
677 // an ARGHEAD, we already know where to terminate. Otherwise we'll have
678 // to loop until the terminating 0.
679 const WORD* arg_stop = with_arghead ? args+args[0] : (WORD*)UINT64_MAX;
680 uint64_t arg_size = 0;
681 if ( with_arghead ) {
682 arg_size = args[0];
683 args += ARGHEAD;
684 }
685
686 // Search for numerical or symbol denominators to create "denpoly".
687 flint::fmpz den_coeff, tmp;
688 fmpz_set_si(den_coeff.d, 1);
689 uint64_t neg_exponent = 0;
690
691 for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
692
693 const WORD* term_stop = term+term[0];
694 const WORD coeff_size = (term_stop)[-1];
695 const WORD* symbol_stop = term_stop - ABS(coeff_size);
696 const WORD* t = term;
697
698 t++; // skip over the total size entry
699 if ( t == symbol_stop ) {
700 // Just a number, no symbols
701 }
702 else {
703 t++; // this entry is SYMBOL
704 t++; // this entry is the size of the symbol array
705 t++; // this is the first (and only) symbol code
706 if ( *t < 0 ) {
707 neg_exponent = MaX(neg_exponent, (uint64_t)(-(*t)) );
708 }
709 }
710
711 // Now check for a denominator in the coefficient:
712 if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
713 flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
714 // Record the LCM of the coefficient denominators:
715 fmpz_lcm(den_coeff.d, den_coeff.d, tmp.d);
716 }
717
718 if ( *term_stop == 0 ) {
719 // + 1 for the terminating 0
720 arg_size = term_stop - args + 1;
721 break;
722 }
723 }
724 // Assemble denpoly.
725 fmpz_poly_set_coeff_fmpz(denpoly, neg_exponent, den_coeff.d);
726
727
728 // For the term coefficients
729 flint::fmpz coeff;
730
731 for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
732
733 const WORD* term_stop = term+term[0];
734 const WORD coeff_size = (term_stop)[-1];
735 const WORD* symbol_stop = term_stop - ABS(coeff_size);
736 const WORD* t = term;
737
738 uint64_t exponent = 0;
739
740 t++; // skip over the total size entry
741 if ( t == symbol_stop ) {
742 // Just a number, no symbols
743 }
744 else {
745 t++; // this entry is SYMBOL
746 t++; // this entry is the size of the symbol array
747 t++; // this is the first (and only) symbol code
748 exponent = *t++;
749 }
750
751 // Now read the coefficient
752 flint::fmpz_set_form(coeff.d, (UWORD*)symbol_stop, coeff_size/2);
753
754 // Multiply by denominator LCM
755 fmpz_mul(coeff.d, coeff.d, den_coeff.d);
756
757 // Shift by neg_exponent
758 exponent += neg_exponent;
759
760 // Read the denominator if there is one, and divide it out of the coefficient
761 if ( *(symbol_stop+ABS(coeff_size/2)) != 1 ) {
762 flint::fmpz_set_form(tmp.d, (UWORD*)(symbol_stop+ABS(coeff_size/2)), ABS(coeff_size/2));
763 // By construction, this is an exact division
764 fmpz_divexact(coeff.d, coeff.d, tmp.d);
765 }
766
767 // Add the term to the poly
768 fmpz_poly_set_coeff_fmpz(poly, exponent, coeff.d);
769
770 if ( *term_stop == 0 ) {
771 break;
772 }
773 }
774
775 return arg_size;
776}
777/*
778 #] flint::from_argument_poly :
779
780 #[ flint::fmpz_get_form :
781*/
782// Write FORM's long integer representation of an fmpz at a, and put the number of WORDs at na.
783// na carries the sign of the integer.
784WORD flint::fmpz_get_form(fmpz_t z, WORD *a) {
785
786 WORD na = 0;
787 const int32_t sgn = fmpz_sgn(z);
788 if ( sgn == -1 ) {
789 fmpz_neg(z, z);
790 }
791 const int64_t nlimbs = fmpz_size(z);
792
793 // This works but is UB?
794 //fmpz_get_ui_array(reinterpret_cast<uint64_t*>(a), nlimbs, z);
795
796 // Use fixed-size functions to get limb data where possible. These probably cover most real
797 // cases.
798 if ( nlimbs == 1 ) {
799 const uint64_t limb = fmpz_get_ui(z);
800 a[0] = (WORD)(limb & 0xFFFFFFFF);
801 na++;
802 a[1] = (WORD)(limb >> BITSINWORD);
803 if ( a[1] != 0 ) {
804 na++;
805 }
806 }
807 else if ( nlimbs == 2 ) {
808 uint64_t limb_hi = 0, limb_lo = 0;
809 fmpz_get_uiui((ulong*)&limb_hi, (ulong*)&limb_lo, z);
810 a[0] = (WORD)(limb_lo & 0xFFFFFFFF);
811 na++;
812 a[1] = (WORD)(limb_lo >> BITSINWORD);
813 na++;
814 a[2] = (WORD)(limb_hi & 0xFFFFFFFF);
815 na++;
816 a[3] = (WORD)(limb_hi >> BITSINWORD);
817 if ( a[3] != 0 ) {
818 na++;
819 }
820 }
821 else {
822 vector<uint64_t> limb_data(nlimbs, 0);
823 fmpz_get_ui_array((ulong*)limb_data.data(), (slong)nlimbs, z);
824 for ( long i = 0; i < nlimbs; i++ ) {
825 a[2*i] = (WORD)(limb_data[i] & 0xFFFFFFFF);
826 na++;
827 a[2*i+1] = (WORD)(limb_data[i] >> BITSINWORD);
828 if ( a[2*i+1] != 0 || i < (nlimbs-1) ) {
829 // The final limb might fit in a single 32bit WORD. Only
830 // increment na if the final WORD is non zero.
831 na++;
832 }
833 }
834 }
835
836 // And now put the sign in the number of limbs
837 if ( sgn == -1 ) {
838 na = -na;
839 }
840
841 return na;
842}
843/*
844 #] flint::fmpz_get_form :
845 #[ flint::fmpz_set_form :
846*/
847// Create an fmpz directly from FORM's long integer representation. fmpz uses 64bit unsigned limbs,
848// but FORM uses 32bit UWORDs on 64bit architectures so we can't use fmpz_set_ui_array directly.
849void flint::fmpz_set_form(fmpz_t z, UWORD *a, WORD na) {
850
851 if ( na == 0 ) {
852 fmpz_zero(z);
853 return;
854 }
855
856 // Negative na represenents a negative number
857 int32_t sgn = 1;
858 if ( na < 0 ) {
859 sgn = -1;
860 na = -na;
861 }
862
863 // Remove padding. FORM stores numerators and denominators with equal numbers of limbs but we
864 // don't need to do this within the fmpz. It is not necessary to do this really, the fmpz
865 // doesn't add zero limbs unnecessarily, but we might be able to avoid creating the limb_data
866 // array below.
867 while ( a[na-1] == 0 ) {
868 na--;
869 }
870
871 // If the number fits in fixed-size fmpz_set functions, we don't need to use additional memory
872 // to convert to uint64_t. These probably cover most real cases.
873 if ( na == 1 ) {
874 fmpz_set_ui(z, (uint64_t)a[0]);
875 }
876 else if ( na == 2 ) {
877 fmpz_set_ui(z, (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
878 }
879 else if ( na == 3 ) {
880 fmpz_set_uiui(z, (uint64_t)a[2], (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
881 }
882 else if ( na == 4 ) {
883 fmpz_set_uiui(z, (((uint64_t)a[3])<<BITSINWORD) + (uint64_t)a[2],
884 (((uint64_t)a[1])<<BITSINWORD) + (uint64_t)a[0]);
885 }
886 else {
887 const int32_t nlimbs = (na+1)/2;
888 vector<uint64_t> limb_data(nlimbs, 0);
889 for ( int32_t i = 0; i < nlimbs; i++ ) {
890 if ( 2*i+1 <= na-1 ) {
891 limb_data[i] = (uint64_t)a[2*i] + (((uint64_t)a[2*i+1])<<BITSINWORD);
892 }
893 else {
894 limb_data[i] = (uint64_t)a[2*i];
895 }
896 }
897 fmpz_set_ui_array(z, (ulong*)limb_data.data(), nlimbs);
898 }
899
900 // Finally set the sign.
901 if ( sgn == -1 ) {
902 fmpz_neg(z, z);
903 }
904
905 return;
906}
907/*
908 #] flint::fmpz_set_form :
909
910 #[ flint::gcd_mpoly :
911*/
912// Return a pointer to a buffer containing the GCD of the 0-terminated term lists at a and b.
913// If must_fit_term, this should be a TermMalloc buffer. Otherwise Malloc1 the buffer.
914// For multi-variate cases.
915WORD* flint::gcd_mpoly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_term,
916 const var_map_t &var_map) {
917
918 flint::mpoly_ctx ctx(var_map.size());
919 flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d), gcd(ctx.d);
920
921 flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
922 flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
923
924 // denpa, denpb must be 1:
925 if ( fmpz_mpoly_is_one(denpa.d, ctx.d) != 1 ) {
926 MLOCK(ErrorMessageLock);
927 MesPrint("flint::gcd_mpoly: error: denpa != 1");
928 MUNLOCK(ErrorMessageLock);
929 Terminate(-1);
930 }
931 if ( fmpz_mpoly_is_one(denpb.d, ctx.d) != 1 ) {
932 MLOCK(ErrorMessageLock);
933 MesPrint("flint::gcd_mpoly: error: denpb != 1");
934 MUNLOCK(ErrorMessageLock);
935 Terminate(-1);
936 }
937
938 // poly returns pa if pa == pb, regardless of the lcoeff sign
939 if ( fmpz_mpoly_equal(pa.d, pb.d, ctx.d) ) {
940 fmpz_mpoly_set(gcd.d, pa.d, ctx.d);
941 }
942 else {
943 // We need some gymnastics to have the same sign conventions as the poly class. It takes the
944 // integer or univar content out of a,b, with the convention that the content sign matches
945 // the lcoeff sign. Since FORM has already taken out the content, we are left with +-1. In
946 // Flint, the content always has a positive sign so here we should find +1. Check this:
947 flint::mpoly tmp(ctx.d);
948 fmpz_mpoly_term_content(tmp.d, pa.d, ctx.d);
949 if ( fmpz_mpoly_is_one(tmp.d, ctx.d) != 1 ) {
950 MLOCK(ErrorMessageLock);
951 MesPrint("flint::gcd_mpoly: error: content of 1st arg != 1");
952 MUNLOCK(ErrorMessageLock);
953 Terminate(-1);
954 }
955 fmpz_mpoly_term_content(tmp.d, pb.d, ctx.d);
956 if ( fmpz_mpoly_is_one(tmp.d, ctx.d) != 1 ) {
957 MLOCK(ErrorMessageLock);
958 MesPrint("flint::gcd_mpoly: error: content of 2nd arg != 1");
959 MUNLOCK(ErrorMessageLock);
960 Terminate(-1);
961 }
962
963 // The poly class now divides the content out of a,b so that they have a positive lcoeff.
964 // Then it multiplies the final gcd (which is given a positive lcoeff also) by
965 // gcd(cont a, cont b). There it has gcd(1,1) = gcd(-1,1) = gcd(1,-1) = 1, and
966 // gcd(-1,-1) = -1 (because of the pa==pb early return). So: if both input polys have a
967 // negative lcoeff, we will flip the sign in the final result.
968 bool flip_sign = 0;
969 if ( ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(pa.d, 0, ctx.d)) == -1 ) &&
970 ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(pb.d, 0, ctx.d)) == -1 ) ) {
971 flip_sign = 1;
972 }
973
974 fmpz_mpoly_gcd(gcd.d, pa.d, pb.d, ctx.d);
975 if ( flip_sign ) {
976 fmpz_mpoly_neg(gcd.d, gcd.d, ctx.d);
977 }
978 }
979
980 // This is freed by the caller
981 WORD *res;
982 if ( must_fit_term ) {
983 res = TermMalloc("flint::gcd_mpoly");
984 }
985 else {
986 // Determine the size of the GCD by passing write = false.
987 const bool with_arghead = false;
988 const bool write = false;
989 const uint64_t prev_size = 0;
990 const uint64_t gcd_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
991 with_arghead, must_fit_term, write, prev_size, gcd.d, var_map, ctx.d);
992
993 res = (WORD *)Malloc1(sizeof(WORD)*gcd_size, "flint::gcd_mpoly");
994 }
995
996 const bool with_arghead = false;
997 const bool write = true;
998 const uint64_t prev_size = 0;
999 flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, gcd.d,
1000 var_map, ctx.d);
1001
1002 return res;
1003}
1004/*
1005 #] flint::gcd_mpoly :
1006 #[ flint::gcd_poly :
1007*/
1008// Return a pointer to a buffer containing the GCD of the 0-terminated term lists at a and b.
1009// If must_fit_term, this should be a TermMalloc buffer. Otherwise Malloc1 the buffer.
1010// For uni-variate cases.
1011WORD* flint::gcd_poly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_term,
1012 const var_map_t &var_map) {
1013
1014 flint::poly pa, pb, denpa, denpb, gcd;
1015
1016 flint::from_argument_poly(pa.d, denpa.d, a, false);
1017 flint::from_argument_poly(pb.d, denpb.d, b, false);
1018
1019 // denpa, denpb must be 1:
1020 if ( fmpz_poly_is_one(denpa.d) != 1 ) {
1021 MLOCK(ErrorMessageLock);
1022 MesPrint("flint::gcd_poly: error: denpa != 1");
1023 MUNLOCK(ErrorMessageLock);
1024 Terminate(-1);
1025 }
1026 if ( fmpz_poly_is_one(denpb.d) != 1 ) {
1027 MLOCK(ErrorMessageLock);
1028 MesPrint("flint::gcd_poly: error: denpb != 1");
1029 MUNLOCK(ErrorMessageLock);
1030 Terminate(-1);
1031 }
1032
1033 // poly returns pa if pa == pb, regardless of the lcoeff sign
1034 if ( fmpz_poly_equal(pa.d, pb.d) ) {
1035 fmpz_poly_set(gcd.d, pa.d);
1036 }
1037 else {
1038 // Here, we don't have to make any sign flips like the mpoly case, because poly's
1039 // integer_gcd(1,1) = integer_gcd(-1,1) = integer_gcd(1,-1) = integer_gcd(-1,-1) = +1.
1040 // Still, verify that the content is 1:
1041 flint::fmpz tmp;
1042 fmpz_poly_content(tmp.d, pa.d);
1043 if ( fmpz_is_one(tmp.d) != 1 ) {
1044 MLOCK(ErrorMessageLock);
1045 MesPrint("flint::gcd_poly: error: content of 1st arg != 1");
1046 MUNLOCK(ErrorMessageLock);
1047 Terminate(-1);
1048 }
1049 fmpz_poly_content(tmp.d, pb.d);
1050 if ( fmpz_is_one(tmp.d) != 1 ) {
1051 MLOCK(ErrorMessageLock);
1052 MesPrint("flint::gcd_poly: error: content of 2nd arg != 1");
1053 MUNLOCK(ErrorMessageLock);
1054 Terminate(-1);
1055 }
1056
1057 fmpz_poly_gcd(gcd.d, pa.d, pb.d);
1058 }
1059
1060 // This is freed by the caller
1061 WORD *res;
1062 if ( must_fit_term ) {
1063 res = TermMalloc("flint::gcd_poly");
1064 }
1065 else {
1066 // Determine the size of the GCD by passing write = false.
1067 const bool with_arghead = false;
1068 const bool write = false;
1069 const uint64_t prev_size = 0;
1070 const uint64_t gcd_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
1071 with_arghead, must_fit_term, write, prev_size, gcd.d, var_map);
1072
1073 res = (WORD *)Malloc1(sizeof(WORD)*gcd_size, "flint::gcd_poly");
1074 }
1075
1076 const bool with_arghead = false;
1077 const uint64_t prev_size = 0;
1078 const bool write = true;
1079 flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, gcd.d,
1080 var_map);
1081
1082 return res;
1083}
1084/*
1085 #] flint::gcd_poly :
1086
1087 #[ flint::get_variables :
1088*/
1089// Get the list of symbols which appear in the vector of expressions. These are polyratfun
1090// numerators, denominators or expressions from calls to gcd_ etc. Return this list as a map
1091// between indices and symbol codes.
1092// TODO FACTORSYMBOL last?
1093flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with_arghead,
1094 const bool sort_vars) {
1095
1096 int32_t num_vars = 0;
1097 // We count the total number of terms to determine "density".
1098 uint32_t num_terms = 0;
1099 // To be used if we sort by highest degree, as the poly code does.
1100 vector<int32_t> degrees;
1101 var_map_t var_map;
1102
1103 // extract all variables
1104 for ( size_t ei = 0; ei < es.size(); ei++ ) {
1105 WORD *e = es[ei];
1106
1107 // fast notation
1108 if ( *e == -SNUMBER ) {
1109 num_terms++;
1110 }
1111 else if ( *e == -SYMBOL ) {
1112 num_terms++;
1113 if ( !var_map.count(e[1]) ) {
1114 var_map[e[1]] = num_vars++;
1115 degrees.push_back(1);
1116 }
1117 }
1118 // JD: Here we need to check for non-symbol/number terms in fast notation.
1119 else if ( *e < 0 ) {
1120 MLOCK(ErrorMessageLock);
1121 MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
1122 MUNLOCK(ErrorMessageLock);
1123 Terminate(1);
1124 }
1125 else {
1126 for ( WORD i = with_arghead ? ARGHEAD:0; with_arghead ? i < e[0]:e[i] != 0; i += e[i] ) {
1127 num_terms++;
1128 if ( i+1 < i+e[i]-ABS(e[i+e[i]-1]) && e[i+1] != SYMBOL ) {
1129 MLOCK(ErrorMessageLock);
1130 MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
1131 MUNLOCK(ErrorMessageLock);
1132 Terminate(1);
1133 }
1134
1135 for ( WORD j = i+3; j<i+e[i]-ABS(e[i+e[i]-1]); j += 2 ) {
1136 if ( !var_map.count(e[j]) ) {
1137 var_map[e[j]] = num_vars++;
1138 degrees.push_back(e[j+1]);
1139 }
1140 else {
1141 degrees[var_map[e[j]]] = MaX(degrees[var_map[e[j]]], e[j+1]);
1142 }
1143 }
1144 }
1145 }
1146 }
1147
1148 if ( sort_vars ) {
1149 // bubble sort variables in decreasing order of degree
1150 // (this seems better for factorization)
1151 for ( size_t i = 0; i < var_map.size(); i++ ) {
1152 for ( size_t j = 0; j+1 < var_map.size(); j++ ) {
1153 if ( degrees[j] < degrees[j+1] ) {
1154 swap(degrees[j], degrees[j+1]);
1155
1156 // Find the map keys associated with the values we want to swap
1157 uint32_t j0 = 0;
1158 uint32_t j1 = 0;
1159 for ( auto x: var_map ) {
1160 if ( x.second == j ) {
1161 j0 = x.first;
1162 }
1163 else if ( x.second == j+1 ) {
1164 j1 = x.first;
1165 }
1166 }
1167 swap(var_map.at(j0), var_map.at(j1));
1168 }
1169 }
1170 }
1171 }
1172 // Otherwise, sort lexicographically in FORM's ordering
1173 else {
1174 for ( size_t i = 0; i < var_map.size(); i++ ) {
1175 for ( size_t j = 0; j+1 < var_map.size(); j++ ) {
1176 uint32_t j0 = 0;
1177 uint32_t j1 = 0;
1178 for ( auto x: var_map ) {
1179 if ( x.second == j ) {
1180 j0 = x.first;
1181 }
1182 else if ( x.second == j+1 ) {
1183 j1 = x.first;
1184 }
1185 }
1186 if ( j0 > j1 ) {
1187 swap(var_map.at(j0), var_map.at(j1));
1188 }
1189 }
1190 }
1191 }
1192
1193 if ( var_map.size() == 1 ) {
1194 // In the univariate case, if the polynomials are sufficiently sparse force the use of the
1195 // multivariate routines, which use a sparse representation, by adding a dummy map entry.
1196 if ( (float)num_terms <= UNIVARIATE_DENSITY_THR * (float)degrees[0] ) {
1197 // -1 will never be a symbol code. Built-in symbols from 0 to 19, and 20 is the first
1198 // user symbol.
1199 var_map[-1] = num_vars;
1200 }
1201 }
1202
1203 return var_map;
1204}
1205/*
1206 #] flint::get_variables :
1207
1208 #[ flint::inverse_poly :
1209*/
1210WORD* flint::inverse_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
1211
1212 flint::poly pa, pb, denpa, denpb;
1213
1214 flint::from_argument_poly(pa.d, denpa.d, a, false);
1215 flint::from_argument_poly(pb.d, denpb.d, b, false);
1216
1217 // fmpz_poly_xgcd is undefined if the content of pa and pb are not 1. Take the content out.
1218 // fmpz_poly_content gives the non-negative content and fmpz_poly_primitive normalizes to a
1219 // non-negative lcoeff, so we need to add the sign to the content if the polys have a negative
1220 // lcoeff. We don't need to keep the content of pb, it is a numerical multiple of the modulus.
1221 flint::fmpz content_a, resultant;
1222 flint::poly inverse, tmp;
1223 fmpz_poly_content(content_a.d, pa.d);
1224 if ( fmpz_sgn(fmpz_poly_lead(pa.d)) == -1 ) {
1225 fmpz_neg(content_a.d, content_a.d);
1226 }
1227 fmpz_poly_primitive_part(pa.d, pa.d);
1228 fmpz_poly_primitive_part(pb.d, pb.d);
1229
1230 // Special cases:
1231 // Possibly strange that we give 1 for inverse_(x1,1) but here we take MMA's convention.
1232 if ( fmpz_poly_is_one(pa.d) && fmpz_poly_is_one(pb.d) ) {
1233 fmpz_poly_one(inverse.d);
1234 fmpz_one(resultant.d);
1235 }
1236 else if ( fmpz_poly_is_one(pb.d) ) {
1237 fmpz_poly_zero(inverse.d);
1238 fmpz_one(resultant.d);
1239 }
1240 else {
1241 // Now use the extended Euclidean algorithm to find inverse, resultant, tmp of the Bezout
1242 // identity: inverse*pa + tmp*pb = resultant. Then inverse/resultant is the multiplicative
1243 // inverse of pa mod pb. We'll divide by resultant in the denscale argument of to_argument_poly.
1244 fmpz_poly_xgcd(resultant.d, inverse.d, tmp.d, pa.d, pb.d);
1245
1246 // If the resultant is zero, the inverse does not exist:
1247 if ( fmpz_is_zero(resultant.d) ) {
1248 MLOCK(ErrorMessageLock);
1249 MesPrint("flint::inverse_poly error: inverse does not exist");
1250 MUNLOCK(ErrorMessageLock);
1251 Terminate(-1);
1252 }
1253 }
1254
1255 // Multiply inverse by denpa. denpb is a numerical multiple of the modulus, so doesn't matter.
1256 // We also need to divide by content_a, which we do in the denscale argument of to_argument_poly
1257 // by multiplying resultant by content_a here.
1258 fmpz_poly_mul(inverse.d, inverse.d, denpa.d);
1259 fmpz_mul(resultant.d, resultant.d, content_a.d);
1260
1261
1262 WORD* res;
1263 // First determine the result size, and malloc. The result should have no arghead. Here we use
1264 // the "scale" argument of to_argument_poly since resultant might not be 1.
1265 const bool with_arghead = false;
1266 bool write = false;
1267 const bool must_fit_term = false;
1268 const uint64_t prev_size = 0;
1269 const uint64_t res_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
1270 with_arghead, must_fit_term, write, prev_size, inverse.d, var_map, resultant.d);
1271 res = (WORD*)Malloc1(sizeof(WORD)*res_size, "flint::inverse_poly");
1272
1273 write = true;
1274 flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, inverse.d,
1275 var_map, resultant.d);
1276
1277 return res;
1278}
1279/*
1280 #] flint::inverse_poly :
1281
1282 #[ flint::mul_mpoly :
1283*/
1284// Return a pointer to a buffer containing the product of the 0-terminated term lists at a and b.
1285// For multi-variate cases.
1286WORD* flint::mul_mpoly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
1287
1288 flint::mpoly_ctx ctx(var_map.size());
1289 flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
1290
1291 flint::from_argument_mpoly(pa.d, denpa.d, a, false, var_map, ctx.d);
1292 flint::from_argument_mpoly(pb.d, denpb.d, b, false, var_map, ctx.d);
1293
1294 // denpa, denpb must be integers. Negative symbol powers have been converted to extra symbols.
1295 if ( fmpz_mpoly_is_fmpz(denpa.d, ctx.d) != 1 ) {
1296 MLOCK(ErrorMessageLock);
1297 MesPrint("flint::mul_mpoly: error: denpa is non-constant");
1298 MUNLOCK(ErrorMessageLock);
1299 Terminate(-1);
1300 }
1301 if ( fmpz_mpoly_is_fmpz(denpb.d, ctx.d) != 1 ) {
1302 MLOCK(ErrorMessageLock);
1303 MesPrint("flint::mul_mpoly: error: denpb is non-constant");
1304 MUNLOCK(ErrorMessageLock);
1305 Terminate(-1);
1306 }
1307
1308 // Multiply numerators, store result in pa
1309 fmpz_mpoly_mul(pa.d, pa.d, pb.d, ctx.d);
1310 // Multiply denominators, store result in denpa, and convert to an fmpz:
1311 fmpz_mpoly_mul(denpa.d, denpa.d, denpb.d, ctx.d);
1312 flint::fmpz den;
1313 fmpz_mpoly_get_fmpz(den.d, denpa.d, ctx.d);
1314
1315 WORD* res;
1316 // First determine the result size, and malloc. The result should have no arghead. Here we use
1317 // the "scale" argument of to_argument_mpoly since den might not be 1.
1318 const bool with_arghead = false;
1319 bool write = false;
1320 const bool must_fit_term = false;
1321 const uint64_t prev_size = 0;
1322 const uint64_t mul_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
1323 with_arghead, must_fit_term, write, prev_size, pa.d, var_map, ctx.d, den.d);
1324 res = (WORD*)Malloc1(sizeof(WORD)*mul_size, "flint::mul_mpoly");
1325
1326 write = true;
1327 flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, pa.d,
1328 var_map, ctx.d, den.d);
1329
1330 return res;
1331}
1332/*
1333 #] flint::mul_mpoly :
1334 #[ flint::mul_poly :
1335*/
1336// Return a pointer to a buffer containing the product of the 0-terminated term lists at a and b.
1337// For uni-variate cases.
1338WORD* flint::mul_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
1339
1340 flint::poly pa, pb, denpa, denpb;
1341
1342 flint::from_argument_poly(pa.d, denpa.d, a, false);
1343 flint::from_argument_poly(pb.d, denpb.d, b, false);
1344
1345 // denpa, denpb must be integers. Negative symbol powers have been converted to extra symbols.
1346 if ( fmpz_poly_degree(denpa.d) != 0 ) {
1347 MLOCK(ErrorMessageLock);
1348 MesPrint("flint::mul_poly: error: denpa is non-constant");
1349 MUNLOCK(ErrorMessageLock);
1350 Terminate(-1);
1351 }
1352 if ( fmpz_poly_degree(denpb.d) != 0 ) {
1353 MLOCK(ErrorMessageLock);
1354 MesPrint("flint::mul_poly: error: denpb is non-constant");
1355 MUNLOCK(ErrorMessageLock);
1356 Terminate(-1);
1357 }
1358
1359 // Multiply numerators, store result in pa
1360 fmpz_poly_mul(pa.d, pa.d, pb.d);
1361 // Multiply denominators, store result in denpa, and convert to an fmpz:
1362 fmpz_poly_mul(denpa.d, denpa.d, denpb.d);
1363 flint::fmpz den;
1364 fmpz_poly_get_coeff_fmpz(den.d, denpa.d, 0);
1365
1366 WORD* res;
1367 // First determine the result size, and malloc. The result should have no arghead. Here we use
1368 // the "scale" argument of to_argument_poly since den might not be 1.
1369 const bool with_arghead = false;
1370 bool write = false;
1371 const bool must_fit_term = false;
1372 const uint64_t prev_size = 0;
1373 const uint64_t mul_size = (uint64_t)flint::to_argument_poly(BHEAD NULL,
1374 with_arghead, must_fit_term, write, prev_size, pa.d, var_map, den.d);
1375 res = (WORD*)Malloc1(sizeof(WORD)*mul_size, "flint::mul_poly");
1376
1377 write = true;
1378 flint::to_argument_poly(BHEAD res, with_arghead, must_fit_term, write, prev_size, pa.d,
1379 var_map, den.d);
1380
1381 return res;
1382}
1383/*
1384 #] flint::mul_poly :
1385
1386 #[ flint::ratfun_add_mpoly :
1387*/
1388// Add the multi-variate FORM rational polynomials at t1 and t2. The result is written at out.
1389void flint::ratfun_add_mpoly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
1390 const var_map_t &var_map) {
1391
1392 flint::mpoly_ctx ctx(var_map.size());
1393 flint::mpoly gcd(ctx.d), num1(ctx.d), den1(ctx.d), num2(ctx.d), den2(ctx.d);
1394
1395 flint::ratfun_read_mpoly(t1, num1.d, den1.d, var_map, ctx.d);
1396 flint::ratfun_read_mpoly(t2, num2.d, den2.d, var_map, ctx.d);
1397
1398 if ( fmpz_mpoly_cmp(den1.d, den2.d, ctx.d) != 0 ) {
1399 fmpz_mpoly_gcd_cofactors(gcd.d, den1.d, den2.d, den1.d, den2.d, ctx.d);
1400
1401 fmpz_mpoly_mul(num1.d, num1.d, den2.d, ctx.d);
1402 fmpz_mpoly_mul(num2.d, num2.d, den1.d, ctx.d);
1403
1404 fmpz_mpoly_add(num1.d, num1.d, num2.d, ctx.d);
1405 fmpz_mpoly_mul(den1.d, den1.d, den2.d, ctx.d);
1406 fmpz_mpoly_mul(den1.d, den1.d, gcd.d, ctx.d);
1407 }
1408 else {
1409 fmpz_mpoly_add(num1.d, num1.d, num2.d, ctx.d);
1410 }
1411
1412 // Finally divide out any common factors between the resulting num1, den1:
1413 fmpz_mpoly_gcd_cofactors(gcd.d, num1.d, den1.d, num1.d, den1.d, ctx.d);
1414
1415 flint::util::fix_sign_fmpz_mpoly_ratfun(num1.d, den1.d, ctx.d);
1416
1417 // Result in FORM notation:
1418 *out++ = AR.PolyFun;
1419 WORD* args_size = out++;
1420 WORD* args_flag = out++;
1421 *args_flag = 0; // clean prf
1422 FILLFUN3(out); // Remainder of funhead, if it is larger than 3
1423
1424 const bool with_arghead = true;
1425 const bool must_fit_term = true;
1426 const bool write = true;
1427 // prev_size + 4, to account for final term size and coeff of "1/1"
1428 out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
1429 num1.d, var_map, ctx.d);
1430 out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
1431 den1.d, var_map, ctx.d);
1432
1433 *args_size = out - args_size + 1; // The +1 is to include the function ID
1434 AT.WorkPointer = out;
1435}
1436/*
1437 #] flint::ratfun_add_mpoly :
1438 #[ flint::ratfun_add_poly :
1439*/
1440// Add the uni-variate FORM rational polynomials at t1 and t2. The result is written at out.
1441void flint::ratfun_add_poly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
1442 const var_map_t &var_map) {
1443
1444 flint::poly gcd, num1, den1, num2, den2;
1445
1446 flint::ratfun_read_poly(t1, num1.d, den1.d);
1447 flint::ratfun_read_poly(t2, num2.d, den2.d);
1448
1449 if ( fmpz_poly_equal(den1.d, den2.d) == 0 ) {
1450 flint::util::simplify_fmpz_poly(den1.d, den2.d, gcd.d);
1451
1452 fmpz_poly_mul(num1.d, num1.d, den2.d);
1453 fmpz_poly_mul(num2.d, num2.d, den1.d);
1454
1455 fmpz_poly_add(num1.d, num1.d, num2.d);
1456 fmpz_poly_mul(den1.d, den1.d, den2.d);
1457 fmpz_poly_mul(den1.d, den1.d, gcd.d);
1458 }
1459 else {
1460 fmpz_poly_add(num1.d, num1.d, num2.d);
1461 }
1462
1463 // Finally divide out any common factors between the resulting num1, den1:
1464 flint::util::simplify_fmpz_poly(num1.d, den1.d, gcd.d);
1465
1466 flint::util::fix_sign_fmpz_poly_ratfun(num1.d, den1.d);
1467
1468 // Result in FORM notation:
1469 *out++ = AR.PolyFun;
1470 WORD* args_size = out++;
1471 WORD* args_flag = out++;
1472 *args_flag = 0; // clean prf
1473 FILLFUN3(out); // Remainder of funhead, if it is larger than 3
1474
1475 const bool with_arghead = true;
1476 const bool must_fit_term = true;
1477 const bool write = true;
1478 // prev_size + 4, to account for final term size and coeff of "1/1"
1479 out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
1480 num1.d, var_map);
1481 out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
1482 den1.d, var_map);
1483
1484 *args_size = out - args_size + 1; // The +1 is to include the function ID
1485 AT.WorkPointer = out;
1486}
1487/*
1488 #] flint::ratfun_add_poly :
1489
1490 #[ flint::ratfun_normalize_mpoly :
1491*/
1492// Multiply and simplify occurrences of the multi-variate FORM rational polynomials found in term.
1493// The final term is written in place, with the rational polynomial at the end.
1494void flint::ratfun_normalize_mpoly(PHEAD WORD *term, const var_map_t &var_map) {
1495
1496 // The length of the coefficient
1497 const WORD ncoeff = (term + *term)[-1];
1498 // The end of the term data, before the coefficient:
1499 const WORD *tstop = term + *term - ABS(ncoeff);
1500
1501 flint::mpoly_ctx ctx(var_map.size());
1502 flint::mpoly num1(ctx.d), den1(ctx.d), num2(ctx.d), den2(ctx.d), gcd(ctx.d);
1503
1504 // Start with "trivial" polynomials, and multiply in the num and den of the prf which appear.
1505 flint::fmpz tmpNum, tmpDen;
1506 flint::fmpz_set_form(tmpNum.d, (UWORD*)tstop, ncoeff/2);
1507 flint::fmpz_set_form(tmpDen.d, (UWORD*)tstop+ABS(ncoeff/2), ABS(ncoeff/2));
1508 fmpz_mpoly_set_fmpz(num1.d, tmpNum.d, ctx.d);
1509 fmpz_mpoly_set_fmpz(den1.d, tmpDen.d, ctx.d);
1510
1511 // Loop over the occurrences of PolyFun in the term, and multiply in to num1, den1.
1512 // s tracks where we are writing the non-PolyFun term data. The final PolyFun will
1513 // go at the end.
1514 WORD* term_size = term;
1515 WORD* s = term + 1;
1516 for ( WORD *t = term + 1; t < tstop; ) {
1517 if ( *t == AR.PolyFun ) {
1518 flint::ratfun_read_mpoly(t, num2.d, den2.d, var_map, ctx.d);
1519
1520 // get gcd of num1,den2 and num2,den1 and then assemble
1521 fmpz_mpoly_gcd_cofactors(gcd.d, num1.d, den2.d, num1.d, den2.d, ctx.d);
1522 fmpz_mpoly_gcd_cofactors(gcd.d, num2.d, den1.d, num2.d, den1.d, ctx.d);
1523
1524 fmpz_mpoly_mul(num1.d, num1.d, num2.d, ctx.d);
1525 fmpz_mpoly_mul(den1.d, den1.d, den2.d, ctx.d);
1526
1527 t += t[1];
1528 }
1529
1530 else {
1531 // Not a PolyFun, just copy or skip over
1532 WORD i = t[1];
1533 if ( s != t ) { NCOPY(s,t,i); }
1534 else { t += i; s += i; }
1535 }
1536 }
1537
1538 flint::util::fix_sign_fmpz_mpoly_ratfun(num1.d, den1.d, ctx.d);
1539
1540 // Result in FORM notation:
1541 WORD* out = s;
1542 *out++ = AR.PolyFun;
1543 WORD* args_size = out++;
1544 WORD* args_flag = out++;
1545 *args_flag &= ~MUSTCLEANPRF;
1546
1547 const bool with_arghead = true;
1548 const bool must_fit_term = true;
1549 const bool write = true;
1550 out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
1551 num1.d, var_map, ctx.d);
1552 out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
1553 den1.d, var_map, ctx.d);
1554
1555 *args_size = out - args_size + 1; // The +1 is to include the function ID
1556
1557 // +3 for the coefficient of 1/1, which is added after the check
1558 if ( sizeof(WORD)*(out-term_size+3) > (size_t)AM.MaxTer ) {
1559 MLOCK(ErrorMessageLock);
1560 MesPrint("flint::ratfun_normalize: output exceeds MaxTermSize");
1561 MUNLOCK(ErrorMessageLock);
1562 Terminate(-1);
1563 }
1564
1565 *out++ = 1;
1566 *out++ = 1;
1567 *out++ = 3; // the term's coefficient is now 1/1
1568
1569 *term_size = out - term_size;
1570}
1571/*
1572 #] flint::ratfun_normalize_mpoly :
1573 #[ flint::ratfun_normalize_poly :
1574*/
1575// Multiply and simplify occurrences of the uni-variate FORM rational polynomials found in term.
1576// The final term is written in place, with the rational polynomial at the end.
1577void flint::ratfun_normalize_poly(PHEAD WORD *term, const var_map_t &var_map) {
1578
1579 // The length of the coefficient
1580 const WORD ncoeff = (term + *term)[-1];
1581 // The end of the term data, before the coefficient:
1582 const WORD *tstop = term + *term - ABS(ncoeff);
1583
1584 flint::poly num1, den1, num2, den2, gcd;
1585
1586 // Start with "trivial" polynomials, and multiply in the num and den of the prf which appear.
1587 flint::fmpz tmpNum, tmpDen;
1588 flint::fmpz_set_form(tmpNum.d, (UWORD*)tstop, ncoeff/2);
1589 flint::fmpz_set_form(tmpDen.d, (UWORD*)tstop+ABS(ncoeff/2), ABS(ncoeff/2));
1590 fmpz_poly_set_fmpz(num1.d, tmpNum.d);
1591 fmpz_poly_set_fmpz(den1.d, tmpDen.d);
1592
1593 // Loop over the occurrences of PolyFun in the term, and multiply in to num1, den1.
1594 // s tracks where we are writing the non-PolyFun term data. The final PolyFun will
1595 // go at the end.
1596 WORD* term_size = term;
1597 WORD* s = term + 1;
1598 for ( WORD *t = term + 1; t < tstop; ) {
1599 if ( *t == AR.PolyFun ) {
1600 flint::ratfun_read_poly(t, num2.d, den2.d);
1601
1602 // get gcd of num1,den2 and num2,den1 and then assemble
1603 flint::util::simplify_fmpz_poly(num1.d, den2.d, gcd.d);
1604 flint::util::simplify_fmpz_poly(num2.d, den1.d, gcd.d);
1605
1606 fmpz_poly_mul(num1.d, num1.d, num2.d);
1607 fmpz_poly_mul(den1.d, den1.d, den2.d);
1608
1609 t += t[1];
1610 }
1611
1612 else {
1613 // Not a PolyFun, just copy or skip over
1614 WORD i = t[1];
1615 if ( s != t ) { NCOPY(s,t,i); }
1616 else { t += i; s += i; }
1617 }
1618 }
1619
1620 flint::util::fix_sign_fmpz_poly_ratfun(num1.d, den1.d);
1621
1622 // Result in FORM notation:
1623 WORD* out = s;
1624 *out++ = AR.PolyFun;
1625 WORD* args_size = out++;
1626 WORD* args_flag = out++;
1627 *args_flag &= ~MUSTCLEANPRF;
1628
1629 const bool with_arghead = true;
1630 const bool must_fit_term = true;
1631 const bool write = true;
1632 out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
1633 num1.d, var_map);
1634 out += flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
1635 den1.d, var_map);
1636
1637 *args_size = out - args_size + 1; // The +1 is to include the function ID
1638
1639 // +3 for the coefficient of 1/1, which is added after the check
1640 if ( sizeof(WORD)*(out-term_size+3) > (size_t)AM.MaxTer ) {
1641 MLOCK(ErrorMessageLock);
1642 MesPrint("flint::ratfun_normalize: output exceeds MaxTermSize");
1643 MUNLOCK(ErrorMessageLock);
1644 Terminate(-1);
1645 }
1646
1647 *out++ = 1;
1648 *out++ = 1;
1649 *out++ = 3; // the term's coefficient is now 1/1
1650
1651 *term_size = out - term_size;
1652}
1653/*
1654 #] flint::ratfun_normalize_poly :
1655
1656 #[ flint::ratfun_read_mpoly :
1657*/
1658// Read the multi-variate FORM rational polynomial at a and create fmpz_mpoly_t numerator and
1659// denominator.
1660void flint::ratfun_read_mpoly(const WORD *a, fmpz_mpoly_t num, fmpz_mpoly_t den,
1661 const var_map_t &var_map, fmpz_mpoly_ctx_t ctx) {
1662
1663 // The end of the arguments:
1664 const WORD* arg_stop = a+a[1];
1665
1666 const bool must_normalize = (a[2] & MUSTCLEANPRF) != 0;
1667
1668 a += FUNHEAD;
1669 if ( a >= arg_stop ) {
1670 MLOCK(ErrorMessageLock);
1671 MesPrint("ERROR: PolyRatFun cannot have zero arguments");
1672 MUNLOCK(ErrorMessageLock);
1673 Terminate(-1);
1674 }
1675
1676 // Polys to collect the "den of the num" and "den of the den".
1677 // Input can arrive like this when enabling the PolyRatFun or moving things into it.
1678 flint::mpoly den_num(ctx), den_den(ctx);
1679
1680 // Read the numerator
1681 flint::from_argument_mpoly(num, den_num.d, a, true, var_map, ctx);
1682 NEXTARG(a);
1683
1684 if ( a < arg_stop ) {
1685 // Read the denominator
1686 flint::from_argument_mpoly(den, den_den.d, a, true, var_map, ctx);
1687 NEXTARG(a);
1688 }
1689 else {
1690 // The denominator is 1
1691 MLOCK(ErrorMessageLock);
1692 MesPrint("implement this");
1693 MUNLOCK(ErrorMessageLock);
1694 Terminate(-1);
1695 }
1696 if ( a < arg_stop ) {
1697 MLOCK(ErrorMessageLock);
1698 MesPrint("ERROR: PolyRatFun cannot have more than two arguments");
1699 MUNLOCK(ErrorMessageLock);
1700 Terminate(-1);
1701 }
1702
1703 // Multiply the num by den_den and den by den_num:
1704 fmpz_mpoly_mul(num, num, den_den.d, ctx);
1705 fmpz_mpoly_mul(den, den, den_num.d, ctx);
1706
1707 if ( must_normalize ) {
1708 flint::mpoly gcd(ctx);
1709 fmpz_mpoly_gcd_cofactors(gcd.d, num, den, num, den, ctx);
1710 }
1711}
1712/*
1713 #] flint::ratfun_read_mpoly :
1714 #[ flint::ratfun_read_poly :
1715*/
1716// Read the uni-variate FORM rational polynomial at a and create fmpz_mpoly_t numerator and
1717// denominator.
1718void flint::ratfun_read_poly(const WORD *a, fmpz_poly_t num, fmpz_poly_t den) {
1719
1720 // The end of the arguments:
1721 const WORD* arg_stop = a+a[1];
1722
1723 const bool must_normalize = (a[2] & MUSTCLEANPRF) != 0;
1724
1725 a += FUNHEAD;
1726 if ( a >= arg_stop ) {
1727 MLOCK(ErrorMessageLock);
1728 MesPrint("ERROR: PolyRatFun cannot have zero arguments");
1729 MUNLOCK(ErrorMessageLock);
1730 Terminate(-1);
1731 }
1732
1733 // Polys to collect the "den of the num" and "den of the den".
1734 // Input can arrive like this when enabling the PolyRatFun or moving things into it.
1735 flint::poly den_num, den_den;
1736
1737 // Read the numerator
1738 flint::from_argument_poly(num, den_num.d, a, true);
1739 NEXTARG(a);
1740
1741 if ( a < arg_stop ) {
1742 // Read the denominator
1743 flint::from_argument_poly(den, den_den.d, a, true);
1744 NEXTARG(a);
1745 }
1746 else {
1747 // The denominator is 1
1748 MLOCK(ErrorMessageLock);
1749 MesPrint("implement this");
1750 MUNLOCK(ErrorMessageLock);
1751 Terminate(-1);
1752 }
1753 if ( a < arg_stop ) {
1754 MLOCK(ErrorMessageLock);
1755 MesPrint("ERROR: PolyRatFun cannot have more than two arguments");
1756 MUNLOCK(ErrorMessageLock);
1757 Terminate(-1);
1758 }
1759
1760 // Multiply the num by den_den and den by den_num:
1761 fmpz_poly_mul(num, num, den_den.d);
1762 fmpz_poly_mul(den, den, den_num.d);
1763
1764 if ( must_normalize ) {
1765 flint::poly gcd;
1766 flint::util::simplify_fmpz_poly(num, den, gcd.d);
1767 }
1768}
1769/*
1770 #] flint::ratfun_read_poly :
1771 #[ flint::to_argument_mpoly :
1772*/
1773// Convert a fmpz_mpoly_t to a FORM argument (or 0-terminated list of terms: with_arghead==false).
1774// If the caller is building an output term, prev_size contains the size of the term so far, to
1775// check that the output fits in AM.MaxTer if must_fit_term.
1776// All coefficients will be divided by denscale (which might just be 1).
1777// If write is false, we never write to out but only track the total would-be size. This lets this
1778// function be repurposed as a "size of FORM notation" function without duplicating the code.
1779#define IFW(x) { if ( write ) {x;} }
1780uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
1781 const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
1782 const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const fmpz_t denscale) {
1783
1784 // out is modified later, keep the pointer at entry
1785 const WORD* out_entry = out;
1786
1787 // Track the total size written. We could do this with pointer differences, but if
1788 // write == false we don't write to or move out to be able to find the size that way.
1789 uint64_t ws = 0;
1790
1791 // Check there is at least space for ARGHEAD WORDs (the arghead or fast-notation number/symbol)
1792 if ( write && must_fit_term && (sizeof(WORD)*(prev_size + ARGHEAD) > (size_t)AM.MaxTer) ) {
1793 MLOCK(ErrorMessageLock);
1794 MesPrint("flint::to_argument_mpoly: output exceeds MaxTermSize");
1795 MUNLOCK(ErrorMessageLock);
1796 Terminate(-1);
1797 }
1798
1799 // Create the inverse of var_map, so we don't have to search it for each symbol written
1800 var_map_t var_map_inv;
1801 for ( auto x: var_map ) {
1802 var_map_inv[x.second] = x.first;
1803 }
1804
1805 vector<int64_t> exponents(var_map.size());
1806 const int64_t n_terms = fmpz_mpoly_length(poly, ctx);
1807
1808 if ( n_terms == 0 ) {
1809 if ( with_arghead ) {
1810 IFW(*out++ = -SNUMBER); ws++;
1811 IFW(*out++ = 0); ws++;
1812 return ws;
1813 }
1814 else {
1815 IFW(*out++ = 0); ws++;
1816 return ws;
1817 }
1818 }
1819
1820 // For dividing out denscale
1821 flint::fmpz coeff, den, gcd;
1822
1823 // The mpoly might be constant or a single symbol with coeff 1. Use fast notation if possible.
1824 if ( with_arghead && n_terms == 1 ) {
1825
1826 if ( fmpz_mpoly_is_fmpz(poly, ctx) ) {
1827 // The mpoly is constant. Use fast notation if the number is an integer and small enough:
1828
1829 fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, 0, ctx);
1830 fmpz_set(den.d, denscale);
1831 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
1832
1833 if ( fmpz_is_one(den.d) && fmpz_fits_si(coeff.d) ) {
1834 const int64_t fast_coeff = fmpz_get_si(coeff.d);
1835 // While ">=", could work here, FORM does not use fast notation for INT_MIN
1836 if ( fast_coeff > INT32_MIN && fast_coeff <= INT32_MAX ) {
1837 IFW(*out++ = -SNUMBER); ws++;
1838 IFW(*out++ = (WORD)fast_coeff); ws++;
1839 return ws;
1840 }
1841 }
1842 }
1843
1844 else {
1845 fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, 0, ctx);
1846 fmpz_set(den.d, denscale);
1847 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
1848
1849 if ( fmpz_is_one(coeff.d) && fmpz_is_one(den.d) ) {
1850 // The coefficient is one. Now check the symbol powers:
1851
1852 fmpz_mpoly_get_term_exp_si((slong*)exponents.data(), poly, 0, ctx);
1853 int64_t use_fast = 0;
1854 uint32_t fast_symbol = 0;
1855
1856 for ( size_t i = 0; i < var_map.size(); i++ ) {
1857 if ( exponents[i] == 1 ) fast_symbol = var_map_inv[i];
1858 use_fast += exponents[i];
1859 }
1860
1861 // use_fast has collected the total degree. If it is 1, then fast_symbol holds the code
1862 if ( use_fast == 1 ) {
1863 IFW(*out++ = -SYMBOL); ws++;
1864 IFW(*out++ = fast_symbol); ws++;
1865 return ws;
1866 }
1867 }
1868 }
1869 }
1870
1871
1872 WORD *tmp_coeff = (WORD *)NumberMalloc("flint::to_argument_mpoly");
1873 WORD *tmp_den = (WORD *)NumberMalloc("flint::to_argument_mpoly");
1874
1875 WORD* arg_size = 0;
1876 WORD* arg_flag = 0;
1877 if ( with_arghead ) {
1878 IFW(arg_size = out++); ws++; // total arg size
1879 IFW(arg_flag = out++); ws++;
1880 IFW(*arg_flag = 0); // clean argument
1881 }
1882
1883 for ( int64_t i = 0; i < n_terms; i++ ) {
1884
1885 fmpz_mpoly_get_term_exp_si((slong*)exponents.data(), poly, i, ctx);
1886
1887 fmpz_mpoly_get_term_coeff_fmpz(coeff.d, poly, i, ctx);
1888 fmpz_set(den.d, denscale);
1889 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
1890
1891 uint32_t num_symbols = 0;
1892 for ( size_t j = 0; j < var_map.size(); j++ ) {
1893 if ( exponents[j] != 0 ) { num_symbols += 1; }
1894 }
1895
1896 // Convert the coefficient, write in temporary space
1897 const WORD num_size = flint::fmpz_get_form(coeff.d, tmp_coeff);
1898 const WORD den_size = flint::fmpz_get_form(den.d, tmp_den);
1899 const WORD coeff_size = [num_size, den_size] () -> WORD {
1900 WORD size = ABS(num_size) > ABS(den_size) ? ABS(num_size) : ABS(den_size);
1901 return size * SGN(num_size) * SGN(den_size);
1902 }();
1903
1904 // Now we have the number of symbols and the coeff size, we can determine the output size.
1905 // Check it fits if necessary: term size, num,den of "coeff_size", +1 for total coeff size
1906 uint64_t current_size = prev_size + ws + 1 + 2*ABS(coeff_size) + 1;
1907 if ( num_symbols ) {
1908 // symbols header, code,power of each symbol:
1909 current_size += 2 + 2*num_symbols;
1910 }
1911 if ( write && must_fit_term && (sizeof(WORD)*current_size > (size_t)AM.MaxTer) ) {
1912 MLOCK(ErrorMessageLock);
1913 MesPrint("flint::to_argument_mpoly: output exceeds MaxTermSize");
1914 MUNLOCK(ErrorMessageLock);
1915 Terminate(-1);
1916 }
1917
1918 WORD* term_size = 0;
1919 IFW(term_size = out++); ws++;
1920 if ( num_symbols ) {
1921 IFW(*out++ = SYMBOL); ws++;
1922 WORD* symbol_size = 0;
1923 IFW(symbol_size = out++); ws++;
1924 IFW(*symbol_size = 2);
1925
1926 for ( size_t j = 0; j < var_map.size(); j++ ) {
1927 if ( exponents[j] != 0 ) {
1928 IFW(*out++ = var_map_inv[j]); ws++;
1929 IFW(*out++ = exponents[j]); ws++;
1930 IFW(*symbol_size += 2);
1931 }
1932 }
1933 }
1934
1935 // Copy numerator
1936 for ( WORD j = 0; j < ABS(num_size); j++ ) {
1937 IFW(*out++ = tmp_coeff[j]); ws++;
1938 }
1939 for ( WORD j = ABS(num_size); j < ABS(coeff_size); j++ ) {
1940 IFW(*out++ = 0); ws++;
1941 }
1942 // Copy denominator
1943 for ( WORD j = 0; j < ABS(den_size); j++ ) {
1944 IFW(*out++ = tmp_den[j]); ws++;
1945 }
1946 for ( WORD j = ABS(den_size); j < ABS(coeff_size); j++ ) {
1947 IFW(*out++ = 0); ws++;
1948 }
1949
1950 IFW(*out = 2*ABS(coeff_size) + 1); // the size of the coefficient
1951 IFW(if ( coeff_size < 0 ) { *out = -(*out); });
1952 IFW(out++); ws++;
1953
1954 IFW(*term_size = out - term_size);
1955 }
1956
1957 if ( with_arghead ) {
1958 IFW(*arg_size = out - arg_size);
1959 if ( write ) {
1960 // Sort into form highfirst ordering
1961 flint::form_sort(BHEAD (WORD*)(out_entry));
1962 }
1963 }
1964 else {
1965 // with no arghead, we write a terminating zero
1966 IFW(*out++ = 0); ws++;
1967 }
1968
1969 NumberFree(tmp_coeff, "flint::to_argument_mpoly");
1970 NumberFree(tmp_den, "flint::to_argument_mpoly");
1971
1972 return ws;
1973}
1974
1975// If no denscale argument is supplied, just set it to 1 and call the usual function
1976uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
1977 const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
1978 const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx) {
1979
1980 flint::fmpz tmp;
1981 fmpz_set_ui(tmp.d, 1);
1982
1983 uint64_t ret = flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write,
1984 prev_size, poly, var_map, ctx, tmp.d);
1985
1986 return ret;
1987}
1988/*
1989 #] flint::to_argument_mpoly :
1990 #[ flint::to_argument_poly :
1991*/
1992// Convert a fmpz_poly_t to a FORM argument (or 0-terminated list of terms: with_arghead==false).
1993// If the caller is building an output term, prev_size contains the size of the term so far, to
1994// check that the output fits in AM.MaxTer if must_fit_term.
1995// All coefficients will be divided by denscale (which might just be 1).
1996// If write is false, we never write to out but only track the total would-be size. This lets this
1997// function be repurposed as a "size of FORM notation" function without duplicating the code.
1998uint64_t flint::to_argument_poly(PHEAD WORD *out, const bool with_arghead,
1999 const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_poly_t poly,
2000 const var_map_t &var_map, const fmpz_t denscale) {
2001
2002 // Track the total size written. We could do this with pointer differences, but if
2003 // write == false we don't write to or move out to be able to find the size that way.
2004 uint64_t ws = 0;
2005
2006 // Check there is at least space for ARGHEAD WORDs (the arghead or fast-notation number/symbol)
2007 if ( write && must_fit_term && (sizeof(WORD)*(prev_size + ARGHEAD) > (size_t)AM.MaxTer) ) {
2008 MLOCK(ErrorMessageLock);
2009 MesPrint("flint::to_argument_poly: output exceeds MaxTermSize");
2010 MUNLOCK(ErrorMessageLock);
2011 Terminate(-1);
2012 }
2013
2014 // Create the inverse of var_map, so we don't have to search it for each symbol written
2015 var_map_t var_map_inv;
2016 for ( auto x: var_map ) {
2017 var_map_inv[x.second] = x.first;
2018 }
2019
2020 const int64_t n_terms = fmpz_poly_length(poly);
2021
2022 // The poly is zero
2023 if ( n_terms == 0 ) {
2024 if ( with_arghead ) {
2025 IFW(*out++ = -SNUMBER); ws++;
2026 IFW(*out++ = 0); ws++;
2027 return ws;
2028 }
2029 else {
2030 IFW(*out++ = 0); ws++;
2031 return ws;
2032 }
2033 }
2034
2035 // For dividing out denscale
2036 flint::fmpz coeff, den, gcd;
2037
2038 // The poly is constant, use fast notation if the coefficient is integer and small enough
2039 if ( with_arghead && n_terms == 1 ) {
2040
2041 fmpz_poly_get_coeff_fmpz(coeff.d, poly, 0);
2042 fmpz_set(den.d, denscale);
2043 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
2044
2045 if ( fmpz_is_one(den.d) && fmpz_fits_si(coeff.d) ) {
2046 const long fast_coeff = fmpz_get_si(coeff.d);
2047 // While ">=", could work here, FORM does not use fast notation for INT_MIN
2048 if ( fast_coeff > INT_MIN && fast_coeff <= INT_MAX ) {
2049 IFW(*out++ = -SNUMBER); ws++;
2050 IFW(*out++ = (WORD)fast_coeff); ws++;
2051 return ws;
2052 }
2053 }
2054 }
2055
2056 // The poly might be a single symbol with coeff 1, use fast notation if so.
2057 if ( with_arghead && n_terms == 2 ) {
2058 if ( fmpz_is_zero(fmpz_poly_get_coeff_ptr(poly, 0)) ) {
2059 // The constant term is zero
2060
2061 fmpz_poly_get_coeff_fmpz(coeff.d, poly, 1);
2062 fmpz_set(den.d, denscale);
2063 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
2064
2065 if ( fmpz_is_one(coeff.d) && fmpz_is_one(den.d) ) {
2066 // Single symbol with coeff 1. Use fast notation:
2067 IFW(*out++ = -SYMBOL); ws++;
2068 IFW(*out++ = var_map_inv[0]); ws++;
2069 return ws;
2070 }
2071 }
2072 }
2073
2074 WORD *tmp_coeff = (WORD *)NumberMalloc("flint::to_argument_poly");
2075 WORD *tmp_den = (WORD *)NumberMalloc("flint::to_argument_mpoly");
2076
2077 WORD* arg_size = 0;
2078 WORD* arg_flag = 0;
2079 if ( with_arghead ) {
2080 IFW(arg_size = out++); ws++; // total arg size
2081 IFW(arg_flag = out++); ws++;
2082 IFW(*arg_flag = 0); // clean argument
2083 }
2084
2085 // In reverse, since we want a "highfirst" output
2086 for ( int64_t i = n_terms-1; i >= 0; i-- ) {
2087
2088 // fmpz_poly is dense, there might be many zero coefficients:
2089 if ( !fmpz_is_zero(fmpz_poly_get_coeff_ptr(poly, i)) ) {
2090
2091 fmpz_poly_get_coeff_fmpz(coeff.d, poly, i);
2092 fmpz_set(den.d, denscale);
2093 flint::util::simplify_fmpz(coeff.d, den.d, gcd.d);
2094
2095 // Convert the coefficient, write in temporary space
2096 const WORD num_size = flint::fmpz_get_form(coeff.d, tmp_coeff);
2097 const WORD den_size = flint::fmpz_get_form(den.d, tmp_den);
2098 const WORD coeff_size = [num_size, den_size] () -> WORD {
2099 WORD size = ABS(num_size) > ABS(den_size) ? ABS(num_size) : ABS(den_size);
2100 return size * SGN(num_size) * SGN(den_size);
2101 }();
2102
2103 // Now we have the coeff size, we can determine the output size
2104 // Check it fits if necessary: symbol code,power, num,den of "coeff_size",
2105 // +1 for total coeff size
2106 uint64_t current_size = prev_size + ws + 1 + 2*ABS(coeff_size) + 1;
2107 if ( i > 0 ) {
2108 // and also symbols header, code,power of the symbol
2109 current_size += 4;
2110 }
2111 if ( write && must_fit_term && (sizeof(WORD)*current_size > (size_t)AM.MaxTer) ) {
2112 MLOCK(ErrorMessageLock);
2113 MesPrint("flint::to_argument_poly: output exceeds MaxTermSize");
2114 MUNLOCK(ErrorMessageLock);
2115 Terminate(-1);
2116 }
2117
2118 WORD* term_size = 0;
2119 IFW(term_size = out++); ws++;
2120
2121 if ( i > 0 ) {
2122 IFW(*out++ = SYMBOL); ws++;
2123 IFW(*out++ = 4); ws++; // The symbol array size, it is univariate
2124 IFW(*out++ = var_map_inv[0]); ws++;
2125 IFW(*out++ = i); ws++;
2126 }
2127
2128 // Copy numerator
2129 for ( WORD j = 0; j < ABS(num_size); j++ ) {
2130 IFW(*out++ = tmp_coeff[j]); ws++;
2131 }
2132 for ( WORD j = ABS(num_size); j < ABS(coeff_size); j++ ) {
2133 IFW(*out++ = 0); ws++;
2134 }
2135 // Copy denominator
2136 for ( WORD j = 0; j < ABS(den_size); j++ ) {
2137 IFW(*out++ = tmp_den[j]); ws++;
2138 }
2139 for ( WORD j = ABS(den_size); j < ABS(coeff_size); j++ ) {
2140 IFW(*out++ = 0); ws++;
2141 }
2142
2143 IFW(*out = 2*ABS(coeff_size) + 1); // the size of the coefficient
2144 IFW(if ( coeff_size < 0 ) { *out = -(*out); });
2145 IFW(out++); ws++;
2146
2147 IFW(*term_size = out - term_size);
2148 }
2149
2150 }
2151
2152 if ( with_arghead ) {
2153 IFW(*arg_size = out - arg_size);
2154 }
2155 else {
2156 // with no arghead, we write a terminating zero
2157 IFW(*out++ = 0); ws++;
2158 }
2159
2160 NumberFree(tmp_coeff, "flint::to_argument_poly");
2161 NumberFree(tmp_den, "flint::to_argument_poly");
2162
2163 return ws;
2164}
2165
2166// If no denscale argument is supplied, just set it to 1 and call the usual function
2167uint64_t flint::to_argument_poly(PHEAD WORD *out, const bool with_arghead,
2168 const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_poly_t poly,
2169 const var_map_t &var_map) {
2170
2171 flint::fmpz tmp;
2172 fmpz_set_ui(tmp.d, 1);
2173
2174 uint64_t ret = flint::to_argument_poly(BHEAD out, with_arghead, must_fit_term, write, prev_size,
2175 poly, var_map, tmp.d);
2176
2177 return ret;
2178}
2179/*
2180 #] flint::to_argument_poly :
2181*/
2182
2183// Utility functions
2184/*
2185 #[ flint::util::simplify_fmpz :
2186*/
2187// Divide the GCD out of num and den
2188inline void flint::util::simplify_fmpz(fmpz_t num, fmpz_t den, fmpz_t gcd) {
2189 fmpz_gcd(gcd, num, den);
2190 if ( !fmpz_is_one(gcd) ) {
2191 fmpz_divexact(num, num, gcd);
2192 fmpz_divexact(den, den, gcd);
2193 }
2194}
2195/*
2196 #] flint::util::simplify_fmpz :
2197 #[ flint::util::simplify_fmpz_poly :
2198*/
2199// Divide the GCD out of num and den
2200inline void flint::util::simplify_fmpz_poly(fmpz_poly_t num, fmpz_poly_t den, fmpz_poly_t gcd) {
2201 fmpz_poly_gcd(gcd, num, den);
2202 if ( !fmpz_poly_is_one(gcd) ) {
2203#if __FLINT_RELEASE >= 30100
2204 // This should be faster than fmpz_poly_div, see https://github.com/flintlib/flint/pull/1766
2205 fmpz_poly_divexact(num, num, gcd);
2206 fmpz_poly_divexact(den, den, gcd);
2207#else
2208 fmpz_poly_div(num, num, gcd);
2209 fmpz_poly_div(den, den, gcd);
2210#endif
2211 }
2212}
2213/*
2214 #] flint::util::simplify_fmpz_poly :
2215 #[ flint::util::fix_sign_fmpz_mpoly_ratfun :
2216*/
2217inline void flint::util::fix_sign_fmpz_mpoly_ratfun(fmpz_mpoly_t num, fmpz_mpoly_t den,
2218 const fmpz_mpoly_ctx_t ctx) {
2219 // Fix sign to align with poly: the leading denominator term should have a positive coeff
2220 if ( fmpz_sgn(fmpz_mpoly_term_coeff_ref(den, 0, ctx)) == -1 ) {
2221 fmpz_mpoly_neg(num, num, ctx);
2222 fmpz_mpoly_neg(den, den, ctx);
2223 }
2224}
2225/*
2226 #] flint::util::fix_sign_fmpz_mpoly_ratfun :
2227 #[ flint::util::fix_sign_fmpz_poly_ratfun :
2228*/
2229inline void flint::util::fix_sign_fmpz_poly_ratfun(fmpz_poly_t num, fmpz_poly_t den) {
2230 // Fix sign to align with poly: the leading denominator term should have a positive coeff
2231 if ( fmpz_sgn(fmpz_poly_get_coeff_ptr(den, fmpz_poly_degree(den))) == -1 ) {
2232 fmpz_poly_neg(num, num);
2233 fmpz_poly_neg(den, den);
2234 }
2235}
2236/*
2237 #] flint::util::fix_sign_fmpz_poly_ratfun :
2238*/
flint::mpoly_ctx
Definition flintinterface.h:125
flint::mpoly_factor
Definition flintinterface.h:115
flint::poly_factor
Definition flintinterface.h:96
poly
Definition poly.h:53
EndSort
LONG EndSort(PHEAD WORD *, int)
Definition sort.c:454
LowerSortLevel
void LowerSortLevel(void)
Definition sort.c:4661
StoreTerm
int StoreTerm(PHEAD WORD *)
Definition sort.c:4244
NewSort
int NewSort(PHEAD0)
Definition sort.c:359
Compare1
WORD Compare1(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2341
CompareSymbols
WORD CompareSymbols(PHEAD WORD *, WORD *, WORD)
Definition sort.c:2804
SymbolNormalize
int SymbolNormalize(WORD *)
Definition normal.c:5195
Generated by doxygen 1.9.8