Initial commitv1.0

50 files changed, 4402 insertions, 0 deletions

diff --git a/src/array.c b/src/array.c
new file mode 100644
index 0000000..772cb39
--- /dev/null
+++ b/src/array.c
@@ -0,0 +1,227 @@
+/*
+Copyright 2016 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 "array.h"
+#include <stdlib.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include "errors.h"
+#include "polynomial.h"
+#include "utils.h"
+
+//--------------------------------------------------
+//
+//  integer array (array_int)
+//
+//--------------------------------------------------
+#include "array_int.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  unsigned integer array (array_uint)
+//
+//--------------------------------------------------
+#include "array_uint.h"
+#include "array_base.c"
+
+//--------------------------------------------------
+//
+//  double array (array_double)
+//
+//--------------------------------------------------
+#include "array_double.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  long double array (array_ldouble)
+//
+//--------------------------------------------------
+#include "array_ldouble.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  MPFR array (array_mpfr)
+//
+//--------------------------------------------------
+#include "array_mpfr.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  array of MPFR arrays (array_2_mpfr)
+//
+//--------------------------------------------------
+#include "array_2_mpfr.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  character array (array_char)
+//
+//--------------------------------------------------
+#include "array_char.h"
+#include "array_base.c"
+
+// append a char*
+int array_char_append_str(char* str, array_char* output){
+  char* ptr;
+  for (ptr=str;*ptr!='\0';ptr++){
+    array_char_append(*ptr, output);
+  }
+  return(0);
+}
+
+// convert to char*
+int array_char_to_str(array_char input, char** output){
+  unsigned int i;
+  (*output)=calloc(input.length+1,sizeof(char));
+  for(i=0;i<input.length;i++){
+    (*output)[i]=input.values[i];
+  }
+  if((*output)[input.length-1]!='\0'){
+    (*output)[input.length]='\0';
+  }
+  return(0);
+}
+// noinit (changes the size of input if needed)
+char* array_char_to_str_noinit(array_char* input){
+  if(input->values[input->length-1]!='\0'){
+    if(input->length==input->memory){
+      array_char_resize(input,input->length+1);
+    }
+    // add final '\0'
+    input->values[input->length]='\0';
+  }
+  return(input->values);
+}
+
+// convert from char*
+int str_to_array_char(char* str, array_char* output){
+  char* ptr;
+  unsigned int str_len=0;
+  for(ptr=str;*ptr!='\0';ptr++){
+    str_len++;
+  }
+  array_char_init(output, str_len);
+  for(ptr=str;*ptr!='\0';ptr++){
+    array_char_append(*ptr,output);
+  }
+  return(0);
+}
+
+// compare an array_char and a char*
+int array_char_cmp_str(array_char array_char, char* str){
+  unsigned int j;
+  for(j=0;j<array_char.length && str[j]!='\0';j++){
+    if(array_char.values[j]!=str[j]){
+      return(0);
+    }
+  }
+  if(j==array_char.length && str[j]=='\0'){
+    return(1);
+  }
+  return(0);
+}
+
+// format strings
+int array_char_snprintf(array_char* output, char* fmt, ...){
+  size_t size=100;
+  unsigned int extra_size;
+  char* out_str=calloc(size,sizeof(char));
+  char* ptr;
+  va_list vaptr;
+
+  // initialize argument list starting after fmt
+  va_start(vaptr, fmt);
+  // print format
+  extra_size=vsnprintf(out_str, size, fmt, vaptr);
+  va_end(vaptr);
+
+  // if too large
+  if(extra_size>size){
+    // resize
+    free(out_str);
+    // +1 for '\0'
+    size=extra_size+1;
+    out_str=calloc(size,sizeof(char));
+    // read format again
+    va_start(vaptr, fmt);
+    vsnprintf(out_str,size,fmt,vaptr);
+    va_end(vaptr);
+  }
+
+  // write to char array
+  for(ptr=out_str;*ptr!='\0';ptr++){
+    array_char_append(*ptr, output);
+  }
+
+  free(out_str);
+
+  return(0);
+}
+
+
+//--------------------------------------------------
+//
+//  string array (array_str)
+//
+//--------------------------------------------------
+#include "array_str.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  mpfr polynomial array (array_polynomial_mpfr)
+//
+//--------------------------------------------------
+#include "array_polynomial_mpfr.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  double polynomial array (array_polynomial_double)
+//
+//--------------------------------------------------
+#include "array_polynomial_double.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  long double polynomial array (array_polynomial_ldouble)
+//
+//--------------------------------------------------
+#include "array_polynomial_ldouble.h"
+#include "array_base.c"
+
+
+//--------------------------------------------------
+//
+//  pthread_t array (array_pthread_t)
+//
+//--------------------------------------------------
+#include "array_pthread_t.h"
+#include "array_base.c"
diff --git a/src/array.h b/src/array.h
new file mode 100644
index 0000000..63639a7
--- /dev/null
+++ b/src/array.h
@@ -0,0 +1,150 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Mutable arrays
+*/
+
+#ifndef LIBINUM_ARRAY_H
+#define LIBINUM_ARRAY_H
+
+#include "types.h"
+
+//--------------------------------------------------
+//
+//  integer array (array_int)
+//
+//--------------------------------------------------
+#include "array_int.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  unsigned integer array (array_uint)
+//
+//--------------------------------------------------
+#include "array_uint.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  double array (array_double)
+//
+//--------------------------------------------------
+#include "array_double.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  long double array (array_ldouble)
+//
+//--------------------------------------------------
+#include "array_ldouble.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  MPFR array (array_mpfr)
+//
+//--------------------------------------------------
+#include "array_mpfr.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  array of MPFR arrays (array_2_mpfr)
+//
+//--------------------------------------------------
+#include "array_2_mpfr.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  character array (array_char)
+//
+//--------------------------------------------------
+#include "array_char.h"
+#include "array_base.h"
+
+// append a char*
+int array_char_append_str(char* str, array_char* output);
+
+// convert to char*
+int array_char_to_str(array_char input, char** output);
+// noinit (changes the size of input if needed)
+char* array_char_to_str_noinit(array_char* input);
+
+// convert from char*
+int str_to_array_char(char* str, array_char* output);
+
+// compare an array_char and a char*
+int array_char_cmp_str(array_char array_char, char* str);
+
+// format strings
+int array_char_snprintf(array_char* output, char* fmt, ...);
+
+
+//--------------------------------------------------
+//
+//  string array (array_str)
+//
+//--------------------------------------------------
+#include "array_str.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  mpfr polynomial array (array_polynomial_mpfr)
+//
+//--------------------------------------------------
+#include "array_polynomial_mpfr.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  double polynomial array (array_polynomial_double)
+//
+//--------------------------------------------------
+#include "array_polynomial_double.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  long double polynomial array (array_polynomial_ldouble)
+//
+//--------------------------------------------------
+#include "array_polynomial_ldouble.h"
+#include "array_base.h"
+
+
+//--------------------------------------------------
+//
+//  pthread_t array (array_pthread_t)
+//
+//--------------------------------------------------
+#include "array_pthread_t.h"
+#include "array_base.h"
+
+#endif
diff --git a/src/array_2_mpfr.h b/src/array_2_mpfr.h
new file mode 100644
index 0000000..afd99d9
--- /dev/null
+++ b/src/array_2_mpfr.h
@@ -0,0 +1,54 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for arrays of MPFR float arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_2_mpfr
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_2_mpfr_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE array_mpfr
+// free
+#define ARRAY_VAL_FREE(VAL) array_mpfr_free(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL=NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) array_mpfr_cpy(FROMVAL, &TOVAL);
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) array_mpfr_cmp(VAL1, VAL2)==0
+#define ARRAY_VAL_IFLT(VAL1, VAL2) array_mpfr_cmp(VAL1, VAL2)<0
+#define ARRAY_VAL_IFGT(VAL1, VAL2) array_mpfr_cmp(VAL1, VAL2)>0
+// print values
+#define ARRAY_VAL_PRINT(VAL) array_mpfr_print(VAL)
+
diff --git a/src/array_base.c b/src/array_base.c
new file mode 100644
index 0000000..a70bd86
--- /dev/null
+++ b/src/array_base.c
@@ -0,0 +1,249 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for arrays
+
+  see polynomial_*.h for the values taken by ARRAY_FUNC, etc...
+*/
+
+
+// init
+int ARRAY_FUNC(init) (ARRAY_TYPENAME* array, unsigned int memory){
+  array->values=calloc(memory,sizeof(ARRAY_VAL_TYPE));
+  array->memory=memory;
+  array->length=0;
+  return(0);
+}
+int ARRAY_FUNC(free) (ARRAY_TYPENAME array){
+  #ifdef ARRAY_VAL_FREE
+    unsigned int i;
+    for(i=0;i<array.length;i++){
+      ARRAY_VAL_FREE(array.values[i]);
+    }
+  #endif
+  free(array.values);
+  return(0);
+}
+// do not free values, only free calloc'ed memory
+int ARRAY_FUNC(free_vects) (ARRAY_TYPENAME array){
+  free(array.values);
+  return(0);
+}
+
+// resize memory
+int ARRAY_FUNC(resize) (ARRAY_TYPENAME* array, unsigned int newsize){
+  unsigned int i;
+  ARRAY_TYPENAME new_array;
+  ARRAY_FUNC(init) (&new_array, newsize);
+  for(i=0;i<array->length;i++){
+    ARRAY_VAL_SET(new_array.values[i], array->values[i]);
+  }
+  free(array->values);
+  *array=new_array;
+  return(0);
+}
+
+// add a value
+int ARRAY_FUNC(append) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output){
+  if(output->length >= output->memory){
+    ARRAY_FUNC(resize) (output,2*output->memory+1);
+  }
+  ARRAY_VAL_CPY(val, output->values[output->length]);
+  output->length++;
+  return(0);
+}
+// do not copy the value, instead let the last element of the array point to 'val'
+// only if the value requires initializing
+#ifdef ARRAY_VAL_FREE
+int ARRAY_FUNC(append_noinit) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output){
+  if(output->length >= output->memory){
+    ARRAY_FUNC(resize) (output,2*output->memory+1);
+  }
+  ARRAY_VAL_SET(output->values[output->length], val);
+  output->length++;
+  return(0);
+}
+#endif
+// add a value only if it is not already present
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(append_unique) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output){
+  if(ARRAY_FUNC(find) (val, *output)<0){
+    ARRAY_FUNC(append) (val, output);
+  }
+  return(0);
+}
+#endif
+
+// copy
+int ARRAY_FUNC(cpy) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output){
+  ARRAY_FUNC(init) (output, input.length);
+  ARRAY_FUNC(cpy_noinit) (input, output);
+  return(0);
+}
+// do not init output array
+int ARRAY_FUNC(cpy_noinit) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output){
+  unsigned int i;
+  if(output->memory<input.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<input.length;i++){
+    ARRAY_VAL_CPY(input.values[i], output->values[i]);
+  }
+  output->length=input.length;
+  return(0);
+}
+
+
+// concatenate
+int ARRAY_FUNC(concat) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output){
+  unsigned int i;
+  for(i=0;i<input.length;i++){
+    ARRAY_FUNC(append) (input.values[i], output);
+  }
+  return(0);
+}
+
+// concat but only add values that are not already present in the array
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(concat_unique) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output){
+  unsigned int i;
+  for(i=0;i<input.length;i++){
+    ARRAY_FUNC(append_unique) (input.values[i], output);
+  }
+  return(0);
+}
+#endif
+
+// sub-array
+int ARRAY_FUNC(subarray) (ARRAY_TYPENAME array, unsigned int begin, unsigned int end, ARRAY_TYPENAME* subarray){
+  unsigned int i;
+  if(begin>end || end>=array.length){
+    return(LIBINUM_ERROR_ILLEGAL_MEMORY_ACCESS);
+  }
+  ARRAY_FUNC(init) (subarray,end-begin);
+  for(i=begin;i<=end;i++){
+    ARRAY_FUNC(append) (array.values[i], subarray);
+  }
+  return(0);
+}
+
+// find (does not assume the array is sorted)
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(find) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME array){
+  unsigned int i;
+  for(i=0;i<array.length;i++){
+    if(ARRAY_VAL_IFEQ(array.values[i], val)){
+      return(i);
+    }
+  }
+  return(-1);
+}
+#endif
+
+// sort (quicksort)
+#ifdef ARRAY_VAL_IFLT
+int ARRAY_FUNC(sort) (ARRAY_TYPENAME array){
+  if(array.length>1){
+    ARRAY_FUNC(sort_sub) (array, 0, array.length-1);
+  }
+  return(0);
+}
+// sort a sub-array
+int ARRAY_FUNC(sort_sub) (ARRAY_TYPENAME array, unsigned int begin, unsigned int end){
+  unsigned int i;
+  unsigned int index;
+  ARRAY_VAL_TYPE tmp;
+  // the pivot: middle of the array
+  unsigned int pivot=(begin+end)/2;
+
+  // if the array is non trivial
+  if(begin<end){
+    // send pivot to the end
+    ARRAY_VAL_SET(tmp, array.values[pivot]);
+    ARRAY_VAL_SET(array.values[pivot], array.values[end]);
+    ARRAY_VAL_SET(array.values[end], tmp);
+
+    // loop over the others
+    for(i=begin, index=begin;i<end;i++){
+      // compare with pivot
+      if(ARRAY_VAL_IFLT(array.values[i], array.values[end])){
+	// if smaller, exchange with reference index
+	ARRAY_VAL_SET(tmp, array.values[index]);
+	ARRAY_VAL_SET(array.values[index], array.values[i]);
+	ARRAY_VAL_SET(array.values[i], tmp);
+	// move reference index
+	index++;
+      }
+    }
+    // put pivot (which we had sent to the end) in the right place
+    ARRAY_VAL_SET(tmp, array.values[end]);
+    ARRAY_VAL_SET(array.values[end], array.values[index]);
+    ARRAY_VAL_SET(array.values[index], tmp);
+
+    // recurse
+    if(index>0){
+      ARRAY_FUNC(sort_sub) (array, begin, index-1);
+    }
+    ARRAY_FUNC(sort_sub) (array, index+1, end);
+  }
+  return(0);
+}
+#endif
+
+// compare arrays
+#if defined ARRAY_VAL_IFLT && defined ARRAY_VAL_IFGT
+int ARRAY_FUNC(cmp) (ARRAY_TYPENAME array1, ARRAY_TYPENAME array2){
+  unsigned int i;
+
+  // compare lengths
+  if(array1.length<array2.length){
+    return(-1);
+  }
+  if(array1.length>array2.length){
+    return(1);
+  }
+
+  // compare terms
+  for(i=0;i<array1.length;i++){
+    if(ARRAY_VAL_IFLT(array1.values[i], array2.values[i])){
+      return(-1);
+    }
+    if(ARRAY_VAL_IFGT(array1.values[i], array2.values[i])){
+      return(1);
+    }
+  }
+
+  // if equal
+  return(0);
+}
+#endif
+
+// print array
+#ifdef ARRAY_VAL_PRINT
+int ARRAY_FUNC(print) (ARRAY_TYPENAME array){
+  unsigned int i;
+  printf("(");
+  for(i=0;i<array.length;i++){
+    ARRAY_VAL_PRINT(array.values[i]);
+    if(i<array.length-1){
+      printf(", ");
+    }
+  }
+  printf(")");
+  return(0);
+}
+#endif
diff --git a/src/array_base.h b/src/array_base.h
new file mode 100644
index 0000000..b9de5d2
--- /dev/null
+++ b/src/array_base.h
@@ -0,0 +1,78 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for arrays
+
+  see polynomial_*.h for the values taken by ARRAY_FUNC, etc...
+*/
+
+// init
+int ARRAY_FUNC(init) (ARRAY_TYPENAME* array, unsigned int memory);
+int ARRAY_FUNC(free) (ARRAY_TYPENAME array);
+// do not free values, only free calloc'ed memory
+int ARRAY_FUNC(free_vects) (ARRAY_TYPENAME array);
+
+// resize memory
+int ARRAY_FUNC(resize) (ARRAY_TYPENAME* array, unsigned int newsize);
+
+// add a value
+int ARRAY_FUNC(append) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output);
+#ifdef ARRAY_VAL_FREE
+// do not copy the value, instead let the last element of the array point to 'val'
+int ARRAY_FUNC(append_noinit) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output);
+#endif
+// add a value only if it is not already present
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(append_unique) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME* output);
+#endif
+
+// copy
+int ARRAY_FUNC(cpy) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output);
+int ARRAY_FUNC(cpy_noinit) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output);
+
+// concatenate
+int ARRAY_FUNC(concat) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output);
+// concat but only add values that are not already present in the array
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(concat_unique) (ARRAY_TYPENAME input, ARRAY_TYPENAME* output);
+#endif
+
+// sub-array
+int ARRAY_FUNC(subarray) (ARRAY_TYPENAME array, unsigned int begin, unsigned int end, ARRAY_TYPENAME* subarray);
+
+// find (does not assume the array is sorted)
+#ifdef ARRAY_VAL_IFEQ
+int ARRAY_FUNC(find) (ARRAY_VAL_TYPE val, ARRAY_TYPENAME array);
+#endif
+
+// sort (quicksort)
+#ifdef ARRAY_VAL_IFLT
+int ARRAY_FUNC(sort) (ARRAY_TYPENAME array);
+// sort a sub-array
+int ARRAY_FUNC(sort_sub) (ARRAY_TYPENAME array, unsigned int begin, unsigned int end);
+#endif
+
+// compare arrays
+#if defined ARRAY_VAL_IFLT && defined ARRAY_VAL_IFGT
+int ARRAY_FUNC(cmp) (ARRAY_TYPENAME array1, ARRAY_TYPENAME array2);
+#endif
+
+// print array
+#ifdef ARRAY_VAL_PRINT
+int ARRAY_FUNC(print) (ARRAY_TYPENAME array);
+#endif
+
diff --git a/src/array_char.h b/src/array_char.h
new file mode 100644
index 0000000..6c95174
--- /dev/null
+++ b/src/array_char.h
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for char arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_char
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_char_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE char
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+// print values
+#define ARRAY_VAL_PRINT(VAL) printf("%c",VAL)
+
diff --git a/src/array_double.h b/src/array_double.h
new file mode 100644
index 0000000..8d3c187
--- /dev/null
+++ b/src/array_double.h
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for double arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_double
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_double_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE double
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+// print values
+#define ARRAY_VAL_PRINT(VAL) fprint_double(stdout, VAL)
+
diff --git a/src/array_int.h b/src/array_int.h
new file mode 100644
index 0000000..9496b89
--- /dev/null
+++ b/src/array_int.h
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for integer arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_int
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_int_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE int
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+// print values
+#define ARRAY_VAL_PRINT(VAL) printf("%d",VAL)
+
diff --git a/src/array_ldouble.h b/src/array_ldouble.h
new file mode 100644
index 0000000..d314906
--- /dev/null
+++ b/src/array_ldouble.h
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for long double arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_ldouble
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_ldouble_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE long double
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+// print values
+#define ARRAY_VAL_PRINT(VAL) fprint_ldouble(stdout, VAL)
+
diff --git a/src/array_mpfr.h b/src/array_mpfr.h
new file mode 100644
index 0000000..46bc54b
--- /dev/null
+++ b/src/array_mpfr.h
@@ -0,0 +1,54 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for MPFR float arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_mpfr
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_mpfr_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE mpfr_t
+// free
+#define ARRAY_VAL_FREE(VAL) mpfr_clear(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL[0]=NEWVAL[0]
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) mpfr_init(TOVAL); mpfr_set(TOVAL, FROMVAL, MPFR_RNDN)
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) mpfr_cmp(VAL1, VAL2)==0
+#define ARRAY_VAL_IFLT(VAL1, VAL2) mpfr_cmp(VAL1, VAL2)<0
+#define ARRAY_VAL_IFGT(VAL1, VAL2) mpfr_cmp(VAL1, VAL2)>0
+// print values
+#define ARRAY_VAL_PRINT(VAL) fprint_mpfr(stdout, VAL)
+
diff --git a/src/array_polynomial_double.h b/src/array_polynomial_double.h
new file mode 100644
index 0000000..d25e689
--- /dev/null
+++ b/src/array_polynomial_double.h
@@ -0,0 +1,50 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for arrays of polynomials with double coefficients
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_polynomial_double
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_polynomial_double_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE polynomial_double
+// free
+#define ARRAY_VAL_FREE(VAL) polynomial_double_free(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) polynomial_double_cpy(FROMVAL, &TOVAL)
+// print values
+#define ARRAY_VAL_PRINT(VAL) polynomial_double_print(VAL)
+
diff --git a/src/array_polynomial_ldouble.h b/src/array_polynomial_ldouble.h
new file mode 100644
index 0000000..49c81cc
--- /dev/null
+++ b/src/array_polynomial_ldouble.h
@@ -0,0 +1,50 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for arrays of polynomials with long double coefficients
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_polynomial_ldouble
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_polynomial_ldouble_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE polynomial_ldouble
+// free
+#define ARRAY_VAL_FREE(VAL) polynomial_ldouble_free(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) polynomial_ldouble_cpy(FROMVAL, &TOVAL)
+// print values
+#define ARRAY_VAL_PRINT(VAL) polynomial_ldouble_print(VAL)
+
diff --git a/src/array_polynomial_mpfr.h b/src/array_polynomial_mpfr.h
new file mode 100644
index 0000000..8cfe6e8
--- /dev/null
+++ b/src/array_polynomial_mpfr.h
@@ -0,0 +1,50 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for arrays of polynomials with mpfr coefficients
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_polynomial_mpfr
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_polynomial_mpfr_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE polynomial_mpfr
+// free
+#define ARRAY_VAL_FREE(VAL) polynomial_mpfr_free(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) polynomial_mpfr_cpy(FROMVAL, &TOVAL)
+// print values
+#define ARRAY_VAL_PRINT(VAL) polynomial_mpfr_print(VAL)
+
diff --git a/src/array_pthread_t.h b/src/array_pthread_t.h
new file mode 100644
index 0000000..e200c5d
--- /dev/null
+++ b/src/array_pthread_t.h
@@ -0,0 +1,50 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for pthread_t arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_pthread_t
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_pthread_t_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE pthread_t
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+
diff --git a/src/array_str.h b/src/array_str.h
new file mode 100644
index 0000000..2d35d1b
--- /dev/null
+++ b/src/array_str.h
@@ -0,0 +1,54 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for string arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_str
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_str_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE array_char
+// free
+#define ARRAY_VAL_FREE(VAL) array_char_free(VAL)
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) array_char_cpy(FROMVAL, &TOVAL)
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) array_char_cmp(VAL1, VAL2)==0
+#define ARRAY_VAL_IFLT(VAL1, VAL2) array_char_cmp(VAL1, VAL2)<0
+#define ARRAY_VAL_IFGT(VAL1, VAL2) array_char_cmp(VAL1, VAL2)>0
+// print values
+#define ARRAY_VAL_PRINT(VAL) array_char_print(VAL)
+
diff --git a/src/array_type.h b/src/array_type.h
new file mode 100644
index 0000000..4aa63fc
--- /dev/null
+++ b/src/array_type.h
@@ -0,0 +1,28 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Structure definition for arrays
+
+  see array_*.h for the values taken by ARRAY_TYPENAME, etc...
+*/
+
+typedef struct ARRAY_TYPENAME {
+  ARRAY_VAL_TYPE* values;
+  unsigned int length;
+  unsigned int memory;
+} ARRAY_TYPENAME;
+
diff --git a/src/array_uint.h b/src/array_uint.h
new file mode 100644
index 0000000..b3cec85
--- /dev/null
+++ b/src/array_uint.h
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for unsigned int arrays
+*/
+
+
+// reset CPP macros
+#undef ARRAY_TYPENAME
+#undef ARRAY_FUNC
+#undef ARRAY_VAL_TYPE
+#undef ARRAY_VAL_FREE
+#undef ARRAY_VAL_SET
+#undef ARRAY_VAL_CPY
+#undef ARRAY_VAL_IFEQ
+#undef ARRAY_VAL_IFLT
+#undef ARRAY_VAL_IFGT
+#undef ARRAY_VAL_PRINT
+
+
+// name of the polynomial type
+#define ARRAY_TYPENAME array_uint
+// prefix of function names
+#define ARRAY_FUNC(NAME) array_uint_ ## NAME
+
+// type of values
+#define ARRAY_VAL_TYPE unsigned int
+// set values
+#define ARRAY_VAL_SET(VAL, NEWVAL) VAL = NEWVAL
+// copy values
+#define ARRAY_VAL_CPY(FROMVAL, TOVAL) TOVAL = FROMVAL
+// compare values
+#define ARRAY_VAL_IFEQ(VAL1, VAL2) VAL1==VAL2
+#define ARRAY_VAL_IFLT(VAL1, VAL2) VAL1<VAL2
+#define ARRAY_VAL_IFGT(VAL1, VAL2) VAL1>VAL2
+// print values
+#define ARRAY_VAL_PRINT(VAL) printf("%u",VAL)
+
diff --git a/src/errors.h b/src/errors.h
new file mode 100644
index 0000000..ee26dd2
--- /dev/null
+++ b/src/errors.h
@@ -0,0 +1,29 @@
+/*
+Copyright 2016 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.
+*/
+
+#ifndef LIBINUM_ERRORS_H
+#define LIBINUM_ERRORS_H
+
+// maximal number of iterations reached
+#define LIBINUM_ERROR_MAXITER -1
+// illegal memoru access
+#define LIBINUM_ERROR_ILLEGAL_MEMORY_ACCESS -2
+// size incompatibilities
+#define LIBINUM_ERROR_SIZE_MISMATCH -3
+// NaN or Inf
+#define LIBINUM_ERROR_NAN -4
+
+#endif
diff --git a/src/integral.c b/src/integral.c
new file mode 100644
index 0000000..43ccd82
--- /dev/null
+++ b/src/integral.c
@@ -0,0 +1,52 @@
+/*
+Copyright 2016 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 "integral.h"
+
+#include <mpfr.h>
+#include <math.h>
+#include <stdlib.h>
+#include "root.h"
+#include "array.h"
+#include "polynomial.h"
+#include "errors.h"
+
+
+//--------------------------------------------------
+//
+//  using doubles
+//
+//--------------------------------------------------
+#include "integral_double.h"
+#include "integral_base.c"
+
+
+//--------------------------------------------------
+//
+//  using long doubles
+//
+//--------------------------------------------------
+#include "integral_ldouble.h"
+#include "integral_base.c"
+
+
+//--------------------------------------------------
+//
+//  using mpfr floats
+//
+//--------------------------------------------------
+#include "integral_mpfr.h"
+#include "integral_base.c"
diff --git a/src/integral.h b/src/integral.h
new file mode 100644
index 0000000..01c1fe4
--- /dev/null
+++ b/src/integral.h
@@ -0,0 +1,53 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Integration of real functions
+*/
+
+#ifndef LIBINUM_INTEGRAL_H
+#define LIBINUM_INTEGRAL_H
+
+#include "types.h"
+
+
+//--------------------------------------------------
+//
+//  using doubles
+//
+//--------------------------------------------------
+#include "integral_double.h"
+#include "integral_base.h"
+
+
+//--------------------------------------------------
+//
+//  using long doubles
+//
+//--------------------------------------------------
+#include "integral_ldouble.h"
+#include "integral_base.h"
+
+//--------------------------------------------------
+//
+//  using mpfr floats
+//
+//--------------------------------------------------
+
+#include "integral_mpfr.h"
+#include "integral_base.h"
+
+#endif
diff --git a/src/integral_base.c b/src/integral_base.c
new file mode 100644
index 0000000..da3f05a
--- /dev/null
+++ b/src/integral_base.c
@@ -0,0 +1,506 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for integration
+
+  see integral_*.h for the values taken by INTEGRAL_FUNC, etc...
+*/
+
+
+// compute the integral of a real function of 1 variable using the Gauss-Legendre scheme
+int INTEGRAL_FUNC(integrate_gauss) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, void* extra_args){
+  unsigned int i;
+  INTEGRAL_FLOAT_TYPE valf, delta, avg, x;
+  int ret;
+
+  // check arguments
+  if(abcissa.length != weights.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  #ifdef INTEGRAL_FLOAT_INIT
+    INTEGRAL_FLOAT_INIT(valf);
+    INTEGRAL_FLOAT_INIT(delta);
+    INTEGRAL_FLOAT_INIT(avg);
+    INTEGRAL_FLOAT_INIT(x);
+  #endif
+
+  // init to 0
+  INTEGRAL_FLOAT_SET_UI(*out, 0);
+
+  // half length of interval
+  INTEGRAL_FLOAT_SUB(delta, upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(delta, delta, 2);
+  // average of interval
+  INTEGRAL_FLOAT_ADD(avg, upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(avg, avg, 2);
+
+  for(i=0;i<abcissa.length;i++){
+    // evaluate at x
+    INTEGRAL_FLOAT_MUL(x, delta, abcissa.values[i]);
+    INTEGRAL_FLOAT_ADD(x, x, avg);
+
+    ret=(*func)(&valf, x, extra_args);
+    if(ret<0){
+      #ifdef INTEGRAL_FLOAT_FREE
+	INTEGRAL_FLOAT_FREE(valf);
+	INTEGRAL_FLOAT_FREE(delta);
+	INTEGRAL_FLOAT_FREE(avg);
+	INTEGRAL_FLOAT_FREE(x);
+      #endif
+      return(ret);
+    }
+    // check whether valf is a number
+    if(! INTEGRAL_FLOAT_ISNUMBER(valf)){
+      #ifdef INTEGRAL_FLOAT_FREE
+	INTEGRAL_FLOAT_FREE(valf);
+	INTEGRAL_FLOAT_FREE(delta);
+	INTEGRAL_FLOAT_FREE(avg);
+	INTEGRAL_FLOAT_FREE(x);
+      #endif
+      return(LIBINUM_ERROR_NAN);
+    }
+      
+    INTEGRAL_FLOAT_MUL(valf, valf, weights.values[i]);
+    INTEGRAL_FLOAT_ADD(*out, *out, valf);
+  }
+
+  INTEGRAL_FLOAT_MUL(*out, *out, delta);
+
+  #ifdef INTEGRAL_FLOAT_FREE
+    INTEGRAL_FLOAT_FREE(valf);
+    INTEGRAL_FLOAT_FREE(delta);
+    INTEGRAL_FLOAT_FREE(avg);
+    INTEGRAL_FLOAT_FREE(x);
+  #endif
+
+  return(0);
+}
+
+// multithreaded version
+// arguments to pass to each thread
+struct INTEGRAL_FUNC(pthread_integrate_gauss_args) {
+  // partial sum of the values assigned to the thread
+  INTEGRAL_FLOAT_TYPE out;
+  // abcissa
+  INTEGRAL_FLOATARRAY_TYPE x;
+  // weights
+  INTEGRAL_FLOATARRAY_TYPE weights;
+  // pointer to the function to evaluate
+  int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*);
+  // extra arguments passed to func
+  void* extra_args;
+  // return value
+  int ret;
+};
+int INTEGRAL_FUNC(integrate_gauss_multithread) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, void* extra_args, unsigned int threads, array_pthread_t* thread_ids){
+  unsigned int i;
+  INTEGRAL_FLOAT_TYPE x, delta, avg;
+  struct INTEGRAL_FUNC(pthread_integrate_gauss_args) args[threads];
+  int ret=0;
+  unsigned int thread_nr;
+
+  array_pthread_t_init(thread_ids, threads);
+  thread_ids->length=threads;
+
+  // check arguments
+  if(abcissa.length != weights.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  #ifdef INTEGRAL_FLOAT_INIT
+    INTEGRAL_FLOAT_INIT(delta);
+    INTEGRAL_FLOAT_INIT(avg);
+    INTEGRAL_FLOAT_INIT(x);
+  #endif
+
+  // init to 0
+  INTEGRAL_FLOAT_SET_UI(*out, 0);
+
+  // half length of interval
+  INTEGRAL_FLOAT_SUB(delta, upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(delta, delta, 2);
+  // average of interval
+  INTEGRAL_FLOAT_ADD(avg, upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(avg, avg, 2);
+
+  // inits
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    INTEGRAL_FLOATARRAY_FUNC(init) (&(args[thread_nr].x),abcissa.length/threads+1);
+    INTEGRAL_FLOATARRAY_FUNC(init) (&(args[thread_nr].weights),abcissa.length/threads+1);
+    #ifdef INTEGRAL_FLOAT_INIT
+      INTEGRAL_FLOAT_INIT(args[thread_nr].out);
+    #endif
+  }
+
+  // set abcissa and weights
+  for(i=0,thread_nr=0;i<abcissa.length;i++,thread_nr=(thread_nr+1)%threads){
+    INTEGRAL_FLOAT_MUL(x, delta, abcissa.values[i]);
+    INTEGRAL_FLOAT_ADD(x, x, avg);
+
+    INTEGRAL_FLOATARRAY_FUNC(append) (x, &(args[thread_nr].x));
+    INTEGRAL_FLOATARRAY_FUNC(append) (weights.values[i], &(args[thread_nr].weights));
+  }
+
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    // set func
+    args[thread_nr].func=func;
+    // set extra_args
+    args[thread_nr].extra_args=extra_args;
+    // init ret
+    args[thread_nr].ret=0;
+    // run threads
+    pthread_create(thread_ids->values+thread_nr, NULL, INTEGRAL_FUNC(integrate_gauss_thread), (void*)(args+thread_nr));
+  }
+
+  // wait for completion and join
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    pthread_join(thread_ids->values[thread_nr], NULL);
+
+    if(args[thread_nr].ret<0){
+      ret=args[thread_nr].ret;
+    }
+    else{
+      INTEGRAL_FLOAT_ADD(*out, *out, args[thread_nr].out);
+    }
+  }
+
+  // multiply by size of interval
+  INTEGRAL_FLOAT_MUL(*out, *out, delta);
+
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    INTEGRAL_FLOATARRAY_FUNC(free) (args[thread_nr].x);
+    INTEGRAL_FLOATARRAY_FUNC(free) (args[thread_nr].weights);
+    #ifdef INTEGRAL_FLOAT_FREE
+      INTEGRAL_FLOAT_FREE(args[thread_nr].out);
+    #endif
+  }
+
+  #ifdef INTEGRAL_FLOAT_FREE
+    INTEGRAL_FLOAT_FREE(delta);
+    INTEGRAL_FLOAT_FREE(avg);
+    INTEGRAL_FLOAT_FREE(x);
+  #endif
+
+
+  return(ret);
+}
+// per-thread function
+void* INTEGRAL_FUNC(integrate_gauss_thread) (void* args){
+  unsigned int i;
+  INTEGRAL_FLOAT_TYPE valf;
+  int ret;
+  struct INTEGRAL_FUNC(pthread_integrate_gauss_args)* argument=((struct INTEGRAL_FUNC(pthread_integrate_gauss_args)*)args);
+  #ifdef INTEGRAL_FLOAT_INIT
+    INTEGRAL_FLOAT_INIT(valf);
+  #endif
+
+  INTEGRAL_FLOAT_SET_UI(argument->out, 0);
+
+  for(i=0;i<argument->x.length;i++){
+    // evaluate
+    ret=(*(argument->func))(&valf, argument->x.values[i], argument->extra_args);
+
+    if(ret<0){
+      #ifdef INTEGRAL_FLOAT_FREE
+	INTEGRAL_FLOAT_FREE(valf);
+      #endif
+      argument->ret=ret;
+      return(NULL);
+    }
+    // check whether valf is a number
+    if(! INTEGRAL_FLOAT_ISNUMBER(valf)){
+      #ifdef INTEGRAL_FLOAT_FREE
+	INTEGRAL_FLOAT_FREE(valf);
+      #endif
+      argument->ret=LIBINUM_ERROR_NAN;
+      return(NULL);
+    }
+
+    INTEGRAL_FLOAT_MUL(valf, valf, argument->weights.values[i]);
+    INTEGRAL_FLOAT_ADD(argument->out, argument->out, valf);
+  }
+  #ifdef INTEGRAL_FLOAT_FREE
+    INTEGRAL_FLOAT_FREE(valf);
+  #endif
+
+  return(NULL);
+}
+
+
+
+// smart management of temporary variables: initialize as many as needed but allow them to be re-used instead of freeing them
+#ifdef INTEGRAL_FLOAT_INIT
+int INTEGRAL_FUNC(integrate_gauss_smarttmp) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, INTEGRAL_FLOATARRAY_TYPE* tmps, void* extra_args){
+  unsigned int i;
+  int ret;
+
+  // check arguments
+  if(abcissa.length != weights.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  // allocate tmps if needed
+  if(tmps->memory<4){
+    // no need to resize since the values should not be kept
+    INTEGRAL_FLOATARRAY_FUNC(free)(*tmps);
+    INTEGRAL_FLOATARRAY_FUNC(init)(tmps, 4);
+  }
+  for (i=tmps->length;i<4;i++){
+    INTEGRAL_FLOAT_INIT(tmps->values[i]);
+    (tmps->length)++;
+  }
+
+  // init to 0
+  INTEGRAL_FLOAT_SET_UI(*out, 0);
+
+  // half length of interval
+  INTEGRAL_FLOAT_SUB(tmps->values[0], upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(tmps->values[0], tmps->values[0], 2);
+  // average of interval
+  INTEGRAL_FLOAT_ADD(tmps->values[1], upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(tmps->values[1], tmps->values[1], 2);
+
+  for(i=0;i<abcissa.length;i++){
+    // evaluation point
+    INTEGRAL_FLOAT_MUL(tmps->values[2], tmps->values[0], abcissa.values[i]);
+    INTEGRAL_FLOAT_ADD(tmps->values[2], tmps->values[2], tmps->values[1]);
+
+    ret=(*func)(tmps->values+3, tmps->values[2], extra_args);
+    if(ret<0){
+      return(ret);
+    }
+    // check whether tmps->values[3] is a number
+    if(! INTEGRAL_FLOAT_ISNUMBER(tmps->values[3])){
+      return(LIBINUM_ERROR_NAN);
+    }
+      
+    INTEGRAL_FLOAT_MUL(tmps->values[3], tmps->values[3], weights.values[i]);
+    INTEGRAL_FLOAT_ADD(*out, *out, tmps->values[3]);
+  }
+
+  INTEGRAL_FLOAT_MUL(*out, *out, tmps->values[0]);
+
+  return(0);
+}
+#endif
+
+// multithreaded version with smart management of temporary variables (only for datatypes that need to be initialized)
+#ifdef INTEGRAL_FLOAT_INIT
+// arguments to pass to each thread
+struct INTEGRAL_FUNC(pthread_integrate_gauss_smarttmp_args) {
+  // partial sum of the values assigned to the thread
+  INTEGRAL_FLOAT_TYPE* out;
+  // abcissa
+  INTEGRAL_FLOATARRAY_TYPE x;
+  // weights
+  INTEGRAL_FLOATARRAY_TYPE weights;
+  // tmp
+  INTEGRAL_FLOAT_TYPE* tmp;
+  // pointer to the function to evaluate
+  int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*);
+  // extra arguments passed to func
+  void* extra_args;
+  // return value
+  int ret;
+};
+int INTEGRAL_FUNC(integrate_gauss_smarttmp_multithread) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, INTEGRAL_FLOATARRAY_TYPE* tmps, void* extra_args, unsigned int threads, array_pthread_t* thread_ids){
+  unsigned int i;
+  struct INTEGRAL_FUNC(pthread_integrate_gauss_smarttmp_args) args[threads];
+  int ret=0;
+  unsigned int thread_nr;
+
+  array_pthread_t_init(thread_ids, threads);
+  thread_ids->length=threads;
+
+  // check arguments
+  if(abcissa.length != weights.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  // allocate tmps if needed
+  if(tmps->memory<2+2*threads+abcissa.length){
+    // no need to resize since the values should not be kept
+    INTEGRAL_FLOATARRAY_FUNC(free)(*tmps);
+    INTEGRAL_FLOATARRAY_FUNC(init)(tmps, 2+2*threads+abcissa.length);
+  }
+  for (i=tmps->length;i<2+2*threads+abcissa.length;i++){
+    INTEGRAL_FLOAT_INIT(tmps->values[i]);
+    (tmps->length)++;
+  }
+
+
+  // init to 0
+  INTEGRAL_FLOAT_SET_UI(*out, 0);
+
+  // half length of interval
+  INTEGRAL_FLOAT_SUB(tmps->values[0], upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(tmps->values[0], tmps->values[0], 2);
+  // average of interval
+  INTEGRAL_FLOAT_ADD(tmps->values[1], upper, lower);
+  INTEGRAL_FLOAT_DIV_UI(tmps->values[1], tmps->values[1], 2);
+
+  // inits
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    INTEGRAL_FLOATARRAY_FUNC(init) (&(args[thread_nr].x),abcissa.length/threads+1);
+    INTEGRAL_FLOATARRAY_FUNC(init) (&(args[thread_nr].weights),abcissa.length/threads+1);
+    args[thread_nr].out=tmps->values+2+thread_nr;
+    args[thread_nr].tmp=tmps->values+2+threads+thread_nr;
+  }
+
+  // set abcissa and weights
+  for(i=0,thread_nr=0;i<abcissa.length;i++,thread_nr=(thread_nr+1)%threads){
+    INTEGRAL_FLOAT_MUL(tmps->values[2+2*threads+i], tmps->values[0], abcissa.values[i]);
+    INTEGRAL_FLOAT_ADD(tmps->values[2+2*threads+i], tmps->values[2+2*threads+i], tmps->values[1]);
+
+    INTEGRAL_FLOATARRAY_FUNC(append_noinit) (tmps->values[2+2*threads+i], &(args[thread_nr].x));
+    INTEGRAL_FLOATARRAY_FUNC(append_noinit) (weights.values[i], &(args[thread_nr].weights));
+  }
+
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    // set func
+    args[thread_nr].func=func;
+    // set extra_args
+    args[thread_nr].extra_args=extra_args;
+    // init ret
+    args[thread_nr].ret=0;
+    // run threads
+    pthread_create(thread_ids->values+thread_nr, NULL, INTEGRAL_FUNC(integrate_gauss_smarttmp_thread), (void*)(args+thread_nr));
+  }
+
+  // wait for completion and join
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    pthread_join(thread_ids->values[thread_nr], NULL);
+
+    if(args[thread_nr].ret<0){
+      ret=args[thread_nr].ret;
+    }
+    else{
+      INTEGRAL_FLOAT_ADD(*out, *out, *(args[thread_nr].out));
+    }
+  }
+
+  // multiply by size of interval
+  INTEGRAL_FLOAT_MUL(*out, *out, tmps->values[0]);
+
+  // free x and weights
+  for(thread_nr=0;thread_nr<threads;thread_nr++){
+    INTEGRAL_FLOATARRAY_FUNC(free_vects)(args[thread_nr].x);
+    INTEGRAL_FLOATARRAY_FUNC(free_vects)(args[thread_nr].weights);
+  }
+  return(ret);
+}
+// per-thread function
+void* INTEGRAL_FUNC(integrate_gauss_smarttmp_thread) (void* args){
+  unsigned int i;
+  int ret;
+  struct INTEGRAL_FUNC(pthread_integrate_gauss_smarttmp_args)* argument=((struct INTEGRAL_FUNC(pthread_integrate_gauss_smarttmp_args)*)args);
+
+  INTEGRAL_FLOAT_SET_UI(*(argument->out), 0);
+
+  for(i=0;i<argument->x.length;i++){
+    // evaluate
+    ret=(*(argument->func))(argument->tmp, argument->x.values[i], argument->extra_args);
+
+    if(ret<0){
+      argument->ret=ret;
+      return(NULL);
+    }
+    // check whether argument->tmp is a number
+    if(! INTEGRAL_FLOAT_ISNUMBER(*(argument->tmp))){
+      argument->ret=LIBINUM_ERROR_NAN;
+      return(NULL);
+    }
+
+    INTEGRAL_FLOAT_MUL(*(argument->tmp), *(argument->tmp), argument->weights.values[i]);
+    INTEGRAL_FLOAT_ADD(*(argument->out), *(argument->out), *(argument->tmp));
+  }
+
+  return(NULL);
+}
+#endif
+
+
+// compute the abcissa and weights for the Gauss-Legendre numerical integration scheme
+int INTEGRAL_FUNC(gauss_legendre_weights) (unsigned int order, INTEGRAL_FLOAT_TYPE tolerance, unsigned int maxiter, INTEGRAL_FLOATARRAY_TYPE* abcissa, INTEGRAL_FLOATARRAY_TYPE* weights){
+  unsigned int i;
+  INTEGRAL_FLOAT_TYPE x, tmp;
+  INTEGRAL_POLYNOMIALARRAY_TYPE legendre;
+  int ret;
+
+  INTEGRAL_FLOATARRAY_FUNC(init) (abcissa, order);
+  INTEGRAL_FLOATARRAY_FUNC(init) (weights, order);
+
+  #ifdef INTEGRAL_FLOAT_INIT
+    INTEGRAL_FLOAT_INIT(x);
+    INTEGRAL_FLOAT_INIT(tmp);
+  #endif
+
+  INTEGRAL_POLYNOMIALARRAY_FUNC(init) (&legendre, 2);
+  // compute the roots of the 'order'-th Legendre polynomial
+  INTEGRAL_POLYNOMIAL_FUNC(legendre) (order, legendre.values);
+  INTEGRAL_POLYNOMIAL_FUNC(derive) (legendre.values[0], legendre.values+1);
+  legendre.length=2;
+
+  for(i=0;i<order;i++){
+    // initial guess
+    INTEGRAL_FLOAT_SET_D(x, cos(3.1415926*(i+0.75)/(order+0.5)));
+    ret=INTEGRAL_FUNC(root_newton_inplace) (&x, &INTEGRAL_FUNC(legendre_wrapper), &INTEGRAL_FUNC(deriv_legendre_wrapper), tolerance, maxiter, &legendre);
+
+    if(ret<0){
+      #ifdef INTEGRAL_FLOAT_FREE
+	INTEGRAL_FLOAT_FREE(x);
+	INTEGRAL_FLOAT_FREE(tmp);
+      #endif
+      INTEGRAL_POLYNOMIALARRAY_FUNC(free) (legendre);
+      INTEGRAL_FLOATARRAY_FUNC(free) (*abcissa);
+      INTEGRAL_FLOATARRAY_FUNC(free) (*weights);
+      return(ret);
+    }
+
+    INTEGRAL_FLOATARRAY_FUNC(append) (x, abcissa);
+
+    // weight: 2/((1-x^2)*(deriv_Legendre(x)))^2
+    INTEGRAL_POLYNOMIAL_FUNC(evaluate) (&tmp, x, legendre.values[1]);
+    INTEGRAL_FLOAT_POW_UI(tmp, tmp, 2);
+    INTEGRAL_FLOAT_POW_UI(x, x, 2);
+    INTEGRAL_FLOAT_UI_SUB(x, 1, x);
+    INTEGRAL_FLOAT_MUL(x, x, tmp);
+    INTEGRAL_FLOAT_SET_UI(tmp, 2);
+    INTEGRAL_FLOAT_DIV(x, tmp, x);
+
+    INTEGRAL_FLOATARRAY_FUNC(append) (x, weights);
+  }
+
+  #ifdef INTEGRAL_FLOAT_FREE
+    INTEGRAL_FLOAT_FREE(x);
+    INTEGRAL_FLOAT_FREE(tmp);
+  #endif
+  INTEGRAL_POLYNOMIALARRAY_FUNC(free) (legendre);
+
+  return(0);
+}
+
+// wrapper functions to evaluate the legendre polynomial and its derivative
+int INTEGRAL_FUNC(legendre_wrapper) (INTEGRAL_FLOAT_TYPE* out, INTEGRAL_FLOAT_TYPE in, void* legendre){
+  INTEGRAL_POLYNOMIAL_FUNC(evaluate) (out, in, ((INTEGRAL_POLYNOMIALARRAY_TYPE*)legendre)->values[0]);
+  return(0);
+}
+int INTEGRAL_FUNC(deriv_legendre_wrapper) (INTEGRAL_FLOAT_TYPE* out, INTEGRAL_FLOAT_TYPE in, void* legendre){
+  INTEGRAL_POLYNOMIAL_FUNC(evaluate) (out, in, ((INTEGRAL_POLYNOMIALARRAY_TYPE*)legendre)->values[1]);
+  return(0);
+}
diff --git a/src/integral_base.h b/src/integral_base.h
new file mode 100644
index 0000000..4440a52
--- /dev/null
+++ b/src/integral_base.h
@@ -0,0 +1,48 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for integration
+
+  see integral_*.h for the values taken by INTEGRAL_FUNC, etc...
+*/
+
+#include <pthread.h>
+
+// compute the integral of a real function of 1 variable using the Gauss-Legendre scheme
+int INTEGRAL_FUNC(integrate_gauss) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, void* extra_args);
+
+// multithreaded version
+int INTEGRAL_FUNC(integrate_gauss_multithread) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, void* extra_args, unsigned int threads, array_pthread_t* thread_ids);
+// per-thread function
+void* INTEGRAL_FUNC(integrate_gauss_thread) (void* args);
+
+// smart management of temporary variables: initialize as many as needed but allow them to be re-used instead of freeing them
+#ifdef INTEGRAL_FLOAT_INIT
+int INTEGRAL_FUNC(integrate_gauss_smarttmp) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, INTEGRAL_FLOATARRAY_TYPE* tmps, void* extra_args);
+// multithreaded version with smart management of temporary variables (only for datatypes that need to be initialized)
+int INTEGRAL_FUNC(integrate_gauss_smarttmp_multithread) (INTEGRAL_FLOAT_TYPE* out, int (*func)(INTEGRAL_FLOAT_TYPE*, INTEGRAL_FLOAT_TYPE, void*), INTEGRAL_FLOAT_TYPE lower, INTEGRAL_FLOAT_TYPE upper, INTEGRAL_FLOATARRAY_TYPE abcissa, INTEGRAL_FLOATARRAY_TYPE weights, INTEGRAL_FLOATARRAY_TYPE* tmps, void* extra_args, unsigned int threads, array_pthread_t* thread_ids);
+// per-thread function
+void* INTEGRAL_FUNC(integrate_gauss_smarttmp_thread) (void* args);
+#endif
+
+// compute the abcissa and weights for the Gauss-Legendre numerical integration scheme
+int INTEGRAL_FUNC(gauss_legendre_weights) (unsigned int order, INTEGRAL_FLOAT_TYPE tolerance, unsigned int maxiter, INTEGRAL_FLOATARRAY_TYPE* abcissa, INTEGRAL_FLOATARRAY_TYPE* weights);
+
+// wrapper functions to evaluate the legendre polynomial and its derivative
+int INTEGRAL_FUNC(legendre_wrapper) (INTEGRAL_FLOAT_TYPE* out, INTEGRAL_FLOAT_TYPE in, void* legendre);
+int INTEGRAL_FUNC(deriv_legendre_wrapper) (INTEGRAL_FLOAT_TYPE* out, INTEGRAL_FLOAT_TYPE in, void* legendre);
+
diff --git a/src/integral_double.h b/src/integral_double.h
new file mode 100644
index 0000000..de179fb
--- /dev/null
+++ b/src/integral_double.h
@@ -0,0 +1,81 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for integration using with doubles
+*/
+
+
+// reset CPP macros
+#undef INTEGRAL_FUNC
+#undef INTEGRAL_FLOAT_TYPE
+#undef INTEGRAL_FLOAT_INIT
+#undef INTEGRAL_FLOAT_FREE
+#undef INTEGRAL_FLOAT_SET_D
+#undef INTEGRAL_FLOAT_SET_UI
+#undef INTEGRAL_FLOAT_ADD
+#undef INTEGRAL_FLOAT_SUB
+#undef INTEGRAL_FLOAT_UI_SUB
+#undef INTEGRAL_FLOAT_MUL
+#undef INTEGRAL_FLOAT_POW_UI
+#undef INTEGRAL_FLOAT_DIV
+#undef INTEGRAL_FLOAT_DIV_UI
+#undef INTEGRAL_FLOAT_ISNUMBER
+#undef INTEGRAL_FLOATARRAY_TYPE
+#undef INTEGRAL_FLOATARRAY_FUNC
+#undef INTEGRAL_POLYNOMIAL_FUNC
+#undef INTEGRAL_POLYNOMIALARRAY_TYPE
+#undef INTEGRAL_POLYNOMIALARRAY_FUNC
+
+// suffix of function names
+#define INTEGRAL_FUNC(NAME) NAME ## _double
+
+// type of floats
+#define INTEGRAL_FLOAT_TYPE double
+// set float from double
+#define INTEGRAL_FLOAT_SET_D(FLOAT, VAL) FLOAT=VAL
+// set float from unsigned int
+#define INTEGRAL_FLOAT_SET_UI(FLOAT, VAL) FLOAT=(double)VAL
+// add floats
+#define INTEGRAL_FLOAT_ADD(FLOAT, VAL1, VAL2) FLOAT=VAL1+VAL2
+// subtract floats
+#define INTEGRAL_FLOAT_SUB(FLOAT, VAL1, VAL2) FLOAT=VAL1-VAL2
+// subtract floats in which the first is specified as an unsigned int
+#define INTEGRAL_FLOAT_UI_SUB(FLOAT, VAL1, VAL2) FLOAT=(double)VAL1-VAL2
+// multiply floats
+#define INTEGRAL_FLOAT_MUL(FLOAT, VAL1, VAL2) FLOAT=VAL1*VAL2
+// power of a float, specified as an unsigned int
+#define INTEGRAL_FLOAT_POW_UI(FLOAT, VAL1, VAL2) FLOAT=pow(VAL1, (double)VAL2)
+// divide floats
+#define INTEGRAL_FLOAT_DIV(FLOAT, VAL1, VAL2) FLOAT=VAL1/VAL2
+// divide floats, one of which is specified as an unsigned int
+#define INTEGRAL_FLOAT_DIV_UI(FLOAT, VAL1, VAL2) FLOAT=VAL1/(double)VAL2
+// check whether a float is a regular number
+#define INTEGRAL_FLOAT_ISNUMBER(FLOAT) (fpclassify(FLOAT)>=FP_ZERO)
+
+// type of float arrays
+#define INTEGRAL_FLOATARRAY_TYPE array_double
+// prefix of float array function names
+#define INTEGRAL_FLOATARRAY_FUNC(NAME) array_double_ ## NAME
+
+// prefix of polynomial function names
+#define INTEGRAL_POLYNOMIAL_FUNC(NAME) polynomial_double_ ## NAME
+
+// type of polynomial arrays
+#define INTEGRAL_POLYNOMIALARRAY_TYPE array_polynomial_double
+// prefix of polynomial array function names
+#define INTEGRAL_POLYNOMIALARRAY_FUNC(NAME) array_polynomial_double_ ## NAME
+
diff --git a/src/integral_ldouble.h b/src/integral_ldouble.h
new file mode 100644
index 0000000..5b1d26f
--- /dev/null
+++ b/src/integral_ldouble.h
@@ -0,0 +1,81 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for integration using with long doubles
+*/
+
+
+// reset CPP macros
+#undef INTEGRAL_FUNC
+#undef INTEGRAL_FLOAT_TYPE
+#undef INTEGRAL_FLOAT_INIT
+#undef INTEGRAL_FLOAT_FREE
+#undef INTEGRAL_FLOAT_SET_D
+#undef INTEGRAL_FLOAT_SET_UI
+#undef INTEGRAL_FLOAT_ADD
+#undef INTEGRAL_FLOAT_SUB
+#undef INTEGRAL_FLOAT_UI_SUB
+#undef INTEGRAL_FLOAT_MUL
+#undef INTEGRAL_FLOAT_POW_UI
+#undef INTEGRAL_FLOAT_DIV
+#undef INTEGRAL_FLOAT_DIV_UI
+#undef INTEGRAL_FLOAT_ISNUMBER
+#undef INTEGRAL_FLOATARRAY_TYPE
+#undef INTEGRAL_FLOATARRAY_FUNC
+#undef INTEGRAL_POLYNOMIAL_FUNC
+#undef INTEGRAL_POLYNOMIALARRAY_TYPE
+#undef INTEGRAL_POLYNOMIALARRAY_FUNC
+
+// suffix of function names
+#define INTEGRAL_FUNC(NAME) NAME ## _ldouble
+
+// type of floats
+#define INTEGRAL_FLOAT_TYPE long double
+// set float from double
+#define INTEGRAL_FLOAT_SET_D(FLOAT, VAL) FLOAT=(long double)VAL
+// set float from unsigned int
+#define INTEGRAL_FLOAT_SET_UI(FLOAT, VAL) FLOAT=(long double)VAL
+// add floats
+#define INTEGRAL_FLOAT_ADD(FLOAT, VAL1, VAL2) FLOAT=VAL1+VAL2
+// subtract floats
+#define INTEGRAL_FLOAT_SUB(FLOAT, VAL1, VAL2) FLOAT=VAL1-VAL2
+// subtract floats in which the first is specified as an unsigned int
+#define INTEGRAL_FLOAT_UI_SUB(FLOAT, VAL1, VAL2) FLOAT=(long double)VAL1-VAL2
+// multiply floats
+#define INTEGRAL_FLOAT_MUL(FLOAT, VAL1, VAL2) FLOAT=VAL1*VAL2
+// power of a float, specified as an unsigned int
+#define INTEGRAL_FLOAT_POW_UI(FLOAT, VAL1, VAL2) FLOAT=powl(VAL1, (long double)VAL2)
+// divide floats
+#define INTEGRAL_FLOAT_DIV(FLOAT, VAL1, VAL2) FLOAT=VAL1/VAL2
+// divide floats, one of which is specified as an unsigned int
+#define INTEGRAL_FLOAT_DIV_UI(FLOAT, VAL1, VAL2) FLOAT=VAL1/(long double)VAL2
+// check whether a float is a regular number
+#define INTEGRAL_FLOAT_ISNUMBER(FLOAT) (fpclassify(FLOAT)>=FP_ZERO)
+
+// type of float arrays
+#define INTEGRAL_FLOATARRAY_TYPE array_ldouble
+// prefix of float array function names
+#define INTEGRAL_FLOATARRAY_FUNC(NAME) array_ldouble_ ## NAME
+
+// prefix of polynomial function names
+#define INTEGRAL_POLYNOMIAL_FUNC(NAME) polynomial_ldouble_ ## NAME
+
+// type of polynomial arrays
+#define INTEGRAL_POLYNOMIALARRAY_TYPE array_polynomial_ldouble
+// prefix of polynomial array function names
+#define INTEGRAL_POLYNOMIALARRAY_FUNC(NAME) array_polynomial_ldouble_ ## NAME
+
diff --git a/src/integral_mpfr.h b/src/integral_mpfr.h
new file mode 100644
index 0000000..30c9b1a
--- /dev/null
+++ b/src/integral_mpfr.h
@@ -0,0 +1,84 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for integration using with multi-precision floats (mpfr)
+*/
+
+
+// reset CPP macros
+#undef INTEGRAL_FUNC
+#undef INTEGRAL_FLOAT_TYPE
+#undef INTEGRAL_FLOAT_INIT
+#undef INTEGRAL_FLOAT_FREE
+#undef INTEGRAL_FLOAT_SET_D
+#undef INTEGRAL_FLOAT_SET_UI
+#undef INTEGRAL_FLOAT_ADD
+#undef INTEGRAL_FLOAT_SUB
+#undef INTEGRAL_FLOAT_UI_SUB
+#undef INTEGRAL_FLOAT_MUL
+#undef INTEGRAL_FLOAT_POW_UI
+#undef INTEGRAL_FLOAT_DIV
+#undef INTEGRAL_FLOAT_DIV_UI
+#undef INTEGRAL_FLOAT_ISNUMBER
+#undef INTEGRAL_FLOATARRAY_TYPE
+#undef INTEGRAL_FLOATARRAY_FUNC
+#undef INTEGRAL_POLYNOMIAL_FUNC
+#undef INTEGRAL_POLYNOMIALARRAY_TYPE
+#undef INTEGRAL_POLYNOMIALARRAY_FUNC
+
+// suffix of function names
+#define INTEGRAL_FUNC(NAME) NAME ## _mpfr
+
+// type of floats
+#define INTEGRAL_FLOAT_TYPE mpfr_t
+// init float
+#define INTEGRAL_FLOAT_INIT(VAR) mpfr_init(VAR)
+// free float
+#define INTEGRAL_FLOAT_FREE(VAR) mpfr_clear(VAR)
+// set float from double
+#define INTEGRAL_FLOAT_SET_D(FLOAT, VAL) mpfr_set_d(FLOAT, VAL, MPFR_RNDN)
+// set float from unsigned int
+#define INTEGRAL_FLOAT_SET_UI(FLOAT, VAL) mpfr_set_ui(FLOAT, VAL, MPFR_RNDN)
+// add floats
+#define INTEGRAL_FLOAT_ADD(FLOAT, VAL1, VAL2) mpfr_add(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// subtract floats
+#define INTEGRAL_FLOAT_SUB(FLOAT, VAL1, VAL2) mpfr_sub(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// subtract floats in which the first is specified as an unsigned int
+#define INTEGRAL_FLOAT_UI_SUB(FLOAT, VAL1, VAL2) mpfr_ui_sub(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// multiply floats
+#define INTEGRAL_FLOAT_MUL(FLOAT, VAL1, VAL2) mpfr_mul(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// power of a float, specified as an unsigned int
+#define INTEGRAL_FLOAT_POW_UI(FLOAT, VAL1, VAL2) mpfr_pow_ui(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// divide floats
+#define INTEGRAL_FLOAT_DIV(FLOAT, VAL1, VAL2) mpfr_div(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// divide floats, one of which is specified as an unsigned int
+#define INTEGRAL_FLOAT_DIV_UI(FLOAT, VAL1, VAL2) mpfr_div_ui(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// check whether a float is a regular number
+#define INTEGRAL_FLOAT_ISNUMBER(FLOAT) mpfr_number_p(FLOAT)!=0
+
+// type of float arrays
+#define INTEGRAL_FLOATARRAY_TYPE array_mpfr
+// prefix of float array function names
+#define INTEGRAL_FLOATARRAY_FUNC(NAME) array_mpfr_ ## NAME
+
+// prefix of polynomial function names
+#define INTEGRAL_POLYNOMIAL_FUNC(NAME) polynomial_mpfr_ ## NAME
+
+// type of polynomial arrays
+#define INTEGRAL_POLYNOMIALARRAY_TYPE array_polynomial_mpfr
+// prefix of polynomial array function names
+#define INTEGRAL_POLYNOMIALARRAY_FUNC(NAME) array_polynomial_mpfr_ ## NAME
diff --git a/src/polynomial.c b/src/polynomial.c
new file mode 100644
index 0000000..9221241
--- /dev/null
+++ b/src/polynomial.c
@@ -0,0 +1,55 @@
+/*
+Copyright 2016 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 "polynomial.h"
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <mpfr.h>
+#include <math.h>
+#include "array.h"
+#include "utils.h"
+#include "errors.h"
+
+
+//--------------------------------------------------
+//
+//  mpfr coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+
+#include "polynomial_mpfr.h"
+#include "polynomial_base.c"
+
+
+//--------------------------------------------------
+//
+//  double coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+
+#include "polynomial_double.h"
+#include "polynomial_base.c"
+
+
+//--------------------------------------------------
+//
+//  long double coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+
+#include "polynomial_ldouble.h"
+#include "polynomial_base.c"
diff --git a/src/polynomial.h b/src/polynomial.h
new file mode 100644
index 0000000..376149b
--- /dev/null
+++ b/src/polynomial.h
@@ -0,0 +1,55 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Real polynomial
+*/
+
+
+#ifndef LIBINUM_POLYNOMIAL_H
+#define LIBINUM_POLYNOMIAL_H
+
+#include "types.h"
+
+//--------------------------------------------------
+//
+//  mpfr coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+#include "polynomial_mpfr.h"
+#include "polynomial_base.h"
+
+
+//--------------------------------------------------
+//
+//  double coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+
+#include "polynomial_double.h"
+#include "polynomial_base.h"
+
+
+//--------------------------------------------------
+//
+//  long double coefficients, unsigned int order (polynomial_mpfr)
+//
+//--------------------------------------------------
+
+#include "polynomial_ldouble.h"
+#include "polynomial_base.h"
+
+#endif
diff --git a/src/polynomialMV.c b/src/polynomialMV.c
new file mode 100644
index 0000000..2cafcf6
--- /dev/null
+++ b/src/polynomialMV.c
@@ -0,0 +1,30 @@
+/*
+Copyright 2016 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 "polynomialMV.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include "array.h"
+#include "errors.h"
+
+// multi-precision integer coefficients
+#include "polynomialMV_mpz.h"
+#include "polynomialMV_base.c"
+
+// integer coefficients
+#include "polynomialMV_int.h"
+#include "polynomialMV_base.c"
diff --git a/src/polynomialMV.h b/src/polynomialMV.h
new file mode 100644
index 0000000..e368d5b
--- /dev/null
+++ b/src/polynomialMV.h
@@ -0,0 +1,34 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Multi-variable polynomials
+*/
+
+#ifndef LIBINUM_POLYNOMIAL_MV_H
+#define LIBINUM_POLYNOMIAL_MV_H
+
+#include "types.h"
+
+// multi-precision integer coefficients
+#include "polynomialMV_mpz.h"
+#include "polynomialMV_base.h"
+
+// integer coefficients
+#include "polynomialMV_int.h"
+#include "polynomialMV_base.h"
+
+#endif
diff --git a/src/polynomialMV_base.c b/src/polynomialMV_base.c
new file mode 100644
index 0000000..814c499
--- /dev/null
+++ b/src/polynomialMV_base.c
@@ -0,0 +1,293 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for multivariable polynomials
+
+  see polynomialMV_*.h for the values taken by POLYNOMIALMV_FUNC, etc...
+*/
+
+// init
+int POLYNOMIALMV_FUNC(init) (POLYNOMIALMV_TYPENAME* polynomial, int size){
+  polynomial->factors=calloc(size, sizeof(POLYNOMIALMV_FACTOR_TYPE));
+  polynomial->coefficients=calloc(size, sizeof(POLYNOMIALMV_COEF_TYPE));
+  polynomial->length=0;
+  polynomial->memory=size;
+  return(0);
+}
+
+// free memory
+int POLYNOMIALMV_FUNC(free) (POLYNOMIALMV_TYPENAME polynomial){
+  unsigned int i;
+  for(i=0;i<polynomial.length;i++){
+    POLYNOMIALMV_FACTOR_FUNC(free) (polynomial.factors[i]);
+    #ifdef POLYNOMIALMV_COEF_FREE
+      POLYNOMIALMV_COEF_FREE(polynomial.coefficients[i]);
+    #endif
+  }
+  free(polynomial.factors);
+  free(polynomial.coefficients);
+  return(0);
+}
+
+// resize the memory allocated to a polynomial
+int POLYNOMIALMV_FUNC(resize) (POLYNOMIALMV_TYPENAME* polynomial, int new_size){
+  POLYNOMIALMV_TYPENAME new_poly;
+  unsigned int i;
+
+  POLYNOMIALMV_FUNC(init) (&new_poly, new_size);
+  for(i=0;i<polynomial->length;i++){
+    new_poly.factors[i]=polynomial->factors[i];
+    POLYNOMIALMV_COEF_SET(new_poly.coefficients[i], polynomial->coefficients[i]);
+  }
+  new_poly.length=polynomial->length;
+
+  free(polynomial->factors);
+  free(polynomial->coefficients);
+
+  *polynomial=new_poly;
+  return(0);
+}
+
+// copy a polynomial
+int POLYNOMIALMV_FUNC(cpy) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* output){
+  POLYNOMIALMV_FUNC(init) (output, input.length);
+  POLYNOMIALMV_FUNC(cpy_noinit) (input, output);
+  return(0);
+}
+int POLYNOMIALMV_FUNC(cpy_noinit) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* output){
+  unsigned int i;
+  if(output->memory<input.length){
+    return(LIBINUM_ERROR_ILLEGAL_MEMORY_ACCESS);
+  }
+  for(i=0;i<input.length;i++){
+    POLYNOMIALMV_FACTOR_FUNC(cpy) (input.factors[i],output->factors+i);
+    POLYNOMIALMV_COEF_CPY(output->coefficients[i], input.coefficients[i]);
+  }
+  output->length=input.length;
+  
+  return(0);
+}
+
+// append an element to a polynomial
+int POLYNOMIALMV_FUNC(append) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  unsigned int offset=output->length;
+
+  if(output->length>=output->memory){
+    POLYNOMIALMV_FUNC(resize) (output,2*output->memory+1);
+  }
+
+  // copy and allocate
+  POLYNOMIALMV_FACTOR_FUNC(cpy) (factor, output->factors+offset);
+  POLYNOMIALMV_COEF_CPY(output->coefficients[offset], coef);
+  //increment length
+  output->length++;
+
+  return(0);
+}
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  unsigned int i;
+  unsigned short int foundit=0;
+  for(i=0;i<output->length;i++){
+    if(POLYNOMIALMV_FACTOR_FUNC(cmp) (factor, output->factors[i])==0){
+      POLYNOMIALMV_COEF_ADD(output->coefficients[i], output->coefficients[i], coef);
+      foundit=1;
+      break;
+    }
+  }
+  if(foundit==0){
+    POLYNOMIALMV_FUNC(append) (factor, coef, output);
+  }
+  return(0);
+}
+// do not allocate memory for the factor or the coefficient
+int POLYNOMIALMV_FUNC(append_noinit) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  int offset=output->length;
+
+  if(output->length>=output->memory){
+    POLYNOMIALMV_FUNC(resize) (output,2*output->memory+1);
+  }
+
+  // copy without allocating
+  output->factors[offset]=factor;
+  POLYNOMIALMV_COEF_SET(output->coefficients[offset], coef);
+  // increment length
+  output->length++;
+  return(0);
+}
+// noinit factor but init coefficient
+int POLYNOMIALMV_FUNC(append_noinitfactor) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  int offset=output->length;
+
+  if(output->length>=output->memory){
+    POLYNOMIALMV_FUNC(resize) (output,2*output->memory+1);
+  }
+
+  // copy without allocating
+  output->factors[offset]=factor;
+  POLYNOMIALMV_COEF_CPY(output->coefficients[offset], coef);
+  // increment length
+  output->length++;
+  return(0);
+}
+// noinit
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_noinit_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  unsigned int i;
+  unsigned short int foundit=0;
+  for(i=0;i<output->length;i++){
+    if(POLYNOMIALMV_FACTOR_FUNC(cmp) (factor, output->factors[i])==0){
+      POLYNOMIALMV_COEF_ADD(output->coefficients[i], output->coefficients[i], coef);
+      foundit=1;
+      // free
+      POLYNOMIALMV_FACTOR_FUNC(free) (factor);
+      #ifdef POLYNOMIALMV_COEF_FREE
+	POLYNOMIALMV_COEF_FREE(coef);
+      #endif
+      break;
+    }
+  }
+  if(foundit==0){
+    POLYNOMIALMV_FUNC(append_noinit) (factor, coef, output);
+  }
+  return(0);
+}
+// noinit factor but init coefficient
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_noinitfactor_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output){
+  unsigned int i;
+  unsigned short int foundit=0;
+  for(i=0;i<output->length;i++){
+    if(POLYNOMIALMV_FACTOR_FUNC(cmp) (factor, output->factors[i])==0){
+      POLYNOMIALMV_COEF_ADD(output->coefficients[i], output->coefficients[i], coef);
+      foundit=1;
+      // free factor
+      POLYNOMIALMV_FACTOR_FUNC(free) (factor);
+      break;
+    }
+  }
+  if(foundit==0){
+    POLYNOMIALMV_FUNC(append_noinitfactor) (factor, coef, output);
+  }
+  return(0);
+}
+
+// add polynomials (inplace)
+int POLYNOMIALMV_FUNC(add_inplace) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* inout){
+  unsigned int i;
+  for(i=0;i<input.length;i++){
+    POLYNOMIALMV_FUNC(append_inplace) (input.factors[i], input.coefficients[i], inout);
+  }
+  return(0);
+}
+// not inplace
+int POLYNOMIALMV_FUNC(add) (POLYNOMIALMV_TYPENAME input1, POLYNOMIALMV_TYPENAME input2, POLYNOMIALMV_TYPENAME* output){
+  POLYNOMIALMV_FUNC(cpy) (input1, output);
+  POLYNOMIALMV_FUNC(add_inplace) (input2, output);
+  return(0);
+}
+
+// multiply a polynomial by a scalar
+int POLYNOMIALMV_FUNC(multiply_scalar) (POLYNOMIALMV_TYPENAME polynomial, POLYNOMIALMV_COEF_TYPE num){
+  unsigned int i;
+  for(i=0;i<polynomial.length;i++){
+    POLYNOMIALMV_COEF_MUL(polynomial.coefficients[i], polynomial.coefficients[i], num);
+  }
+  return(0);
+}
+
+// multiply polynomials
+int POLYNOMIALMV_FUNC(prod) (POLYNOMIALMV_TYPENAME input1, POLYNOMIALMV_TYPENAME input2, POLYNOMIALMV_TYPENAME* output){
+  // position in polys
+  unsigned int pos1, pos2;
+  POLYNOMIALMV_FACTOR_TYPE out_factor;
+  POLYNOMIALMV_COEF_TYPE out_num;
+
+  POLYNOMIALMV_FUNC(init) (output, input1.length);
+  #ifdef POLYNOMIALMV_COEF_INIT
+    POLYNOMIALMV_COEF_INIT(out_num);
+  #endif
+
+  // loop over terms
+  for(pos1=0;pos1<input1.length;pos1++){
+    for(pos2=0;pos2<input2.length;pos2++){
+      // allocate
+      POLYNOMIALMV_FACTOR_FUNC(init) (&out_factor, input1.factors[pos1].length + input2.factors[pos2].length);
+
+      // concatenate factor
+      POLYNOMIALMV_FACTOR_FUNC(concat) (input1.factors[pos1],&out_factor);
+      POLYNOMIALMV_FACTOR_FUNC(concat) (input2.factors[pos2],&out_factor);
+      // sort factor
+      POLYNOMIALMV_FACTOR_FUNC(sort) (out_factor);
+      // multiply coefficient
+      POLYNOMIALMV_COEF_MUL(out_num, input1.coefficients[pos1], input2.coefficients[pos2]);
+
+      // write factor and coef
+      POLYNOMIALMV_FUNC(append_noinit_inplace) (out_factor, out_num, output);
+    }
+  }
+
+  #ifdef POLYNOMIALMV_COEF_FREE
+    POLYNOMIALMV_COEF_FREE(out_num);
+  #endif
+  return(0);
+}
+
+// order factors
+int POLYNOMIALMV_FUNC(order) (POLYNOMIALMV_TYPENAME polynomial){
+  unsigned int i,j;
+  for(i=0;i<polynomial.length;i++){
+    for(j=0;j<polynomial.factors[i].length-1;j++){
+      if(polynomial.factors[i].values[j] > polynomial.factors[i].values[j+1]){
+	POLYNOMIALMV_FACTOR_FUNC(sort) (polynomial.factors[i]);
+	break;
+      }
+    }
+  }
+  return(0);
+}
+
+// print
+int POLYNOMIALMV_FUNC(print) (POLYNOMIALMV_TYPENAME polynomial){
+  unsigned int i,j;
+
+  // for each monomial
+  for(i=0;i<polynomial.length;i++){
+    if(i==0){
+      printf("  ");
+    }
+    else{
+      printf("+ ");
+    }
+
+    // print num
+    POLYNOMIALMV_COEF_PRINT(polynomial.coefficients[i]);
+
+    // print factors
+    for(j=0;j<polynomial.factors[i].length;j++){
+      POLYNOMIALMV_FACTOR_ELT_PRINT(polynomial.factors[i].values[j]);
+    }
+
+    printf("\n");
+  }
+  return(0);
+}
+
+
diff --git a/src/polynomialMV_base.h b/src/polynomialMV_base.h
new file mode 100644
index 0000000..dde20ac
--- /dev/null
+++ b/src/polynomialMV_base.h
@@ -0,0 +1,67 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for multivariable polynomials
+
+  polynomialMV_*.h for the values taken by POLYNOMIALMV_FUNC, etc...
+*/
+
+// init
+int POLYNOMIALMV_FUNC(init) (POLYNOMIALMV_TYPENAME* polynomial, int size);
+int POLYNOMIALMV_FUNC(free) (POLYNOMIALMV_TYPENAME polynomial);
+
+// resize the memory allocated to a polynomial
+int POLYNOMIALMV_FUNC(resize) (POLYNOMIALMV_TYPENAME* polynomial, int new_size);
+
+// copy a polynomial
+int POLYNOMIALMV_FUNC(cpy) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* output);
+int POLYNOMIALMV_FUNC(cpy_noinit) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* output);
+
+// append an element to a polynomial
+int POLYNOMIALMV_FUNC(append) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+// do not allocate memory for factor or coefficient
+int POLYNOMIALMV_FUNC(append_noinit) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+// noinit factor but init coefficient
+int POLYNOMIALMV_FUNC(append_noinitfactor) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+// noinit
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_noinit_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+// noinit factors but init coefficient
+// if there already exists an element with the same factor, then just add coefficients
+// requires the factors to be ordered
+int POLYNOMIALMV_FUNC(append_noinitfactor_inplace) (POLYNOMIALMV_FACTOR_TYPE factor, POLYNOMIALMV_COEF_TYPE coef, POLYNOMIALMV_TYPENAME* output);
+
+// add polynomials (inplace)
+int POLYNOMIALMV_FUNC(add_inplace) (POLYNOMIALMV_TYPENAME input, POLYNOMIALMV_TYPENAME* inout);
+// not inplace
+int POLYNOMIALMV_FUNC(add) (POLYNOMIALMV_TYPENAME input1, POLYNOMIALMV_TYPENAME input2, POLYNOMIALMV_TYPENAME* output);
+
+// multiply a polynomial by a scalar
+int POLYNOMIALMV_FUNC(multiply_scalar) (POLYNOMIALMV_TYPENAME polynomial, POLYNOMIALMV_COEF_TYPE num);
+
+// multiply polynomials
+int POLYNOMIALMV_FUNC(prod) (POLYNOMIALMV_TYPENAME input1, POLYNOMIALMV_TYPENAME input2, POLYNOMIALMV_TYPENAME* output);
+
+// order factors
+int POLYNOMIALMV_FUNC(order) (POLYNOMIALMV_TYPENAME polynomial);
+
+// print
+int POLYNOMIALMV_FUNC(print) (POLYNOMIALMV_TYPENAME polynomial);
diff --git a/src/polynomialMV_int.h b/src/polynomialMV_int.h
new file mode 100644
index 0000000..cb25866
--- /dev/null
+++ b/src/polynomialMV_int.h
@@ -0,0 +1,62 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for multivariable polynomials with integer coefficients and integer-indexed variables
+*/
+
+
+// reset CPP macros
+#undef POLYNOMIALMV_TYPENAME
+#undef POLYNOMIALMV_FUNC
+#undef POLYNOMIALMV_COEF_TYPE
+#undef POLYNOMIALMV_COEF_INIT
+#undef POLYNOMIALMV_COEF_FREE
+#undef POLYNOMIALMV_COEF_SET
+#undef POLYNOMIALMV_COEF_CPY
+#undef POLYNOMIALMV_COEF_ADD
+#undef POLYNOMIALMV_COEF_MUL
+#undef POLYNOMIALMV_COEF_PRINT
+#undef POLYNOMIALMV_FACTOR_TYPE
+#undef POLYNOMIALMV_FACTOR_FUNC
+#undef POLYNOMIALMV_FACTOR_ELT_PRINT
+
+
+// name of the polynomial type
+#define POLYNOMIALMV_TYPENAME polynomialMV_int
+// prefix of function names
+#define POLYNOMIALMV_FUNC(NAME) polynomialMV_int_ ## NAME
+
+// type of the coefficient
+#define POLYNOMIALMV_COEF_TYPE int
+// set coefficient
+#define POLYNOMIALMV_COEF_SET(COEF, VAL) COEF = VAL
+// copy coefficient
+#define POLYNOMIALMV_COEF_CPY(COEF, VAL) COEF = VAL
+// add coefficients
+#define POLYNOMIALMV_COEF_ADD(COEF, VAL1, VAL2) COEF = VAL1 + VAL2
+// multiply coefficients
+#define POLYNOMIALMV_COEF_MUL(COEF, VAL1, VAL2) COEF = VAL1 * VAL2
+// print a coefficient
+#define POLYNOMIALMV_COEF_PRINT(COEF) printf("%d", COEF)
+
+// type of the factor (must be an array)
+#define POLYNOMIALMV_FACTOR_TYPE array_int
+// prefix of factor function names
+#define POLYNOMIALMV_FACTOR_FUNC(NAME) array_int_ ## NAME
+// print an element of a factor
+#define POLYNOMIALMV_FACTOR_ELT_PRINT(ELT) printf("[x%d]", ELT)
+
diff --git a/src/polynomialMV_mpz.h b/src/polynomialMV_mpz.h
new file mode 100644
index 0000000..4b99aa6
--- /dev/null
+++ b/src/polynomialMV_mpz.h
@@ -0,0 +1,66 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for multivariable polynomials with multi-precision integer (mpz) coefficients and integer-indexed variables
+*/
+
+
+// reset CPP macros
+#undef POLYNOMIALMV_TYPENAME
+#undef POLYNOMIALMV_FUNC
+#undef POLYNOMIALMV_COEF_TYPE
+#undef POLYNOMIALMV_COEF_INIT
+#undef POLYNOMIALMV_COEF_FREE
+#undef POLYNOMIALMV_COEF_SET
+#undef POLYNOMIALMV_COEF_CPY
+#undef POLYNOMIALMV_COEF_ADD
+#undef POLYNOMIALMV_COEF_MUL
+#undef POLYNOMIALMV_COEF_PRINT
+#undef POLYNOMIALMV_FACTOR_TYPE
+#undef POLYNOMIALMV_FACTOR_FUNC
+#undef POLYNOMIALMV_FACTOR_ELT_PRINT
+
+
+// name of the polynomial type
+#define POLYNOMIALMV_TYPENAME polynomialMV_mpz
+// prefix of function names
+#define POLYNOMIALMV_FUNC(NAME) polynomialMV_mpz_ ## NAME
+
+// type of the coefficient
+#define POLYNOMIALMV_COEF_TYPE mpz_t
+// init coefficient
+#define POLYNOMIALMV_COEF_INIT(VAR) mpz_init(VAR)
+// free coefficient
+#define POLYNOMIALMV_COEF_FREE(VAR) mpz_clear(VAR)
+// set coefficient
+#define POLYNOMIALMV_COEF_SET(COEF, VAL) COEF[0]=VAL[0]
+// copy coefficient
+#define POLYNOMIALMV_COEF_CPY(COEF, VAL) mpz_init(COEF); mpz_set(COEF, VAL)
+// add coefficients
+#define POLYNOMIALMV_COEF_ADD(COEF, VAL1, VAL2) mpz_add(COEF, VAL1, VAL2)
+// multiply coefficients
+#define POLYNOMIALMV_COEF_MUL(COEF, VAL1, VAL2) mpz_mul(COEF, VAL1, VAL2)
+// print a coefficient
+#define POLYNOMIALMV_COEF_PRINT(COEF) gmp_printf("%Zd", COEF)
+
+// type of the factor (must be an array)
+#define POLYNOMIALMV_FACTOR_TYPE array_int
+// prefix of factor function names
+#define POLYNOMIALMV_FACTOR_FUNC(NAME) array_int_ ## NAME
+// print an element of a factor
+#define POLYNOMIALMV_FACTOR_ELT_PRINT(ELT) printf("[x%d]", ELT)
+
diff --git a/src/polynomialMV_type.h b/src/polynomialMV_type.h
new file mode 100644
index 0000000..6fbf778
--- /dev/null
+++ b/src/polynomialMV_type.h
@@ -0,0 +1,29 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Structure definition for multivariable polynomials
+
+  see polynomialMV_*.h for the values taken by POLYNOMIALMV_TYPENAME, etc...
+*/
+
+typedef struct POLYNOMIALMV_TYPENAME {
+  POLYNOMIALMV_COEF_TYPE* coefficients;
+  POLYNOMIALMV_FACTOR_TYPE* factors;
+  unsigned int length;
+  unsigned int memory;
+} POLYNOMIALMV_TYPENAME;
+
diff --git a/src/polynomial_base.c b/src/polynomial_base.c
new file mode 100644
index 0000000..1b4137d
--- /dev/null
+++ b/src/polynomial_base.c
@@ -0,0 +1,337 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for real polynomials
+
+  see polynomial_*.h for the values taken by POLYNOMIAL_FUNC, etc...
+*/
+
+
+// init
+int POLYNOMIAL_FUNC(init) (POLYNOMIAL_TYPENAME* poly, unsigned int memory){
+  POLYNOMIAL_COEFSARRAY_FUNC(init) (&(poly->coefficients), memory);
+  POLYNOMIAL_ORDERSARRAY_FUNC(init) (&(poly->orders), memory);
+  return(0);
+}
+int POLYNOMIAL_FUNC(free) (POLYNOMIAL_TYPENAME poly){
+  POLYNOMIAL_COEFSARRAY_FUNC(free) (poly.coefficients);
+  POLYNOMIAL_ORDERSARRAY_FUNC(free) (poly.orders);
+  return(0);
+}
+
+// resize memory
+int POLYNOMIAL_FUNC(resize) (POLYNOMIAL_TYPENAME* poly, unsigned int newsize){
+  POLYNOMIAL_COEFSARRAY_FUNC(resize) (&(poly->coefficients), newsize);
+  POLYNOMIAL_ORDERSARRAY_FUNC(resize) (&(poly->orders), newsize);
+  return(0);
+}
+
+// add a monomial
+int POLYNOMIAL_FUNC(add_monomial) (POLYNOMIAL_COEF_TYPE val, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME* output){
+  unsigned int i;
+  if((output->coefficients.length != output->orders.length) || (output->coefficients.memory != output->orders.memory)){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  // check whether the order already exists in the polynomial
+  for(i=0;i< output->coefficients.length;i++){
+    if(POLYNOMIAL_ORDER_CMP(order, output->orders.values[i])==0){
+      POLYNOMIAL_COEF_ADD(output->coefficients.values[i], output->coefficients.values[i], val);
+      return(0);
+    }
+  }
+  // if the order does not already exist
+  if(output->coefficients.length >= output->coefficients.memory){
+    POLYNOMIAL_FUNC(resize) (output,2*output->coefficients.memory+1);
+  }
+  POLYNOMIAL_COEF_CPY(val, output->coefficients.values[output->coefficients.length]);
+  POLYNOMIAL_ORDER_CPY(order, output->orders.values[output->coefficients.length]);
+  output->coefficients.length++;
+  output->orders.length++;
+  return(0);
+}
+// from a double-valued coefficient and an unsigned int order
+int POLYNOMIAL_FUNC(add_monomial_dui) (double val, unsigned int order, POLYNOMIAL_TYPENAME* output){
+  unsigned int i;
+  if((output->coefficients.length != output->orders.length) || (output->coefficients.memory != output->orders.memory)){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  // check whether the order already exists in the polynomial
+  for(i=0;i< output->coefficients.length;i++){
+    if(POLYNOMIAL_ORDER_CMP_UI(output->orders.values[i], order)==0){
+      POLYNOMIAL_COEF_ADD_D(output->coefficients.values[i], output->coefficients.values[i], val);
+      return(0);
+    }
+  }
+  // if the order does not already exist
+  if(output->coefficients.length >= output->coefficients.memory){
+    POLYNOMIAL_FUNC(resize) (output,2*output->coefficients.memory+1);
+  }
+  POLYNOMIAL_COEF_CPY_D(val, output->coefficients.values[output->coefficients.length]);
+  POLYNOMIAL_ORDER_CPY_UI(order, output->orders.values[output->coefficients.length]);
+  output->coefficients.length++;
+  output->orders.length++;
+  return(0);
+}
+
+// copy
+int POLYNOMIAL_FUNC(cpy) (POLYNOMIAL_TYPENAME input, POLYNOMIAL_TYPENAME* output){
+  POLYNOMIAL_COEFSARRAY_FUNC(cpy) (input.coefficients, &(output->coefficients));
+  POLYNOMIAL_ORDERSARRAY_FUNC(cpy) (input.orders, &(output->orders));
+  return(0);
+}
+int POLYNOMIAL_FUNC(cpy_noinit) (POLYNOMIAL_TYPENAME input, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  ret=POLYNOMIAL_COEFSARRAY_FUNC(cpy_noinit) (input.coefficients, &(output->coefficients));
+  if(ret<0){
+    return(ret);
+  }
+  POLYNOMIAL_ORDERSARRAY_FUNC(cpy_noinit) (input.orders, &(output->orders));
+  return(ret);
+}
+
+// add
+int POLYNOMIAL_FUNC(add_inplace) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME* poly2){
+  unsigned int i;
+  if(poly1.coefficients.length != poly1.orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<poly1.coefficients.length;i++){
+    POLYNOMIAL_FUNC(add_monomial) (poly1.coefficients.values[i], poly1.orders.values[i], poly2);
+  }
+  return(0);
+}
+int POLYNOMIAL_FUNC(add) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME poly2, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  POLYNOMIAL_FUNC(cpy) (poly2, output);
+  ret=POLYNOMIAL_FUNC(add_inplace) (poly1, output);
+  return(ret);
+}
+
+// multiply by a scalar
+int POLYNOMIAL_FUNC(mul_scalar_inplace) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME* poly){
+  unsigned int i;
+  if(poly->coefficients.length != poly->orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<poly->coefficients.length;i++){
+    POLYNOMIAL_COEF_MUL(poly->coefficients.values[i], poly->coefficients.values[i], x);
+  }
+  return(0);
+}
+int POLYNOMIAL_FUNC(mul_scalar) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  POLYNOMIAL_FUNC(cpy) (poly, output);
+  ret=POLYNOMIAL_FUNC(mul_scalar_inplace) (x, output);
+  return(ret);
+}
+
+// multiply by a monomial
+int POLYNOMIAL_FUNC(mul_monomial_inplace) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME* poly){
+  unsigned int i;
+  if(poly->coefficients.length != poly->orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<poly->coefficients.length;i++){
+    POLYNOMIAL_COEF_MUL(poly->coefficients.values[i], poly->coefficients.values[i], x);
+    POLYNOMIAL_ORDER_ADD(poly->orders.values[i], poly->orders.values[i], order);
+  }
+  return(0);
+}
+int POLYNOMIAL_FUNC(mul_monomial) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  POLYNOMIAL_FUNC(cpy) (poly,output);
+  ret=POLYNOMIAL_FUNC(mul_monomial_inplace) (x, order, output);
+  return(ret);
+}
+
+// multiply two polynomials
+int POLYNOMIAL_FUNC(mul_inplace) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME* poly2){
+  unsigned int i;
+  int ret;
+  if(poly1.coefficients.length != poly1.orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<poly1.coefficients.length;i++){
+    ret=POLYNOMIAL_FUNC(mul_monomial_inplace) (poly1.coefficients.values[i], poly1.orders.values[i], poly2);
+    if(ret<0){
+      return(ret);
+    }
+  }
+  return(0);
+}
+int POLYNOMIAL_FUNC(mul) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME poly2, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  POLYNOMIAL_FUNC(cpy) (poly2, output);
+  ret=POLYNOMIAL_FUNC(mul_inplace) (poly1, output);
+  return(ret);
+}
+
+// derive
+int POLYNOMIAL_FUNC(derive_inplace) (POLYNOMIAL_TYPENAME* poly){
+  unsigned int i;
+  if(poly->coefficients.length != poly->orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+  for(i=0;i<poly->coefficients.length;i++){
+    if(POLYNOMIAL_ORDER_CMP_UI(poly->orders.values[i], 0)>0){
+      POLYNOMIAL_COEF_MUL_ORDER(poly->coefficients.values[i], poly->coefficients.values[i], poly->orders.values[i]);
+      POLYNOMIAL_ORDER_SUB_UI(poly->orders.values[i], poly->orders.values[i], 1);
+    }
+    else{
+      // remove the term by setting it to the last term and reducing the length
+      #ifdef POLYNOMIAL_COEF_FREE
+	POLYNOMIAL_COEF_FREE(poly->coefficients.values[i]);
+      #endif
+      #ifdef POLYNOMIAL_ORDER_FREE
+	POLYNOMIAL_ORDER_FREE(poly->orders.values[i]);
+      #endif
+      if(i<poly->coefficients.length-1){
+	POLYNOMIAL_COEF_SET(poly->coefficients.values[i], poly->coefficients.values[poly->coefficients.length-1]);
+	POLYNOMIAL_ORDER_SET(poly->orders.values[i], poly->orders.values[poly->coefficients.length-1]);
+	i--;
+      }
+      poly->coefficients.length--;
+      poly->orders.length--;
+    }
+  }
+  return(0);
+}
+int POLYNOMIAL_FUNC(derive) (POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output){
+  int ret;
+  POLYNOMIAL_FUNC(cpy) (poly, output);
+  ret=POLYNOMIAL_FUNC(derive_inplace) (output);
+  return(ret);
+}
+
+// evaluate
+int POLYNOMIAL_FUNC(evaluate) (POLYNOMIAL_COEF_TYPE* out, POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME poly){
+  unsigned int i;
+  POLYNOMIAL_COEF_TYPE tmp;
+
+  #ifdef POLYNOMIAL_COEF_INIT
+    POLYNOMIAL_COEF_INIT(tmp);
+  #endif
+
+  if(poly.coefficients.length != poly.orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  POLYNOMIAL_COEF_SET_UI(*out, 0);
+  for(i=0;i<poly.coefficients.length;i++){
+    POLYNOMIAL_COEF_POW_UI(tmp, x, poly.orders.values[i]);
+    POLYNOMIAL_COEF_MUL(tmp, tmp, poly.coefficients.values[i]);
+    POLYNOMIAL_COEF_ADD(*out, *out, tmp);
+  }
+  #ifdef POLYNOMIAL_COEF_FREE
+    POLYNOMIAL_COEF_FREE(tmp);
+  #endif
+  return(0);
+}
+
+// print
+int POLYNOMIAL_FUNC(print) (POLYNOMIAL_TYPENAME poly){
+  unsigned int j;
+
+  if(poly.coefficients.length != poly.orders.length){
+    return(LIBINUM_ERROR_SIZE_MISMATCH);
+  }
+
+  printf("  ");
+  for(j=0;j<poly.coefficients.length;j++){
+    if(j>0){
+      printf("+ ");
+    }
+    POLYNOMIAL_COEF_PRINT(poly.coefficients.values[j]);
+    printf(" x^");
+    POLYNOMIAL_ORDER_PRINT(poly.orders.values[j]);
+    printf("\n");
+  }
+
+  return(0);
+}
+
+
+// n-th Legendre polynomial
+int POLYNOMIAL_FUNC(legendre) (unsigned int n, POLYNOMIAL_TYPENAME* output){
+  POLYNOMIAL_TYPENAME prevpoly, tmppoly;
+  unsigned int i;
+  int j;
+  POLYNOMIAL_COEF_TYPE tmp;
+  POLYNOMIAL_ORDER_TYPE tmp_o;
+  
+  if(n==0){
+    POLYNOMIAL_FUNC(init) (output, 1);
+    POLYNOMIAL_FUNC(add_monomial_dui) (1., 0, output);
+    return(0);
+  }
+  else if(n==1){
+    POLYNOMIAL_FUNC(init) (output, n);
+    POLYNOMIAL_FUNC(add_monomial_dui) (1., 1, output);
+    return(0);
+  }
+
+  #ifdef POLYNOMIAL_COEF_INIT
+    POLYNOMIAL_COEF_INIT(tmp);
+  #endif
+  #ifdef POLYNOMIAL_ORDER_INIT
+    POLYNOMIAL_ORDER_INIT(tmp_o);
+  #endif
+
+  // init: prevpoly=1
+  POLYNOMIAL_FUNC(init) (&prevpoly, n-1);
+  POLYNOMIAL_FUNC(add_monomial_dui) (1., 0, &prevpoly);
+  // init: output=x
+  POLYNOMIAL_FUNC(init) (output, n);
+  POLYNOMIAL_FUNC(add_monomial_dui) (1., 1, output);
+
+  // implement i*p_i=(2*i-1)*p_{i-1}-(i-1)*p_{i-2}
+  for(i=2;i<=n;i++){
+    // save current poly to copy it to prevpoly at the end of the loop
+    POLYNOMIAL_FUNC(cpy) (*output, &tmppoly);
+
+    // (2*i-1)/i*p_{i-1}
+    POLYNOMIAL_COEF_SET_UI(tmp, 2*i-1);
+    POLYNOMIAL_COEF_DIV_UI(tmp, tmp, i);
+    POLYNOMIAL_ORDER_SET_UI(tmp_o, 1);
+    POLYNOMIAL_FUNC(mul_monomial_inplace) (tmp, tmp_o, output);
+
+    // -(i-1)/i*p_{i-2}
+    // copy i to a signed int
+    j=i;
+    POLYNOMIAL_COEF_SET_SI(tmp, -j+1);
+    POLYNOMIAL_COEF_DIV_UI(tmp, tmp, i);
+    POLYNOMIAL_FUNC(mul_scalar_inplace) (tmp, &prevpoly);
+
+    // add it to p_n
+    POLYNOMIAL_FUNC(add_inplace) (prevpoly, output);
+
+    // replace prevpoly
+    POLYNOMIAL_FUNC(free) (prevpoly);
+    prevpoly=tmppoly;
+  }
+
+  POLYNOMIAL_FUNC(free) (prevpoly);
+  #ifdef POLYNOMIAL_COEF_FREE
+    POLYNOMIAL_COEF_FREE(tmp);
+  #endif
+  #ifdef POLYNOMIAL_ORDER_FREE
+    POLYNOMIAL_ORDER_FREE(tmp_o);
+  #endif
+
+  return(0);
+}
+
diff --git a/src/polynomial_base.h b/src/polynomial_base.h
new file mode 100644
index 0000000..3165a48
--- /dev/null
+++ b/src/polynomial_base.h
@@ -0,0 +1,67 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for real polynomials
+
+  see polynomial_*.h for the values taken by POLYNOMIAL_FUNC, etc...
+*/
+
+// init
+int POLYNOMIAL_FUNC(init) (POLYNOMIAL_TYPENAME* poly, unsigned int memory);
+int POLYNOMIAL_FUNC(free) (POLYNOMIAL_TYPENAME poly);
+
+// resize memory
+int POLYNOMIAL_FUNC(resize) (POLYNOMIAL_TYPENAME* poly, unsigned int newsize);
+
+// add a monomial
+int POLYNOMIAL_FUNC(add_monomial) (POLYNOMIAL_COEF_TYPE val, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME* output);
+// from a double-valued coefficient and an unsigned int order
+int POLYNOMIAL_FUNC(add_monomial_dui) (double val, unsigned int order, POLYNOMIAL_TYPENAME* output);
+
+// copy
+int POLYNOMIAL_FUNC(cpy) (POLYNOMIAL_TYPENAME input, POLYNOMIAL_TYPENAME* output);
+int POLYNOMIAL_FUNC(cpy_noinit) (POLYNOMIAL_TYPENAME input, POLYNOMIAL_TYPENAME* output);
+
+// add
+int POLYNOMIAL_FUNC(add_inplace) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME* poly2);
+int POLYNOMIAL_FUNC(add) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME poly2, POLYNOMIAL_TYPENAME* output);
+
+// multiply by a scalar
+int POLYNOMIAL_FUNC(mul_scalar_inplace) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME* poly);
+int POLYNOMIAL_FUNC(mul_scalar) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output);
+
+// multiply by a monomial
+int POLYNOMIAL_FUNC(mul_monomial_inplace) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME* poly);
+int POLYNOMIAL_FUNC(mul_monomial) (POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_ORDER_TYPE order, POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output);
+
+// multiply two polynomials
+int POLYNOMIAL_FUNC(mul_inplace) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME* poly2);
+int POLYNOMIAL_FUNC(mul) (POLYNOMIAL_TYPENAME poly1, POLYNOMIAL_TYPENAME poly2, POLYNOMIAL_TYPENAME* output);
+
+// derive
+int POLYNOMIAL_FUNC(derive_inplace) (POLYNOMIAL_TYPENAME* poly);
+int POLYNOMIAL_FUNC(derive) (POLYNOMIAL_TYPENAME poly, POLYNOMIAL_TYPENAME* output);
+
+// evaluate
+int POLYNOMIAL_FUNC(evaluate) (POLYNOMIAL_COEF_TYPE* out, POLYNOMIAL_COEF_TYPE x, POLYNOMIAL_TYPENAME poly);
+
+// print
+int POLYNOMIAL_FUNC(print) (POLYNOMIAL_TYPENAME poly);
+
+
+// n-th Legendre polynomial
+int POLYNOMIAL_FUNC(legendre) (unsigned int n, POLYNOMIAL_TYPENAME* output);
diff --git a/src/polynomial_double.h b/src/polynomial_double.h
new file mode 100644
index 0000000..f876119
--- /dev/null
+++ b/src/polynomial_double.h
@@ -0,0 +1,117 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for real polynomials with double coefficients and unsigned integer powers
+*/
+
+
+// reset CPP macros
+#undef POLYNOMIAL_TYPENAME
+#undef POLYNOMIAL_FUNC
+#undef POLYNOMIAL_COEF_TYPE
+#undef POLYNOMIAL_COEFSARRAY_TYPE
+#undef POLYNOMIAL_COEFSARRAY_FUNC
+#undef POLYNOMIAL_COEF_INIT
+#undef POLYNOMIAL_COEF_FREE
+#undef POLYNOMIAL_COEF_SET
+#undef POLYNOMIAL_COEF_SET_SI
+#undef POLYNOMIAL_COEF_SET_UI
+#undef POLYNOMIAL_COEF_CPY
+#undef POLYNOMIAL_COEF_CPY_D
+#undef POLYNOMIAL_COEF_ADD
+#undef POLYNOMIAL_COEF_ADD_D
+#undef POLYNOMIAL_COEF_MUL
+#undef POLYNOMIAL_COEF_MUL_ORDER
+#undef POLYNOMIAL_COEF_DIV_UI
+#undef POLYNOMIAL_COEF_POW_UI
+#undef POLYNOMIAL_COEF_PRINT
+#undef POLYNOMIAL_ORDER_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_FUNC
+#undef POLYNOMIAL_ORDER_INIT
+#undef POLYNOMIAL_ORDER_FREE
+#undef POLYNOMIAL_ORDER_SET
+#undef POLYNOMIAL_ORDER_SET_UI
+#undef POLYNOMIAL_ORDER_CPY
+#undef POLYNOMIAL_ORDER_CPY_UI
+#undef POLYNOMIAL_ORDER_ADD
+#undef POLYNOMIAL_ORDER_SUB_UI
+#undef POLYNOMIAL_ORDER_CMP
+#undef POLYNOMIAL_ORDER_CMP_UI
+#undef POLYNOMIAL_ORDER_PRINT
+
+// name of the polynomial type
+#define POLYNOMIAL_TYPENAME polynomial_double
+// prefix of function names
+#define POLYNOMIAL_FUNC(NAME) polynomial_double_ ## NAME
+
+// type of the coefficient
+#define POLYNOMIAL_COEF_TYPE double
+// type of coefficient arrays
+#define POLYNOMIAL_COEFSARRAY_TYPE array_double
+// prefix of coefficient array function names
+#define POLYNOMIAL_COEFSARRAY_FUNC(NAME) array_double_ ## NAME
+// set coefficient
+#define POLYNOMIAL_COEF_SET(COEF, VAL) COEF=VAL
+// set coefficient from signed int
+#define POLYNOMIAL_COEF_SET_SI(COEF, VAL) COEF=(double)VAL
+// set coefficient from unsigned int
+#define POLYNOMIAL_COEF_SET_UI(COEF, VAL) COEF=(double)VAL
+// copy coefficient
+#define POLYNOMIAL_COEF_CPY(VAL, COEF) COEF=VAL
+// copy coefficient from double
+#define POLYNOMIAL_COEF_CPY_D(VAL, COEF) COEF=VAL
+// add coefficients
+#define POLYNOMIAL_COEF_ADD(COEF, VAL1, VAL2) COEF=VAL1+VAL2
+// add coefficients, one of which is specified as a double
+#define POLYNOMIAL_COEF_ADD_D(COEF, VAL1, VAL2) COEF=VAL1+VAL2
+// multiply coefficients
+#define POLYNOMIAL_COEF_MUL(COEF, VAL1, VAL2) COEF=VAL1*VAL2
+// multiply coefficients, one of which is specified as a POLYNOMIAL_ORDER_TYPE
+#define POLYNOMIAL_COEF_MUL_ORDER(COEF, VAL1, VAL2) COEF=VAL1*(double)VAL2
+// divide coefficients, one of which is specified as an unsigned int
+#define POLYNOMIAL_COEF_DIV_UI(COEF, VAL1, VAL2) COEF=VAL1/(double)VAL2
+// power of a coefficient, specified as an unsigned int
+#define POLYNOMIAL_COEF_POW_UI(COEF, VAL1, VAL2) COEF=pow(VAL1, (double)VAL2)
+// print a coefficient
+#define POLYNOMIAL_COEF_PRINT(COEF) fprint_double(stdout, COEF)
+
+// type of the order
+#define POLYNOMIAL_ORDER_TYPE unsigned int
+// type of order arrays
+#define POLYNOMIAL_ORDERSARRAY_TYPE array_uint
+// prefix of order array function names
+#define POLYNOMIAL_ORDERSARRAY_FUNC(NAME) array_uint_ ## NAME
+// set order
+#define POLYNOMIAL_ORDER_SET(ORDER, VAL) ORDER=VAL
+// set order from unsigned int
+#define POLYNOMIAL_ORDER_SET_UI(ORDER, VAL) ORDER=VAL
+// copy order
+#define POLYNOMIAL_ORDER_CPY(VAL, ORDER) ORDER=VAL
+// copy order from unsigned int
+#define POLYNOMIAL_ORDER_CPY_UI(VAL, ORDER) ORDER=VAL
+// add orders
+#define POLYNOMIAL_ORDER_ADD(ORDER, VAL1, VAL2) ORDER=VAL1+VAL2
+// subtract an order and an unsigned integer
+#define POLYNOMIAL_ORDER_SUB_UI(ORDER, VAL1, VAL2) ORDER=VAL1-VAL2
+// compare orders (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP(VAL1, VAL2) (int)VAL1-(int)VAL2
+// compare orders, one of which is specified as an unsigned int (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP_UI(VAL1, VAL2) (int)VAL1-(int)VAL2
+// print orders
+#define POLYNOMIAL_ORDER_PRINT(ORDER) printf("%u",ORDER)
+
diff --git a/src/polynomial_ldouble.h b/src/polynomial_ldouble.h
new file mode 100644
index 0000000..54e5501
--- /dev/null
+++ b/src/polynomial_ldouble.h
@@ -0,0 +1,117 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for real polynomials with long double coefficients and unsigned integer powers
+*/
+
+
+// reset CPP macros
+#undef POLYNOMIAL_TYPENAME
+#undef POLYNOMIAL_FUNC
+#undef POLYNOMIAL_COEF_TYPE
+#undef POLYNOMIAL_COEFSARRAY_TYPE
+#undef POLYNOMIAL_COEFSARRAY_FUNC
+#undef POLYNOMIAL_COEF_INIT
+#undef POLYNOMIAL_COEF_FREE
+#undef POLYNOMIAL_COEF_SET
+#undef POLYNOMIAL_COEF_SET_SI
+#undef POLYNOMIAL_COEF_SET_UI
+#undef POLYNOMIAL_COEF_CPY
+#undef POLYNOMIAL_COEF_CPY_D
+#undef POLYNOMIAL_COEF_ADD
+#undef POLYNOMIAL_COEF_ADD_D
+#undef POLYNOMIAL_COEF_MUL
+#undef POLYNOMIAL_COEF_MUL_ORDER
+#undef POLYNOMIAL_COEF_DIV_UI
+#undef POLYNOMIAL_COEF_POW_UI
+#undef POLYNOMIAL_COEF_PRINT
+#undef POLYNOMIAL_ORDER_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_FUNC
+#undef POLYNOMIAL_ORDER_INIT
+#undef POLYNOMIAL_ORDER_FREE
+#undef POLYNOMIAL_ORDER_SET
+#undef POLYNOMIAL_ORDER_SET_UI
+#undef POLYNOMIAL_ORDER_CPY
+#undef POLYNOMIAL_ORDER_CPY_UI
+#undef POLYNOMIAL_ORDER_ADD
+#undef POLYNOMIAL_ORDER_SUB_UI
+#undef POLYNOMIAL_ORDER_CMP
+#undef POLYNOMIAL_ORDER_CMP_UI
+#undef POLYNOMIAL_ORDER_PRINT
+
+// name of the polynomial type
+#define POLYNOMIAL_TYPENAME polynomial_ldouble
+// prefix of function names
+#define POLYNOMIAL_FUNC(NAME) polynomial_ldouble_ ## NAME
+
+// type of the coefficient
+#define POLYNOMIAL_COEF_TYPE long double
+// type of coefficient arrays
+#define POLYNOMIAL_COEFSARRAY_TYPE array_ldouble
+// prefix of coefficient array function names
+#define POLYNOMIAL_COEFSARRAY_FUNC(NAME) array_ldouble_ ## NAME
+// set coefficient
+#define POLYNOMIAL_COEF_SET(COEF, VAL) COEF=VAL
+// set coefficient from signed int
+#define POLYNOMIAL_COEF_SET_SI(COEF, VAL) COEF=(long double)VAL
+// set coefficient from unsigned int
+#define POLYNOMIAL_COEF_SET_UI(COEF, VAL) COEF=(long double)VAL
+// copy coefficient
+#define POLYNOMIAL_COEF_CPY(VAL, COEF) COEF=VAL
+// copy coefficient from double
+#define POLYNOMIAL_COEF_CPY_D(VAL, COEF) COEF=(long double)VAL
+// add coefficients
+#define POLYNOMIAL_COEF_ADD(COEF, VAL1, VAL2) COEF=VAL1+VAL2
+// add coefficients, one of which is specified as a double
+#define POLYNOMIAL_COEF_ADD_D(COEF, VAL1, VAL2) COEF=VAL1+(long double)VAL2
+// multiply coefficients
+#define POLYNOMIAL_COEF_MUL(COEF, VAL1, VAL2) COEF=VAL1*VAL2
+// multiply coefficients, one of which is specified as a POLYNOMIAL_ORDER_TYPE
+#define POLYNOMIAL_COEF_MUL_ORDER(COEF, VAL1, VAL2) COEF=VAL1*(long double)VAL2
+// divide coefficients, one of which is specified as an unsigned int
+#define POLYNOMIAL_COEF_DIV_UI(COEF, VAL1, VAL2) COEF=VAL1/(long double)VAL2
+// power of a coefficient, specified as an unsigned int
+#define POLYNOMIAL_COEF_POW_UI(COEF, VAL1, VAL2) COEF=powl(VAL1, (long double)VAL2)
+// print a coefficient
+#define POLYNOMIAL_COEF_PRINT(COEF) fprint_ldouble(stdout, COEF)
+
+// type of the order
+#define POLYNOMIAL_ORDER_TYPE unsigned int
+// type of order arrays
+#define POLYNOMIAL_ORDERSARRAY_TYPE array_uint
+// prefix of order array function names
+#define POLYNOMIAL_ORDERSARRAY_FUNC(NAME) array_uint_ ## NAME
+// set order
+#define POLYNOMIAL_ORDER_SET(ORDER, VAL) ORDER=VAL
+// set order from unsigned int
+#define POLYNOMIAL_ORDER_SET_UI(ORDER, VAL) ORDER=VAL
+// copy order
+#define POLYNOMIAL_ORDER_CPY(VAL, ORDER) ORDER=VAL
+// copy order from unsigned int
+#define POLYNOMIAL_ORDER_CPY_UI(VAL, ORDER) ORDER=VAL
+// add orders
+#define POLYNOMIAL_ORDER_ADD(ORDER, VAL1, VAL2) ORDER=VAL1+VAL2
+// subtract an order and an unsigned integer
+#define POLYNOMIAL_ORDER_SUB_UI(ORDER, VAL1, VAL2) ORDER=VAL1-VAL2
+// compare orders (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP(VAL1, VAL2) (int)VAL1-(int)VAL2
+// compare orders, one of which is specified as an unsigned int (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP_UI(VAL1, VAL2) (int)VAL1-(int)VAL2
+// print orders
+#define POLYNOMIAL_ORDER_PRINT(ORDER) printf("%u",ORDER)
+
diff --git a/src/polynomial_mpfr.h b/src/polynomial_mpfr.h
new file mode 100644
index 0000000..f8bc983
--- /dev/null
+++ b/src/polynomial_mpfr.h
@@ -0,0 +1,121 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for real polynomials with multi-precision float (mpfr) coefficients and unsigned integer powers
+*/
+
+
+// reset CPP macros
+#undef POLYNOMIAL_TYPENAME
+#undef POLYNOMIAL_FUNC
+#undef POLYNOMIAL_COEF_TYPE
+#undef POLYNOMIAL_COEFSARRAY_TYPE
+#undef POLYNOMIAL_COEFSARRAY_FUNC
+#undef POLYNOMIAL_COEF_INIT
+#undef POLYNOMIAL_COEF_FREE
+#undef POLYNOMIAL_COEF_SET
+#undef POLYNOMIAL_COEF_SET_SI
+#undef POLYNOMIAL_COEF_SET_UI
+#undef POLYNOMIAL_COEF_CPY
+#undef POLYNOMIAL_COEF_CPY_D
+#undef POLYNOMIAL_COEF_ADD
+#undef POLYNOMIAL_COEF_ADD_D
+#undef POLYNOMIAL_COEF_MUL
+#undef POLYNOMIAL_COEF_MUL_ORDER
+#undef POLYNOMIAL_COEF_DIV_UI
+#undef POLYNOMIAL_COEF_POW_UI
+#undef POLYNOMIAL_COEF_PRINT
+#undef POLYNOMIAL_ORDER_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_TYPE
+#undef POLYNOMIAL_ORDERSARRAY_FUNC
+#undef POLYNOMIAL_ORDER_INIT
+#undef POLYNOMIAL_ORDER_FREE
+#undef POLYNOMIAL_ORDER_SET
+#undef POLYNOMIAL_ORDER_SET_UI
+#undef POLYNOMIAL_ORDER_CPY
+#undef POLYNOMIAL_ORDER_CPY_UI
+#undef POLYNOMIAL_ORDER_ADD
+#undef POLYNOMIAL_ORDER_SUB_UI
+#undef POLYNOMIAL_ORDER_CMP
+#undef POLYNOMIAL_ORDER_CMP_UI
+#undef POLYNOMIAL_ORDER_PRINT
+
+// name of the polynomial type
+#define POLYNOMIAL_TYPENAME polynomial_mpfr
+// prefix of function names
+#define POLYNOMIAL_FUNC(NAME) polynomial_mpfr_ ## NAME
+
+// type of the coefficient
+#define POLYNOMIAL_COEF_TYPE mpfr_t
+// type of coefficient arrays
+#define POLYNOMIAL_COEFSARRAY_TYPE array_mpfr
+// prefix of coefficient array function names
+#define POLYNOMIAL_COEFSARRAY_FUNC(NAME) array_mpfr_ ## NAME
+// init coefficient
+#define POLYNOMIAL_COEF_INIT(VAR) mpfr_init(VAR)
+// free coefficient
+#define POLYNOMIAL_COEF_FREE(VAR) mpfr_clear(VAR)
+// set coefficient
+#define POLYNOMIAL_COEF_SET(COEF, VAL) COEF[0]=VAL[0]
+// set coefficient from signed int
+#define POLYNOMIAL_COEF_SET_SI(COEF, VAL) mpfr_set_si(COEF, VAL, MPFR_RNDN)
+// set coefficient from unsigned int
+#define POLYNOMIAL_COEF_SET_UI(COEF, VAL) mpfr_set_ui(COEF, VAL, MPFR_RNDN)
+// copy coefficient
+#define POLYNOMIAL_COEF_CPY(VAL, COEF) mpfr_init(COEF); mpfr_set(COEF, VAL, MPFR_RNDN)
+// copy coefficient from double
+#define POLYNOMIAL_COEF_CPY_D(VAL, COEF) mpfr_init(COEF); mpfr_set_d(COEF, VAL, MPFR_RNDN)
+// add coefficients
+#define POLYNOMIAL_COEF_ADD(COEF, VAL1, VAL2) mpfr_add(COEF, VAL1, VAL2, MPFR_RNDN)
+// add coefficients, one of which is specified as a double
+#define POLYNOMIAL_COEF_ADD_D(COEF, VAL1, VAL2) mpfr_add_d(COEF, VAL1, VAL2, MPFR_RNDN)
+// multiply coefficients
+#define POLYNOMIAL_COEF_MUL(COEF, VAL1, VAL2) mpfr_mul(COEF, VAL1, VAL2, MPFR_RNDN)
+// multiply coefficients, one of which is specified as a POLYNOMIAL_ORDER_TYPE
+#define POLYNOMIAL_COEF_MUL_ORDER(COEF, VAL1, VAL2) mpfr_mul_ui(COEF, VAL1, VAL2, MPFR_RNDN)
+// divide coefficients, one of which is specified as an unsigned int
+#define POLYNOMIAL_COEF_DIV_UI(COEF, VAL1, VAL2) mpfr_div_ui(COEF, VAL1, VAL2, MPFR_RNDN)
+// power of a coefficient, specified as an unsigned int
+#define POLYNOMIAL_COEF_POW_UI(COEF, VAL1, VAL2) mpfr_pow_ui(COEF, VAL1, VAL2, MPFR_RNDN)
+// print a coefficient
+#define POLYNOMIAL_COEF_PRINT(COEF) fprint_mpfr(stdout, COEF)
+
+// type of the order
+#define POLYNOMIAL_ORDER_TYPE unsigned int
+// type of order arrays
+#define POLYNOMIAL_ORDERSARRAY_TYPE array_uint
+// prefix of order array function names
+#define POLYNOMIAL_ORDERSARRAY_FUNC(NAME) array_uint_ ## NAME
+// set order
+#define POLYNOMIAL_ORDER_SET(ORDER, VAL) ORDER=VAL
+// set order from unsigned int
+#define POLYNOMIAL_ORDER_SET_UI(ORDER, VAL) ORDER=VAL
+// copy order
+#define POLYNOMIAL_ORDER_CPY(VAL, ORDER) ORDER=VAL
+// copy order from unsigned int
+#define POLYNOMIAL_ORDER_CPY_UI(VAL, ORDER) ORDER=VAL
+// add orders
+#define POLYNOMIAL_ORDER_ADD(ORDER, VAL1, VAL2) ORDER=VAL1+VAL2
+// subtract an order and an unsigned integer
+#define POLYNOMIAL_ORDER_SUB_UI(ORDER, VAL1, VAL2) ORDER=VAL1-VAL2
+// compare orders (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP(VAL1, VAL2) (int)VAL1-(int)VAL2
+// compare orders, one of which is specified as an unsigned int (0 if equal, -1 if <, +1 if >)
+#define POLYNOMIAL_ORDER_CMP_UI(VAL1, VAL2) (int)VAL1-(int)VAL2
+// print orders
+#define POLYNOMIAL_ORDER_PRINT(ORDER) printf("%u",ORDER)
+
diff --git a/src/polynomial_type.h b/src/polynomial_type.h
new file mode 100644
index 0000000..3361d67
--- /dev/null
+++ b/src/polynomial_type.h
@@ -0,0 +1,27 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Structure definition for real polynomials
+
+  see polynomial_*.h for the values taken by POLYNOMIAL_TYPENAME, etc...
+*/
+
+typedef struct POLYNOMIAL_TYPENAME {
+  POLYNOMIAL_COEFSARRAY_TYPE coefficients;
+  POLYNOMIAL_ORDERSARRAY_TYPE orders;
+} POLYNOMIAL_TYPENAME;
+
diff --git a/src/root.c b/src/root.c
new file mode 100644
index 0000000..c79745f
--- /dev/null
+++ b/src/root.c
@@ -0,0 +1,51 @@
+/*
+Copyright 2016 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 "root.h"
+
+#include <mpfr.h>
+#include <math.h>
+#include "errors.h"
+
+
+//--------------------------------------------------
+//
+//  using doubles
+//
+//--------------------------------------------------
+
+#include "root_double.h"
+#include "root_base.c"
+
+
+//--------------------------------------------------
+//
+//  using long doubles
+//
+//--------------------------------------------------
+
+#include "root_ldouble.h"
+#include "root_base.c"
+
+
+//--------------------------------------------------
+//
+//  using mpfr floats
+//
+//--------------------------------------------------
+
+#include "root_mpfr.h"
+#include "root_base.c"
diff --git a/src/root.h b/src/root.h
new file mode 100644
index 0000000..7106389
--- /dev/null
+++ b/src/root.h
@@ -0,0 +1,57 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Root finding
+*/
+
+#ifndef LIBINUM_ROOT_H
+#define LIBINUM_ROOT_H
+
+#include "types.h"
+
+
+//--------------------------------------------------
+//
+//  using doubles
+//
+//--------------------------------------------------
+
+#include "root_double.h"
+#include "root_base.h"
+
+
+//--------------------------------------------------
+//
+//  using long doubles
+//
+//--------------------------------------------------
+
+#include "root_ldouble.h"
+#include "root_base.h"
+
+
+//--------------------------------------------------
+//
+//  using mpfr floats
+//
+//--------------------------------------------------
+
+#include "root_mpfr.h"
+#include "root_base.h"
+
+#endif
+
diff --git a/src/root_base.c b/src/root_base.c
new file mode 100644
index 0000000..b13d9db
--- /dev/null
+++ b/src/root_base.c
@@ -0,0 +1,134 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for root finding
+
+  see integral_*.h for the values taken by ROOT_FUNC, etc...
+*/
+
+
+// compute the root of a real function of 1 variable using the Newton method
+int ROOT_FUNC(root_newton_inplace) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args){
+  unsigned int count=0;
+  ROOT_FLOAT_TYPE valf, valdf, delta, tmp;
+  int ret;
+
+  #ifdef ROOT_FLOAT_INIT
+    ROOT_FLOAT_INIT(valf);
+    ROOT_FLOAT_INIT(valdf);
+    ROOT_FLOAT_INIT(delta);
+    ROOT_FLOAT_INIT(tmp);
+  #endif
+
+  // evaluate at guess
+  ret=(*func)(&valf, *out, extra_args);
+  if(ret<0){
+    #ifdef ROOT_FLOAT_FREE
+      ROOT_FLOAT_FREE(valf);
+      ROOT_FLOAT_FREE(valdf);
+      ROOT_FLOAT_FREE(delta);
+      ROOT_FLOAT_FREE(tmp);
+    #endif
+    return(ret);
+  }
+  // check that valf is a number
+  if(! ROOT_FLOAT_ISNUMBER(valf)){
+    #ifdef ROOT_FLOAT_FREE
+      ROOT_FLOAT_FREE(valf);
+      ROOT_FLOAT_FREE(valdf);
+      ROOT_FLOAT_FREE(delta);
+      ROOT_FLOAT_FREE(tmp);
+    #endif
+    return(LIBINUM_ERROR_NAN);
+  }
+
+  ROOT_FLOAT_ABS(delta, valf);
+
+  // loop until tolerance is reached
+  while(ROOT_FLOAT_CMP(delta, tolerance) > 0 && count < maxiter){
+    // new guess
+    (*deriv_func)(&valdf, *out, extra_args);
+    if(ret<0){
+      #ifdef ROOT_FLOAT_FREE
+	ROOT_FLOAT_FREE(valf);
+	ROOT_FLOAT_FREE(valdf);
+	ROOT_FLOAT_FREE(delta);
+	ROOT_FLOAT_FREE(tmp);
+      #endif
+      return(ret);
+    }
+    // check that valdf is a number that is not 0
+    if(! ROOT_FLOAT_ISNUMBER(valdf) || ROOT_FLOAT_ISZERO(valdf)){
+      #ifdef ROOT_FLOAT_FREE
+	ROOT_FLOAT_FREE(valf);
+	ROOT_FLOAT_FREE(valdf);
+	ROOT_FLOAT_FREE(delta);
+	ROOT_FLOAT_FREE(tmp);
+      #endif
+      return(LIBINUM_ERROR_NAN);
+    }
+
+    ROOT_FLOAT_DIV(tmp, valf, valdf);
+    ROOT_FLOAT_SUB(*out, *out, tmp);
+    
+    // evaluate at new guess
+    (*func)(&valf, *out, extra_args);
+    if(ret<0){
+      #ifdef ROOT_FLOAT_FREE
+	ROOT_FLOAT_FREE(valf);
+	ROOT_FLOAT_FREE(valdf);
+	ROOT_FLOAT_FREE(delta);
+	ROOT_FLOAT_FREE(tmp);
+      #endif
+      return(ret);
+    }
+    // check that valf is a number
+    if(! ROOT_FLOAT_ISNUMBER(valf)){
+      #ifdef ROOT_FLOAT_FREE
+	ROOT_FLOAT_FREE(valf);
+	ROOT_FLOAT_FREE(valdf);
+	ROOT_FLOAT_FREE(delta);
+	ROOT_FLOAT_FREE(tmp);
+      #endif
+      return(LIBINUM_ERROR_NAN);
+    }
+
+    ROOT_FLOAT_ABS(delta, valf);
+
+    // increase counter
+    count++;
+  }
+
+  // free variables
+  #ifdef ROOT_FLOAT_FREE
+    ROOT_FLOAT_FREE(valf);
+    ROOT_FLOAT_FREE(valdf);
+    ROOT_FLOAT_FREE(delta);
+    ROOT_FLOAT_FREE(tmp);
+  #endif
+
+  // fail if maxiter was reached
+  if(count==maxiter){
+    return(LIBINUM_ERROR_MAXITER);
+  }
+  return(0);
+}
+int ROOT_FUNC(root_newton) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE init, ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args){
+  ROOT_FLOAT_SET(*out, init);
+  ROOT_FUNC(root_newton_inplace) (out, func, deriv_func, tolerance, maxiter, extra_args);
+  return(0);
+}
diff --git a/src/root_base.h b/src/root_base.h
new file mode 100644
index 0000000..7a69616
--- /dev/null
+++ b/src/root_base.h
@@ -0,0 +1,27 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Base functions for root finding
+
+  see integral_*.h for the values taken by ROOT_FUNC, etc...
+*/
+
+
+// compute the root of a real function of 1 variable using the Newton method
+int ROOT_FUNC(root_newton_inplace) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args);
+int ROOT_FUNC(root_newton) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE init, ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args);
+
diff --git a/src/root_double.h b/src/root_double.h
new file mode 100644
index 0000000..f564998
--- /dev/null
+++ b/src/root_double.h
@@ -0,0 +1,55 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for root finding using doubles
+*/
+
+
+// reset CPP macros
+#undef ROOT_FUNC
+#undef ROOT_FLOAT_TYPE
+#undef ROOT_FLOAT_INIT
+#undef ROOT_FLOAT_FREE
+#undef ROOT_FLOAT_SET
+#undef ROOT_FLOAT_SUB
+#undef ROOT_FLOAT_DIV
+#undef ROOT_FLOAT_CMP
+#undef ROOT_FLOAT_ABS
+#undef ROOT_FLOAT_ISNUMBER
+#undef ROOT_FLOAT_ISZERO
+
+
+// suffix of function names
+#define ROOT_FUNC(NAME) NAME ## _double
+
+// type of floats
+#define ROOT_FLOAT_TYPE double
+// set float
+#define ROOT_FLOAT_SET(FLOAT, VAL) FLOAT=VAL
+// subtract floats
+#define ROOT_FLOAT_SUB(FLOAT, VAL1, VAL2) FLOAT=VAL1-VAL2
+// divide floats
+#define ROOT_FLOAT_DIV(FLOAT, VAL1, VAL2) FLOAT=VAL1/VAL2
+// compare floats (0 if equal, -1 if <, +1 if >)
+#define ROOT_FLOAT_CMP(VAL1, VAL2) VAL1-VAL2
+// abs of float
+#define ROOT_FLOAT_ABS(FLOAT, VAL) FLOAT=fabs(VAL)
+// check whether a float is a regular number
+#define ROOT_FLOAT_ISNUMBER(FLOAT) (fpclassify(FLOAT)>=FP_ZERO)
+// whether float is 0
+#define ROOT_FLOAT_ISZERO(VAL) (VAL==0.)
+
diff --git a/src/root_ldouble.h b/src/root_ldouble.h
new file mode 100644
index 0000000..a20643b
--- /dev/null
+++ b/src/root_ldouble.h
@@ -0,0 +1,55 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for root finding using long doubles
+*/
+
+
+// reset CPP macros
+#undef ROOT_FUNC
+#undef ROOT_FLOAT_TYPE
+#undef ROOT_FLOAT_INIT
+#undef ROOT_FLOAT_FREE
+#undef ROOT_FLOAT_SET
+#undef ROOT_FLOAT_SUB
+#undef ROOT_FLOAT_DIV
+#undef ROOT_FLOAT_CMP
+#undef ROOT_FLOAT_ABS
+#undef ROOT_FLOAT_ISNUMBER
+#undef ROOT_FLOAT_ISZERO
+
+
+// suffix of function names
+#define ROOT_FUNC(NAME) NAME ## _ldouble
+
+// type of floats
+#define ROOT_FLOAT_TYPE long double
+// set float
+#define ROOT_FLOAT_SET(FLOAT, VAL) FLOAT=VAL
+// subtract floats
+#define ROOT_FLOAT_SUB(FLOAT, VAL1, VAL2) FLOAT=VAL1-VAL2
+// divide floats
+#define ROOT_FLOAT_DIV(FLOAT, VAL1, VAL2) FLOAT=VAL1/VAL2
+// compare floats (0 if equal, -1 if <, +1 if >)
+#define ROOT_FLOAT_CMP(VAL1, VAL2) VAL1-VAL2
+// abs of float
+#define ROOT_FLOAT_ABS(FLOAT, VAL) FLOAT=fabsl(VAL)
+// check whether a float is a regular number
+#define ROOT_FLOAT_ISNUMBER(FLOAT) (fpclassify(FLOAT)>=FP_ZERO)
+// whether float is 0
+#define ROOT_FLOAT_ISZERO(VAL) (VAL==0.)
+
diff --git a/src/root_mpfr.h b/src/root_mpfr.h
new file mode 100644
index 0000000..6535e0b
--- /dev/null
+++ b/src/root_mpfr.h
@@ -0,0 +1,59 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Preprocessor macros for root finding using multi-precision floats (mpfr)
+*/
+
+
+// reset CPP macros
+#undef ROOT_FUNC
+#undef ROOT_FLOAT_TYPE
+#undef ROOT_FLOAT_INIT
+#undef ROOT_FLOAT_FREE
+#undef ROOT_FLOAT_SET
+#undef ROOT_FLOAT_SUB
+#undef ROOT_FLOAT_DIV
+#undef ROOT_FLOAT_CMP
+#undef ROOT_FLOAT_ABS
+#undef ROOT_FLOAT_ISNUMBER
+#undef ROOT_FLOAT_ISZERO
+
+
+// suffix of function names
+#define ROOT_FUNC(NAME) NAME ## _mpfr
+
+// type of floats
+#define ROOT_FLOAT_TYPE mpfr_t
+// init float
+#define ROOT_FLOAT_INIT(VAR) mpfr_init(VAR)
+// free float
+#define ROOT_FLOAT_FREE(VAR) mpfr_clear(VAR)
+// set float
+#define ROOT_FLOAT_SET(FLOAT, VAL) mpfr_set(FLOAT, VAL, MPFR_RNDN)
+// subtract floats
+#define ROOT_FLOAT_SUB(FLOAT, VAL1, VAL2) mpfr_sub(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// divide floats
+#define ROOT_FLOAT_DIV(FLOAT, VAL1, VAL2) mpfr_div(FLOAT, VAL1, VAL2, MPFR_RNDN)
+// compare floats (0 if equal, -1 if <, +1 if >)
+#define ROOT_FLOAT_CMP(VAL1, VAL2) mpfr_cmp(VAL1, VAL2)
+// abs of float
+#define ROOT_FLOAT_ABS(FLOAT, VAL) mpfr_abs(FLOAT, VAL, MPFR_RNDN)
+// check whether a float is a regular number
+#define ROOT_FLOAT_ISNUMBER(FLOAT) mpfr_number_p(FLOAT)!=0
+// whether float is 0
+#define ROOT_FLOAT_ISZERO(VAL) mpfr_zero_p(VAL)!=0
+
diff --git a/src/types.h b/src/types.h
new file mode 100644
index 0000000..e438133
--- /dev/null
+++ b/src/types.h
@@ -0,0 +1,84 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Types and structures
+*/
+
+#ifndef LIBINUM_TYPES_H
+#define LIBINUM_TYPES_H
+
+#include <mpfr.h>
+#include <pthread.h>
+
+// arrays
+#include "array_int.h"
+#include "array_type.h"
+
+#include "array_uint.h"
+#include "array_type.h"
+
+#include "array_double.h"
+#include "array_type.h"
+
+#include "array_ldouble.h"
+#include "array_type.h"
+
+#include "array_mpfr.h"
+#include "array_type.h"
+
+#include "array_2_mpfr.h"
+#include "array_type.h"
+
+#include "array_char.h"
+#include "array_type.h"
+
+#include "array_str.h"
+#include "array_type.h"
+
+#include "array_pthread_t.h"
+#include "array_type.h"
+
+
+// 1-variable polynomial
+#include "polynomial_mpfr.h"
+#include "polynomial_type.h"
+
+#include "polynomial_double.h"
+#include "polynomial_type.h"
+
+#include "polynomial_ldouble.h"
+#include "polynomial_type.h"
+
+// polynomial array
+#include "array_polynomial_mpfr.h"
+#include "array_type.h"
+
+#include "array_polynomial_double.h"
+#include "array_type.h"
+
+#include "array_polynomial_ldouble.h"
+#include "array_type.h"
+
+
+// multivariable polynomials
+#include "polynomialMV_mpz.h"
+#include "polynomialMV_type.h"
+
+#include "polynomialMV_int.h"
+#include "polynomialMV_type.h"
+
+#endif
diff --git a/src/utils.c b/src/utils.c
new file mode 100644
index 0000000..d1fad7e
--- /dev/null
+++ b/src/utils.c
@@ -0,0 +1,79 @@
+/*
+Copyright 2016 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 "utils.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+// define MPFR_USE_FILE to enable the use of mpfr_printf
+#define MPFR_USE_FILE
+#include <mpfr.h>
+#include <math.h>
+#include "array.h"
+
+// stringify macros
+#define STR(S) XSTR(S)
+#define XSTR(S) #S
+
+// print a double at maximal precision
+int fprint_double(FILE* file, double x){
+  fprintf(file, "% ." STR(__DBL_DIG__) "le",x);
+  return(0);
+}
+
+// print a long double at maximal precision
+int fprint_ldouble(FILE* file, long double x){
+  fprintf(file, "% ." STR(__LDBL_DIG__) "Le",x);
+  return(0);
+}
+
+// log10(2)
+#define LOG2 0.3010299956639812
+// print an mpfr at maximal precision
+int fprint_mpfr(FILE* file, mpfr_t x){
+  // the printf format
+  array_char printf_format;
+  array_char_init(&printf_format,12);
+  array_char_snprintf(&printf_format,"%% .%dRe\0", (int)((mpfr_get_default_prec()-1)*LOG2));
+  mpfr_fprintf(file, printf_format.values, x);
+  array_char_free(printf_format);
+  return(0);
+}
+
+
+// print information about data types
+int print_datatype_info(FILE* file){
+
+  fprintf(file, "       int8_t = " STR(__INT8_TYPE__) "\n");
+  fprintf(file, "      int16_t = " STR(__INT16_TYPE__) "\n");
+  fprintf(file, "      int32_t = " STR(__INT32_TYPE__) "\n");
+  fprintf(file, "      int64_t = " STR(__INT64_TYPE__) "\n");
+
+  fprintf(file, "\n");
+
+  fprintf(file, "     double: precision: " STR(__DBL_MANT_DIG__) ", emax: " STR(__DBL_MAX_EXP__) ", emin: " STR(__DBL_MIN_EXP__) "\n");
+  fprintf(file, "long double: precision: " STR(__LDBL_MANT_DIG__) ", emax: " STR(__LDBL_MAX_EXP__) ", emin: " STR(__LDBL_MIN_EXP__) "\n");
+
+  fprintf(file, "\n");
+
+  #if _MPFR_PREC_FORMAT == 2
+    fprintf(file, "mpfr precision and emax: int\n");
+  #elif _MPFR_PREC_FORMAT == 3
+    fprintf(file, "mpfr precision and emax: long int\n");
+  #endif
+
+  return(0);
+}
diff --git a/src/utils.h b/src/utils.h
new file mode 100644
index 0000000..c3954b4
--- /dev/null
+++ b/src/utils.h
@@ -0,0 +1,36 @@
+/*
+Copyright 2016 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.
+*/
+
+/*
+  Useful functions
+*/
+
+#ifndef LIBINUM_UTILS_H
+#define LIBINUM_UTILS_H
+
+#include <mpfr.h>
+#include <stdio.h>
+
+// print a double at maximal precision
+int fprint_double(FILE* file, double x);
+// print a long double at maximal precision
+int fprint_ldouble(FILE* file, long double x);
+// print an mpfr at maximal precision
+int fprint_mpfr(FILE* file, mpfr_t x);
+
+// print information about data types
+int print_datatype_info();
+#endif