Statistiques
| Révision :

root / pobysoC-4.0 / src / pobyso.c @ 133

Historique | Voir | Annoter | Télécharger (25,12 ko)

1 26 storres
/** @file pobyso.c
2 26 storres
 * Name & purpose
3 26 storres
 *
4 26 storres
 * @author S.T.
5 26 storres
 * @date 2011-10-12
6 26 storres
 *
7 27 storres
 *
8 26 storres
 */
9 50 storres
10 26 storres
/******************************************************************************/
11 26 storres
/* Headers, applying the "particular to general" convention.*/
12 26 storres
13 26 storres
#include "pobyso.h"
14 26 storres
15 26 storres
/* includes of local headers */
16 26 storres
17 26 storres
/* includes of project headers */
18 26 storres
19 26 storres
/* includes of system headers */
20 128 storres
#include <string.h>
21 128 storres
#include <stdlib.h>
22 128 storres
#include <stdio.h>
23 26 storres
24 26 storres
/* Other declarations */
25 26 storres
26 26 storres
/* Internal prototypes */
27 26 storres
void
28 26 storres
pobyso_error_message(char *functionName, char *messageName, char* message);
29 26 storres
/* Types, constants and macros definitions */
30 26 storres
31 26 storres
/* Global variables */
32 26 storres
33 26 storres
/* Functions */
34 127 storres
35 127 storres
/* @see pobyso.h#pobyso_autprint */
36 26 storres
void
37 26 storres
pobyso_autoprint(sollya_obj_t solObj, ...)
38 26 storres
{
39 26 storres
  va_list va;
40 26 storres
  va_start(va, solObj);
41 26 storres
  sollya_lib_v_autoprint(solObj, va);
42 26 storres
  va_end(va);
43 26 storres
} /* End pobyso_autoprint. */
44 26 storres
45 133 storres
/* @see pobyso.h#pobyso_is_constant_expression*/
46 133 storres
int
47 133 storres
pobyso_is_constant_expression(sollya_obj_t obj_to_test)
48 133 storres
{
49 133 storres
  mpfr_t dummy;
50 133 storres
  int test;
51 133 storres
  sollya_obj_t verbosity_level      = NULL;
52 133 storres
  /* Test argument. */
53 133 storres
  if (obj_to_test == NULL)
54 133 storres
  {
55 133 storres
    pobyso_error_message("pobyso_parse_string",
56 133 storres
                        "NULL_POINTER_ARGUMENT",
57 133 storres
                        "The expression is a NULL pointer");
58 133 storres
    return 0;
59 133 storres
  }
60 133 storres
  verbosity_level = pobyso_set_verbosity_off();
61 128 storres
62 133 storres
  if (! sollya_lib_obj_is_function(obj_to_test))
63 133 storres
  {
64 133 storres
    pobyso_set_verbosity_to(verbosity_level);
65 133 storres
    sollya_lib_clear_obj(verbosity_level);
66 133 storres
    return 0;
67 133 storres
  }
68 133 storres
  mpfr_init2(dummy,64);
69 133 storres
  /* Call to previous Sollya function resets verbosity. */
70 133 storres
  verbosity_level = pobyso_set_verbosity_off();
71 133 storres
  test = sollya_lib_get_constant(dummy, obj_to_test);
72 133 storres
  mpfr_clear(dummy);
73 133 storres
  pobyso_set_verbosity_to(verbosity_level);
74 133 storres
  sollya_lib_clear_obj(verbosity_level);
75 133 storres
  return test;
76 133 storres
} /* End pobyso_is_constant_expression. */
77 133 storres
78 133 storres
/** @see pobyso.h#pobyso_new_monomial. */
79 133 storres
pobyso_func_exp_t
80 133 storres
pobyso_new_monomial(pobyso_func_exp_t coefficientSa, long degree)
81 133 storres
{
82 133 storres
  sollya_obj_t degreeSa   = NULL;
83 133 storres
  sollya_obj_t varToPowSa = NULL;
84 133 storres
  sollya_obj_t monomialSa = NULL;
85 133 storres
  if (coefficientSa == NULL)
86 133 storres
  {
87 133 storres
    pobyso_error_message("pobyso_parse_string",
88 133 storres
                        "NULL_POINTER_ARGUMENT",
89 133 storres
                        "The expression is a NULL pointer");
90 133 storres
    return sollya_lib_error();
91 133 storres
  }
92 133 storres
  if (! pobyso_is_constant_expression(coefficientSa))
93 133 storres
  {
94 133 storres
    return sollya_lib_error();
95 133 storres
  }
96 133 storres
  if (degree < 0)
97 133 storres
  {
98 133 storres
    return sollya_lib_error();
99 133 storres
  }
100 133 storres
  degreeSa    = sollya_lib_constant_from_int64(degree);
101 133 storres
  varToPowSa  = sollya_lib_build_function_pow(sollya_lib_free_variable(),
102 133 storres
                                              degreeSa);
103 133 storres
  monomialSa = sollya_lib_mul(coefficientSa,varToPowSa);
104 133 storres
  sollya_lib_clear_obj(varToPowSa);
105 133 storres
  return monomialSa;
106 133 storres
} /* End pobyso_new_monomial. */
107 133 storres
108 127 storres
/* @see pobyso.h#pobyso_parse_string*/
109 130 storres
pobyso_func_exp_t
110 26 storres
pobyso_parse_string(const char* expression)
111 26 storres
{
112 128 storres
  int expressionLength, i;
113 127 storres
  char *expressionWithSemiCo;
114 127 storres
  sollya_obj_t expressionSa;
115 128 storres
116 127 storres
  /* Arguments check. */
117 26 storres
  if (expression == NULL)
118 26 storres
  {
119 26 storres
    pobyso_error_message("pobyso_parse_string",
120 26 storres
                        "NULL_POINTER_ARGUMENT",
121 26 storres
                        "The expression is a NULL pointer");
122 127 storres
    return(sollya_lib_error());
123 26 storres
  }
124 127 storres
  expressionLength = strlen(expression);
125 127 storres
  if (expressionLength == 0)
126 127 storres
  {
127 127 storres
    pobyso_error_message("pobyso_parse_string",
128 127 storres
                        "EMPTY_STRING_ARGUMENT",
129 127 storres
                        "The expression is an empty string");
130 132 storres
    return sollya_lib_error();
131 127 storres
  }
132 128 storres
  /* Search from the last char of the expression until, whichever happens
133 128 storres
   * first:
134 128 storres
   * a ";" is found;
135 128 storres
   * a char  != ';' is found the the ";" is appended.
136 128 storres
   * If the head of the string is reached before any of these two events happens
137 128 storres
   * return an error.
138 128 storres
   */
139 128 storres
  for (i = expressionLength - 1 ; i >= 0 ; i--)
140 127 storres
  {
141 128 storres
    if (expression[i] == ';') /* Nothing special to do:
142 128 storres
                                 try to parse the string*/
143 127 storres
    {
144 128 storres
      expressionSa = sollya_lib_parse_string(expression);
145 132 storres
      return expressionSa;
146 127 storres
    }
147 128 storres
    else
148 128 storres
    {
149 128 storres
      if (expression[i] == ' ' || expression[i] == '\t') /* A blank,
150 128 storres
                                                           just continue. */
151 128 storres
      {
152 128 storres
        continue;
153 128 storres
      }
154 128 storres
      else /* a character != ';' and from a blank: create the ';'ed string. */
155 128 storres
      {
156 128 storres
        /* Create a new string for the expression, add the ";" and
157 128 storres
         * and call sollya_lib_parse_string. */
158 128 storres
        expressionWithSemiCo = calloc(i + 3, sizeof(char));
159 128 storres
        if (expressionWithSemiCo == NULL)
160 128 storres
        {
161 128 storres
          pobyso_error_message("pobyso_parse_string",
162 128 storres
                                "MEMORY_ALLOCATION_ERROR",
163 128 storres
                                "Could not allocate the expression string");
164 132 storres
          return sollya_lib_error();
165 128 storres
        }
166 132 storres
        strncpy(expressionWithSemiCo, expression, i+1);
167 128 storres
        expressionWithSemiCo[i + 1] = ';';
168 128 storres
        expressionWithSemiCo[i + 2] = '\0';
169 128 storres
        expressionSa = sollya_lib_parse_string(expressionWithSemiCo);
170 128 storres
        free(expressionWithSemiCo);
171 132 storres
        return expressionSa;
172 128 storres
      } /* End character != ';' and from a blank. */
173 128 storres
    /* Create a new string for the expression, add the ";" and
174 128 storres
     * and call sollya_lib_parse_string. */
175 128 storres
    } /* End else. */
176 128 storres
  } /* End for. */
177 128 storres
  /* We get here, it is because we did not find any char == anything different
178 128 storres
   * from ' ' or '\t': the string is empty.
179 128 storres
   */
180 128 storres
  pobyso_error_message("pobyso_parse_string",
181 128 storres
                       "ONLY_BLANK_ARGUMENT",
182 128 storres
                        "The expression is only made of blanks");
183 128 storres
  return(sollya_lib_error());
184 128 storres
185 26 storres
} /* pobyso_parse_string */
186 26 storres
187 132 storres
pobyso_func_exp_t
188 132 storres
pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function,
189 132 storres
                                      long int degree,
