diff options
Diffstat (limited to 'src')
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 |