Statistiques
| Révision :

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

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

1 26 storres
/** @file pobyso.c
2 137 storres
 * Integration of Sollya to C programs
3 26 storres
 *
4 26 storres
 * @author S.T.
5 26 storres
 * @date 2011-10-12
6 26 storres
 *
7 137 storres
 * @todo write pobyso_is_monomial function <br>
8 137 storres
 *       write pobyso_is_free_var_int_poson_power function
9 26 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 134 storres
/* @see pobyso.h#pobyso_autoprint */
36 26 storres
void
37 137 storres
pobyso_autoprint(sollya_obj_t objSo)
38 26 storres
{
39 137 storres
  sollya_lib_autoprint(objSo, NULL);
40 26 storres
} /* End pobyso_autoprint. */
41 26 storres
42 134 storres
/* @see pobyso.h#pobyso_get_verbosity */
43 134 storres
int
44 134 storres
pobyso_get_verbosity()
45 134 storres
{
46 134 storres
  sollya_obj_t verbositySo = NULL;
47 134 storres
  int verbosity            = 0;
48 134 storres
49 134 storres
  verbositySo = sollya_lib_get_verbosity();
50 134 storres
  sollya_lib_get_constant_as_int(&verbosity,verbositySo);
51 134 storres
  sollya_lib_clear_obj(verbositySo);
52 134 storres
  return verbosity;
53 134 storres
} /* End pobyso_get_verbosity. */
54 134 storres
55 137 storres
/** @see pobyso.h#pobyso_is_constant_expression
56 137 storres
 * Strategy: rely on sollya_lib_get_constant. It return 1, when the
57 137 storres
 * expression is constant.
58 137 storres
 */
59 133 storres
int
60 133 storres
pobyso_is_constant_expression(sollya_obj_t obj_to_test)
61 133 storres
{
62 133 storres
  mpfr_t dummy;
63 133 storres
  int test;
64 134 storres
  int initial_verbosity_level      = 0;
65 134 storres
66 133 storres
  /* Test argument. */
67 133 storres
  if (obj_to_test == NULL)
68 133 storres
  {
69 137 storres
    pobyso_error_message("pobyso_is_constant_expression",
70 133 storres
                        "NULL_POINTER_ARGUMENT",
71 133 storres
                        "The expression is a NULL pointer");
72 133 storres
    return 0;
73 133 storres
  }
74 134 storres
  initial_verbosity_level = pobyso_set_verbosity_off();
75 128 storres
76 133 storres
  if (! sollya_lib_obj_is_function(obj_to_test))
77 133 storres
  {
78 134 storres
    pobyso_set_verbosity_to(initial_verbosity_level);
79 133 storres
    return 0;
80 133 storres
  }
81 133 storres
  mpfr_init2(dummy,64);
82 133 storres
  /* Call to previous Sollya function resets verbosity. */
83 137 storres
  /* Todo: change verbosity suppression strategy with a message call back. */
84 134 storres
  pobyso_set_verbosity_off();
85 133 storres
  test = sollya_lib_get_constant(dummy, obj_to_test);
86 134 storres
  pobyso_set_verbosity_to(initial_verbosity_level);
87 138 storres
  if (test)
88 137 storres
  {
89 138 storres
    if(mpfr_number_p(dummy))
90 138 storres
    {
91 138 storres
      mpfr_clear(dummy);
92 138 storres
      return test;
93 138 storres
    }
94 138 storres
    else
95 138 storres
    {
96 138 storres
      mpfr_clear(dummy);
97 138 storres
      return 0;
98 138 storres
    }
99 137 storres
  }
100 138 storres
  else
101 137 storres
  {
102 137 storres
    return 0;
103 137 storres
  }
104 138 storres
} /* End pobyso_is_constant_expression. */
105 137 storres
106 137 storres
/** @see pobyso.h#pobyso_is_monomial. */
107 137 storres
int
108 137 storres
pobyso_is_int(pobyso_func_exp_t exprSo)
109 137 storres
{
110 137 storres
  mpfr_t float1M;
111 137 storres
  mpfr_t float2M;
112 137 storres
  mpfr_t tempFloat1M;
113 137 storres
  mpfr_t tempFloat2M;
114 137 storres
  mpfr_prec_t prec;
115 137 storres
  int64_t asInt;
116 137 storres
  sollya_obj_t newConstantSo = NULL;
117 137 storres
  /* Arguments check. */
118 137 storres
  if (exprSo == NULL)
119 137 storres
  {
120 137 storres
    pobyso_error_message("pobyso_is_free_var_posze_int_power",
121 137 storres
                        "NULL_POINTER_ARGUMENT",
122 137 storres
                        "The expression is a NULL pointer");
123 137 storres
    return 0;
124 137 storres
  }
125 137 storres
  //fprintf(stdout, "Not NULL.\n"); pobyso_autoprint(exprSo);
126 137 storres
  if (! pobyso_is_constant_expression(exprSo)) return 0;
127 137 storres
  if (! sollya_lib_get_constant_as_int64(&asInt, exprSo)) return 0;
128 137 storres
  if (asInt == INT64_MIN || asInt == INT64_MAX) return 0;
129 137 storres
  /* Some constant integer expression can't have their precision computed
130 137 storres
   * (e.g. cos(pi). */
131 137 storres
  if (! sollya_lib_get_prec_of_constant(&prec, exprSo))
132 137 storres
  {
133 137 storres
    mpfr_init2(tempFloat1M, 165);
134 137 storres
    mpfr_init2(tempFloat2M, 165);
135 137 storres
    mpfr_abs(tempFloat1M, tempFloat1M, MPFR_RNDN);
136 137 storres
    mpfr_log2(tempFloat2M, tempFloat1M, MPFR_RNDU);
137 137 storres
    mpfr_rint_ceil(tempFloat1M, tempFloat2M, MPFR_RNDU);
138 137 storres
    prec = mpfr_get_si(tempFloat1M, MPFR_RNDN) + 10;
139 137 storres
    if (prec < 1024) prec = 1024;
140 137 storres
    mpfr_clear(tempFloat1M);
141 137 storres
    mpfr_clear(tempFloat2M);
142 137 storres
    mpfr_init2(float1M, prec);
143 137 storres
    if (!sollya_lib_get_constant(float1M, exprSo))
144 137 storres
    {
145 137 storres
      mpfr_clear(float1M);
146 137 storres
      return 0;
147 137 storres
    }
148 137 storres
  }
149 137 storres
  else /* Precision could be given. */
150 137 storres
  {
151 137 storres
    mpfr_init2(float1M, prec);
152 137 storres
    if (! sollya_lib_get_constant(float1M, exprSo))
153 137 storres
    {
154 137 storres
      mpfr_clear(float1M);
155 137 storres
      return 0;
156 137 storres
    }
157 137 storres
  }
158 137 storres
  if (mpfr_nan_p(float1M) || mpfr_inf_p(float1M))
159 137 storres
  {
160 137 storres
    mpfr_clear(float1M);
161 137 storres
    return 0;
162 137 storres
  }
163 137 storres
  if ((newConstantSo = sollya_lib_constant_from_int64(asInt)) == NULL)
164 137 storres
  {
165 137 storres
    mpfr_clear(float1M);
166 137 storres
    return 0;
167 137 storres
  }
168 137 storres
  sollya_lib_get_prec_of_constant(&prec, newConstantSo);
169 137 storres
  mpfr_init2(float2M, prec);
170 137 storres
  sollya_lib_get_constant(float2M, newConstantSo);
171 137 storres
  if (mpfr_cmp(float1M, float2M) == 0)
172 137 storres
  {
173 137 storres
    mpfr_clear(float1M);
174 137 storres
    mpfr_clear(float2M);
175 137 storres
    sollya_lib_clear_obj(newConstantSo);
176 137 storres
    return 1;
177 137 storres
  }
178 137 storres
  else
179 137 storres
  {
180 137 storres
    mpfr_clear(float1M);
181 137 storres
    mpfr_clear(float2M);
182 137 storres
    sollya_lib_clear_obj(newConstantSo);
183 137 storres
    return 0;
184 137 storres
  }
185 137 storres
} /* End pobyso_is_int. */
186 137 storres
187 137 storres
/** @see pobyso.h#pobyso_is_monomial.
188 137 storres
 * Strategy: check that the object is a functional expression and
189 137 storres
 * if so check that it is made of cte * free_var ^ some_power where :
190 137 storres
 * - cte is a constant expression (a per pobyso_is_constant_experession;
191 137 storres
 * - some_power is a positive or null power. t*/
