Laboratoire de l'Informatique et du Parallélisme » Pobyso

/** @file pobyso.c

 * Name & purpose

 *

 * @author S.T.

 * @date 2011-10-12

 *

 */

/******************************************************************************/

/* Headers, applying the "particular to general" convention.*/

#include "pobyso.h"

/* includes of local headers */

/* includes of project headers */

/* includes of system headers */

/* Other declarations */

/* Internal prototypes */

void

pobyso_error_message(char *functionName, char *messageName, char* message);

/* Types, constants and macros definitions */

/* Global variables */

/* Functions */

void

pobyso_autoprint(sollya_obj_t solObj, ...)

{

  va_list va;

  va_start(va, solObj);

  sollya_lib_v_autoprint(solObj, va);

  va_end(va);

} /* End pobyso_autoprint. */

sollya_obj_t

pobyso_parse_string(const char* expression)

{

  if (expression == NULL)

  {

    pobyso_error_message("pobyso_parse_string",

                        "NULL_POINTER_ARGUMENT",

                        "The expression is a NULL pointer");

    return(NULL);

  }

  return(sollya_lib_parse_string(expression));

} /* pobyso_parse_string */

void

pobyso_set_verbosity_off(sollya_obj_t* currentVerbosityLevel)

{

  sollya_obj_t verbosityLevelZero;

  if (currentVerbosityLevel != NULL)

  {

     *currentVerbosityLevel = sollya_lib_get_verbosity();

  }

  verbosityLevelZero = sollya_lib_constant_from_int(0);

  sollya_lib_set_verbosity(verbosityLevelZero);

  sollya_lib_clear_obj(verbosityLevelZero);

} /* End of pobyso_set_verbosity_off. */

void

pobyso_set_verbosity_to(sollya_obj_t newVerbosityLevel)

{

  if (newVerbosityLevel == NULL)

  {

    pobyso_error_message("pobyso_set_verbosity_to",

                        "NULL_POINTER_ARGUMENT",

                        "The new verbosity level is a NULL pointer");

    return;

  }

  sollya_lib_set_verbosity(newVerbosityLevel);

} /* End of pobyso_set_verbosity_to. */

/* Attic from the sollya_lib < 4. */

#if 0

chain*

pobyso_create_canonical_monomials_base(const unsigned int degree)

{

  int i    = 0;

  chain *monomials  = NULL;

  node  *monomial   = NULL;

  for(i = degree ; i >= 0  ; i--)

  {

     monomial   = makePow(makeVariable(), makeConstantDouble((double)i));

     monomials  = addElement(monomials, monomial);

     fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);

  }

  if (monomials == NULL)

  {

    pobyso_error_message("pobyso_create_canonical_monomial_base",

                        "CHAIN_CREATION_ERROR",

                        "Could not create the monomials chain");

    return(NULL);

  }

  return(monomials);

} /* End pobyso_create_canonical_monomials_base. */

chain*

pobyso_create_chain_from_int_array(int* intArray,

                                    const unsigned int arrayLength)

{

  int i = 0;

  chain *newChain = NULL;

  int *currentInt;

  if (arrayLength == 0) return(NULL);

  if (intArray == NULL)

  {

    pobyso_error_message("pobyso_create_chain_from_int_array",

                        "NULL_POINTER_ARGUMENT",

                        "The array is a NULL pointer");

    return(NULL);

  }

  for (i = arrayLength - 1 ; i >= 0 ; i--)

  {

    currentInt = malloc(sizeof(int));

    if (currentInt == NULL)

    {

      pobyso_error_message("pobyso_create_chain_from_int_array",

                          "MEMORY_ALLOCATION_ERROR",

                          "Could not allocate one of the integers");

      freeChain(newChain, free);

      return(NULL);

    }

    *currentInt = intArray[i];

    newChain = addElement(newChain, currentInt);

  }

  return(newChain);

} // End pobyso_create_chain_from_int_array. */

chain*

pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,

                                        const unsigned int arrayLength)

{

  int i = 0;

  chain *newChain = NULL;

  unsigned int *currentInt;

  /* Argument checking. */

  if (arrayLength == 0) return(NULL);

  if (intArray == NULL)

  {

    pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",

                        "NULL_POINTER_ARGUMENT",

                        "The array is a NULL pointer");

    return(NULL);

  }

  for (i = arrayLength - 1 ; i >= 0 ; i--)

  {

    currentInt = malloc(sizeof(unsigned int));

    if (currentInt == NULL)

    {

      pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",

                          "MEMORY_ALLOCATION_ERROR",

                          "Could not allocate one of the integers");

      freeChain(newChain, free);

      return(NULL);

    }

    *currentInt = intArray[i];

    newChain = addElement(newChain, currentInt);

  }

  return(newChain);

} // End pobyso_create_chain_from_unsigned_int_array. */

