Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'src/coefficient.c')
-rw-r--r--src/coefficient.c245
1 files changed, 237 insertions, 8 deletions
diff --git a/src/coefficient.c b/src/coefficient.c
index d4cb6eb..1efe8fd 100644
--- a/src/coefficient.c
+++ b/src/coefficient.c
@@ -1,5 +1,5 @@
/*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
@@ -202,6 +202,235 @@ int coefficient_simplify(Coefficient* coefficient){
return(0);
}
+// put all terms under a common denominator and simplify the resulting fraction
+int coefficient_simplify_rational(Coefficient constant, Coefficient* coefficient){
+ int ret;
+ Coefficient remainder;
+ Coefficient quotient;
+ Coefficient quotient_prev;
+ Coefficient out;
+ int power;
+ int max_power;
+
+ // common denominator
+ coefficient_common_denominator(constant, coefficient);
+
+ // init
+ init_Coefficient(&out, COEF_SIZE);
+
+ // simplify, one power at a time
+ // largest power (larger powers are at the end)
+ max_power=(*coefficient).denoms[(*coefficient).length-1].power;
+
+ quotient_prev=*coefficient;
+ for(power=max_power;power>=1;power--){
+ ret=coefficient_simplify_fraction(constant, quotient_prev, &remainder, &quotient);
+
+ // if fail to simplify, stop
+ if(ret<0){
+ if(power<max_power){
+ coefficient_concat_noinit(quotient_prev, &out);
+ }
+ else{
+ coefficient_concat(quotient_prev, &out);
+ }
+ break;
+ }
+
+ // add to output
+ coefficient_concat_noinit(remainder, &out);
+ }
+ // if the factorization always succeeded
+ if(max_power>=1 && power==0){
+ coefficient_concat_noinit(quotient, &out);
+ }
+
+ coefficient_simplify(&out);
+
+ // set coefficient to out
+ free_Coefficient(*coefficient);
+ *coefficient=out;
+
+ return 0;
+}
+
+// put all terms under a common denominator
+// only supports coefficients with only one constant
+int coefficient_common_denominator(Coefficient constant, Coefficient* coefficient){
+ int max_power;
+ int i,j;
+ Coefficient tmp;
+ Coefficient out;
+ Coefficient* C_n;
+
+ init_Coefficient(&out, COEF_SIZE);
+
+ // largest power (larger powers are at the end)
+ max_power=(*coefficient).denoms[(*coefficient).length-1].power;
+
+ // store powers of the constant
+ C_n=calloc(sizeof(Coefficient), max_power-1);
+ for(i=0;i<max_power-1;i++){
+ // start from previous product
+ if(i==0){
+ coefficient_cpy(constant, C_n+i);
+ }
+ else{
+ coefficient_cpy(C_n[i-1], C_n+i);
+ }
+ // multiply by constant
+ coefficient_prod_chain(constant, C_n+i);
+ }
+
+ // multiply each term
+ for (i=0;i<(*coefficient).length;i++){
+ init_Coefficient(&tmp, COEF_SIZE);
+ // start with numerator
+ coefficient_append_noinit((*coefficient).factors[i], (*coefficient).nums[i], (*coefficient).denoms[i], &tmp);
+ // multiply
+ if((*coefficient).denoms[i].power<max_power){
+ if((*coefficient).denoms[i].power==max_power-1){
+ coefficient_prod_chain(constant, &tmp);
+ }
+ else{
+ coefficient_prod_chain(C_n[max_power-(*coefficient).denoms[i].power-2], &tmp);
+ }
+ }
+
+ // set denom
+ for(j=0;j<tmp.length;j++){
+ tmp.denoms[j].power=max_power;
+ }
+
+ // add to out
+ coefficient_concat_noinit(tmp, &out);
+ }
+
+ // free C_n
+ for(i=0;i<max_power-1;i++){
+ free_Coefficient(C_n[i]);
+ }
+ free(C_n);
+
+ // free coefficient vectors
+ free((*coefficient).factors);
+ free((*coefficient).nums);
+ free((*coefficient).denoms);
+
+ // set output
+ *coefficient=out;
+
+ coefficient_simplify(coefficient);
+
+ return(0);
+}
+
+// simplify coefficient / constant
+// returns both the remainder and the quotient
+// assumes both coefficient and constant are ordered with the highest order terms last
+int coefficient_simplify_fraction(Coefficient constant, Coefficient coefficient, Coefficient* remainder, Coefficient* out){
+ Coefficient tmp;
+ int step_counter=0;
+ int max_order;
+ int i,j,k;
+ Int_Array rfactors;
+
+ if(constant.length==0){
+ // nothing to do
+ return 0;
+ }
+
+ coefficient_cpy(coefficient, remainder);
+ init_Coefficient(out, COEF_SIZE);
+
+ // continue until (*remainder) is of lower order than constant
+ while((*remainder).length>0 && (*remainder).factors[(*remainder).length-1].length>=constant.factors[constant.length-1].length){
+ step_counter++;
+
+ // interrupt if too long
+ if(step_counter>=coefficient.length*100){
+ free_Coefficient(*remainder);
+ free_Coefficient(*out);
+ return -1;
+ }
+
+ // try to find a term in the constant that divides the last term of the (*remainder)
+ rfactors=(*remainder).factors[(*remainder).length-1];
+
+ // highest order in constant
+ max_order=constant.factors[constant.length-1].length;
+
+ // start from one of the highest order term and look for a common factor
+ for(i=constant.length-1; i>=0; i--){
+ // fail: no highest order terms have been matched
+ if(constant.factors[i].length<max_order){
+ free_Coefficient(*remainder);
+ free_Coefficient(*out);
+ return -2;
+ }
+
+ // check whether the term can be a factor of the last term of the (*remainder)
+ if(int_array_is_subarray_ordered(constant.factors[i], rfactors)==1){
+ // extract the factors that are not in constant
+ init_Coefficient(&tmp, constant.length);
+
+ // init with one term
+ tmp.length=1;
+ init_Int_Array(tmp.factors,MONOMIAL_SIZE);
+
+ for(j=0,k=0;j<rfactors.length;j++){
+ // check that index is not in constant
+ if(k<constant.factors[i].length){
+ if(rfactors.values[j]!=constant.factors[i].values[k]){
+ int_array_append(rfactors.values[j],tmp.factors);
+ }
+ else{
+ // move to next term in constant
+ k++;
+ }
+ }
+ }
+
+ // numerical prefactor: term in the (*remainder) / term in the constant
+ number_quot((*remainder).nums[(*remainder).length-1], constant.nums[i], tmp.nums);
+ // denominator (dummy)
+ tmp.denoms[0]=(*remainder).denoms[(*remainder).length-1];
+
+ // add to out
+ coefficient_concat(tmp, out);
+
+ // multiply by -1
+ Q minus_1;
+ minus_1.numerator=-1;
+ minus_1.denominator=1;
+ number_Qprod_chain(minus_1, tmp.nums);
+
+ // multiply by constant
+ coefficient_prod_chain(constant, &tmp);
+
+ // add to remainder
+ coefficient_concat(tmp, remainder);
+
+ // free memory
+ free_Coefficient(tmp);
+
+ // simplify
+ coefficient_simplify(remainder);
+
+ break;
+ }
+ }
+ }
+
+ // success!
+ // decrease power of constant
+ for(i=0;i<(*out).length;i++){
+ (*out).denoms[i].power=(*out).denoms[i].power-1;
+ }
+
+ return(0);
+}
+
// sort the terms in an equation (quicksort algorithm)
int sort_coefficient(Coefficient* coefficient, int begin, int end){
int i;
@@ -251,7 +480,7 @@ int exchange_coefficient_terms(int i, int j, Coefficient* coefficient){
return(0);
}
-// derive a coefficient with respect to an index
+// differentiate a coefficient with respect to an index
int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
int i,j;
// temp list of indices
@@ -325,7 +554,7 @@ int coefficient_deriv(Coefficient input, int index, Coefficient* output){
}
/*
-// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p )
+// differentiate a coefficient with respect to an index (as a polynomial) (does not differentiate the 1/(1+C)^p )
int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
int i;
// temp list of indices
@@ -365,7 +594,7 @@ int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){
return(0);
}
-// derive a monomial with respect to an index
+// differentiate a monomial with respect to an index
int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){
int j;
@@ -414,14 +643,14 @@ int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output){
int_array_concat(coef2.factors[j],&factor);
// don't throw an error if the power is 0
- if(coef2.denoms[i].power==0){
- coef2.denoms[i].index=coef1.denoms[i].index;
+ if(coef2.denoms[j].power==0){
+ coef2.denoms[j].index=coef1.denoms[i].index;
}
else if(coef1.denoms[i].power==0){
- coef1.denoms[i].index=coef2.denoms[i].index;
+ coef1.denoms[i].index=coef2.denoms[j].index;
}
if(coef1.denoms[i].index!=coef2.denoms[j].index){
- fprintf(stderr,"error: cannot multiply flow equations with different constants\n");
+ fprintf(stderr,"error: cannot multiply flow equations with different constants: got %d and %d\n", coef1.denoms[i].index, coef2.denoms[j].index);
exit(-1);
}
denom=coef1.denoms[i];