@@ -202,6 +202,20 @@ rat_fit_length(rat_t rat, ulong len)
202202 }
203203}
204204
205+ void
206+ rat_set (rat_t rat, const rat_t src)
207+ {
208+ rat_one (rat);
209+ fmpq_set (rat->numfactor , src->numfactor );
210+ rat_fit_length (rat, src->num );
211+ for (ulong i = 0 ; i < src->num ; i++) {
212+ fmpz_mpoly_init (&rat->factors [i], ctx);
213+ fmpz_mpoly_set (&rat->factors [i], &src->factors [i], ctx);
214+ rat->powers [i] = src->powers [i];
215+ }
216+ rat->num = src->num ;
217+ }
218+
205219void
206220rat_mul_fmpq (rat_t rat, const fmpq_t x, const int power)
207221{
@@ -386,11 +400,8 @@ rat_cofactorize(rat_t rat)
386400{
387401 LOGME;
388402 logd (" Cofactorizing %r" , rat);
389- rat_sort (rat);
390- fmpz_mpoly_t gcd, aprime, bprime;
403+ fmpz_mpoly_t gcd;
391404 fmpz_mpoly_init (gcd, ctx);
392- fmpz_mpoly_init (aprime, ctx);
393- fmpz_mpoly_init (bprime, ctx);
394405 for (ulong i = 0 ; i < rat->num ; i++) {
395406 #define A (&rat->factors[i])
396407 #define Apower rat->powers[i]
@@ -403,11 +414,9 @@ rat_cofactorize(rat_t rat)
403414 if (fmpz_mpoly_is_one (gcd, ctx)) {
404415 continue ;
405416 }
406- logd (" Split factor: %p^%d from %p^%d %p^%d" , gcd, Apower + Bpower, A, Apower, B, Bpower);
407- if (fmpz_mpoly_divides (aprime, A, gcd, ctx) != 1 ) { assert (0 ); }
408- if (fmpz_mpoly_divides (bprime, B, gcd, ctx) != 1 ) { assert (0 ); }
409- fmpz_mpoly_swap (A, aprime, ctx);
410- fmpz_mpoly_swap (B, bprime, ctx);
417+ logd (" Split factor %p^%d from %p^%d %p^%d" , gcd, Apower + Bpower, A, Apower, B, Bpower);
418+ if (fmpz_mpoly_divides (A, A, gcd, ctx) != 1 ) { assert (0 ); }
419+ if (fmpz_mpoly_divides (B, B, gcd, ctx) != 1 ) { assert (0 ); }
411420 rat_mul_fmpz_mpoly_setx (rat, gcd, Apower + Bpower);
412421 if (fmpz_mpoly_is_fmpz (B, ctx)) {
413422 rat_swapoff_fmpz_factor (rat, j);
@@ -425,9 +434,6 @@ rat_cofactorize(rat_t rat)
425434 #undef Apower
426435 }
427436 fmpz_mpoly_clear (gcd, ctx);
428- fmpz_mpoly_clear (aprime, ctx);
429- fmpz_mpoly_clear (bprime, ctx);
430- rat_sort (rat);
431437}
432438
433439void
@@ -559,7 +565,9 @@ rat_add_setx(rat_t res, rat_t rat1, rat_t rat2)
559565 #undef Apower
560566 #undef B
561567 #undef Bpower
568+ rat_sort (res);
562569 rat_cofactorize (res);
570+ rat_sort (res);
563571 logd (" Result is %r" , res);
564572}
565573
@@ -655,24 +663,11 @@ ratsum_sum_setx(rat_t rat, ratsum_t sum)
655663 }
656664 assert ((idx1 >= 0 ) && (idx2 >= 0 ) && (idx1 != idx2));
657665 logd (" Adding pair %ld of %ld" , n+1 , sum->num -1 );
658- // rat_set(tmp1, &sum->terms[idx1]);
659- // rat_set(tmp2, &sum->terms[idx2]);
660666 rat_add_setx (rat, &sum->terms [idx1], &sum->terms [idx2]);
661- /*
662- if (1) {
663- logd("Same pair, but with combined den");
664- rat_combine_denominators(tmp1);
665- rat_combine_denominators(tmp2);
666- rat_add_setx(tmp3, tmp1, tmp2);
667- }
668- */
669667 rat_swap (&sum->terms [idx1], rat);
670668 rat_swap (&sum->terms [idx2], &sum->terms [sum->num - n - 1 ]);
671669 }
672670 rat_swap (&sum->terms [0 ], rat);
673- // rat_clear(tmp1);
674- // rat_clear(tmp2);
675- // rat_clear(tmp3);
676671}
677672
678673/* GiNaC conversion
@@ -791,7 +786,9 @@ ratsum_of_ginac(ratsum_t sum, const GiNaC::ex &expr)
791786 rat_t rat;
792787 rat_init (rat);
793788 rat_of_ginac (rat, term);
789+ rat_sort (rat);
794790 rat_cofactorize (rat);
791+ rat_sort (rat);
795792 ratsum_add_rat_setx (sum, rat);
796793 rat_clear (rat);
797794 });
@@ -940,10 +937,11 @@ main(int argc, char *argv[])
940937 int nthreads = 1 ;
941938 const char *inputfile = " -" ;
942939 const char *outputfile = " -" ;
943- for (int opt; (opt = getopt (argc, (char *const *)argv, " j:hC " )) != -1 ;) {
940+ for (int opt; (opt = getopt (argc, (char *const *)argv, " j:hCV " )) != -1 ;) {
944941 switch (opt) {
945942 case ' j' : nthreads = atoi (optarg); break ;
946943 case ' h' : usage (stdout); return 0 ;
944+ case ' V' : printf (" %s" , VERSION); return 0 ;
947945 case ' C' : COLORS = true ; break ;
948946 default : return 1 ;
949947 }
0 commit comments