190 132 storres
                                      pobyso_range_t interval,
191 132 storres
                                      pobyso_func_exp_t weight,
192 132 storres
                                      double quality,
193 132 storres
                                      pobyso_range_t bounds)
194 132 storres
{
195 132 storres
  sollya_obj_t  degreeSa = NULL;
196 132 storres
  sollya_obj_t qualitySa = NULL;
197 132 storres
198 132 storres
  degreeSa = sollya_lib_constant_from_int(degree);
199 132 storres
  qualitySa = sollya_lib_constant_from_double(quality);
200 132 storres
201 132 storres
  return NULL;
202 132 storres
} /* End pobyso_remez_canonical_monomials_base. */
203 132 storres
204 132 storres
int
205 132 storres
pobyso_set_canonical_on()
206 132 storres
{
207 132 storres
  pobyso_on_off_t currentCanonicalModeSo;
208 132 storres
  int currentCanonicalMode;
209 132 storres
  pobyso_on_off_t on;
210 132 storres
211 132 storres
  on = sollya_lib_on();
212 132 storres
  currentCanonicalModeSo = sollya_lib_get_canonical();
213 132 storres
  if (sollya_lib_is_on(currentCanonicalModeSo))
214 132 storres
  {
215 132 storres
    currentCanonicalMode = POBYSO_ON;
216 132 storres
  }
217 132 storres
  else
218 132 storres
  {
219 132 storres
    currentCanonicalMode = POBYSO_OFF;
220 132 storres
  }
221 132 storres
  sollya_lib_set_canonical(on);
222 132 storres
  sollya_lib_clear_obj(on);
223 132 storres
  sollya_lib_clear_obj(currentCanonicalModeSo);
224 132 storres
  return currentCanonicalMode;
225 132 storres
} /* End pobyso_set_canonical_on. */
226 132 storres
227 132 storres
pobyso_sollya_verbosity_t
228 128 storres
pobyso_set_verbosity_off()
229 26 storres
{
230 26 storres
  sollya_obj_t verbosityLevelZero;
231 128 storres
  sollya_obj_t currentVerbosityLevel = NULL;
232 128 storres
233 128 storres
  currentVerbosityLevel = sollya_lib_get_verbosity();
234 26 storres
  verbosityLevelZero = sollya_lib_constant_from_int(0);
235 26 storres
  sollya_lib_set_verbosity(verbosityLevelZero);
236 26 storres
  sollya_lib_clear_obj(verbosityLevelZero);
237 128 storres
  return(currentVerbosityLevel);
238 26 storres
} /* End of pobyso_set_verbosity_off. */
239 26 storres
240 132 storres
pobyso_sollya_verbosity_t
241 132 storres
pobyso_set_verbosity_to(pobyso_sollya_verbosity_t newVerbosityLevel)
242 26 storres
{
243 131 storres
  sollya_obj_t currentVerbosityLevel = NULL;
244 26 storres
  if (newVerbosityLevel == NULL)
245 26 storres
  {
246 26 storres
    pobyso_error_message("pobyso_set_verbosity_to",
247 26 storres
                        "NULL_POINTER_ARGUMENT",
248 132 storres
                        "The new verbosity level is a null pointer");
249 128 storres
    return(sollya_lib_error());
250 26 storres
  }
251 131 storres
  currentVerbosityLevel = sollya_lib_get_verbosity();
252 26 storres
  sollya_lib_set_verbosity(newVerbosityLevel);
253 131 storres
  return(currentVerbosityLevel);
254 26 storres
} /* End of pobyso_set_verbosity_to. */
255 26 storres
256 132 storres
pobyso_sollya_verbosity_t
257 132 storres
pobyso_set_verbosity_to_int(int intVerbosityLevel)
258 132 storres
{
259 132 storres
  sollya_obj_t currentVerbosityLevel = NULL;
260 132 storres
  sollya_obj_t newVerbosityLevel     = NULL;
261 132 storres
  if (intVerbosityLevel < 0)
262 132 storres
  {
263 132 storres
    pobyso_error_message("pobyso_set_verbosity_to",
264 132 storres
                        "INVALID_VALUE",
265 132 storres
                        "The new verbosity level is a negative number");
266 132 storres
    return(sollya_lib_error());
267 132 storres
  }
268 132 storres
  newVerbosityLevel = sollya_lib_constant_from_int(intVerbosityLevel);
269 132 storres
  currentVerbosityLevel = sollya_lib_get_verbosity();
270 132 storres
  sollya_lib_set_verbosity(newVerbosityLevel);
271 132 storres
  return(currentVerbosityLevel);
272 132 storres
} /* End of pobyso_set_verbosity_to_int. */
273 132 storres
274 26 storres
/* Attic from the sollya_lib < 4. */
275 26 storres
#if 0
276 26 storres
chain*
277 26 storres
pobyso_create_canonical_monomials_base(const unsigned int degree)
278 26 storres
{
279 26 storres
  int i    = 0;
280 26 storres
  chain *monomials  = NULL;
281 26 storres
  node  *monomial   = NULL;
282 26 storres

283 26 storres
  for(i = degree ; i >= 0  ; i--)
284 26 storres
  {
285 26 storres
     monomial   = makePow(makeVariable(), makeConstantDouble((double)i));
286 26 storres
     monomials  = addElement(monomials, monomial);
287 26 storres
     fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
288 26 storres
  }
289 26 storres
  if (monomials == NULL)
290 26 storres
  {
291 26 storres
    pobyso_error_message("pobyso_create_canonical_monomial_base",
292 26 storres
                        "CHAIN_CREATION_ERROR",
293 26 storres
                        "Could not create the monomials chain");
294 26 storres
    return(NULL);
295 26 storres
  }
296 26 storres
  return(monomials);
297 26 storres
} /* End pobyso_create_canonical_monomials_base. */
298 26 storres

