path: root/src/coefficient.c

/*
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 <stdarg.h>
// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
#define MPFR_USE_VA_LIST
#include <mpfr.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);
}

// evaluate a coefficient on a vector (using mpfr floats)
int evalcoef_mpfr(RCC_mpfr rccs, Coefficient coef, mpfr_t out){
  int i,j;
  mpfr_t num_factor;
  // tmp number (do not initialize Z)
  mpfr_t x, y, Z;

  // init numbers
  mpfr_inits(num_factor, x, y, (mpfr_ptr) NULL);

  mpfr_init(out);
  mpfr_set_zero(out, 1);

  // for each monomial
  for(i=0;i<coef.length;i++){
    // product of factors
    mpfr_set_flt(num_factor, 1., MPFR_RNDN);
    for(j=0;j<coef.factors[i].length;j++){
      mpfr_mul(x,num_factor,rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.factors[i].values[j])], MPFR_RNDN);
      mpfr_set(num_factor,x, MPFR_RNDN);
    }
    // denominator
    if(coef.denoms[i].power>0){
      mpfr_pow_si(y, rccs.values[intlist_find_err(rccs.indices,rccs.length,coef.denoms[i].index)], coef.denoms[i].power, MPFR_RNDN);
      mpfr_div(x, num_factor, y, MPFR_RNDN);
      mpfr_set(num_factor, x, MPFR_RNDN);
    }

    number_mpfr_val(Z, coef.nums[i]);
    mpfr_mul(x, num_factor, Z, MPFR_RNDN);
    mpfr_add(y, x, out, MPFR_RNDN);
    mpfr_set(out, y, MPFR_RNDN);
    mpfr_clear(Z);
  }

  // free numbers
  mpfr_clears(num_factor, x, y, (mpfr_ptr)NULL);

  return(0);
}