192 137 storres
int
193 137 storres
pobyso_is_monomial(sollya_obj_t objSo)
194 137 storres
{
195 138 storres
  int arity;
196 138 storres
  sollya_obj_t subFun1So = NULL;
197 138 storres
  sollya_obj_t subFun2So = NULL;
198 138 storres
  sollya_obj_t subFun3So = NULL;
199 138 storres
  sollya_base_function_t head = 0;
200 138 storres
  long int exponent = 0;
201 138 storres
  long int exprIntValue = 0;
202 138 storres
203 138 storres
  /* Arguments check. */
204 138 storres
  if (objSo == NULL)
205 138 storres
  {
206 138 storres
    pobyso_error_message("pobyso_is_monomial",
207 138 storres
                        "NULL_POINTER_ARGUMENT",
208 138 storres
                        "The expression is a NULL pointer");
209 138 storres
    return 0;
210 138 storres
  }
211 138 storres
  /* The object must be a function. */
212 138 storres
  if (! sollya_lib_obj_is_function(objSo)) return 0;
213 138 storres
  /* Check if it is the 1 constant. */
214 138 storres
  if (pobyso_is_int(objSo))
215 138 storres
  {
216 138 storres
    if (! sollya_lib_get_constant_as_int64(&exprIntValue, objSo))
217 138 storres
    {
218 138 storres
      return 0;
219 138 storres
    }
220 138 storres
    else
221 138 storres
    {
222 138 storres
      if (exprIntValue == 1) return 1;
223 138 storres
      else return 0;
224 138 storres
    }
225 138 storres
  }
226 138 storres
  if (! sollya_lib_decompose_function(objSo,
227 138 storres
                                      &head,
228 138 storres
                                      &arity,
229 138 storres
                                      &subFun1So,
230 138 storres
                                      &subFun2So,
231 138 storres
                                      NULL)) return 0;
232 138 storres
  if (arity > 2)
233 138 storres
  {
234 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
235 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
236 138 storres
    return 0;
237 138 storres
  }
238 138 storres
  /* Arity == 1 must be free variable by itself. */
239 138 storres
  if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE)
240 138 storres
  {
241 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
242 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
243 138 storres
    return 1;
244 138 storres
  }
245 138 storres
  else
246 138 storres
  {
247 138 storres
    /* Another expression. Must be: free variable  ^ poze integer. */
248 138 storres
    if (arity == 2)
249 138 storres
    {
250 138 storres
      if (head != SOLLYA_BASE_FUNC_POW)
251 138 storres
      {
252 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
253 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
254 138 storres
        return 0;
255 138 storres
      }
256 138 storres
      if (! pobyso_is_int(subFun2So))
257 138 storres
      {
258 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
259 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
260 138 storres
        return 0;
261 138 storres
      }
262 138 storres
      if (! sollya_lib_get_constant_as_int64(&exponent, subFun2So))
263 138 storres
      {
264 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
265 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
266 138 storres
        return 0;
267 138 storres
      }
268 138 storres
      if (exponent < 0)
269 138 storres
      {
270 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
271 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
272 138 storres
        return 0;
273 138 storres
      }
274 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
275 138 storres
      /* Check that the first subfunction is the free variable. */
276 138 storres
      if (! sollya_lib_decompose_function(subFun1So,
277 138 storres
                                          &head,
278 138 storres
                                          &arity,
279 138 storres
                                          &subFun2So,
280 138 storres
                                          &subFun3So,
281 138 storres
                                          NULL))
282 138 storres
      {
283 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
284 138 storres
        return 0;
285 138 storres
      }
286 138 storres
      if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE)
287 138 storres
      {
288 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
289 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
290 138 storres
        if (subFun3So != NULL) sollya_lib_clear_obj(subFun3So);
291 138 storres
        return 1;
292 138 storres
      }
293 138 storres
      else
294 138 storres
      {
295 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
296 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
297 138 storres
        return 0;
298 138 storres
      }
299 138 storres
    } /* End if arity == 2. */
300 138 storres
  } /* End else if arity == 1. */
301 137 storres
  return 0;
302 137 storres
} /* End pobyso_is_monomial. */
303 137 storres
304 138 storres
/** @see pobyso.h#pobyso_is_polynomial_term.
305 138 storres
 * Strategy: check that the object is a functional expression and
306 138 storres
 * if so check that it is made of cte * monomial.
307 138 storres
 */
