diff options
Diffstat (limited to 'src/coefficient.c')
| -rw-r--r-- | src/coefficient.c | 739 |
1 files changed, 739 insertions, 0 deletions
diff --git a/src/coefficient.c b/src/coefficient.c new file mode 100644 index 0000000..c26b668 --- /dev/null +++ b/src/coefficient.c @@ -0,0 +1,739 @@ +/* +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 "coefficient.h" + +#include <stdio.h> +#include <stdlib.h> +#include "definitions.cpp" +#include "rational.h" +#include "istring.h" +#include "array.h" +#include "number.h" +#include "tools.h" + + +// allocate memory +int init_Coefficient(Coefficient* coef,int size){ + (*coef).factors=calloc(size,sizeof(Int_Array)); + (*coef).nums=calloc(size,sizeof(Number)); + (*coef).denoms=calloc(size,sizeof(coef_denom)); + (*coef).length=0; + (*coef).memory=size; + + return(0); +} + +// free memory +int free_Coefficient(Coefficient coef){ + int i; + for(i=0;i<coef.length;i++){ + free_Int_Array(coef.factors[i]); + free_Number(coef.nums[i]); + } + free(coef.factors); + free(coef.nums); + free(coef.denoms); + + return(0); +} + +// copy a coefficient +int coefficient_cpy(Coefficient input, Coefficient* output){ + init_Coefficient(output,input.length); + coefficient_cpy_noinit(input,output); + return(0); +} +int coefficient_cpy_noinit(Coefficient input, Coefficient* output){ + int i; + + // if output does not have enough memory allocated to it + if(input.length>(*output).memory){ + resize_Coefficient(output,input.length); + } + + for(i=0;i<input.length;i++){ + int_array_cpy(input.factors[i],(*output).factors+i); + number_cpy(input.nums[i],(*output).nums+i); + (*output).denoms[i]=input.denoms[i]; + } + (*output).length=input.length; + return(0); +} + + +// resize the memory allocated to a coefficient +int resize_Coefficient(Coefficient* coefficient,int new_size){ + Coefficient new_coef; + int i; + + init_Coefficient(&new_coef,new_size); + for(i=0;i<(*coefficient).length;i++){ + new_coef.factors[i]=(*coefficient).factors[i]; + new_coef.nums[i]=(*coefficient).nums[i]; + new_coef.denoms[i]=(*coefficient).denoms[i]; + } + new_coef.length=(*coefficient).length; + + free((*coefficient).factors); + free((*coefficient).nums); + free((*coefficient).denoms); + + *coefficient=new_coef; + return(0); +} + + +// append an element to a coefficient +int coefficient_append(Int_Array factor,Number num, coef_denom denom, Coefficient* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy and allocate + int_array_cpy(factor,(*output).factors+offset); + number_cpy(num,(*output).nums+offset); + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} +// append an element to a coefficient without allocating memory +int coefficient_append_noinit(Int_Array factor, Number num, coef_denom denom, Coefficient* output){ + int offset=(*output).length; + + if((*output).length>=(*output).memory){ + resize_Coefficient(output,2*(*output).memory+1); + } + // copy without allocating + (*output).factors[offset]=factor; + (*output).nums[offset]=num; + (*output).denoms[offset]=denom; + // increment length + (*output).length++; + return(0); +} + +// concatenate coefficients and simplify result +int coefficient_concat(Coefficient input, Coefficient* output){ + int i; + for(i=0;i<input.length;i++){ + coefficient_append(input.factors[i],input.nums[i],input.denoms[i],output); + } + + coefficient_simplify(output); + return(0); +} +int coefficient_concat_noinit(Coefficient input, Coefficient* output){ + int i; + for(i=0;i<input.length;i++){ + coefficient_append_noinit(input.factors[i],input.nums[i],input.denoms[i],output); + } + + coefficient_simplify(output); + + // free input arrays + free(input.factors); + free(input.nums); + free(input.denoms); + return(0); +} + + +// simplify a Coefficient +int coefficient_simplify(Coefficient* coefficient){ + int i; + Coefficient output; + init_Coefficient(&output,(*coefficient).length); + // the combination of numerical factors + Number new_num; + init_Number(&new_num,NUMBER_SIZE); + + // sort the factors in the coefficient + for(i=0;i<(*coefficient).length;i++){ + int_array_sort((*coefficient).factors[i],0,(*coefficient).factors[i].length-1); + } + // in order to perform a simplification, the list of terms must be + // sorted (so that terms that are proportional are next to each other) + sort_coefficient(coefficient,0,(*coefficient).length-1); + + for(i=0;i<(*coefficient).length;i++){ + // if the term actually exists + if(number_is_zero((*coefficient).nums[i])!=1){ + // combine numerical factors + number_add_chain((*coefficient).nums[i],&new_num); + } + // if the number is 0, the previous terms that may have the same factors should still be added, hence the 'if' ends here + + // if the factor is different from the next then add term + if(i==(*coefficient).length-1 || (int_array_cmp((*coefficient).factors[i],(*coefficient).factors[i+1])!=0) || coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[i+1])!=0){ + // check that the coefficient is not trivial + if(number_is_zero(new_num)!=1){ + coefficient_append((*coefficient).factors[i],new_num,(*coefficient).denoms[i],&output); + } + + // reset new numerical factor + free_Number(new_num); + init_Number(&new_num,NUMBER_SIZE); + } + } + + free_Number(new_num); + free_Coefficient(*coefficient); + *coefficient=output; + return(0); +} + +// sort the terms in an equation (quicksort algorithm) +int sort_coefficient(Coefficient* coefficient, int begin, int end){ + int i; + int index; + // the pivot: middle of the array + int pivot=(begin+end)/2; + // if the array is non trivial + if(begin<end){ + // send pivot to the end + exchange_coefficient_terms(pivot,end,coefficient); + // loop over the others + for(i=begin, index=begin;i<end;i++){ + // compare with pivot + if(coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[end])<0 || ( coef_denom_cmp((*coefficient).denoms[i],(*coefficient).denoms[end])==0 && (int_array_cmp((*coefficient).factors[i],(*coefficient).factors[end])<0)) ){ + // if smaller, exchange with reference index + exchange_coefficient_terms(i,index,coefficient); + // move reference index + index++; + } + } + // put pivot (which we had sent to the end) in the right place + exchange_coefficient_terms(index,end,coefficient); + // recurse + sort_coefficient(coefficient, begin, index-1); + sort_coefficient(coefficient, index+1, end); + } + return(0); +} + +// exchange two terms (for the sorting algorithm) +int exchange_coefficient_terms(int i, int j, Coefficient* coefficient){ + Int_Array ptmp; + Number tmpq; + coef_denom tmpc; + + ptmp=(*coefficient).factors[j]; + (*coefficient).factors[j]=(*coefficient).factors[i]; + (*coefficient).factors[i]=ptmp; + + tmpq=(*coefficient).nums[j]; + (*coefficient).nums[j]=(*coefficient).nums[i]; + (*coefficient).nums[i]=tmpq; + + tmpc=(*coefficient).denoms[j]; + (*coefficient).denoms[j]=(*coefficient).denoms[i]; + (*coefficient).denoms[i]=tmpc; + return(0); +} + +// derive a coefficient with respect to an index +int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){ + int i,j; + // temp list of indices + Int_Array factor; + // number of times index was found + int match_count; + // whether the output contains at least one factor + int at_least_one=0; + coef_denom denom; + + // loop over monomials + for(i=0;i<input.length;i++){ + init_Int_Array(&factor,input.factors[i].length); + // init match count + match_count=0; + // loop over indices + for(j=0;j<input.factors[i].length;j++){ + // if found + if(input.factors[i].values[j]==index){ + // if it's the first match, don't add it + if(match_count!=0){ + int_array_append(index,&factor); + } + match_count++; + } + // add the index + else{ + int_array_append(input.factors[i].values[j],&factor); + } + } + + denom=input.denoms[i]; + // if the index is that of 1/C + if(index==input.denoms[i].index){ + // if no C in the numerator (which is normal behavior) + if(match_count==0){ + denom.power++; + } + match_count-=input.denoms[i].power; + } + + + // if the derivative doesn't vanish, add it to the coefficient + if(match_count!=0){ + at_least_one=1; + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),input.nums[i]), denom, output); + } + else{ + free_Int_Array(factor); + } + } + + if(at_least_one==1){ + coefficient_simplify(output); + } + 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); + } + + return(0); +} + +int coefficient_deriv(Coefficient input, int index, Coefficient* output){ + init_Coefficient(output, COEF_SIZE); + coefficient_deriv_noinit(input, index, output); + return(0); +} + +/* +// derive a coefficient with respect to an index (as a polynomial) (does not derive the 1/(1+C)^p ) +int coefficient_deriv_noinit(Coefficient input, int index, Coefficient* output){ + int i; + // temp list of indices + Int_Array factor; + // number of times index was found + int match_count; + // whether the output contains at least one factor + int at_least_one=0; + coef_denom denom; + + // loop over monomials + for(i=0;i<input.length;i++){ + // derivative of monomial + monomial_deriv(input.factors[i], index, &factor, &match_count); + + // if the derivative doesn't vanish, add it to the coefficient + if(match_count>0){ + at_least_one=1; + coefficient_append_noinit(factor,number_Qprod_ret(quot(match_count,1),input.nums[i]), input.denoms[i],output); + } + else{ + free_Int_Array(factor); + } + } + + if(at_least_one>0){ + coefficient_simplify(output); + } + 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); + } + + return(0); +} + +// derive a monomial with respect to an index +int monomial_deriv(Int_Array factor, int index, Int_Array* out_factor, int* match_count){ + int j; + + init_Int_Array(out_factor,factor.length); + // init match count + *match_count=0; + + // loop over indices + for(j=0;j<factor.length;j++){ + // if found + if(factor.values[j]==index){ + // if it's the first match, don't add it + if(*match_count!=0){ + int_array_append(index,out_factor); + } + (*match_count)++; + } + // add the index + else{ + int_array_append(factor.values[j],out_factor); + } + } + + return(0); +} +*/ + + + +// product of two coefficients +int coefficient_prod(Coefficient coef1, Coefficient coef2, Coefficient* output){ + int i,j; + // tmp factor + Int_Array factor; + coef_denom denom; + + // init + init_Coefficient(output,COEF_SIZE); + + // loop over coef1 + for(i=0;i<coef1.length;i++){ + // loop over coef2 + for(j=0;j<coef2.length;j++){ + init_Int_Array(&factor,coef1.factors[i].length+coef2.factors[j].length); + int_array_concat(coef1.factors[i],&factor); + 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; + } + else if(coef1.denoms[i].power==0){ + coef1.denoms[i].index=coef2.denoms[i].index; + } + if(coef1.denoms[i].index!=coef2.denoms[j].index){ + fprintf(stderr,"error: cannot multiply flow equations with different constants\n"); + exit(-1); + } + denom=coef1.denoms[i]; + denom.power+=coef2.denoms[j].power; + coefficient_append_noinit(factor,number_prod_ret(coef1.nums[i],coef2.nums[j]), denom, output); + } + } + + // simplify output + coefficient_simplify(output); + return(0); +} + +// product of coefficients, output replaces the second coefficient +int coefficient_prod_chain(Coefficient in, Coefficient* inout){ + Coefficient tmp_coef; + coefficient_prod(in,*inout,&tmp_coef); + free_Coefficient(*inout); + *inout=tmp_coef; + return(0); +} + + +// print coefficient +// offset specifies the amount of blank space to the left of the terms after the first +// prepend indices by ind_pre +int coefficient_sprint(Coefficient coef, Char_Array* output, int offset, char index_pre){ + Char_Array buffer; + int i,j,k; + int dcount; + + if(coef.length==0){ + char_array_snprintf(output, " (0)\n"); + } + + for(i=0;i<coef.length;i++){ + if(i==0){ + char_array_snprintf(output, " "); + } + else{ + for(j=0;j<=offset;j++){ + char_array_append(' ',output); + } + char_array_append('+',output); + } + + // print numerical coefficient + char_array_append('(',output); + init_Char_Array(&buffer, STR_SIZE); + number_sprint(coef.nums[i], &buffer); + char_array_concat(buffer, output); + free_Char_Array(buffer); + char_array_append(')',output); + + // print factors + for(j=0;j<coef.factors[i].length;j++){ + // constant indices + if(coef.factors[i].values[j]<0){ + // count derivatives + dcount=-coef.factors[i].values[j]/DOFFSET; + char_array_append('[',output); + for(k=0;k<dcount;k++){ + char_array_append('d',output); + } + char_array_snprintf(output,"C%d]",-coef.factors[i].values[j]-dcount*DOFFSET); + } + else{ + // count derivatives + dcount=coef.factors[i].values[j]/DOFFSET; + char_array_append('[',output); + for(k=0;k<dcount;k++){ + char_array_append('d',output); + } + char_array_snprintf(output,"%c%d]",index_pre,coef.factors[i].values[j]-dcount*DOFFSET); + } + } + + // print constant denominators + if(coef.denoms[i].power!=0){ + char_array_snprintf(output,"[/C%d^%d]",-coef.denoms[i].index,coef.denoms[i].power); + } + + char_array_append('\n',output); + } + + return(0); +} + + +// read from a string +#define PP_NULL_MODE 0 +#define PP_BRACKET_MODE 1 +#define PP_INDICES_MODE 2 +#define PP_POWER_MODE 3 +#define PP_COMMENT_MODE 4 +#define PP_NUMBER_MODE 5 +#define PP_CONSTANT_MODE 6 +#define PP_CONSTANT_DENOM_MODE 7 +int char_array_to_Coefficient(Char_Array str_coef, Coefficient* output){ + // buffer + char* buffer=calloc(str_coef.length+1,sizeof(char)); + char* buffer_ptr=buffer; + Int_Array indices; + coef_denom denom; + Number num, tmp1_num; + int mode; + int i,j; + int parenthesis_count=0; + int dcount=0; + + // allocate memory + init_Coefficient(output,COEF_SIZE); + + // init + init_Int_Array(&indices, MONOMIAL_SIZE); + num=number_one(); + denom.index=-1; + denom.power=0; + + *buffer_ptr='\0'; + // loop over the input polynomial + // start in null mode + mode=PP_NULL_MODE; + for(j=0;j<str_coef.length;j++){ + if(mode==PP_COMMENT_MODE){ + if(str_coef.str[j]=='\n'){ + mode=PP_NULL_MODE; + } + } + else{ + switch(str_coef.str[j]){ + // new indices + case '+': + if(mode==PP_NULL_MODE){ + coefficient_append_noinit(indices, num, denom, output); + // reset indices, num + init_Int_Array(&indices, MONOMIAL_SIZE); + num=number_one(); + denom.index=-1; + denom.power=0; + } + break; + + // enter indices or power mode + case '[': + if(mode==PP_NULL_MODE){ + mode=PP_BRACKET_MODE; + // reset derivatives count + dcount=0; + } + break; + // indices mode + case '%': + if(mode==PP_BRACKET_MODE){ + mode=PP_INDICES_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + break; + case 'C': + if(mode==PP_BRACKET_MODE){ + mode=PP_CONSTANT_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + break; + case '/': + if(mode==PP_BRACKET_MODE){ + mode=PP_CONSTANT_DENOM_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + else if(mode!=PP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]); + } + break; + // derivatives + case 'd': + if(mode==PP_BRACKET_MODE || mode==PP_INDICES_MODE || mode==PP_CONSTANT_MODE){ + dcount++; + } + break; + // power mode + case '^': + if(mode==PP_CONSTANT_DENOM_MODE){ + sscanf(buffer,"%d",&i); + denom.index=-i; + mode=PP_POWER_MODE; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + else{ + buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]); + } + break; + // read indices or power + case ']': + sscanf(buffer,"%d",&i); + if(mode==PP_INDICES_MODE){ + int_array_append(i+dcount*DOFFSET,&indices); + } + else if(mode==PP_CONSTANT_MODE){ + int_array_append(-i-dcount*DOFFSET,&indices); + } + else if(mode==PP_POWER_MODE){ + denom.power=i; + } + // switch back to null mode + mode=PP_NULL_MODE; + break; + + // numerical pre-factor + case '(': + if(mode==PP_NULL_MODE){ + mode=PP_NUMBER_MODE; + parenthesis_count=0; + buffer_ptr=buffer; + *buffer_ptr='\0'; + } + else if(mode==PP_NUMBER_MODE){ + // match parentheses + parenthesis_count++; + buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]); + } + break; + case ')': + if(mode==PP_NUMBER_MODE){ + if(parenthesis_count==0){ + // write num + str_to_Number(buffer,&tmp1_num); + number_prod_chain(tmp1_num,&num); + free_Number(tmp1_num); + // back to null mode + mode=PP_NULL_MODE; + } + else{ + parenthesis_count--; + buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]); + } + } + break; + + // characters to ignore + case ' ':break; + case '&':break; + case '\n':break; + + // comments + case '#': + mode=PP_COMMENT_MODE; + break; + + default: + if(mode!=PP_NULL_MODE){ + // write to buffer + buffer_ptr=str_addchar(buffer_ptr,str_coef.str[j]); + } + break; + } + } + } + + // last term + coefficient_append_noinit(indices, num, denom, output); + + free(buffer); + return(0); +} + +int str_to_Coefficient(char* str, Coefficient* output){ + Char_Array array; + array.length=str_len(str); + array.str=str; + char_array_to_Coefficient(array, output); + return(0); +} + +// compare coefficient denominators +int coef_denom_cmp(coef_denom denom1, coef_denom denom2){ + if(denom1.index<denom2.index){ + return(1); + } + else if(denom1.index>denom2.index){ + return(-1); + } + + if(denom1.power<denom2.power){ + return(-1); + } + else if(denom1.power>denom2.power){ + return(1); + } + + return(0); +} + + +// evaluate a coefficient on a vector +int evalcoef(RCC rccs, Coefficient coef, long double* out){ + int i,j; + long double num_factor; + + *out=0; + + // for each monomial + for(i=0;i<coef.length;i++){ + // product of factors + for(j=0, num_factor=1.;j<coef.factors[i].length;j++){ + num_factor*=rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.factors[i].values[j])]; + } + // denominator + if(coef.denoms[i].power>0){ + num_factor/=dpower(rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.denoms[i].index)],coef.denoms[i].power); + } + *out+=num_factor*number_double_val(coef.nums[i]); + } + return(0); +} |