node*

pobyso_differentiate(node *functionNode)

{

  /* Argument checking. */

  node *differentialNode;

  if (functionNode == NULL)

  {

    pobyso_error_message("pobyso_differentiate",

                        "NULL_POINTER_ARGUMENT",

                        "The function to differentiate is a NULL pointer");

    return(NULL);

  }

  differentialNode = differentiate(functionNode);

  if (differentialNode == NULL)

  {

    pobyso_error_message("pobyso_differentiate",

                        "INTERNAL ERROR",

                        "Sollya could not differentiate the function");

  }

  return(differentialNode);

} // End pobyso_differentiate

int

pobyso_dirty_infnorm(mpfr_t infNorm,

                      node *functionNode,

                      mpfr_t lowerBound,

                      mpfr_t upperBound,

                      mp_prec_t precision)

{

  int functionCallResult;

  /* Arguments checking. */

  if (functionNode == NULL)

  {

    pobyso_error_message("pobyso_dirty_infnorm",

                        "NULL_POINTER_ARGUMENT",

                        "The function to compute is a NULL pointer");

    return(1);

  }

  if (mpfr_cmp(lowerBound, upperBound) > 0)

  {

    pobyso_error_message("pobyso_dirty_infnorm",

                        "INCOHERENT_INPUT_DATA",

                        "The lower bond is greater than the upper bound");

    return(1);

  }

  /* Particular cases. */

  if (mpfr_cmp(lowerBound, upperBound) == 0)

  {

    functionCallResult = pobyso_evaluate_faithful(infNorm,

                                                  functionNode,

                                                  lowerBound,

                                                  precision);

    return(functionCallResult);

  }

  if (isConstant(functionNode))

  {

    functionCallResult = pobyso_evaluate_faithful(infNorm,

                                                  functionNode,

                                                  lowerBound,

                                                  precision);

    if (!functionCallResult)

    {

      mpfr_abs(infNorm, infNorm, MPFR_RNDN);

    }

    return(functionCallResult);

  }

  uncertifiedInfnorm(infNorm,

                      functionNode,

                      lowerBound,

                      upperBound,

                      POBYSO_DEFAULT_POINTS,

                      precision);

  return(0);

} /* End pobyso_dirty_infnorm. */

int

pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,

                          node *nodeToEvaluate,

                          mpfr_t argument,

                          mpfr_prec_t precision)

{

  /* Check input arguments. */

  if (nodeToEvaluate == NULL)

  {

    pobyso_error_message("pobyso_evaluate_faithful",

                        "NULL_POINTER_ARGUMENT",

                        "nodeToEvaluate is a NULL pointer");

    return(1);

  }

  evaluateFaithful(faithfulEvaluation,

                   nodeToEvaluate,

                   argument,

                   precision);

  return(0);

} /* End pobyso_evaluate_faithfull. */

chain*

pobyso_find_zeros(node *function,

                  mpfr_t *lower_bound,

                  mpfr_t *upper_bound)

{

  mp_prec_t currentPrecision;

  mpfr_t currentDiameter;

  rangetype bounds;

  currentPrecision = getToolPrecision();

  mpfr_init2(currentDiameter, currentPrecision);

  bounds.a = lower_bound;

  bounds.b = upper_bound;

  if (bounds.a == NULL || bounds.b == NULL)

  {

    pobyso_error_message("pobyso_find_zeros",

                        "MEMORY_ALLOCATION_ERROR",

                        "Could not allocate one of the bounds");

    return(NULL);

  }

  return(findZerosFunction(function,

                            bounds,

                            currentPrecision,

                            currentDiameter));

} /* End pobyso_find_zeros. */

void

pobyso_free_chain_of_nodes(chain *theChainOfNodes)

{

  node *currentNode           = NULL;

  chain *currentChainElement  = NULL;

  chain *nextChainElement     = NULL;

  nextChainElement = theChainOfNodes;

  while(nextChainElement != NULL)

  {

    currentChainElement = nextChainElement;

    currentNode = (node*)(currentChainElement->value);

    nextChainElement = nextChainElement->next;

    free_memory(currentNode);

    free((void*)(currentChainElement));

  }

} /* End pobyso_free_chain_of_nodes. */