308 138 storres
int
309 138 storres
pobyso_is_polynomial_term(sollya_obj_t objSo)
310 138 storres
{
311 138 storres
  int arity;
312 138 storres
  sollya_obj_t subFun1So = NULL;
313 138 storres
  sollya_obj_t subFun2So = NULL;
314 138 storres
  sollya_base_function_t head = 0;
315 138 storres
316 138 storres
  /* Arguments check. */
317 138 storres
  if (objSo == NULL)
318 138 storres
  {
319 138 storres
    pobyso_error_message("pobyso_is_polynomial_term",
320 138 storres
                        "NULL_POINTER_ARGUMENT",
321 138 storres
                        "The expression is a NULL pointer");
322 138 storres
    return 0;
323 138 storres
  }
324 138 storres
  /* The object must be a function. */
325 138 storres
  if (! sollya_lib_obj_is_function(objSo)) return 0;
326 138 storres
  /* A constant expression is degree 0 polynomial term. */
327 138 storres
  if (pobyso_is_constant_expression(objSo)) return 1;
328 138 storres
  /* A monomial is a polynomial term. */
329 138 storres
  if (pobyso_is_monomial(objSo)) return 1;
330 138 storres
  /* Decompose the functional expression and study the elements. */
331 138 storres
  if (! sollya_lib_decompose_function(objSo,
332 138 storres
                                      &head,
333 138 storres
                                      &arity,
334 138 storres
                                      &subFun1So,
335 138 storres
                                      &subFun2So,
336 138 storres
                                      NULL)) return 0;
337 138 storres
  /* Monomial case has been dealt with abobe. */
338 138 storres
  if (arity != 2)
339 138 storres
  {
340 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
341 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
342 138 storres
    return 0;
343 138 storres
  }
344 138 storres
  /* The expression must be: cte  * monomial or monomial * cte. */
345 138 storres
  if (head != SOLLYA_BASE_FUNC_MUL)
346 138 storres
  {
347 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
348 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
349 138 storres
    return 0;
350 138 storres
  }
351 138 storres
  if (! pobyso_is_monomial(subFun2So))
352 138 storres
  {
353 138 storres
    if (! pobyso_is_constant_expression(subFun2So) ||
354 138 storres
        ! pobyso_is_monomial(subFun1So))
355 138 storres
    {
356 138 storres
      if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
357 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
358 138 storres
      return 0;
359 138 storres
    }
360 138 storres
  }
361 138 storres
  else
362 138 storres
  {
363 138 storres
    if (! pobyso_is_constant_expression(subFun1So))
364 138 storres
    {
365 138 storres
      if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
366 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
367 138 storres
      return 0;
368 138 storres
    }
369 138 storres
  }
370 138 storres
  return 1;
371 138 storres
} /* End pobyso_is_polynomial_term. */
372 133 storres
/** @see pobyso.h#pobyso_new_monomial. */
373 133 storres
pobyso_func_exp_t
374 134 storres
pobyso_new_monomial(pobyso_func_exp_t coefficientSo, long degree)
375 133 storres
{
376 134 storres
  sollya_obj_t degreeSo   = NULL;
377 134 storres
  sollya_obj_t varToPowSo = NULL;
378 134 storres
  sollya_obj_t monomialSo = NULL;
379 134 storres
  int initial_verbosity_level = 0;
380 134 storres
381 134 storres
  /* Arguments check. */
382 134 storres
  if (coefficientSo == NULL)
383 133 storres
  {
384 133 storres
    pobyso_error_message("pobyso_parse_string",
385 133 storres
                        "NULL_POINTER_ARGUMENT",
386 133 storres
                        "The expression is a NULL pointer");
387 134 storres
    return NULL;
388 133 storres
  }
389 134 storres
  if (! pobyso_is_constant_expression(coefficientSo))
390 133 storres
  {
391 134 storres
    return NULL;
392 133 storres
  }
393 133 storres
  if (degree < 0)
394 133 storres
  {
395 134 storres
    pobyso_error_message("pobyso_new_monomial",
396 134 storres
                        "NEGATIVE_DEGREE_ARGUMENT",
397 134 storres
                        "The degree is a negative integer");
398 134 storres
    return NULL;
399 133 storres
  }
400 134 storres
  /* If degree == 0, just return a copy of the coefficient. Do not
401 134 storres
   * return the coefficient itself to avoid "double clear" issues. */
402 134 storres
  if (degree == 0)
403 134 storres
  {
404 134 storres
    initial_verbosity_level = pobyso_set_verbosity_off();
405 134 storres
    monomialSo = sollya_lib_copy_obj(coefficientSo);
406 134 storres
    pobyso_set_verbosity_to(initial_verbosity_level);
407 134 storres
  }
408 134 storres
  degreeSo    = sollya_lib_constant_from_int64(degree);
409 134 storres
  varToPowSo  = sollya_lib_build_function_pow(sollya_lib_free_variable(),
410 134 storres
                                              degreeSo);
411 134 storres
  /* Do not use the "build" version to avoid "eating up" the coefficient. */
412 134 storres
  monomialSo = sollya_lib_mul(coefficientSo,varToPowSo);
413 134 storres
  sollya_lib_clear_obj(varToPowSo);
414 134 storres
  /* Do not clear degreeSa: it was "eaten" by sollya-lib_build_function. */
415 134 storres
  return monomialSo;
416 133 storres
} /* End pobyso_new_monomial. */
417 133 storres
418 127 storres
/* @see pobyso.h#pobyso_parse_string*/
419 130 storres
pobyso_func_exp_t
420 26 storres
pobyso_parse_string(const char* expression)
421 26 storres
{
422 128 storres
  int expressionLength, i;
423 127 storres
  char *expressionWithSemiCo;
424 136 storres
  sollya_obj_t expressionSo;
425 128 storres
426 127 storres
  /* Arguments check. */
427 26 storres
  if (expression == NULL)
428 26 storres
  {
429 26 storres
    pobyso_error_message("pobyso_parse_string",
430 26 storres
                        "NULL_POINTER_ARGUMENT",
431 26 storres
                        "The expression is a NULL pointer");
432 134 storres
    return  NULL;
433 26 storres
  }
434 127 storres
  expressionLength = strlen(expression);
435 127 storres
  if (expressionLength == 0)
436 127 storres
  {
437 127 storres
    pobyso_error_message("pobyso_parse_string",
438 127 storres
                        "EMPTY_STRING_ARGUMENT",
439 127 storres
                        "The expression is an empty string");
440 134 storres
    return NULL;
441 127 storres
  }
442 128 storres
  /* Search from the last char of the expression until, whichever happens
443 128 storres
   * first:
444 128 storres
   * a ";" is found;
445 128 storres
   * a char  != ';' is found the the ";" is appended.
446 128 storres
   * If the head of the string is reached before any of these two events happens
447 128 storres
   * return an error.
448 128 storres
   */
449 128 storres
  for (i = expressionLength - 1 ; i >= 0 ; i--)
450 127 storres
  {
451 128 storres
    if (expression[i] == ';') /* Nothing special to do:
452 128 storres
                                 try to parse the string*/
453 127 storres
    {
454 136 storres
      expressionSo = sollya_lib_parse_string(expression);
455 136 storres
      if (sollya_lib_obj_is_error(expressionSo))
456 136 storres
      {
457 136 storres
        sollya_lib_clear_obj(expressionSo);
458 136 storres
        return NULL;
459 136 storres
      }
460 136 storres
      else
461 136 storres
      {
462 136 storres
        return expressionSo;
463 136 storres
      }
464 127 storres
    }
465 128 storres
    else
466 128 storres
    {
467 128 storres
      if (expression[i] == ' ' || expression[i] == '\t') /* A blank,
468 128 storres
                                                           just continue. */
469 128 storres
      {
470 128 storres
        continue;
471 128 storres
      }
472 128 storres
      else /* a character != ';' and from a blank: create the ';'ed string. */
473 128 storres
      {
474 128 storres
        /* Create a new string for the expression, add the ";" and
475 128 storres
         * and call sollya_lib_parse_string. */
476 128 storres
        expressionWithSemiCo = calloc(i + 3, sizeof(char));
477 128 storres
        if (expressionWithSemiCo == NULL)
478 128 storres
        {
479 128 storres
          pobyso_error_message("pobyso_parse_string",
480 128 storres
                                "MEMORY_ALLOCATION_ERROR",
481 128 storres
                                "Could not allocate the expression string");
482 134 storres
          return NULL;
483 128 storres
        }
484 132 storres
        strncpy(expressionWithSemiCo, expression, i+1);
485 128 storres
        expressionWithSemiCo[i + 1] = ';';
486 128 storres
        expressionWithSemiCo[i + 2] = '\0';
487 136 storres
        expressionSo = sollya_lib_parse_string(expressionWithSemiCo);
488 128 storres
        free(expressionWithSemiCo);
489 136 storres
        if (sollya_lib_obj_is_error(expressionSo))
490 136 storres
        {
491 136 storres
          sollya_lib_clear_obj(expressionSo);
492 136 storres
          return NULL;
493 136 storres
        }
494 136 storres
        else
495 136 storres
        {
496 136 storres
          return expressionSo;
497 136 storres
        }
498 128 storres
      } /* End character != ';' and from a blank. */
499 128 storres
    /* Create a new string for the expression, add the ";" and
500 128 storres
     * and call sollya_lib_parse_string. */
501 128 storres
    } /* End else. */
502 128 storres
  } /* End for. */
503 128 storres
  /* We get here, it is because we did not find any char == anything different
504 128 storres
   * from ' ' or '\t': the string is empty.
505 128 storres
   */
506 128 storres
  pobyso_error_message("pobyso_parse_string",
507 128 storres
                       "ONLY_BLANK_ARGUMENT",
508 128 storres
                        "The expression is only made of blanks");
509 134 storres
  return NULL;
510 26 storres
} /* pobyso_parse_string */
511 26 storres
512 132 storres
pobyso_func_exp_t
513 132 storres
pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function,
514 132 storres
                                      long int degree,
