diff options
Diffstat (limited to 'src/grouped_polynomial.c')
| -rw-r--r-- | src/grouped_polynomial.c | 138 |
1 files changed, 114 insertions, 24 deletions
diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c index 0d4d75d..9e245a8 100644 --- a/src/grouped_polynomial.c +++ b/src/grouped_polynomial.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. @@ -205,9 +205,9 @@ int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomi if(index==-2){ fprintf(stderr,"error: monomial ("); - for(j=0;j<polynomial.monomials[i].length;j++){ - fprintf(stderr,"%d", polynomial.monomials[i].values[j]); - if(j<polynomial.monomials[i].length-1){ + for(j=0;j<remainder.monomials[i].length;j++){ + fprintf(stderr,"%d", remainder.monomials[i].values[j]); + if(j<remainder.monomials[i].length-1){ fprintf(stderr,","); } } @@ -436,7 +436,7 @@ int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){ } -// derive a flow equation with respect to an unknown variable +// differentiate a flow equation with respect to an unknown variable // equivalent to DB.dl where dl are symbols for the derivatives of the indices in the flow equation with respect to the unknown variable // indices specifies the list of indices that depend on the variable int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow){ @@ -487,7 +487,7 @@ int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Gr /* -// derive a flow equation with respect to an index +// differentiate a flow equation with respect to an index int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){ int i,k; // temp list of indices @@ -603,7 +603,7 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre init_Char_Array(&buffer, STR_SIZE); coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre); if(buffer.length>0){ - printf("%s",buffer.str); + printf("%s",char_array_to_str_noinit(&buffer)); } free_Char_Array(buffer); @@ -732,29 +732,33 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output) } -// eValuate an equation on a vector -int evaleq(RCC* rccs, Grouped_Polynomial poly){ +// evaluate an equation on a vector +int evaleq(RCC out, RCC in, Grouped_Polynomial poly){ int i; - long double* res=calloc((*rccs).length,sizeof(long double)); + long double* res=calloc(out.length,sizeof(long double)); - if((*rccs).length!=poly.length){ - fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length); + if(in.length!=poly.length){ + fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length); + exit(-1); + } + if(out.length!=poly.length){ + fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length); exit(-1); } - // initialize vectors to 0 - for(i=0;i<(*rccs).length;i++){ + // initialize vectors to 0 in an auxiliary vector (to allow for out=in without interference) + for(i=0;i<in.length;i++){ res[i]=0.; } // for each equation for(i=0;i<poly.length;i++){ - evalcoef(*rccs, poly.coefs[i], res+i); + evalcoef(in, poly.coefs[i], res+i); } // copy res to rccs - for(i=0;i<(*rccs).length;i++){ - (*rccs).values[i]=res[i]; + for(i=0;i<out.length;i++){ + out.values[i]=res[i]; } // free memory @@ -763,25 +767,29 @@ int evaleq(RCC* rccs, Grouped_Polynomial poly){ } // evaluate an equation on a vector (using mpfr floats) -int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){ +int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly){ int i; mpfr_t* res; - if((*rccs).length!=poly.length){ - fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length); + if(in.length!=poly.length){ + fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length); + exit(-1); + } + if(out.length!=poly.length){ + fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length); exit(-1); } - res=calloc((*rccs).length,sizeof(mpfr_t)); + res=calloc(out.length,sizeof(mpfr_t)); // for each equation for(i=0;i<poly.length;i++){ - evalcoef_mpfr(*rccs, poly.coefs[i], res[i]); + evalcoef_mpfr(in, poly.coefs[i], res[i]); } // copy res to rccs - for(i=0;i<(*rccs).length;i++){ - mpfr_set((*rccs).values[i], res[i], MPFR_RNDN); + for(i=0;i<out.length;i++){ + mpfr_set(out.values[i], res[i], MPFR_RNDN); mpfr_clear(res[i]); } @@ -791,3 +799,85 @@ int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){ } +// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2) +int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out){ + if(flow1.length!=flow2.length){ + fprintf(stderr, "error: trying to compose two flow equations of different size\n"); + exit(-1); + } + int i,j,k; + Coefficient constant; + + // init + init_Grouped_Polynomial(out, flow1.length); + (*out).length=flow1.length; + + // init constant (so we can tell when the constant was not found) + constant.length=0; + + // loop over rcc's + for(i=0;i<flow1.length;i++){ + // set indices + (*out).indices[i]=flow1.indices[i]; + + // passthrough constant terms + if((*out).indices[i]<0){ + int index=intlist_find_err(flow2.indices,flow2.length,(*out).indices[i]); + coefficient_cpy(flow2.coefs[index], (*out).coefs+i); + constant=flow2.coefs[index]; + continue; + } + + // init + init_Coefficient((*out).coefs+i, COEF_SIZE); + + // loop over terms + for(j=0;j<flow1.coefs[i].length;j++){ + Coefficient tmp_coef; + + // init + init_Coefficient(&tmp_coef, COEF_SIZE); + // init factor + Int_Array tmp_factor; + init_Int_Array(&tmp_factor, MONOMIAL_SIZE); + // init denom + coef_denom denom; + // index should be that appearing in flow2 + if(flow2.coefs[i].length<1){ + fprintf(stderr,"error: composing two flow equations: the %d-th term in the flow equation is empty\n",flow1.indices[i]); + exit(-1); + } + denom.index=flow2.coefs[i].denoms[0].index; + denom.power=0; + // init num + Number tmp_num; + number_cpy(flow1.coefs[i].nums[j], &tmp_num); + + // init coefficient with numerical prefactor + coefficient_append_noinit(tmp_factor, tmp_num, denom, &tmp_coef); + + // loop over factors + for(k=0;k<flow1.coefs[i].factors[j].length;k++){ + // multiply factors together + coefficient_prod_chain(flow2.coefs[intlist_find_err(flow2.indices,flow2.length,flow1.coefs[i].factors[j].values[k])], &tmp_coef); + } + + // add to out + coefficient_concat_noinit(tmp_coef, (*out).coefs+i); + } + + } + + // simplify fractions + if(constant.length!=0){ + for(i=0;i<(*out).length;i++){ + if((*out).indices[i]>=0){ + // reduce them to a common denominator (not much is gained from trying to simplify them) + coefficient_common_denominator(constant, (*out).coefs+i); + //coefficient_simplify_rational(constant, (*out).coefs+i); + } + } + } + + return(0); +} |