void

pobyso_free_range(rangetype range)

{

  mpfr_clear(*(range.a));

  mpfr_clear(*(range.b));

  free(range.a);

  free(range.b);

} /* End pobyso_free_range. */

node*

pobyso_fp_minimax_canonical_monomials_base(node *function,

                                      int degree,

                                      chain *formats,

                                      chain *points,

                                      mpfr_t lowerBound,

                                      mpfr_t upperBound,

                                      int fpFixedArg,

                                      int absRel,

                                      node *constPart,

                                      node *minimax)

{

  return(NULL);

} /* End pobyso_fp_minimax_canonical_monomials_base. */

node*

pobyso_parse_function(char *functionString,

                      char *freeVariableNameString)

{

  if (functionString == NULL || freeVariableNameString == NULL)

  {

    pobyso_error_message("pobyso_parse_function",

                        "NULL_POINTER_ARGUMENT",

                        "One of the arguments is a NULL pointer");

    return(NULL);

  }

  return(parseString(functionString));

} /* End pobyso_parse_function */

node*

pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,

                                                  unsigned int mode,

                                                  mpfr_t lowerBound,

                                                  mpfr_t upperBound,

                                                  mpfr_t eps)

{

  node *weight              = NULL;

  node *bestApproxPolyNode  = NULL;

  node *bestApproxHorner    = NULL;

  node *errorNode           = NULL;

  rangetype degreeRange;

  mpfr_t quality;

  mpfr_t currentError;

  unsigned int degree;

  /* Check the parameters. */

  if (functionNode == NULL)

  {

    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",

                        "NULL_POINTER_ARGUMENT",

                        "functionNode is a NULL pointer");

    return(NULL);

  }

  if (mpfr_cmp(lowerBound, upperBound) >= 0)

  {

    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",

                        "INCOHERENT_INPUT_DATA",

                        "the lower_bound >= upper_bound");

    return(NULL);

  }

  /* Set the weight. */

  if (mode == POBYSO_ABSOLUTE)

  {

    /* Set the weight to 1 for the ABSOLUTE_MODE. */

    weight = makeConstantDouble(1.0);

  }

  else

  {

    if (mode == POBYSO_RELATIVE)

    {

      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",

                          "NOT_IMPLEMENTED",

                          "the search for relative error is not implemented yet");

      return(NULL);

    }

    else

    {

      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",

                          "INCOHERENT_INPUT_DATA",

                          "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");

      return(NULL);

    }

  }

  //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");

  degreeRange = guessDegree(functionNode,

                            weight,

                            lowerBound,

                            upperBound,

                            eps,

                            POBYSO_GUESS_DEGREE_BOUND);

  degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);

  //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");

  fprintf(stderr, "Guessed degree: ");

  mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);

  fprintf(stderr, " - as int: %u\n", degree);

  /* Reduce the degree by 1 in the foolish hope it could work. */

  if (degree > 0) degree--;

  /* Both elements of degreeRange where "inited" within guessDegree. */

  mpfr_clear(*(degreeRange.a));

  mpfr_clear(*(degreeRange.b));

  free(degreeRange.a);

  free(degreeRange.b);

  /* Mimic the default behavior of interactive Sollya. */

  mpfr_init(quality);

  mpfr_set_d(quality, 1e-5, MPFR_RNDN);

  mpfr_init2(currentError, getToolPrecision());

  mpfr_set_inf(currentError,1);

  /* Try to refine the initial guess: loop with increasing degrees until

   * we find a satisfactory one. */

  while(mpfr_cmp(currentError, eps) > 0)

  {

    /* Get rid of the previous polynomial, if any. */

    if (bestApproxPolyNode != NULL)

    {

      free_memory(bestApproxPolyNode);

    }

    fprintf(stderr, "Degree: %u\n", degree);

    fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");

    /* Try to find a polynomial with the guessed degree. */

    bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,

                                                            weight,

                                                            degree,

                                                            lowerBound,

                                                            upperBound,

                                                            quality);

    if (bestApproxPolyNode == NULL)

    {

      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",

                          "INTERNAL_ERROR",

                          "could not compute the bestApproxPolyNode");

      mpfr_clear(currentError);

      mpfr_clear(quality);

      return(NULL);

    }

    setDisplayMode(DISPLAY_MODE_DECIMAL);

    fprintTree(stderr, bestApproxPolyNode);

    fprintf(stderr, "\n\n");

    errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));

    /* Check the error with the computed polynomial. */

    uncertifiedInfnorm(currentError,

                        errorNode,

                        lowerBound,

                        upperBound,

                        POBYSO_INF_NORM_NUM_POINTS,

                        getToolPrecision());

    fprintf(stderr, "Inf norm: ");

    mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);

    fprintf(stderr, "\n\n");

    /* Free the errorNode but not the bestApproxPolyNode (we need it if

     * we exit the loop at the next iteration). */

    free_memory(errorNode);

    degree++;

  }

  /* Use an intermediate variable, since horner() creates a new node

   * and does not reorder the argument "in place". This allows for the memory

   * reclaim of bestApproxHorner.

   */

  bestApproxHorner = horner(bestApproxPolyNode);

  free_memory(bestApproxPolyNode);

  mpfr_clear(currentError);

  mpfr_clear(quality);

  free_memory(weight);

  return(bestApproxHorner);

} /* End pobyso_remez_approx_canonical_monomials_base_for_error */

