diff options
Diffstat (limited to 'src/coefficient.c')
| -rw-r--r-- | src/coefficient.c | 245 |
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, "ient); + + // 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]; |