515 132 storres
                                      pobyso_range_t interval,
516 132 storres
                                      pobyso_func_exp_t weight,
517 132 storres
                                      double quality,
518 132 storres
                                      pobyso_range_t bounds)
519 132 storres
{
520 134 storres
  sollya_obj_t  degreeSo = NULL;
521 134 storres
  sollya_obj_t qualitySo = NULL;
522 132 storres
523 134 storres
  degreeSo = sollya_lib_constant_from_int(degree);
524 134 storres
  qualitySo = sollya_lib_constant_from_double(quality);
525 132 storres
526 134 storres
  sollya_lib_clear_obj(degreeSo);
527 134 storres
  sollya_lib_clear_obj(qualitySo);
528 132 storres
  return NULL;
529 132 storres
} /* End pobyso_remez_canonical_monomials_base. */
530 132 storres
531 132 storres
int
532 132 storres
pobyso_set_canonical_on()
533 132 storres
{
534 132 storres
  pobyso_on_off_t currentCanonicalModeSo;
535 132 storres
  pobyso_on_off_t on;
536 132 storres
537 132 storres
  currentCanonicalModeSo = sollya_lib_get_canonical();
538 132 storres
  if (sollya_lib_is_on(currentCanonicalModeSo))
539 132 storres
  {
540 134 storres
    sollya_lib_clear_obj(currentCanonicalModeSo);
541 134 storres
    return POBYSO_ON;
542 132 storres
  }
543 132 storres
  else
544 132 storres
  {
545 134 storres
    on = sollya_lib_on();
546 134 storres
    sollya_lib_set_canonical(on);
547 134 storres
    sollya_lib_clear_obj(on);
548 134 storres
    sollya_lib_clear_obj(currentCanonicalModeSo);
549 134 storres
    return POBYSO_OFF;
550 132 storres
  }
551 132 storres
} /* End pobyso_set_canonical_on. */
552 132 storres
553 134 storres
int
554 128 storres
pobyso_set_verbosity_off()
555 26 storres
{
556 134 storres
  sollya_obj_t verbosityLevelZeroSo;
557 134 storres
  sollya_obj_t currentVerbosityLevelSo = NULL;
558 134 storres
  int currentVerbosityLevel = 0;
559 128 storres
560 134 storres
  currentVerbosityLevelSo = sollya_lib_get_verbosity();
561 134 storres
  sollya_lib_get_constant_as_int(&currentVerbosityLevel,
562 134 storres
                                 currentVerbosityLevelSo);
563 134 storres
  verbosityLevelZeroSo = sollya_lib_constant_from_int(0);
564 134 storres
  sollya_lib_set_verbosity(verbosityLevelZeroSo);
565 134 storres
  sollya_lib_clear_obj(verbosityLevelZeroSo);
566 134 storres
  sollya_lib_clear_obj(currentVerbosityLevelSo);
567 134 storres
  return currentVerbosityLevel;
568 26 storres
} /* End of pobyso_set_verbosity_off. */
569 26 storres
570 134 storres
int
571 134 storres
pobyso_set_verbosity_to(int newVerbosityLevel)
572 26 storres
{
573 134 storres
  int initialVerbosityLevel = 0;
574 134 storres
  sollya_obj_t initialVerbosityLevelSo = NULL;
575 134 storres
  sollya_obj_t newVerbosityLevelSo = NULL;
576 134 storres
577 134 storres
  initialVerbosityLevelSo = sollya_lib_get_verbosity();
578 134 storres
  sollya_lib_get_constant_as_int(&initialVerbosityLevel,
579 134 storres
                                 initialVerbosityLevelSo);
580 134 storres
  sollya_lib_clear_obj(initialVerbosityLevelSo);
581 134 storres
  if (newVerbosityLevel < 0)
582 26 storres
  {
583 26 storres
    pobyso_error_message("pobyso_set_verbosity_to",
584 134 storres
                        "INVALID_VALUE",
585 134 storres
                        "The new verbosity level is a negative number");
586 134 storres
    return initialVerbosityLevel;
587 26 storres
  }
588 134 storres
  newVerbosityLevelSo = sollya_lib_constant_from_int(newVerbosityLevel);
589 134 storres
  sollya_lib_set_verbosity(newVerbosityLevelSo);
590 134 storres
  sollya_lib_clear_obj(newVerbosityLevelSo);
591 134 storres
  return initialVerbosityLevel;
592 26 storres
} /* End of pobyso_set_verbosity_to. */
593 26 storres
594 134 storres
/**
595 134 storres
 * @see pobyso.h#pobyso_subpoly
596 134 storres
 */
