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

CC           = gcc

CFLAGS       = -Wall -g  

LDLIB        = #-Wl,--verbose

OBJS         = 

INCS         = 

SOURCES      = $(INCS)

TEST_DRIVER  = # CUnit test driver binary program

ALL_TARGETS = 

SCRATCH_PRG := $(patsubst %.c,%,$(wildcard *.c))

SCRATCH_ALL := $(SCRATCH_PRG)

LAST_ARCH_FILE = lastArch.txt

ARCH := $(shell uname -m)

DUMMY := $(shell touch $(LAST_ARCH_FILE))

LAST_ARCH := $(shell cat $(LAST_ARCH_FILE))

ifneq ($(LAST_ARCH), $(ARCH))

  DUMMY := $(shell echo $(ARCH) > $(LAST_ARCH_FILE))

  DUMMY := $(shell $(MAKE) clean) 

endif

ifeq (x86_64, $(ARCH))

  CFLAGS := -m64 $(CFLAGS)

  # LDLIB :=

endif

ifeq (ppc64, $(ARCH))

  CFLAGS := -m64 $(CFLAGS)

  # LDLIB :=

endif

ifeq (i686, $(ARCH))

  # CFLAGS := $(CFLAGS)

  # LDLIB :=

endif

# -------------------------------------------------------------------------

default: all

# Remove all implicit (suffix) rules. No need to remove pattern rules.

# Just redefine them.

.SUFFIXES:

# Keep the OBJS we want to link with user code.

.PRECIOUS: $(OBJS)

all: $(ALL_TARGETS)

test: $(TEST_DRIVER)

% : %.o $(OBJS)

	$(CC) $(CFLAGS) -o $@ $< $(OBJS) $(LDLIB)

%.o : %.c $(INCS)

	$(CC) $(CFLAGS) -c -o $@ $<

$(TEST_DRIVER) : 

	make -C test

	cp test/$(TEST_DRIVER) .

# -------------------------------------------------------------------------

.PHONY: clean scratch

clean:

	@rm -f $(ALL_TARGETS) $(TEST_DRIVER)

	@rm -f *.o a.out $(OBJS)

	@rm -f *~ *% #*#

scratch: clean

	rm -f $(SCRATCH_ALL)

# Change the documentation configuration and path accordingly.

documentation: pobyso.doxyconf $(SOURCES)

	doxygen pobyso.doxyconf

/** @file pobyso.h

 * Integration of Sollya to C programs

 * @author S.T.

 * @date 2011-10-11

 * Note: pobyso stands for POwered BY SOllya.

 *                         --      -- --

 */

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

/*

 * Add below all the headers needed to get this header work.

 */

/* <stdio.h> is needed *before* <mpfr.h> for all MPFR input/output functions

 * prototypes be defined. */

#include <stdio.h>

#include <sollya.h>

#include <mpfr.h>

#ifndef POBYSO_h

#define POBYSO_ABSOLUTE 1

#define POBYSO_RELATIVE 2

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

#define POBYSO_DEFAULT_POINTS 501

#define POBYSO_INF_NORM_NUM_POINTS (POBYSO_DEFAULT_POINTS)

#define POBYSO_GUESS_DEGREE_BOUND 1024

/**

 * Print object(s) to stdout.

 * A very thin wrapper around the lib_sollya_v_autoprint() function.

 * Should be called with NULL as the final argument.

 */

void

pobyso_autoprint(sollya_obj_t soObj, ...);

/**

 * Print object(s) to stdout: the va_list companion function.

 * A very thin wrapper around the lib_sollya_v_autoprint() function.

 * The last argument in the va_list should be NULL.

 */

void

pobyso_autoprint_v(sollya_obj_t soObj, va_list va);

/**

 * Parse a string to create a Sollya object.

 * A very thin wrapper around the sollya_lib_parse_string() function.

 */

sollya_obj_t

pobyso_parse_string(const char* expression);

/**

 * Set the verbosity mode off (level 0).

 * If currentVerbosity != NULL, currentVerbosity is set

 * to the current verbosity level. It may be used to reset the

 * verbosity level later.

 */

void

pobyso_set_verbosity_off(sollya_obj_t* currentVerbosityLevel);

/**

 * Set the verbosity level to newVerbosityLevel.

 * @param newVerbosityLevel must be a Sollay object corresponding to an

 *        integer constant.

 */

void

pobyso_set_verbosity_to(sollya_obj_t newVerbosityLevel);

#if 0

/**

 * Create the canonical (non sparse) base of monomials for a given degree.

 */

chain*

pobyso_create_canonical_monomials_base(const unsigned int degree);

/**

 * Create a chain from an array of int.

 */

sollya_int_list

pobyso_create_int_list_from_int_Array(int* intArray,

                                    const unsigned int arrayLength);

/**

 * Create a chain from an array of int.

 */

sollya_int_list_t

pobyso_create_int_list_from_unsigned_int_array(unsigned int* intArray,

                                            const unsigned int arrayLength);

/**

 * Differentiation of a function.

 * A slim wrapper around the Sollya differentiate function.

 * @param functionNode - the Sollya node to differentiate;

 * @return a node representing the function differentiated or NULL, if

 *         something goes wrong.

 */

sollya_obj_t

pobyso_diff(sollya_obj_t function);

