Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
committerIan Jauslin <ian.jauslin@roma1.infn.it>2015-06-14 00:52:45 +0000
commitaa0f3ae2988d372b190b9bde2e75a6d17e744e93 (patch)
tree14482245c2fca27fcdad3078e97d0871352d52a7 /src/coefficient.c
Initial commitv1.2
Diffstat (limited to 'src/coefficient.c')
-rw-r--r--src/coefficient.c739
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);
+}