597 134 storres
pobyso_func_exp_t
598 136 storres
pobyso_subpoly(pobyso_func_exp_t polynomialSo, long expsNum, long int* expsList)
599 132 storres
{
600 136 storres
  sollya_obj_t  expsListSo    = NULL;
601 136 storres
  sollya_obj_t* expsList_pso  = NULL;
602 136 storres
  sollya_obj_t subpoly        = NULL;
603 136 storres
  int i, j;
604 134 storres
605 134 storres
  /* Arguments check. */
606 134 storres
  if (polynomialSo == NULL) return NULL;
607 134 storres
  if (expsNum < 0) return NULL;
608 134 storres
  if (expsNum == 0) return sollya_lib_copy_obj(polynomialSo);
609 136 storres
  if (expsList == 0) return NULL;
610 136 storres
  /* Create a list of Sollya constants. */
611 136 storres
  expsList_pso = (sollya_obj_t*) malloc(expsNum * sizeof(sollya_obj_t));
612 136 storres
  if (expsList_pso == NULL)
613 132 storres
  {
614 134 storres
    pobyso_error_message("pobyso_subpoly",
615 134 storres
                          "MEMORY_ALLOCATION_ERROR",
616 136 storres
                          "Could not allocate the Sollya exponents list");
617 134 storres
    return NULL;
618 132 storres
  }
619 136 storres
  /* Fill up the list. */
620 134 storres
  for (i = 0 ; i < expsNum ; i++)
621 134 storres
  {
622 136 storres
    /* Abort if an exponent is negative. */
623 136 storres
    if (expsList[i] < 0 )
624 136 storres
    {
625 136 storres
      for (j = 0 ; j < i ; j++)
626 136 storres
      {
627 136 storres
        sollya_lib_clear_obj(expsList_pso[j]);
628 136 storres
      }
629 136 storres
      free(expsList_pso);
630 136 storres
      return NULL;
631 136 storres
    }
632 136 storres
    expsList_pso[i] = sollya_lib_constant_from_int64(expsList[i]);
633 134 storres
  } /* End for */
634 136 storres
  expsListSo = sollya_lib_list(expsList_pso, expsNum);
635 136 storres
  for (i = 0 ; i < expsNum ; i++)
636 134 storres
  {
637 136 storres
    sollya_lib_clear_obj(expsList_pso[i]);
638 136 storres
  }
639 136 storres
  free(expsList_pso);
640 136 storres
  if (expsListSo == NULL)
641 136 storres
  {
642 134 storres
    pobyso_error_message("pobyso_subpoly",
643 134 storres
                          "LIST_CREATIONERROR",
644 134 storres
                          "Could not create the exponents list");
645 134 storres
    return NULL;
646 134 storres
  }
647 134 storres
  subpoly = sollya_lib_subpoly(polynomialSo, expsListSo);
648 136 storres
  pobyso_autoprint(expsListSo);
649 134 storres
  sollya_lib_clear_obj(expsListSo);
650 134 storres
  return subpoly;
651 134 storres
} /* pobyso_subpoly. */
652 132 storres
653 26 storres
/* Attic from the sollya_lib < 4. */
654 26 storres
#if 0
655 26 storres
chain*
656 26 storres
pobyso_create_canonical_monomials_base(const unsigned int degree)
657 26 storres
{
658 26 storres
  int i    = 0;
659 26 storres
  chain *monomials  = NULL;
660 26 storres
  node  *monomial   = NULL;
661 26 storres

662 26 storres
  for(i = degree ; i >= 0  ; i--)
663 26 storres
  {
664 26 storres
     monomial   = makePow(makeVariable(), makeConstantDouble((double)i));
665 26 storres
     monomials  = addElement(monomials, monomial);
666 26 storres
     fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
667 26 storres
  }
668 26 storres
  if (monomials == NULL)
669 26 storres
  {
670 26 storres
    pobyso_error_message("pobyso_create_canonical_monomial_base",
671 26 storres
                        "CHAIN_CREATION_ERROR",
672 26 storres
                        "Could not create the monomials chain");
673 26 storres
    return(NULL);
674 26 storres
  }
675 26 storres
  return(monomials);
676 26 storres
} /* End pobyso_create_canonical_monomials_base. */
677 26 storres

678 26 storres
chain*
679 26 storres
pobyso_create_chain_from_int_array(int* intArray,
680 26 storres
                                    const unsigned int arrayLength)
681 26 storres
{
682 26 storres
  int i = 0;
683 26 storres
  chain *newChain = NULL;
684 26 storres
  int *currentInt;
685 26 storres

686 26 storres
  if (arrayLength == 0) return(NULL);
687 26 storres
  if (intArray == NULL)
688 26 storres
  {
689 26 storres
    pobyso_error_message("pobyso_create_chain_from_int_array",
690 26 storres
                        "NULL_POINTER_ARGUMENT",
691 26 storres
                        "The array is a NULL pointer");
692 26 storres
    return(NULL);
693 26 storres
  }
694 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
695 26 storres
  {
696 26 storres
    currentInt = malloc(sizeof(int));
697 26 storres
    if (currentInt == NULL)
698 26 storres
    {
699 26 storres
      pobyso_error_message("pobyso_create_chain_from_int_array",
700 26 storres
                          "MEMORY_ALLOCATION_ERROR",
701 26 storres
                          "Could not allocate one of the integers");
702 26 storres
      freeChain(newChain, free);
703 26 storres
      return(NULL);
704 26 storres
    }
705 26 storres
    *currentInt = intArray[i];
706 26 storres
    newChain = addElement(newChain, currentInt);
707 26 storres
  }
708 26 storres
  return(newChain);
709 26 storres
} // End pobyso_create_chain_from_int_array. */
710 26 storres

711 26 storres
chain*
712 26 storres
pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
713 26 storres
                                        const unsigned int arrayLength)
714 26 storres
{
715 26 storres
  int i = 0;
716 26 storres
  chain *newChain = NULL;
717 26 storres
  unsigned int *currentInt;
718 26 storres

719 26 storres
  /* Argument checking. */
720 26 storres
  if (arrayLength == 0) return(NULL);
721 26 storres
  if (intArray == NULL)
722 26 storres
  {
723 26 storres
    pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
724 26 storres
                        "NULL_POINTER_ARGUMENT",
725 26 storres
                        "The array is a NULL pointer");
726 26 storres
    return(NULL);
727 26 storres
  }
728 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
729 26 storres
  {
730 26 storres
    currentInt = malloc(sizeof(unsigned int));
731 26 storres
    if (currentInt == NULL)
732 26 storres
    {
733 26 storres
      pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
734 26 storres
                          "MEMORY_ALLOCATION_ERROR",
735 26 storres
                          "Could not allocate one of the integers");
736 26 storres
      freeChain(newChain, free);
737 26 storres
      return(NULL);
738 26 storres
    }