299 26 storres
chain*
300 26 storres
pobyso_create_chain_from_int_array(int* intArray,
301 26 storres
                                    const unsigned int arrayLength)
302 26 storres
{
303 26 storres
  int i = 0;
304 26 storres
  chain *newChain = NULL;
305 26 storres
  int *currentInt;
306 26 storres

307 26 storres
  if (arrayLength == 0) return(NULL);
308 26 storres
  if (intArray == NULL)
309 26 storres
  {
310 26 storres
    pobyso_error_message("pobyso_create_chain_from_int_array",
311 26 storres
                        "NULL_POINTER_ARGUMENT",
312 26 storres
                        "The array is a NULL pointer");
313 26 storres
    return(NULL);
314 26 storres
  }
315 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
316 26 storres
  {
317 26 storres
    currentInt = malloc(sizeof(int));
318 26 storres
    if (currentInt == NULL)
319 26 storres
    {
320 26 storres
      pobyso_error_message("pobyso_create_chain_from_int_array",
321 26 storres
                          "MEMORY_ALLOCATION_ERROR",
322 26 storres
                          "Could not allocate one of the integers");
323 26 storres
      freeChain(newChain, free);
324 26 storres
      return(NULL);
325 26 storres
    }
326 26 storres
    *currentInt = intArray[i];
327 26 storres
    newChain = addElement(newChain, currentInt);
328 26 storres
  }
329 26 storres
  return(newChain);
330 26 storres
} // End pobyso_create_chain_from_int_array. */
331 26 storres

