diff options
Diffstat (limited to 'src/grouped_polynomial.c')
| -rw-r--r-- | src/grouped_polynomial.c | 765 |
1 files changed, 765 insertions, 0 deletions
diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c new file mode 100644 index 0000000..b7bea42 --- /dev/null +++ b/src/grouped_polynomial.c @@ -0,0 +1,765 @@ +/* +Copyright 2015 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "grouped_polynomial.h" + +#include <stdio.h> +#include <stdlib.h> +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "coefficient.h" +#include "polynomial.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Grouped_Polynomial(Grouped_Polynomial* gpolynomial, int size){ + (*gpolynomial).coefs=calloc(size,sizeof(Coefficient)); + (*gpolynomial).indices=calloc(size,sizeof(int)); + (*gpolynomial).length=0; + (*gpolynomial).memory=size; + + return(0); +} + +// free memory +int free_Grouped_Polynomial(Grouped_Polynomial gpolynomial){ + int i; + for(i=0;i<gpolynomial.length;i++){ + free_Coefficient(gpolynomial.coefs[i]); + } + free(gpolynomial.coefs); + free(gpolynomial.indices); + + return(0); +} + +// resize the memory allocated to a grouped_polynomial +int resize_grouped_polynomial(Grouped_Polynomial* grouped_polynomial,int new_size){ + Grouped_Polynomial new_poly; + int i; + + init_Grouped_Polynomial(&new_poly,new_size); + for(i=0;i<(*grouped_polynomial).length;i++){ + new_poly.indices[i]=(*grouped_polynomial).indices[i]; + new_poly.coefs[i]=(*grouped_polynomial).coefs[i]; + } + new_poly.length=(*grouped_polynomial).length; + + free((*grouped_polynomial).indices); + free((*grouped_polynomial).coefs); + + *grouped_polynomial=new_poly; + return(0); +} + +// copy a grouped_polynomial +int grouped_polynomial_cpy(Grouped_Polynomial input, Grouped_Polynomial* output){ + init_Grouped_Polynomial(output,input.length); + grouped_polynomial_cpy_noinit(input,output); + return(0); +} +int grouped_polynomial_cpy_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){ + int i; + if((*output).memory<input.length){ + fprintf(stderr,"error: trying to copy a grouped polynomial of length %d to another with memory %d\n",input.length,(*output).memory); + exit(-1); + } + for(i=0;i<input.length;i++){ + (*output).indices[i]=input.indices[i]; + coefficient_cpy(input.coefs[i], (*output).coefs+i); + } + (*output).length=input.length; + + return(0); +} + +// append an element to a grouped_polynomial +int grouped_polynomial_append(int index, Coefficient coef, Grouped_Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy and allocate + (*output).indices[offset]=index; + coefficient_cpy(coef, (*output).coefs+offset); + //increment length + (*output).length++; + + return(0); +} +// append an element to a grouped_polynomial without allocating memory +int grouped_polynomial_append_noinit(int index, Coefficient coef, Grouped_Polynomial* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_grouped_polynomial(output,2*(*output).memory+1); + } + + // copy without allocating + (*output).indices[offset]=index; + (*output).coefs[offset]=coef; + // increment length + (*output).length++; + return(0); +} + +// concatenate two grouped_polynomials +int grouped_polynomial_concat(Grouped_Polynomial input, Grouped_Polynomial* output){ + int i; + for(i=0;i<input.length;i++){ + grouped_polynomial_append(input.indices[i],input.coefs[i],output); + } + return(0); +} +int grouped_polynomial_concat_noinit(Grouped_Polynomial input, Grouped_Polynomial* output){ + int i; + for(i=0;i<input.length;i++){ + grouped_polynomial_append_noinit(input.indices[i],input.coefs[i],output); + } + + // free input arrays + free(input.indices); + free(input.coefs); + return(0); +} + + +// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table +// robust algorithm: allows for a term in the polynomial to contribute to several id table entries +int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable, Fields_Table fields){ + int i,j; + // what is left to group + Polynomial remainder; + int index; + Number ratio; + Number num; + Int_Array factor; + int security=0; + int pos; + coef_denom denom; + + // init remainder + polynomial_cpy(polynomial, &remainder); + + // allocate memory + init_Grouped_Polynomial(grouped_polynomial, idtable.length+1); + + // copy indices from idtable and allocate + // the constant term + (*grouped_polynomial).indices[0]=-1; + init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE); + for(i=1;i<=idtable.length;i++){ + (*grouped_polynomial).indices[i]=idtable.indices[i-1]; + init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE); + } + (*grouped_polynomial).length=idtable.length+1; + + // keep on going as long as there are terms in the remainder + while(remainder.length>0){ + // stop if the number of iterations exceeds 100 times the length of the polynomial + if(security >= 100*polynomial.length){ + fprintf(stderr,"error: polynomial could not be grouped in less than %d groupings\n", 100*polynomial.length); + exit(-1); + } + security++; + + // index of the last element + i=remainder.length-1; + + // find entry + if(remainder.monomials[i].length==0){ + // constant + index=-1; + } + else{ + // loop over entries + for(j=0,index=-2;j<idtable.length && index==-2;j++){ + // loop over terms in the polynomial + for(pos=0;pos<idtable.polynomials[j].length;pos++){ + if(int_array_cmp(idtable.polynomials[j].monomials[pos],remainder.monomials[i])==0){ + index=j; + break; + } + } + } + } + + 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){ + fprintf(stderr,","); + } + } + fprintf(stderr,") not found in idtable\n"); + exit(-1); + } + + // if not constant + if(index>=0){ + ratio=number_quot_ret(remainder.nums[i],idtable.polynomials[index].nums[pos]); + factor=remainder.factors[i]; + // add to coefficient + denom.index=-1; + denom.power=1; + coefficient_append(factor, ratio, denom, (*grouped_polynomial).coefs+index+1); + + // remove from remainder + free_Int_Array(remainder.monomials[i]); + // do not free factor yet + free_Number(remainder.nums[i]); + remainder.length--; + + // add terms from idtable with minus sign + for(j=0;j<idtable.polynomials[index].length;j++){ + if(j!=pos){ + num=number_prod_ret(ratio, idtable.polynomials[index].nums[j]); + number_Qprod_chain(quot(-1,1),&num); + polynomial_append(idtable.polynomials[index].monomials[j], factor, num, &remainder); + free_Number(num); + } + } + + free_Int_Array(factor); + free_Number(ratio); + + // simplify remainder + polynomial_simplify(&remainder, fields); + } + // constant + else if(index==-1){ + // add to coefficient + denom.index=-1; + denom.power=0; + coefficient_append(remainder.factors[i], remainder.nums[i], denom, (*grouped_polynomial).coefs); + // remove from remainder + free_Int_Array(remainder.monomials[i]); + free_Int_Array(remainder.factors[i]); + free_Number(remainder.nums[i]); + remainder.length--; + } + } + + // simplify the result + simplify_grouped_polynomial(grouped_polynomial); + + free_Polynomial(remainder); + return(0); +} + + +// construct a grouped polynomial from a polynomial, grouping together the terms specified in the id table. +// identifies sub-polynomials in the polynomial corresponding to the entire rhs of an entry in the id table. +// requires the polynomial and the idtable to be sorted +// can only treat cases in which monomials in different polynomials of the idtable are distinct +int group_polynomial_pickandchoose(Polynomial polynomial, Grouped_Polynomial* grouped_polynomial, Id_Table idtable){ + int i,j,k; + // a mask specifying which terms of the polynomial have already been grouped + int* mask=calloc(polynomial.length, sizeof(int)); + int index; + Number ratio, ratio_check; + // whether ratio was ever allocated + int alloc_ratio=0; + // whether the correct index was found + int found_index; + int start_index_search; + Int_Array mask_tmp_flips; + coef_denom denom; + + // allocate memory + init_Grouped_Polynomial(grouped_polynomial, idtable.length+1); + + // copy indices from idtable and allocate + // the constant term + (*grouped_polynomial).indices[0]=-1; + init_Coefficient((*grouped_polynomial).coefs, COEF_SIZE); + for(i=1;i<=idtable.length;i++){ + (*grouped_polynomial).indices[i]=idtable.indices[i-1]; + init_Coefficient((*grouped_polynomial).coefs+i, COEF_SIZE); + } + (*grouped_polynomial).length=idtable.length+1; + + // loop over monomials + for(i=0;i<polynomial.length;i++){ + // check that the term hasn't already been added + if(mask[i]==0){ + // loop until the correct index is found (the polynomial must contain all the terms in the index and the numerical factors must match) + found_index=0; + start_index_search=0; + while(found_index==0){ + found_index=1; + // find entry + index=find_id(polynomial.monomials[i], idtable,start_index_search); + // easier to debug if the error is here instead of inside find_id + if(index==-2){ + fprintf(stderr,"error: monomial not found in idtable\n"); + exit(-1); + } + + // if not constant + if(index>=0){ + // a vector in which to store the indices that were masked + init_Int_Array(&mask_tmp_flips,idtable.polynomials[index].length); + + // loop over all monomials in that entry of the idtable + for(j=0;j<idtable.polynomials[index].length && found_index==1;j++){ + // find the monomial in the polynomial + for(k=i;k<polynomial.length;k++){ + // only check if mask==0 + // only check if the factors are correct + if(mask[k]==0 && int_array_cmp(polynomial.factors[i],polynomial.factors[k])==0 && int_array_cmp(idtable.polynomials[index].monomials[j],polynomial.monomials[k])==0){ + ratio_check=number_quot_ret(polynomial.nums[k],idtable.polynomials[index].nums[j]); + + // if the factors don't factor + if(alloc_ratio!=0 && number_compare(ratio,ratio_check)==0){ + found_index=0; + break; + } + // check that ratio was allocated + if(alloc_ratio!=0){ + free_Number(ratio); + } + ratio=ratio_check; + alloc_ratio=1; + + // added to polynomial + mask[k]=1; + // keep track of the flips so that they can be undone if the index turns out to be incorrect + int_array_append(k,&mask_tmp_flips); + break; + } + } + + // error if the monomial could not be found + if(k==polynomial.length){ + found_index=0; + } + } + + // if the index was incorrect + if(found_index==0){ + // reset mask + for(j=0;j<mask_tmp_flips.length;j++){ + mask[mask_tmp_flips.values[j]]=0; + } + // start index search at next item + start_index_search=index+1; + } + else{ + // add to grouped polynomial + denom.index=-1; + denom.power=1; + coefficient_append(polynomial.factors[i], ratio, denom, (*grouped_polynomial).coefs+index+1); + } + + if(alloc_ratio==1){ + free_Number(ratio); + alloc_ratio=0; + } + free_Int_Array(mask_tmp_flips); + } + // constant + else if(index==-1){ + mask[i]=1; + denom.index=-1; + denom.power=0; + coefficient_append(polynomial.factors[i], polynomial.nums[i], denom, (*grouped_polynomial).coefs); + } + } + } + } + + // check all the terms were grouped + for(i=0;i<polynomial.length;i++){ + if(mask[i]==0){ + fprintf(stderr,"error: this polynomial could not be grouped: no matches were found for some of the terms\n"); + exit(-1); + } + } + + free(mask); + return(0); +} + + +// find the entry in the idtable containing monomial +// start search at the specified index +int find_id(Int_Array monomial, Id_Table idtable, int start){ + int i,j; + + // constant + if(monomial.length==0){ + return(-1); + } + + // loop over entries + for(i=start;i<idtable.length;i++){ + // loop over terms in the polynomial + for(j=0;j<idtable.polynomials[i].length;j++){ + if(int_array_cmp(idtable.polynomials[i].monomials[j],monomial)==0){ + return(i); + } + } + } + + return(-2); +} + + +// simplify grouped polynomial +int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){ + int i; + for(i=0;i<(*polynomial).length;i++){ + coefficient_simplify((*polynomial).coefs+i); + } + return(0); +} + + +// derive 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){ + int i,j,k; + Coefficient tmp_coef; + + // alloc + init_Grouped_Polynomial(dflow, flow_equation.length); + + // for each equation + for(i=0;i<flow_equation.length;i++){ + // copy indices + if(flow_equation.indices[i]>=0){ + (*dflow).indices[i]=flow_equation.indices[i]+DOFFSET; + } + else{ + (*dflow).indices[i]=flow_equation.indices[i]-DOFFSET; + } + + init_Coefficient((*dflow).coefs+i, COEF_SIZE); + // for each index + for(j=0;j<indices.length;j++){ + coefficient_deriv(flow_equation.coefs[i], indices.values[j], &tmp_coef); + // multiply each coefficient by the appropriate dl[j] + for(k=0;k<tmp_coef.length;k++){ + // only in non-trivial cases + if(number_is_zero(tmp_coef.nums[k])==0){ + // non-constants + if(indices.values[j]>=0){ + int_array_append(DOFFSET + indices.values[j], tmp_coef.factors+k); + } + // constants are offset with -doffset (so that the derivatives of constants also have a negative index) + else{ + int_array_append(-DOFFSET + indices.values[j], tmp_coef.factors+k); + } + } + } + + // add to output + coefficient_concat_noinit(tmp_coef, (*dflow).coefs+i); + } + } + + (*dflow).length=flow_equation.length; + + return(0); +} + + +/* +// derive 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 + Int_Array factor; + // number of times index was found + int match_count; + coef_denom denom; + // store the computation of the derivative of the constant + int previous_constant_index=0; + Coefficient dC; + Coefficient tmp_coef; + + init_Grouped_Polynomial(output, flow_equation.length); + + // loop over equations + for(k=0;k<flow_equation.length;k++){ + init_Coefficient((*output).coefs+k, COEF_SIZE); + + // loop over monomials + for(i=0;i<flow_equation.coefs[k].length;i++){ + + // derivative of the numerator + monomial_deriv(flow_equation.coefs[k].factors[i], index, &factor, &match_count); + + // if the derivative doesn't vanish, add it to the coefficient + if(match_count>0){ + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], (*output).coefs+k); + } + else{ + free_Int_Array(factor); + } + + // derivative of the denominator + if(flow_equation.coefs[k].denoms[i].power>0){ + // check whether the derivative was already computed + if(flow_equation.coefs[k].denoms[i].index!=previous_constant_index){ + // if not first, then free + if(previous_constant_index!=0){ + free_Coefficient(dC); + previous_constant_index=0; + } + init_Coefficient(&dC,COEF_SIZE); + coefficient_deriv_noinit(flow_equation.coefs[intlist_find_err(flow_equation.indices, flow_equation.length, flow_equation.coefs[k].denoms[i].index)], index, &dC); + previous_constant_index=flow_equation.coefs[k].denoms[i].index; + } + + init_Coefficient(&tmp_coef, dC.length); + coefficient_append(flow_equation.coefs[k].factors[i], number_Qprod_ret(quot(-flow_equation.coefs[k].denoms[i].power,1), flow_equation.coefs[k].nums[i]), flow_equation.coefs[k].denoms[i], &tmp_coef); + (tmp_coef.denoms[0].power)++; + + coefficient_prod_chain(dC, &tmp_coef); + + coefficient_concat_noinit(tmp_coef, (*output).coefs+k); + } + } + + // memory safe + if((*output).coefs[k].length>0){ + coefficient_simplify((*output).coefs+k); + } + else{ + // add a trivial element to the coefficient + init_Int_Array(&factor,0); + denom.index=-1; + denom.power=0; + coefficient_append_noinit(factor,number_zero(),denom,(*output).coefs+k); + } + } + + free_Coefficient(dC); + return(0); +} +*/ + + +// print a grouped polynomial +// prepend the indices on the left side with lhs_pre, and those on the right by rhs_pre +int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre, char rhs_pre){ + int i,j; + Char_Array buffer; + int dcount; + + // for each equation + for(i=0;i<grouped_polynomial.length;i++){ + //print lhs + // negative indices are constants + if(grouped_polynomial.indices[i]<0){ + // count derivatives + dcount=-grouped_polynomial.indices[i]/DOFFSET; + for(j=0;j<3-dcount;j++){ + printf(" "); + } + printf("["); + for(j=0;j<dcount;j++){ + printf("d"); + } + printf("C%d] =",-grouped_polynomial.indices[i]-dcount*DOFFSET); + } + else{ + // count derivatives + dcount=grouped_polynomial.indices[i]/DOFFSET; + for(j=0;j<2-dcount;j++){ + printf(" "); + } + printf("["); + for(j=0;j<dcount;j++){ + printf("d"); + } + printf("%c%2d] =",lhs_pre,grouped_polynomial.indices[i]-dcount*DOFFSET); + } + + // rhs + init_Char_Array(&buffer, STR_SIZE); + coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre); + if(buffer.length>0){ + printf("%s",buffer.str); + } + free_Char_Array(buffer); + + if(i<grouped_polynomial.length-1){ + printf(","); + } + // extra \n + printf("\n"); + } + return(0); +} + +// read from string +#define PP_NULL_MODE 0 +#define PP_COEF_MODE 1 +#define PP_INDEX_MODE 3 +#define PP_COMMENT_MODE 4 +#define PP_BRACKET_MODE 5 +#define PP_CONSTANT_MODE 6 +int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output){ + // buffer + char* buffer=calloc(str.length+1,sizeof(char)); + char* buffer_ptr=buffer; + int index=-2; + Coefficient coef; + int i,j; + int mode; + int dcount=0; + + init_Grouped_Polynomial(output, EQUATION_SIZE); + + // loop over input + mode=PP_NULL_MODE; + for(j=0;j<str.length;j++){ + if(mode==PP_COMMENT_MODE){ + if(str.str[j]=='\n'){ + mode=PP_NULL_MODE; + } + } + // stay in polynomial mode until ',' + else if(mode==PP_COEF_MODE){ + if(str.str[j]==','){ + // parse polynomial + str_to_Coefficient(buffer, &coef); + // write index and polynomial + grouped_polynomial_append_noinit(index, coef, output); + mode=PP_NULL_MODE; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,str.str[j]); + } + } + else{ + switch(str.str[j]){ + // index + case '[': + if(mode==PP_NULL_MODE){ + mode=PP_BRACKET_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + // reset derivatives count + dcount=0; + } + break; + case '%': + if(mode==PP_BRACKET_MODE){ + mode=PP_INDEX_MODE; + } + break; + // constant term + case 'C': + if(mode==PP_BRACKET_MODE){ + mode=PP_CONSTANT_MODE; + } + break; + // derivatives + case 'd': + if(mode==PP_BRACKET_MODE || mode==PP_INDEX_MODE || mode==PP_CONSTANT_MODE){ + dcount++; + } + break; + // write index + case ']': + sscanf(buffer,"%d",&i); + if(mode==PP_INDEX_MODE){ + index=i+dcount*DOFFSET; + } + else if(mode==PP_CONSTANT_MODE){ + index=-i-dcount*DOFFSET; + } + mode=PP_NULL_MODE; + + break; + + // coef mode + case '=': + if(mode==PP_NULL_MODE){ + buffer_ptr=buffer; + *buffer_ptr='\0'; + mode=PP_COEF_MODE; + } + break; + + // comment + case '#': + mode=PP_COMMENT_MODE; + break; + + default: + if(mode!=PP_NULL_MODE){ + buffer_ptr=str_addchar(buffer_ptr,str.str[j]); + } + break; + } + } + } + + // last step + if(mode==PP_COEF_MODE){ + str_to_Coefficient(buffer, &coef); + grouped_polynomial_append_noinit(index, coef, output); + } + + free(buffer); + return(0); +} + + +// evaluate an equation on a vector +int evaleq(RCC* rccs, Grouped_Polynomial poly){ + int i; + long double* res=calloc((*rccs).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); + exit(-1); + } + + // initialize vectors to 0 + for(i=0;i<(*rccs).length;i++){ + res[i]=0.; + } + + // for each equation + for(i=0;i<poly.length;i++){ + evalcoef(*rccs, poly.coefs[i], res+i); + } + + // copy res to rccs + for(i=0;i<(*rccs).length;i++){ + (*rccs).values[i]=res[i]; + } + + // free memory + free(res); + return(0); + +} + |