739 26 storres
    *currentInt = intArray[i];
740 26 storres
    newChain = addElement(newChain, currentInt);
741 26 storres
  }
742 26 storres
  return(newChain);
743 26 storres
} // End pobyso_create_chain_from_unsigned_int_array. */
744 26 storres

745 26 storres
node*
746 26 storres
pobyso_differentiate(node *functionNode)
747 26 storres
{
748 26 storres
  /* Argument checking. */
749 26 storres
  node *differentialNode;
750 26 storres
  if (functionNode == NULL)
751 26 storres
  {
752 26 storres
    pobyso_error_message("pobyso_differentiate",
753 26 storres
                        "NULL_POINTER_ARGUMENT",
754 26 storres
                        "The function to differentiate is a NULL pointer");
755 26 storres
    return(NULL);
756 26 storres
  }
757 26 storres
  differentialNode = differentiate(functionNode);
758 26 storres
  if (differentialNode == NULL)
759 26 storres
  {
760 26 storres
    pobyso_error_message("pobyso_differentiate",
761 26 storres
                        "INTERNAL ERROR",
762 26 storres
                        "Sollya could not differentiate the function");
763 26 storres
  }
764 26 storres
  return(differentialNode);
765 26 storres
} // End pobyso_differentiate
766 26 storres

767 26 storres

768 26 storres
int
769 26 storres
pobyso_dirty_infnorm(mpfr_t infNorm,
770 26 storres
                      node *functionNode,
771 26 storres
                      mpfr_t lowerBound,
772 26 storres
                      mpfr_t upperBound,
773 26 storres
                      mp_prec_t precision)
774 26 storres
{
775 26 storres
  int functionCallResult;
776 26 storres
  /* Arguments checking. */
777 26 storres
  if (functionNode == NULL)
778 26 storres
  {
779 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
780 26 storres
                        "NULL_POINTER_ARGUMENT",
781 26 storres
                        "The function to compute is a NULL pointer");
782 26 storres
    return(1);
783 26 storres
  }
784 26 storres
  if (mpfr_cmp(lowerBound, upperBound) > 0)
785 26 storres
  {
786 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
787 26 storres
                        "INCOHERENT_INPUT_DATA",
788 26 storres
                        "The lower bond is greater than the upper bound");
789 26 storres
    return(1);
790 26 storres
  }
791 26 storres
  /* Particular cases. */
792 26 storres
  if (mpfr_cmp(lowerBound, upperBound) == 0)
793 26 storres
  {
794 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
795 26 storres
                                                  functionNode,
796 26 storres
                                                  lowerBound,
797 26 storres
                                                  precision);
798 26 storres
    return(functionCallResult);
799 26 storres
  }
800 26 storres
  if (isConstant(functionNode))
801 26 storres
  {
802 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
803 26 storres
                                                  functionNode,
804 26 storres
                                                  lowerBound,
805 26 storres
                                                  precision);
806 26 storres
    if (!functionCallResult)
807 26 storres
    {
808 26 storres
      mpfr_abs(infNorm, infNorm, MPFR_RNDN);
809 26 storres
    }
810 26 storres
    return(functionCallResult);
811 26 storres
  }
812 26 storres
  uncertifiedInfnorm(infNorm,
813 26 storres
                      functionNode,
814 26 storres
                      lowerBound,
815 26 storres
                      upperBound,
816 26 storres
                      POBYSO_DEFAULT_POINTS,
817 26 storres
                      precision);
818 26 storres
  return(0);
819 26 storres
} /* End pobyso_dirty_infnorm. */
820 26 storres

821 26 storres
int
822 26 storres
pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
823 26 storres
                          node *nodeToEvaluate,
824 26 storres
                          mpfr_t argument,
825 26 storres
                          mpfr_prec_t precision)
826 26 storres
{
827 26 storres
  /* Check input arguments. */
828 26 storres
  if (nodeToEvaluate == NULL)
829 26 storres
  {
830 26 storres
    pobyso_error_message("pobyso_evaluate_faithful",
831 26 storres
                        "NULL_POINTER_ARGUMENT",
832 26 storres
                        "nodeToEvaluate is a NULL pointer");
833 26 storres
    return(1);
834 26 storres
  }
835 26 storres
  evaluateFaithful(faithfulEvaluation,
836 26 storres
                   nodeToEvaluate,
837 26 storres
                   argument,
838 26 storres
                   precision);
839 26 storres
  return(0);
840 26 storres
} /* End pobyso_evaluate_faithfull. */
841 26 storres

842 26 storres
chain*
843 26 storres
pobyso_find_zeros(node *function,
844 26 storres
                  mpfr_t *lower_bound,
845 26 storres
                  mpfr_t *upper_bound)
846 26 storres
{
847 26 storres
  mp_prec_t currentPrecision;
848 26 storres
  mpfr_t currentDiameter;
849 26 storres
  rangetype bounds;
850 26 storres

851 26 storres
  currentPrecision = getToolPrecision();
852 26 storres
  mpfr_init2(currentDiameter, currentPrecision);
853 26 storres

854 26 storres
  bounds.a = lower_bound;
855 26 storres
  bounds.b = upper_bound;
856 26 storres

857 26 storres
  if (bounds.a == NULL || bounds.b == NULL)
858 26 storres
  {
859 26 storres
    pobyso_error_message("pobyso_find_zeros",
860 26 storres
                        "MEMORY_ALLOCATION_ERROR",
861 26 storres
                        "Could not allocate one of the bounds");
862 26 storres
    return(NULL);
863 26 storres
  }
864 26 storres
  return(findZerosFunction(function,
865 26 storres
                            bounds,
866 26 storres
                            currentPrecision,
867 26 storres
                            currentDiameter));
868 26 storres
} /* End pobyso_find_zeros. */
869 26 storres

870 26 storres
void
871 26 storres
pobyso_free_chain_of_nodes(chain *theChainOfNodes)
872 26 storres
{
873 26 storres
  node *currentNode           = NULL;
874 26 storres
  chain *currentChainElement  = NULL;
875 26 storres
  chain *nextChainElement     = NULL;
876 26 storres

877 26 storres
  nextChainElement = theChainOfNodes;
878 26 storres

879 26 storres
  while(nextChainElement != NULL)
880 26 storres
  {
881 26 storres
    currentChainElement = nextChainElement;
882 26 storres
    currentNode = (node*)(currentChainElement->value);
883 26 storres
    nextChainElement = nextChainElement->next;
884 26 storres
    free_memory(currentNode);
885 26 storres
    free((void*)(currentChainElement));
886 26 storres
  }
887 26 storres
} /* End pobyso_free_chain_of_nodes. */
888 26 storres