332 26 storres
chain*
333 26 storres
pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
334 26 storres
                                        const unsigned int arrayLength)
335 26 storres
{
336 26 storres
  int i = 0;
337 26 storres
  chain *newChain = NULL;
338 26 storres
  unsigned int *currentInt;
339 26 storres

340 26 storres
  /* Argument checking. */
341 26 storres
  if (arrayLength == 0) return(NULL);
342 26 storres
  if (intArray == NULL)
343 26 storres
  {
344 26 storres
    pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
345 26 storres
                        "NULL_POINTER_ARGUMENT",
346 26 storres
                        "The array is a NULL pointer");
347 26 storres
    return(NULL);
348 26 storres
  }
349 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
350 26 storres
  {
351 26 storres
    currentInt = malloc(sizeof(unsigned int));
352 26 storres
    if (currentInt == NULL)
353 26 storres
    {
354 26 storres
      pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
355 26 storres
                          "MEMORY_ALLOCATION_ERROR",
356 26 storres
                          "Could not allocate one of the integers");
357 26 storres
      freeChain(newChain, free);
358 26 storres
      return(NULL);
359 26 storres
    }
360 26 storres
    *currentInt = intArray[i];
361 26 storres
    newChain = addElement(newChain, currentInt);
362 26 storres
  }
363 26 storres
  return(newChain);
364 26 storres
} // End pobyso_create_chain_from_unsigned_int_array. */
365 26 storres

366 26 storres
node*
367 26 storres
pobyso_differentiate(node *functionNode)
368 26 storres
{
369 26 storres
  /* Argument checking. */
370 26 storres
  node *differentialNode;
371 26 storres
  if (functionNode == NULL)
372 26 storres
  {
373 26 storres
    pobyso_error_message("pobyso_differentiate",
374 26 storres
                        "NULL_POINTER_ARGUMENT",
375 26 storres
                        "The function to differentiate is a NULL pointer");
376 26 storres
    return(NULL);
377 26 storres
  }
378 26 storres
  differentialNode = differentiate(functionNode);
379 26 storres
  if (differentialNode == NULL)
380 26 storres
  {
381 26 storres
    pobyso_error_message("pobyso_differentiate",
382 26 storres
                        "INTERNAL ERROR",
383 26 storres
                        "Sollya could not differentiate the function");
384 26 storres
  }
385 26 storres
  return(differentialNode);
386 26 storres
} // End pobyso_differentiate
387 26 storres