/**

 * A match to the Sollya dirtyinfnorm.

 * A slim wrapper around the Sollya function.

 * @param infnorm - out parameter to return the result, must be "inited"

 *                  and "cleared" by the caller;

 * @param functionNode - the Sollya node to compute the infinite norm of;

 * @param lowerBound - the lower bound of the interval;

 * @param upperBound - the upper bound of the interval;

 * @param precision  - the internal precision Sollya must use.

 * @return 0 if everything is OK, != 0 if something goes wrong.

 */

int

pobyso_dirty_infnorm(mpfr_t infNorm,

                      node *functionNode,

                      mpfr_t lowerBound,

                      mpfr_t upperBound,

                      mp_prec_t precision);

/**

 * Faithful evaluation of an expression.

 *

 * @param faitufulEvaluation - holds the result, must be "inited" by the

 *                             caller;

 */

int

pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,

                          node *nodeToEvaluate,

                          mpfr_t argument,

                          mpfr_prec_t precision);

/**

 * Find the zeros of a function on a given interval.

 */

chain*

pobyso_find_zeros(node *function,

                  mpfr_t *lowerBound,

                  mpfr_t *upperBound);

/**

 * Free a chain of node.

 * All elements of the chain have to be nodes.

 */

void

pobyso_free_chain_of_nodes(chain *theChain);

/**

 * Free a range.

 * It involves clearing the mpfr_t elements and deallocating the

 * pointers.

 */

void

pobyso_free_range(rangetype range);

/**

 * Computes a good polynomial approximation with fixed-point or floating-point

 * coefficients.

 */

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);

/**

 * Parses a string to build the node representing the function.

 * In fact, does nothing for the moment: the string must be a correct function

 * definition. No error correction takes place here.

 */

node*

pobyso_parse_function(char *functionString,

                      char *freeVariableNameString);

/** Compute a polynomial approximation in the canonical monomials basis for

 *  a function, for a given precision. The returned polynomial has the minimal

 *  degree to achieve the required precision.

 */

node*

pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,

                                                      unsigned int mode,

                                                      mpfr_t lowerBound,

                                                      mpfr_t upperBound,

                                                      mpfr_t eps);

/**

 * Computes a the remez approximation of a function.

 * @param function   - the node holding the function to approximate;

 * @param weight     - the node holding the weight function, can be NULL. In

 *                     this case a default weight of "1" will be provided and

 *                     the approximation is related is with respect to the

 *                     absolute error.

 * @param degree     - the degree of the approximation polynomial;

 * @param lowerBound - the lower bound of the approximation interval;

 * @param upperBound - the upper bound of the approximation interval;

 * @param quality    - quality = (eps - eps*) / eps*; the search stop when the required

 *                     quality is achieved.

 *

 * @return a node holding the approximation polynomial in Horner form.

 */

node*

pobyso_remez_canonical_monomials_base(node *function,

                                       node *weight,

                                       unsigned int degree,

                                       mpfr_t lowerBound,

                                       mpfr_t upperBound,

                                       mpfr_t quality);

#endif

#define POBYSO_h  

#endif

/** @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 */

CC           = gcc

CFLAGS       = -Wall -g -I /warehouse/storres/root/include

LDLIB        = #-Wl,--verbose

OBJS         = pobyso.o

INCS         = pobyso.h 

SOURCES      = $(INCS) pobyso.c

TEST_DRIVER  = # CUnit test driver binary program

TEST_DIR     = test

TEST_SOURCES = $(TEST_DIR)/$(TEST_DRIVER).c \

               # Add the suites definition sources

ALL_TARGETS = pobyso.o

SCRATCH_PRG := $(patsubst %.c,%,$(wildcard *.c))

SCRATCH_ALL := $(SCRATCH_PRG)

LAST_ARCH_FILE = lastArch.txt

ARCH := $(shell uname -m)

DUMMY := $(shell touch $(LAST_ARCH_FILE))

LAST_ARCH := $(shell cat $(LAST_ARCH_FILE))

ifneq ($(LAST_ARCH), $(ARCH))

  DUMMY := $(shell echo $(ARCH) > $(LAST_ARCH_FILE))

  DUMMY := $(shell $(MAKE) clean) 

endif

ifeq (x86_64, $(ARCH))

  CFLAGS := -m64 $(CFLAGS)

  # LDLIB :=

endif

ifeq (ppc64, $(ARCH))

  CFLAGS := -m64 $(CFLAGS)

  # LDLIB :=

endif

ifeq (i686, $(ARCH))

  # CFLAGS := $(CFLAGS)

  # LDLIB :=

endif

# -------------------------------------------------------------------------

default: all

# Remove all implicit (suffix) rules. No need to remove pattern rules.

# Just redefine them.

.SUFFIXES:

# Keep the OBJS we want to link with user code.

.PRECIOUS: $(OBJS)

all: $(ALL_TARGETS)

test: $(TEST_DRIVER)

% : %.o $(OBJS)

	$(CC) $(CFLAGS) -o $@ $< $(OBJS) $(LDLIB)

%.o : %.c $(INCS)

	$(CC) $(CFLAGS) -c -o $@ $<

$(TEST_DRIVER) : 

	make -C test

	cp test/$(TEST_DRIVER) .

# -------------------------------------------------------------------------

.PHONY: clean scratch

clean:

	@rm -f $(ALL_TARGETS) $(TEST_DRIVER)

	@rm -f *.o a.out $(OBJS)

	@rm -f *~ *% #*#

	@make -C $(TEST_DIR) clean

scratch: clean

	rm -f $(SCRATCH_ALL)

# Change the documentation configuration and path accordingly.

documentation: pobyso.doxyconf $(SOURCES)

	doxygen pobyso.doxyconf