889 26 storres
void
890 26 storres
pobyso_free_range(rangetype range)
891 26 storres
{
892 26 storres

893 26 storres
  mpfr_clear(*(range.a));
894 26 storres
  mpfr_clear(*(range.b));
895 26 storres
  free(range.a);
896 26 storres
  free(range.b);
897 26 storres
} /* End pobyso_free_range. */
898 26 storres

899 26 storres
node*
900 26 storres
pobyso_fp_minimax_canonical_monomials_base(node *function,
901 26 storres
                                      int degree,
902 26 storres
                                      chain *formats,
903 26 storres
                                      chain *points,
904 26 storres
                                      mpfr_t lowerBound,
905 26 storres
                                      mpfr_t upperBound,
906 26 storres
                                      int fpFixedArg,
907 26 storres
                                      int absRel,
908 26 storres
                                      node *constPart,
909 26 storres
                                      node *minimax)
910 26 storres
{
911 26 storres
  return(NULL);
912 26 storres
} /* End pobyso_fp_minimax_canonical_monomials_base. */
913 26 storres

914 26 storres
node*
915 26 storres
pobyso_parse_function(char *functionString,
916 26 storres
                      char *freeVariableNameString)
917 26 storres
{
918 26 storres
  if (functionString == NULL || freeVariableNameString == NULL)
919 26 storres
  {
920 26 storres
    pobyso_error_message("pobyso_parse_function",
921 26 storres
                        "NULL_POINTER_ARGUMENT",
922 26 storres
                        "One of the arguments is a NULL pointer");
923 26 storres
    return(NULL);
924 26 storres
  }
925 26 storres
  return(parseString(functionString));
926 26 storres

927 26 storres
} /* End pobyso_parse_function */
928 26 storres

929 26 storres
node*
930 26 storres
pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
931 26 storres
                                                  unsigned int mode,
932 26 storres
                                                  mpfr_t lowerBound,
933 26 storres
                                                  mpfr_t upperBound,
934 26 storres
                                                  mpfr_t eps)
935 26 storres
{
936 26 storres
  node *weight              = NULL;
937 26 storres
  node *bestApproxPolyNode  = NULL;
938 26 storres
  node *bestApproxHorner    = NULL;
939 26 storres
  node *errorNode           = NULL;
940 26 storres
  rangetype degreeRange;
941 26 storres
  mpfr_t quality;
942 26 storres
  mpfr_t currentError;
943 26 storres
  unsigned int degree;
944 26 storres

945 26 storres
  /* Check the parameters. */
946 26 storres
  if (functionNode == NULL)
947 26 storres
  {
948 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
949 26 storres
                        "NULL_POINTER_ARGUMENT",
950 26 storres
                        "functionNode is a NULL pointer");
951 26 storres
    return(NULL);
952 26 storres
  }
953 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
954 26 storres
  {
955 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
956 26 storres
                        "INCOHERENT_INPUT_DATA",
957 26 storres
                        "the lower_bound >= upper_bound");
958 26 storres
    return(NULL);
959 26 storres
  }
960 26 storres
  /* Set the weight. */
961 26 storres
  if (mode == POBYSO_ABSOLUTE)
962 26 storres
  {
963 26 storres
    /* Set the weight to 1 for the ABSOLUTE_MODE. */
964 26 storres
    weight = makeConstantDouble(1.0);
965 26 storres
  }
966 26 storres
  else
967 26 storres
  {
968 26 storres
    if (mode == POBYSO_RELATIVE)
969 26 storres
    {
970 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
971 26 storres
                          "NOT_IMPLEMENTED",
972 26 storres
                          "the search for relative error is not implemented yet");
973 26 storres
      return(NULL);
974 26 storres
    }
975 26 storres
    else
976 26 storres
    {
977 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
978 26 storres
                          "INCOHERENT_INPUT_DATA",
979 26 storres
                          "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
980 26 storres
      return(NULL);
981 26 storres
    }
982 26 storres
  }
983 26 storres
  //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
984 26 storres
  degreeRange = guessDegree(functionNode,
985 26 storres
                            weight,
986 26 storres
                            lowerBound,
987 26 storres
                            upperBound,
988 26 storres
                            eps,
989 26 storres
                            POBYSO_GUESS_DEGREE_BOUND);
990 26 storres
  degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
991 26 storres
  //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
992 26 storres
  fprintf(stderr, "Guessed degree: ");
993 26 storres
  mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
994 26 storres
  fprintf(stderr, " - as int: %u\n", degree);
995 26 storres
  /* Reduce the degree by 1 in the foolish hope it could work. */
996 26 storres
  if (degree > 0) degree--;
997 26 storres
  /* Both elements of degreeRange where "inited" within guessDegree. */
998 26 storres
  mpfr_clear(*(degreeRange.a));
999 26 storres
  mpfr_clear(*(degreeRange.b));
1000 26 storres
  free(degreeRange.a);
1001 26 storres
  free(degreeRange.b);
1002 26 storres
  /* Mimic the default behavior of interactive Sollya. */
1003 26 storres
  mpfr_init(quality);
1004 26 storres
  mpfr_set_d(quality, 1e-5, MPFR_RNDN);
1005 26 storres
  mpfr_init2(currentError, getToolPrecision());
1006 26 storres
  mpfr_set_inf(currentError,1);
1007 26 storres

1008 26 storres
  /* Try to refine the initial guess: loop with increasing degrees until
1009 26 storres
   * we find a satisfactory one. */
1010 26 storres
  while(mpfr_cmp(currentError, eps) > 0)
1011 26 storres
  {
1012 26 storres
    /* Get rid of the previous polynomial, if any. */
1013 26 storres
    if (bestApproxPolyNode != NULL)
1014 26 storres
    {
1015 26 storres
      free_memory(bestApproxPolyNode);
1016 26 storres
    }
1017 26 storres
    fprintf(stderr, "Degree: %u\n", degree);
1018 26 storres
    fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
1019 26 storres
    /* Try to find a polynomial with the guessed degree. */
1020 26 storres
    bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
1021 26 storres
                                                            weight,
1022 26 storres
                                                            degree,
1023 26 storres
                                                            lowerBound,
1024 26 storres
                                                            upperBound,
1025 26 storres
                                                            quality);
1026 26 storres

1027 26 storres
    if (bestApproxPolyNode == NULL)
1028 26 storres
    {
1029 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
1030 26 storres
                          "INTERNAL_ERROR",
1031 26 storres
                          "could not compute the bestApproxPolyNode");
1032 26 storres
      mpfr_clear(currentError);
1033 26 storres
      mpfr_clear(quality);
1034 26 storres
      return(NULL);
1035 26 storres
    }
1036 26 storres

1037 26 storres
    setDisplayMode(DISPLAY_MODE_DECIMAL);
1038 26 storres
    fprintTree(stderr, bestApproxPolyNode);
1039 26 storres
    fprintf(stderr, "\n\n");