388 26 storres

389 26 storres
int
390 26 storres
pobyso_dirty_infnorm(mpfr_t infNorm,
391 26 storres
                      node *functionNode,
392 26 storres
                      mpfr_t lowerBound,
393 26 storres
                      mpfr_t upperBound,
394 26 storres
                      mp_prec_t precision)
395 26 storres
{
396 26 storres
  int functionCallResult;
397 26 storres
  /* Arguments checking. */
398 26 storres
  if (functionNode == NULL)
399 26 storres
  {
400 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
401 26 storres
                        "NULL_POINTER_ARGUMENT",
402 26 storres
                        "The function to compute is a NULL pointer");
403 26 storres
    return(1);
404 26 storres
  }
405 26 storres
  if (mpfr_cmp(lowerBound, upperBound) > 0)
406 26 storres
  {
407 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
408 26 storres
                        "INCOHERENT_INPUT_DATA",
409 26 storres
                        "The lower bond is greater than the upper bound");
410 26 storres
    return(1);
411 26 storres
  }
412 26 storres
  /* Particular cases. */
413 26 storres
  if (mpfr_cmp(lowerBound, upperBound) == 0)
414 26 storres
  {
415 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
416 26 storres
                                                  functionNode,
417 26 storres
                                                  lowerBound,
418 26 storres
                                                  precision);
419 26 storres
    return(functionCallResult);
420 26 storres
  }
421 26 storres
  if (isConstant(functionNode))
422 26 storres
  {
423 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
424 26 storres
                                                  functionNode,
425 26 storres
                                                  lowerBound,
426 26 storres
                                                  precision);
427 26 storres
    if (!functionCallResult)
428 26 storres
    {
429 26 storres
      mpfr_abs(infNorm, infNorm, MPFR_RNDN);
430 26 storres
    }
431 26 storres
    return(functionCallResult);
432 26 storres
  }
433 26 storres
  uncertifiedInfnorm(infNorm,
434 26 storres
                      functionNode,
435 26 storres
                      lowerBound,
436 26 storres
                      upperBound,
437 26 storres
                      POBYSO_DEFAULT_POINTS,
438 26 storres
                      precision);
439 26 storres
  return(0);
440 26 storres
} /* End pobyso_dirty_infnorm. */
441 26 storres

442 26 storres
int
443 26 storres
pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
444 26 storres
                          node *nodeToEvaluate,
445 26 storres
                          mpfr_t argument,
446 26 storres
                          mpfr_prec_t precision)
447 26 storres
{
448 26 storres
  /* Check input arguments. */
449 26 storres
  if (nodeToEvaluate == NULL)
450 26 storres
  {
451 26 storres
    pobyso_error_message("pobyso_evaluate_faithful",
452 26 storres
                        "NULL_POINTER_ARGUMENT",
453 26 storres
                        "nodeToEvaluate is a NULL pointer");
454 26 storres
    return(1);
455 26 storres
  }
456 26 storres
  evaluateFaithful(faithfulEvaluation,
457 26 storres
                   nodeToEvaluate,
458 26 storres
                   argument,
459 26 storres
                   precision);
460 26 storres
  return(0);
461 26 storres
} /* End pobyso_evaluate_faithfull. */
462 26 storres

463 26 storres
chain*
464 26 storres
pobyso_find_zeros(node *function,
465 26 storres
                  mpfr_t *lower_bound,
466 26 storres
                  mpfr_t *upper_bound)
467 26 storres
{
468 26 storres
  mp_prec_t currentPrecision;
469 26 storres
  mpfr_t currentDiameter;
470 26 storres
  rangetype bounds;
471 26 storres

472 26 storres
  currentPrecision = getToolPrecision();
473 26 storres
  mpfr_init2(currentDiameter, currentPrecision);
474 26 storres

475 26 storres
  bounds.a = lower_bound;
476 26 storres
  bounds.b = upper_bound;
477 26 storres

478 26 storres
  if (bounds.a == NULL || bounds.b == NULL)
479 26 storres
  {
480 26 storres
    pobyso_error_message("pobyso_find_zeros",
481 26 storres
                        "MEMORY_ALLOCATION_ERROR",
482 26 storres
                        "Could not allocate one of the bounds");
483 26 storres
    return(NULL);
484 26 storres
  }
485 26 storres
  return(findZerosFunction(function,
486 26 storres
                            bounds,
487 26 storres
                            currentPrecision,
488 26 storres
                            currentDiameter));
489 26 storres
} /* End pobyso_find_zeros. */
490 26 storres

