Initial commitv1.2

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);
+}