1040 26 storres

1041 26 storres
    errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
1042 26 storres
    /* Check the error with the computed polynomial. */
1043 26 storres
    uncertifiedInfnorm(currentError,
1044 26 storres
                        errorNode,
1045 26 storres
                        lowerBound,
1046 26 storres
                        upperBound,
1047 26 storres
                        POBYSO_INF_NORM_NUM_POINTS,
1048 26 storres
                        getToolPrecision());
1049 26 storres
    fprintf(stderr, "Inf norm: ");
1050 26 storres
    mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
1051 26 storres
    fprintf(stderr, "\n\n");
1052 26 storres
    /* Free the errorNode but not the bestApproxPolyNode (we need it if
1053 26 storres
     * we exit the loop at the next iteration). */
1054 26 storres
    free_memory(errorNode);
1055 26 storres
    degree++;
1056 26 storres
  }
1057 26 storres
  /* Use an intermediate variable, since horner() creates a new node
1058 26 storres
   * and does not reorder the argument "in place". This allows for the memory
1059 26 storres
   * reclaim of bestApproxHorner.
1060 26 storres
   */
1061 26 storres
  bestApproxHorner = horner(bestApproxPolyNode);
1062 26 storres
  free_memory(bestApproxPolyNode);
1063 26 storres
  mpfr_clear(currentError);
1064 26 storres
  mpfr_clear(quality);
1065 26 storres
  free_memory(weight);
1066 26 storres
  return(bestApproxHorner);
1067 26 storres
} /* End pobyso_remez_approx_canonical_monomials_base_for_error */
1068 26 storres

1069 26 storres
node*
1070 26 storres
pobyso_remez_canonical_monomials_base(node *function,
1071 26 storres
                                     node *weight,
1072 26 storres
                                     unsigned int degree,
1073 26 storres
                                     mpfr_t lowerBound,
1074 26 storres
                                     mpfr_t upperBound,
1075 26 storres
                                     mpfr_t quality)
1076 26 storres
{
1077 26 storres
  node  *bestApproxPoly = NULL;
1078 26 storres
  chain *monomials      = NULL;
1079 26 storres
  chain *curMonomial    = NULL;
1080 26 storres

1081 26 storres
  mpfr_t satisfying_error;
1082 26 storres
  mpfr_t target_error;
1083 26 storres

1084 26 storres
  /* Argument checking */
1085 26 storres
  /* Function tree. */
1086 26 storres
  if (function == NULL)
1087 26 storres
  {
1088 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1089 26 storres
                        "NULL_POINTER_ARGUMENT",
1090 26 storres
                        "the \"function\" argument is a NULL pointer");
1091 26 storres
    return(NULL);
1092 26 storres
  }
1093 26 storres
  if (weight == NULL)
1094 26 storres
  {
1095 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1096 26 storres
                        "NULL_POINTER_ARGUMENT",
1097 26 storres
                        "the \"weight\" argument is a NULL pointer");
1098 26 storres
    return(NULL);
1099 26 storres
  }
1100 26 storres
  /* Check the bounds. */
1101 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
1102 26 storres
  {
1103 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1104 26 storres
                        "INCOHERENT_IMPUT_DATA",
1105 26 storres
                        "the lower_bound >= upper_bound");
1106 26 storres
    return(NULL);
1107 26 storres
  }
1108 26 storres
  /* The quality must be a non null positive number. */
1109 26 storres
  if (mpfr_sgn(quality) <= 0)
1110 26 storres
  {
1111 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1112 26 storres
                        "INCOHERENT_INPUT_DATA",
1113 26 storres
                        "the quality <= 0");
1114 26 storres
  }
1115 26 storres
  /* End argument checking. */
1116 26 storres
  /* Create the monomials nodes chains. */
1117 26 storres
  monomials = pobyso_create_canonical_monomials_base(degree);
1118 26 storres
  fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
1119 26 storres
  if (monomials == NULL || (lengthChain(monomials) != degree + 1))
1120 26 storres
  {
1121 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1122 26 storres
                        "CHAIN_CREATION_ERROR",
1123 26 storres
                        "could not create the monomials chain");
1124 26 storres
    return(NULL);
1125 26 storres
  }
1126 26 storres
  curMonomial = monomials;
1127 26 storres

1128 26 storres
  while (curMonomial != NULL)
1129 26 storres
  {
1130 26 storres
    fprintf(stderr, "monomial tree: ");
1131 26 storres
    //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
1132 26 storres
    fprintTree(stderr, (node*)(curMonomial->value));
1133 26 storres
    fprintf(stderr, "\n");
1134 26 storres
    curMonomial = curMonomial->next;
1135 26 storres
  }
1136 26 storres

1137 26 storres
  /* Deal with NULL weight. */
1138 26 storres
  if (weight == NULL)
1139 26 storres
  {
1140 26 storres
    weight = makeConstantDouble(1.0);
1141 26 storres
  }
1142 26 storres
  /* Compute the best polynomial with the required quality.
1143 26 storres
     The behavior is as if satisfying error and target_error had
1144 26 storres
     not been used.*/
1145 26 storres
  mpfr_init(satisfying_error);
1146 26 storres
  mpfr_init(target_error);
1147 26 storres
  mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
1148 26 storres
  mpfr_set_inf(target_error, 1);
1149 26 storres

1150 26 storres

1151 26 storres
  fprintf(stderr, "satisfying_error: ");
1152 26 storres
  mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
1153 26 storres
  fprintf(stderr, ".\n");
1154 26 storres
  fprintf(stderr, "target_error: ");
1155 26 storres
  mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
1156 26 storres
  fprintf(stderr, ".\n");
1157 26 storres

1158 26 storres
  fprintf(stderr,
1159 26 storres
          "current precision: %li\n", getToolPrecision());
1160 26 storres
  /* Call the Sollya function. */
1161 26 storres
  bestApproxPoly = remez(function,
1162 26 storres
                          weight,
1163 26 storres
                          monomials,
1164 26 storres
                          lowerBound,
1165 26 storres
                          upperBound,
1166 26 storres
                          quality,
1167 26 storres
                          satisfying_error,
1168 26 storres
                          target_error,
1169 26 storres
                          getToolPrecision());
1170 26 storres

1171 26 storres
  mpfr_clear(satisfying_error);
1172 26 storres
  mpfr_clear(target_error);
1173 26 storres
  pobyso_free_chain_of_nodes(monomials);
1174 26 storres

1175 26 storres
  return(bestApproxPoly);
1176 26 storres
} /* End pobyso_remez_canonical_monomials_base. */
1177 26 storres

1178 26 storres
#endif
1179 26 storres
1180 26 storres
void
1181 26 storres
pobyso_error_message(char *functionName, char *messageName, char* message)
1182 26 storres
{
1183 26 storres
  fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
1184 26 storres
} /* End pobyso_error_message */