491 26 storres
void
492 26 storres
pobyso_free_chain_of_nodes(chain *theChainOfNodes)
493 26 storres
{
494 26 storres
  node *currentNode           = NULL;
495 26 storres
  chain *currentChainElement  = NULL;
496 26 storres
  chain *nextChainElement     = NULL;
497 26 storres

498 26 storres
  nextChainElement = theChainOfNodes;
499 26 storres

500 26 storres
  while(nextChainElement != NULL)
501 26 storres
  {
502 26 storres
    currentChainElement = nextChainElement;
503 26 storres
    currentNode = (node*)(currentChainElement->value);
504 26 storres
    nextChainElement = nextChainElement->next;
505 26 storres
    free_memory(currentNode);
506 26 storres
    free((void*)(currentChainElement));
507 26 storres
  }
508 26 storres
} /* End pobyso_free_chain_of_nodes. */
509 26 storres

510 26 storres
void
511 26 storres
pobyso_free_range(rangetype range)
512 26 storres
{
513 26 storres

514 26 storres
  mpfr_clear(*(range.a));
515 26 storres
  mpfr_clear(*(range.b));
516 26 storres
  free(range.a);
517 26 storres
  free(range.b);
518 26 storres
} /* End pobyso_free_range. */
519 26 storres

520 26 storres
node*
521 26 storres
pobyso_fp_minimax_canonical_monomials_base(node *function,
522 26 storres
                                      int degree,
523 26 storres
                                      chain *formats,
524 26 storres
                                      chain *points,
525 26 storres
                                      mpfr_t lowerBound,
526 26 storres
                                      mpfr_t upperBound,
527 26 storres
                                      int fpFixedArg,
528 26 storres
                                      int absRel,
529 26 storres
                                      node *constPart,
530 26 storres
                                      node *minimax)
531 26 storres
{
532 26 storres
  return(NULL);
533 26 storres
} /* End pobyso_fp_minimax_canonical_monomials_base. */
534 26 storres

535 26 storres
node*
536 26 storres
pobyso_parse_function(char *functionString,
537 26 storres
                      char *freeVariableNameString)
538 26 storres
{
539 26 storres
  if (functionString == NULL || freeVariableNameString == NULL)
540 26 storres
  {
541 26 storres
    pobyso_error_message("pobyso_parse_function",
542 26 storres
                        "NULL_POINTER_ARGUMENT",
543 26 storres
                        "One of the arguments is a NULL pointer");
544 26 storres
    return(NULL);
545 26 storres
  }
546 26 storres
  return(parseString(functionString));
547 26 storres

548 26 storres
} /* End pobyso_parse_function */
549 26 storres

550 26 storres
node*
551 26 storres
pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
552 26 storres
                                                  unsigned int mode,
553 26 storres
                                                  mpfr_t lowerBound,
554 26 storres
                                                  mpfr_t upperBound,
555 26 storres
                                                  mpfr_t eps)
556 26 storres
{
557 26 storres
  node *weight              = NULL;
558 26 storres
  node *bestApproxPolyNode  = NULL;
559 26 storres
  node *bestApproxHorner    = NULL;
560 26 storres
  node *errorNode           = NULL;
561 26 storres
  rangetype degreeRange;
562 26 storres
  mpfr_t quality;
563 26 storres
  mpfr_t currentError;
564 26 storres
  unsigned int degree;
565 26 storres

566 26 storres
  /* Check the parameters. */
567 26 storres
  if (functionNode == NULL)
568 26 storres
  {
569 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
570 26 storres
                        "NULL_POINTER_ARGUMENT",
571 26 storres
                        "functionNode is a NULL pointer");
572 26 storres
    return(NULL);
573 26 storres
  }
574 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
575 26 storres
  {
576 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
577 26 storres
                        "INCOHERENT_INPUT_DATA",
578 26 storres
                        "the lower_bound >= upper_bound");
579 26 storres
    return(NULL);
580 26 storres
  }
581 26 storres
  /* Set the weight. */
582 26 storres
  if (mode == POBYSO_ABSOLUTE)
583 26 storres
  {
584 26 storres
    /* Set the weight to 1 for the ABSOLUTE_MODE. */
585 26 storres
    weight = makeConstantDouble(1.0);
586 26 storres
  }
587 26 storres
  else
588 26 storres
  {
589 26 storres
    if (mode == POBYSO_RELATIVE)
590 26 storres
    {
591 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
592 26 storres
                          "NOT_IMPLEMENTED",
593 26 storres
                          "the search for relative error is not implemented yet");
594 26 storres
      return(NULL);
595 26 storres
    }
596 26 storres
    else
597 26 storres
    {
598 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
599 26 storres
                          "INCOHERENT_INPUT_DATA",
600 26 storres
                          "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
601 26 storres
      return(NULL);
602 26 storres
    }
