diff options
Diffstat (limited to 'src/polynomial_base.c')
| -rw-r--r-- | src/polynomial_base.c | 337 |
1 files changed, 337 insertions, 0 deletions
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); +} + |