node*

pobyso_remez_canonical_monomials_base(node *function,

                                     node *weight,

                                     unsigned int degree,

                                     mpfr_t lowerBound,

                                     mpfr_t upperBound,

                                     mpfr_t quality)

{

  node  *bestApproxPoly = NULL;

  chain *monomials      = NULL;

  chain *curMonomial    = NULL;

  mpfr_t satisfying_error;

  mpfr_t target_error;

  /* Argument checking */

  /* Function tree. */

  if (function == NULL)

  {

    pobyso_error_message("pobyso_remez_canonical_monomials_base",

                        "NULL_POINTER_ARGUMENT",

                        "the \"function\" argument is a NULL pointer");

    return(NULL);

  }

  if (weight == NULL)

  {

    pobyso_error_message("pobyso_remez_canonical_monomials_base",

                        "NULL_POINTER_ARGUMENT",

                        "the \"weight\" argument is a NULL pointer");

    return(NULL);

  }

  /* Check the bounds. */

  if (mpfr_cmp(lowerBound, upperBound) >= 0)

  {

    pobyso_error_message("pobyso_remez_canonical_monomials_base",

                        "INCOHERENT_IMPUT_DATA",

                        "the lower_bound >= upper_bound");

    return(NULL);

  }

  /* The quality must be a non null positive number. */

  if (mpfr_sgn(quality) <= 0)

  {

    pobyso_error_message("pobyso_remez_canonical_monomials_base",

                        "INCOHERENT_INPUT_DATA",

                        "the quality <= 0");

  }

  /* End argument checking. */

  /* Create the monomials nodes chains. */

  monomials = pobyso_create_canonical_monomials_base(degree);

  fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));

  if (monomials == NULL || (lengthChain(monomials) != degree + 1))

  {

    pobyso_error_message("pobyso_remez_canonical_monomials_base",

                        "CHAIN_CREATION_ERROR",

                        "could not create the monomials chain");

    return(NULL);

  }

  curMonomial = monomials;

  while (curMonomial != NULL)

  {

    fprintf(stderr, "monomial tree: ");

    //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);

    fprintTree(stderr, (node*)(curMonomial->value));

    fprintf(stderr, "\n");

    curMonomial = curMonomial->next;

  }

  /* Deal with NULL weight. */

  if (weight == NULL)

  {

    weight = makeConstantDouble(1.0);

  }

  /* Compute the best polynomial with the required quality.

     The behavior is as if satisfying error and target_error had

     not been used.*/

  mpfr_init(satisfying_error);

  mpfr_init(target_error);

  mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);

  mpfr_set_inf(target_error, 1);

  fprintf(stderr, "satisfying_error: ");

  mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);

  fprintf(stderr, ".\n");

  fprintf(stderr, "target_error: ");

  mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);

  fprintf(stderr, ".\n");

  fprintf(stderr,

          "current precision: %li\n", getToolPrecision());

  /* Call the Sollya function. */

  bestApproxPoly = remez(function,

                          weight,

                          monomials,

                          lowerBound,

                          upperBound,

                          quality,

                          satisfying_error,

                          target_error,

                          getToolPrecision());

  mpfr_clear(satisfying_error);

  mpfr_clear(target_error);

  pobyso_free_chain_of_nodes(monomials);

  return(bestApproxPoly);

} /* End pobyso_remez_canonical_monomials_base. */

#endif

void

pobyso_error_message(char *functionName, char *messageName, char* message)

{

  fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);

} /* End pobyso_error_message */