603 26 storres
  }
604 26 storres
  //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
605 26 storres
  degreeRange = guessDegree(functionNode,
606 26 storres
                            weight,
607 26 storres
                            lowerBound,
608 26 storres
                            upperBound,
609 26 storres
                            eps,
610 26 storres
                            POBYSO_GUESS_DEGREE_BOUND);
611 26 storres
  degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
612 26 storres
  //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
613 26 storres
  fprintf(stderr, "Guessed degree: ");
614 26 storres
  mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
615 26 storres
  fprintf(stderr, " - as int: %u\n", degree);
616 26 storres
  /* Reduce the degree by 1 in the foolish hope it could work. */
617 26 storres
  if (degree > 0) degree--;
618 26 storres
  /* Both elements of degreeRange where "inited" within guessDegree. */
619 26 storres
  mpfr_clear(*(degreeRange.a));
620 26 storres
  mpfr_clear(*(degreeRange.b));
621 26 storres
  free(degreeRange.a);
622 26 storres
  free(degreeRange.b);
623 26 storres
  /* Mimic the default behavior of interactive Sollya. */
624 26 storres
  mpfr_init(quality);
625 26 storres
  mpfr_set_d(quality, 1e-5, MPFR_RNDN);
626 26 storres
  mpfr_init2(currentError, getToolPrecision());
627 26 storres
  mpfr_set_inf(currentError,1);
628 26 storres

629 26 storres
  /* Try to refine the initial guess: loop with increasing degrees until
630 26 storres
   * we find a satisfactory one. */
631 26 storres
  while(mpfr_cmp(currentError, eps) > 0)
632 26 storres
  {
633 26 storres
    /* Get rid of the previous polynomial, if any. */
634 26 storres
    if (bestApproxPolyNode != NULL)
635 26 storres
    {
636 26 storres
      free_memory(bestApproxPolyNode);
637 26 storres
    }
638 26 storres
    fprintf(stderr, "Degree: %u\n", degree);
639 26 storres
    fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
640 26 storres
    /* Try to find a polynomial with the guessed degree. */
641 26 storres
    bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
642 26 storres
                                                            weight,
643 26 storres
                                                            degree,
644 26 storres
                                                            lowerBound,
645 26 storres
                                                            upperBound,
646 26 storres
                                                            quality);
647 26 storres

648 26 storres
    if (bestApproxPolyNode == NULL)
649 26 storres
    {
650 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
651 26 storres
                          "INTERNAL_ERROR",
652 26 storres
                          "could not compute the bestApproxPolyNode");
653 26 storres
      mpfr_clear(currentError);
654 26 storres
      mpfr_clear(quality);
655 26 storres
      return(NULL);
656 26 storres
    }
657 26 storres

658 26 storres
    setDisplayMode(DISPLAY_MODE_DECIMAL);
659 26 storres
    fprintTree(stderr, bestApproxPolyNode);
660 26 storres
    fprintf(stderr, "\n\n");
661 26 storres

662 26 storres
    errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
663 26 storres
    /* Check the error with the computed polynomial. */
664 26 storres
    uncertifiedInfnorm(currentError,
665 26 storres
                        errorNode,
666 26 storres
                        lowerBound,
667 26 storres
                        upperBound,
668 26 storres
                        POBYSO_INF_NORM_NUM_POINTS,
669 26 storres
                        getToolPrecision());
670 26 storres
    fprintf(stderr, "Inf norm: ");
671 26 storres
    mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
672 26 storres
    fprintf(stderr, "\n\n");
673 26 storres
    /* Free the errorNode but not the bestApproxPolyNode (we need it if
674 26 storres
     * we exit the loop at the next iteration). */
675 26 storres
    free_memory(errorNode);
676 26 storres
    degree++;
677 26 storres
  }
678 26 storres
  /* Use an intermediate variable, since horner() creates a new node
679 26 storres
   * and does not reorder the argument "in place". This allows for the memory
680 26 storres
   * reclaim of bestApproxHorner.
681 26 storres
   */
682 26 storres
  bestApproxHorner = horner(bestApproxPolyNode);
683 26 storres
  free_memory(bestApproxPolyNode);
684 26 storres
  mpfr_clear(currentError);
685 26 storres
  mpfr_clear(quality);
686 26 storres
  free_memory(weight);
687 26 storres
  return(bestApproxHorner);
688 26 storres
} /* End pobyso_remez_approx_canonical_monomials_base_for_error */
689 26 storres

690 26 storres
node*
691 26 storres
pobyso_remez_canonical_monomials_base(node *function,
692 26 storres
                                     node *weight,
693 26 storres
                                     unsigned int degree,
694 26 storres
                                     mpfr_t lowerBound,
695 26 storres
                                     mpfr_t upperBound,
696 26 storres
                                     mpfr_t quality)
697 26 storres
{
698 26 storres
  node  *bestApproxPoly = NULL;
699 26 storres
  chain *monomials      = NULL;
700 26 storres
  chain *curMonomial    = NULL;
701 26 storres

702 26 storres
  mpfr_t satisfying_error;
703 26 storres
  mpfr_t target_error;
704 26 storres

705 26 storres
  /* Argument checking */
706 26 storres
  /* Function tree. */
707 26 storres
  if (function == NULL)
708 26 storres
  {
709 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
710 26 storres
                        "NULL_POINTER_ARGUMENT",
711 26 storres
                        "the \"function\" argument is a NULL pointer");
712 26 storres
    return(NULL);
713 26 storres
  }
714 26 storres
  if (weight == NULL)
715 26 storres
  {
716 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
717 26 storres
                        "NULL_POINTER_ARGUMENT",
718 26 storres
                        "the \"weight\" argument is a NULL pointer");
719 26 storres
    return(NULL);
720 26 storres
  }
721 26 storres
  /* Check the bounds. */
722 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
723 26 storres
  {
724 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
725 26 storres
                        "INCOHERENT_IMPUT_DATA",
726 26 storres
                        "the lower_bound >= upper_bound");
727 26 storres
    return(NULL);
728 26 storres
  }
729 26 storres
  /* The quality must be a non null positive number. */
730 26 storres
  if (mpfr_sgn(quality) <= 0)
731 26 storres
  {
732 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
733 26 storres
                        "INCOHERENT_INPUT_DATA",
734 26 storres
                        "the quality <= 0");
735 26 storres
  }
736 26 storres
  /* End argument checking. */
737 26 storres
  /* Create the monomials nodes chains. */
738 26 storres
  monomials = pobyso_create_canonical_monomials_base(degree);
739 26 storres
  fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
740 26 storres
  if (monomials == NULL || (lengthChain(monomials) != degree + 1))
741 26 storres
  {
742 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
743 26 storres
                        "CHAIN_CREATION_ERROR",
744 26 storres
                        "could not create the monomials chain");
745 26 storres
    return(NULL);
746 26 storres
  }
747 26 storres
  curMonomial = monomials;
748 26 storres

749 26 storres
  while (curMonomial != NULL)
750 26 storres
  {
751 26 storres
    fprintf(stderr, "monomial tree: ");
752 26 storres
    //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
753 26 storres
    fprintTree(stderr, (node*)(curMonomial->value));
754 26 storres
    fprintf(stderr, "\n");
755 26 storres
    curMonomial = curMonomial->next;
756 26 storres
  }
757 26 storres

758 26 storres
  /* Deal with NULL weight. */
759 26 storres
  if (weight == NULL)
760 26 storres
  {
761 26 storres
    weight = makeConstantDouble(1.0);
762 26 storres
  }
763 26 storres
  /* Compute the best polynomial with the required quality.
764 26 storres
     The behavior is as if satisfying error and target_error had
765 26 storres
     not been used.*/
766 26 storres
  mpfr_init(satisfying_error);
767 26 storres
  mpfr_init(target_error);
768 26 storres
  mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
769 26 storres
  mpfr_set_inf(target_error, 1);
770 26 storres

771 26 storres

772 26 storres
  fprintf(stderr, "satisfying_error: ");
773 26 storres
  mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
774 26 storres
  fprintf(stderr, ".\n");
775 26 storres
  fprintf(stderr, "target_error: ");
776 26 storres
  mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
777 26 storres
  fprintf(stderr, ".\n");
778 26 storres

779 26 storres
  fprintf(stderr,
780 26 storres
          "current precision: %li\n", getToolPrecision());
781 26 storres
  /* Call the Sollya function. */
782 26 storres
  bestApproxPoly = remez(function,
783 26 storres
                          weight,
784 26 storres
                          monomials,
785 26 storres
                          lowerBound,
786 26 storres
                          upperBound,
787 26 storres
                          quality,
788 26 storres
                          satisfying_error,
789 26 storres
                          target_error,
790 26 storres
                          getToolPrecision());
791 26 storres

792 26 storres
  mpfr_clear(satisfying_error);
793 26 storres
  mpfr_clear(target_error);
794 26 storres
  pobyso_free_chain_of_nodes(monomials);
795 26 storres

796 26 storres
  return(bestApproxPoly);
797 26 storres
} /* End pobyso_remez_canonical_monomials_base. */
798 26 storres

799 26 storres
#endif
800 26 storres
801 26 storres
void
802 26 storres
pobyso_error_message(char *functionName, char *messageName, char* message)
803 26 storres
{
804 26 storres
  fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
805 26 storres
} /* End pobyso_error_message */