Statistiques
| Révision :

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

Historique | Voir | Annoter | Télécharger (44,01 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 147 storres
/* @see pobyso.h#pobyso_dirty_find_zeros */
43 147 storres
mpfr_t*
44 147 storres
pobyso_dirty_find_zeros_bounds(pobyso_func_exp_t funcExpSo,
45 147 storres
                              mpfr_t lowerBound,
46 147 storres
                              mpfr_t upperBound,
47 147 storres
                              int* zerosCount)
48 147 storres
{
49 147 storres
  pobyso_range_t rangeSo;
50 147 storres
  sollya_obj_t zerosListSo    = NULL;
51 147 storres
  sollya_obj_t* zerosArraySo  = NULL;
52 147 storres
  mpfr_t* zerosArrayMp        = NULL;
53 147 storres
  pobyso_precision_t prec;
54 147 storres
55 147 storres
  int endEll = 0;
56 147 storres
  int i,j;
57 147 storres
58 147 storres
  /* Arguments check. */
59 147 storres
  *zerosCount = -1;
60 147 storres
  if (funcExpSo == NULL )
61 147 storres
  {
62 147 storres
    pobyso_error_message("pobyso_dirty_find_zeros",
63 147 storres
                        "NULL_POINTER_ARGUMENT",
64 147 storres
                        "The funcExpSo argument is a NULL pointer");
65 147 storres
    return NULL;
66 147 storres
  }
67 147 storres
  if (mpfr_cmp(lowerBound, upperBound) > 0)
68 147 storres
  {
69 147 storres
    pobyso_error_message("pobyso_dirty_find_zeros",
70 147 storres
                        "INVALID_INTERVAL_ARGUMENT",
71 147 storres
                        "The lower bound is larger than the upper bound");
72 147 storres
    return NULL;
73 147 storres
  }
74 147 storres
  /* Make a range out of the bounds. */
75 147 storres
  rangeSo = sollya_lib_range_from_bounds(lowerBound, upperBound);
76 147 storres
  if (rangeSo == NULL)
77 147 storres
  {
78 147 storres
    return NULL;
79 147 storres
  }
80 147 storres
  zerosListSo = sollya_lib_dirtyfindzeros(funcExpSo, rangeSo);
81 147 storres
  if (zerosListSo == NULL)
82 147 storres
  {
83 147 storres
    return NULL;
84 147 storres
  }
85 147 storres
  sollya_lib_clear_obj(rangeSo);
86 147 storres
  /* Transform the Sollya list into an MPFR list. */
87 147 storres
  if (! sollya_lib_get_list_elements(&zerosArraySo,
88 147 storres
                                      zerosCount,
89 147 storres
                                      &endEll,
90 147 storres
                                      zerosListSo))
91 147 storres
  {
92 147 storres
    sollya_lib_clear_obj(zerosListSo);
93 147 storres
    *zerosCount = -1;
94 147 storres
    return NULL;
95 147 storres
  }
96 147 storres
  sollya_lib_clear_obj(zerosListSo);
97 147 storres
  zerosArrayMp = (mpfr_t*) malloc(*zerosCount * sizeof(mpfr_t));
98 147 storres
  if (zerosArrayMp == NULL)
99 147 storres
  {
100 147 storres
    pobyso_error_message("pobyso_dirty_find_zeros",
101 147 storres
                          "MEMORY_ALLOCATION_ERROR",
102 147 storres
                          "Could not allocate zeroes array");
103 147 storres
    *zerosCount = -1;
104 147 storres
    return NULL;
105 147 storres
  }
106 147 storres
  for (i = 0 ; i < *zerosCount ; i++)
107 147 storres
  {
108 147 storres
    if (! sollya_lib_get_prec_of_constant(&prec, zerosArraySo[i]))
109 147 storres
    {
110 147 storres
      /* Clean up the already allocated MPFRs. */
111 147 storres
      for (j = 0 ; j < i ; j++) mpfr_clear(zerosArrayMp[j]);
112 147 storres
      /* Clean up the zerosArrayMp array itself. */
113 147 storres
      free(zerosArrayMp);
114 147 storres
      /* Clean up what is left in the zerosArraySo. */
115 147 storres
      for (j = i ; j < *zerosCount ; j++) sollya_lib_clear_obj(zerosArraySo[j]);
116 147 storres
      /* Clean up the zerosArraySo array itself. */
117 147 storres
      sollya_lib_free(zerosArraySo);
118 147 storres
      *zerosCount = -1;
119 147 storres
      return NULL;
120 147 storres
    }
121 147 storres
    mpfr_init2(zerosArrayMp[i], prec);
122 147 storres
    if (! sollya_lib_get_constant(zerosArrayMp[i], zerosArraySo[i]))
123 147 storres
    {
124 147 storres
      /* Clean up the already allocated MPFRs. */
125 147 storres
      for (j = 0 ; j <= i ; j++) mpfr_clear(zerosArrayMp[j]);
126 147 storres
      /* Clean up the zerosArrayMp array itself. */
127 147 storres
      free(zerosArrayMp);
128 147 storres
      /* Clean up what is left in the zerosArraySo. */
129 147 storres
      for (j = i ; j < *zerosCount ; j++) sollya_lib_clear_obj(zerosArraySo[j]);
130 147 storres
      /* Clean up the zerosArraySo array itself. */
131 147 storres
      sollya_lib_free(zerosArraySo);
132 147 storres
      *zerosCount = -1;
133 147 storres
      return NULL;
134 147 storres
    }
135 147 storres
    sollya_lib_clear_obj(zerosArraySo[i]);
136 147 storres
  } /* End for i. */
137 147 storres
  sollya_lib_free(zerosArraySo);
138 147 storres
  return zerosArrayMp;
139 147 storres
} /* End pobyso_dirty_find_zeros. */
140 147 storres
141 140 storres
/* @see pobyso.h#pobyso_evaluate_constant */
142 139 storres
int
143 139 storres
pobyso_evaluate_constant(pobyso_func_exp_t functionSo,
144 139 storres
                          mpfr_t argumentMp,
145 139 storres
                          mpfr_t evaluationMp)
146 139 storres
{
147 139 storres
  sollya_obj_t argumentSo   = NULL;
148 139 storres
  sollya_obj_t evaluationSo = NULL;
149 140 storres
  mpfr_t evaluationTmpMp1;
150 140 storres
  mpfr_t evaluationTmpMp2;
151 140 storres
  mpfr_prec_t evalPrec      = 0;
152 139 storres
153 139 storres
  /* Test argument. */
154 139 storres
  if (functionSo == NULL || argumentMp == NULL || evaluationMp == NULL)
155 139 storres
  {
156 139 storres
    pobyso_error_message("pobyso_evaluate_constant",
157 139 storres
                        "NULL_POINTER_ARGUMENT",
158 139 storres
                        "One of the arguments is a NULL pointer");
159 139 storres
    return 1;
160 139 storres
  }
161 139 storres
  if (! sollya_lib_obj_is_function(functionSo))
162 139 storres
  {
163 139 storres
    pobyso_error_message("pobyso_evaluate_constant",
164 139 storres
                        "INVALID_TYPE_ARGUMENT",
165 139 storres
                     "The functionSo argument is not a functional expression");
166 139 storres
    return 1;
167 139 storres
  }
168 140 storres
  /* Function evaluation and checks. */
169 139 storres
  argumentSo    = sollya_lib_constant(argumentMp);
170 139 storres
  evaluationSo  = sollya_lib_evaluate(functionSo, argumentSo);
171 147 storres
  /* Not needed any more. */
172 147 storres
  //pobyso_autoprint(evaluationSo);
173 139 storres
  sollya_lib_clear_obj(argumentSo);
174 140 storres
  /* The range case: we return the mean of the bounds. The result
175 140 storres
   * is not faithfully rounded. */
176 140 storres
  if (sollya_lib_obj_is_range(evaluationSo))
177 140 storres
  {
178 147 storres
    //pobyso_autoprint(evaluationSo);
179 140 storres
    if (sollya_lib_get_prec_of_range(&evalPrec, evaluationSo))
180 140 storres
    {
181 140 storres
      mpfr_init2(evaluationTmpMp1, evalPrec);
182 140 storres
      mpfr_init2(evaluationTmpMp2, evalPrec);
183 140 storres
      if (sollya_lib_get_bounds_from_range(evaluationTmpMp1,
184 140 storres
                                       evaluationTmpMp2,
185 140 storres
                                       evaluationSo))
186 140 storres
      {
187 140 storres
        /* Compute the mean of the bounds. */
188 140 storres
        mpfr_clear(evaluationMp);
189 140 storres
        mpfr_init2(evaluationMp, evalPrec);
190 140 storres
        mpfr_add(evaluationMp,evaluationTmpMp1,evaluationTmpMp2,MPFR_RNDN);
191 140 storres
        mpfr_div_2ui(evaluationMp, evaluationMp, 1, MPFR_RNDN);
192 140 storres
        mpfr_clear(evaluationTmpMp1);
193 140 storres
        mpfr_clear(evaluationTmpMp2);
194 140 storres
        sollya_lib_clear_obj(evaluationSo);
195 147 storres
        /* It may happen, in this case, when the bounds are -Infty  and
196 147 storres
         * +Infty, that the average is NaN. */
197 147 storres
        if (mpfr_nan_p(evaluationMp))
198 147 storres
        {
199 147 storres
          return POBYSO_NAN;
200 147 storres
        }
201 147 storres
        else
202 147 storres
        {
203 147 storres
          return POBYSO_UNFAITHFUL;
204 147 storres
        }
205 140 storres
      }
206 140 storres
      else /* Could not get the values of the bounds. */
207 140 storres
      {
208 140 storres
        sollya_lib_clear_obj(evaluationSo);
209 140 storres
        return 1;
210 140 storres
      }
211 140 storres
    }
212 140 storres
    else /* Could not get the precision of the range. */
213 140 storres
    {
214 140 storres
      sollya_lib_clear_obj(evaluationSo);
215 140 storres
      return 1;
216 140 storres
    }
217 140 storres
  } /* End the evaluation is a range. */
218 140 storres
  /* From now on, we assume that the evaluation is constant. */
219 147 storres
  if (sollya_lib_get_prec_of_constant(&evalPrec, evaluationSo))
220 140 storres
  {
221 140 storres
    mpfr_init2(evaluationTmpMp1, evalPrec);
222 140 storres
    if (sollya_lib_get_constant(evaluationTmpMp1, evaluationSo))
223 140 storres
    {
224 140 storres
      /* Deal with the NaN case. */
225 140 storres
      if (mpfr_nan_p(evaluationTmpMp1))
226 140 storres
      {
227 140 storres
        mpfr_clear(evaluationTmpMp1);
228 140 storres
        sollya_lib_clear_obj(evaluationSo);
229 140 storres
        return POBYSO_NAN;
230 140 storres
      }
231 140 storres
      else /* The evaluation is not NaN. */
232 140 storres
      {
233 140 storres
        mpfr_clear(evaluationMp);
234 140 storres
        mpfr_init2(evaluationMp, evalPrec);
235 140 storres
        mpfr_set(evaluationMp, evaluationTmpMp1, MPFR_RNDN);
236 140 storres
        mpfr_clear(evaluationTmpMp1);
237 140 storres
        sollya_lib_clear_obj(evaluationSo);
238 140 storres
        return 0;
239 140 storres
      }
240 140 storres
    }
241 140 storres
    else /* Could not recover the value of the evaluation. */
242 140 storres
    {
243 140 storres
      mpfr_clear(evaluationTmpMp1);
244 140 storres
      sollya_lib_clear_obj(evaluationSo);
245 140 storres
      return 1;
246 140 storres
    }
247 140 storres
  }
248 140 storres
  else /* Could not get the precision of the evaluation. */
249 140 storres
  {
250 140 storres
    sollya_lib_clear_obj(evaluationSo);
251 140 storres
    return 0;
252 140 storres
  }
253 139 storres
}
254 139 storres
/* End pobyso_evaluate_constant. */
255 139 storres
256 134 storres
/* @see pobyso.h#pobyso_get_verbosity */
257 134 storres
int
258 134 storres
pobyso_get_verbosity()
259 134 storres
{
260 134 storres
  sollya_obj_t verbositySo = NULL;
261 134 storres
  int verbosity            = 0;
262 134 storres
263 134 storres
  verbositySo = sollya_lib_get_verbosity();
264 134 storres
  sollya_lib_get_constant_as_int(&verbosity,verbositySo);
265 134 storres
  sollya_lib_clear_obj(verbositySo);
266 134 storres
  return verbosity;
267 134 storres
} /* End pobyso_get_verbosity. */
268 134 storres
269 137 storres
/** @see pobyso.h#pobyso_is_constant_expression
270 137 storres
 * Strategy: rely on sollya_lib_get_constant. It return 1, when the
271 137 storres
 * expression is constant.
272 137 storres
 */
273 133 storres
int
274 133 storres
pobyso_is_constant_expression(sollya_obj_t obj_to_test)
275 133 storres
{
276 133 storres
  mpfr_t dummy;
277 133 storres
  int test;
278 134 storres
  int initial_verbosity_level      = 0;
279 134 storres
280 133 storres
  /* Test argument. */
281 133 storres
  if (obj_to_test == NULL)
282 133 storres
  {
283 137 storres
    pobyso_error_message("pobyso_is_constant_expression",
284 133 storres
                        "NULL_POINTER_ARGUMENT",
285 133 storres
                        "The expression is a NULL pointer");
286 133 storres
    return 0;
287 133 storres
  }
288 134 storres
  initial_verbosity_level = pobyso_set_verbosity_off();
289 139 storres
  /* In Sollya, constant are functional expressions. */
290 133 storres
  if (! sollya_lib_obj_is_function(obj_to_test))
291 133 storres
  {
292 134 storres
    pobyso_set_verbosity_to(initial_verbosity_level);
293 133 storres
    return 0;
294 133 storres
  }
295 133 storres
  mpfr_init2(dummy,64);
296 133 storres
  /* Call to previous Sollya function resets verbosity. */
297 137 storres
  /* Todo: change verbosity suppression strategy with a message call back. */
298 134 storres
  pobyso_set_verbosity_off();
299 139 storres
  /* Try to convert the would be constant into an MPFR number. */
300 139 storres
  /* If OK, we indeed have a constant. If the conversion fails, return 0. */
301 133 storres
  test = sollya_lib_get_constant(dummy, obj_to_test);
302 134 storres
  pobyso_set_verbosity_to(initial_verbosity_level);
303 138 storres
  if (test)
304 137 storres
  {
305 147 storres
    if(mpfr_number_p(dummy) || mpfr_inf_p(dummy))
306 138 storres
    {
307 138 storres
      mpfr_clear(dummy);
308 138 storres
      return test;
309 138 storres
    }
310 147 storres
    else /* We do not consider NaNs as constants. */
311 138 storres
    {
312 138 storres
      mpfr_clear(dummy);
313 138 storres
      return 0;
314 138 storres
    }
315 137 storres
  }
316 138 storres
  else
317 137 storres
  {
318 139 storres
    mpfr_clear(dummy);
319 137 storres
    return 0;
320 137 storres
  }
321 138 storres
} /* End pobyso_is_constant_expression. */
322 137 storres
323 137 storres
/** @see pobyso.h#pobyso_is_monomial. */
324 137 storres
int
325 137 storres
pobyso_is_int(pobyso_func_exp_t exprSo)
326 137 storres
{
327 137 storres
  mpfr_t float1M;
328 137 storres
  mpfr_t float2M;
329 137 storres
  mpfr_t tempFloat1M;
330 137 storres
  mpfr_t tempFloat2M;
331 137 storres
  mpfr_prec_t prec;
332 137 storres
  int64_t asInt;
333 137 storres
  sollya_obj_t newConstantSo = NULL;
334 137 storres
  /* Arguments check. */
335 137 storres
  if (exprSo == NULL)
336 137 storres
  {
337 137 storres
    pobyso_error_message("pobyso_is_free_var_posze_int_power",
338 137 storres
                        "NULL_POINTER_ARGUMENT",
339 137 storres
                        "The expression is a NULL pointer");
340 137 storres
    return 0;
341 137 storres
  }
342 137 storres
  //fprintf(stdout, "Not NULL.\n"); pobyso_autoprint(exprSo);
343 137 storres
  if (! pobyso_is_constant_expression(exprSo)) return 0;
344 137 storres
  if (! sollya_lib_get_constant_as_int64(&asInt, exprSo)) return 0;
345 137 storres
  if (asInt == INT64_MIN || asInt == INT64_MAX) return 0;
346 137 storres
  /* Some constant integer expression can't have their precision computed
347 137 storres
   * (e.g. cos(pi). */
348 137 storres
  if (! sollya_lib_get_prec_of_constant(&prec, exprSo))
349 137 storres
  {
350 137 storres
    mpfr_init2(tempFloat1M, 165);
351 137 storres
    mpfr_init2(tempFloat2M, 165);
352 137 storres
    mpfr_abs(tempFloat1M, tempFloat1M, MPFR_RNDN);
353 137 storres
    mpfr_log2(tempFloat2M, tempFloat1M, MPFR_RNDU);
354 137 storres
    mpfr_rint_ceil(tempFloat1M, tempFloat2M, MPFR_RNDU);
355 137 storres
    prec = mpfr_get_si(tempFloat1M, MPFR_RNDN) + 10;
356 137 storres
    if (prec < 1024) prec = 1024;
357 137 storres
    mpfr_clear(tempFloat1M);
358 137 storres
    mpfr_clear(tempFloat2M);
359 137 storres
    mpfr_init2(float1M, prec);
360 137 storres
    if (!sollya_lib_get_constant(float1M, exprSo))
361 137 storres
    {
362 137 storres
      mpfr_clear(float1M);
363 137 storres
      return 0;
364 137 storres
    }
365 137 storres
  }
366 137 storres
  else /* Precision could be given. */
367 137 storres
  {
368 137 storres
    mpfr_init2(float1M, prec);
369 137 storres
    if (! sollya_lib_get_constant(float1M, exprSo))
370 137 storres
    {
371 137 storres
      mpfr_clear(float1M);
372 137 storres
      return 0;
373 137 storres
    }
374 137 storres
  }
375 137 storres
  if (mpfr_nan_p(float1M) || mpfr_inf_p(float1M))
376 137 storres
  {
377 137 storres
    mpfr_clear(float1M);
378 137 storres
    return 0;
379 137 storres
  }
380 137 storres
  if ((newConstantSo = sollya_lib_constant_from_int64(asInt)) == NULL)
381 137 storres
  {
382 137 storres
    mpfr_clear(float1M);
383 137 storres
    return 0;
384 137 storres
  }
385 137 storres
  sollya_lib_get_prec_of_constant(&prec, newConstantSo);
386 137 storres
  mpfr_init2(float2M, prec);
387 137 storres
  sollya_lib_get_constant(float2M, newConstantSo);
388 137 storres
  if (mpfr_cmp(float1M, float2M) == 0)
389 137 storres
  {
390 137 storres
    mpfr_clear(float1M);
391 137 storres
    mpfr_clear(float2M);
392 137 storres
    sollya_lib_clear_obj(newConstantSo);
393 137 storres
    return 1;
394 137 storres
  }
395 137 storres
  else
396 137 storres
  {
397 137 storres
    mpfr_clear(float1M);
398 137 storres
    mpfr_clear(float2M);
399 137 storres
    sollya_lib_clear_obj(newConstantSo);
400 137 storres
    return 0;
401 137 storres
  }
402 137 storres
} /* End pobyso_is_int. */
403 137 storres
404 137 storres
/** @see pobyso.h#pobyso_is_monomial.
405 137 storres
 * Strategy: check that the object is a functional expression and
406 137 storres
 * if so check that it is made of cte * free_var ^ some_power where :
407 137 storres
 * - cte is a constant expression (a per pobyso_is_constant_experession;
408 137 storres
 * - some_power is a positive or null power. t*/
409 137 storres
int
410 137 storres
pobyso_is_monomial(sollya_obj_t objSo)
411 137 storres
{
412 138 storres
  int arity;
413 138 storres
  sollya_obj_t subFun1So = NULL;
414 138 storres
  sollya_obj_t subFun2So = NULL;
415 138 storres
  sollya_obj_t subFun3So = NULL;
416 138 storres
  sollya_base_function_t head = 0;
417 138 storres
  long int exponent = 0;
418 138 storres
  long int exprIntValue = 0;
419 138 storres
420 138 storres
  /* Arguments check. */
421 138 storres
  if (objSo == NULL)
422 138 storres
  {
423 138 storres
    pobyso_error_message("pobyso_is_monomial",
424 138 storres
                        "NULL_POINTER_ARGUMENT",
425 138 storres
                        "The expression is a NULL pointer");
426 138 storres
    return 0;
427 138 storres
  }
428 138 storres
  /* The object must be a function. */
429 138 storres
  if (! sollya_lib_obj_is_function(objSo)) return 0;
430 138 storres
  /* Check if it is the 1 constant. */
431 138 storres
  if (pobyso_is_int(objSo))
432 138 storres
  {
433 138 storres
    if (! sollya_lib_get_constant_as_int64(&exprIntValue, objSo))
434 138 storres
    {
435 138 storres
      return 0;
436 138 storres
    }
437 138 storres
    else
438 138 storres
    {
439 138 storres
      if (exprIntValue == 1) return 1;
440 138 storres
      else return 0;
441 138 storres
    }
442 138 storres
  }
443 138 storres
  if (! sollya_lib_decompose_function(objSo,
444 138 storres
                                      &head,
445 138 storres
                                      &arity,
446 138 storres
                                      &subFun1So,
447 138 storres
                                      &subFun2So,
448 138 storres
                                      NULL)) return 0;
449 138 storres
  if (arity > 2)
450 138 storres
  {
451 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
452 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
453 138 storres
    return 0;
454 138 storres
  }
455 138 storres
  /* Arity == 1 must be free variable by itself. */
456 138 storres
  if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE)
457 138 storres
  {
458 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
459 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
460 138 storres
    return 1;
461 138 storres
  }
462 138 storres
  else
463 138 storres
  {
464 138 storres
    /* Another expression. Must be: free variable  ^ poze integer. */
465 138 storres
    if (arity == 2)
466 138 storres
    {
467 138 storres
      if (head != SOLLYA_BASE_FUNC_POW)
468 138 storres
      {
469 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
470 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
471 138 storres
        return 0;
472 138 storres
      }
473 138 storres
      if (! pobyso_is_int(subFun2So))
474 138 storres
      {
475 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
476 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
477 138 storres
        return 0;
478 138 storres
      }
479 138 storres
      if (! sollya_lib_get_constant_as_int64(&exponent, subFun2So))
480 138 storres
      {
481 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
482 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
483 138 storres
        return 0;
484 138 storres
      }
485 138 storres
      if (exponent < 0)
486 138 storres
      {
487 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
488 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
489 138 storres
        return 0;
490 138 storres
      }
491 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
492 138 storres
      /* Check that the first subfunction is the free variable. */
493 138 storres
      if (! sollya_lib_decompose_function(subFun1So,
494 138 storres
                                          &head,
495 138 storres
                                          &arity,
496 138 storres
                                          &subFun2So,
497 138 storres
                                          &subFun3So,
498 138 storres
                                          NULL))
499 138 storres
      {
500 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
501 138 storres
        return 0;
502 138 storres
      }
503 138 storres
      if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE)
504 138 storres
      {
505 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
506 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
507 138 storres
        if (subFun3So != NULL) sollya_lib_clear_obj(subFun3So);
508 138 storres
        return 1;
509 138 storres
      }
510 138 storres
      else
511 138 storres
      {
512 138 storres
        if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
513 138 storres
        if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
514 138 storres
        return 0;
515 138 storres
      }
516 138 storres
    } /* End if arity == 2. */
517 138 storres
  } /* End else if arity == 1. */
518 137 storres
  return 0;
519 137 storres
} /* End pobyso_is_monomial. */
520 137 storres
521 138 storres
/** @see pobyso.h#pobyso_is_polynomial_term.
522 138 storres
 * Strategy: check that the object is a functional expression and
523 138 storres
 * if so check that it is made of cte * monomial.
524 138 storres
 */
525 138 storres
int
526 138 storres
pobyso_is_polynomial_term(sollya_obj_t objSo)
527 138 storres
{
528 138 storres
  int arity;
529 138 storres
  sollya_obj_t subFun1So = NULL;
530 138 storres
  sollya_obj_t subFun2So = NULL;
531 138 storres
  sollya_base_function_t head = 0;
532 138 storres
533 138 storres
  /* Arguments check. */
534 138 storres
  if (objSo == NULL)
535 138 storres
  {
536 138 storres
    pobyso_error_message("pobyso_is_polynomial_term",
537 138 storres
                        "NULL_POINTER_ARGUMENT",
538 138 storres
                        "The expression is a NULL pointer");
539 138 storres
    return 0;
540 138 storres
  }
541 138 storres
  /* The object must be a function. */
542 138 storres
  if (! sollya_lib_obj_is_function(objSo)) return 0;
543 138 storres
  /* A constant expression is degree 0 polynomial term. */
544 138 storres
  if (pobyso_is_constant_expression(objSo)) return 1;
545 138 storres
  /* A monomial is a polynomial term. */
546 138 storres
  if (pobyso_is_monomial(objSo)) return 1;
547 138 storres
  /* Decompose the functional expression and study the elements. */
548 138 storres
  if (! sollya_lib_decompose_function(objSo,
549 138 storres
                                      &head,
550 138 storres
                                      &arity,
551 138 storres
                                      &subFun1So,
552 138 storres
                                      &subFun2So,
553 138 storres
                                      NULL)) return 0;
554 138 storres
  /* Monomial case has been dealt with abobe. */
555 138 storres
  if (arity != 2)
556 138 storres
  {
557 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
558 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
559 138 storres
    return 0;
560 138 storres
  }
561 138 storres
  /* The expression must be: cte  * monomial or monomial * cte. */
562 138 storres
  if (head != SOLLYA_BASE_FUNC_MUL)
563 138 storres
  {
564 138 storres
    if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
565 138 storres
    if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
566 138 storres
    return 0;
567 138 storres
  }
568 138 storres
  if (! pobyso_is_monomial(subFun2So))
569 138 storres
  {
570 138 storres
    if (! pobyso_is_constant_expression(subFun2So) ||
571 138 storres
        ! pobyso_is_monomial(subFun1So))
572 138 storres
    {
573 138 storres
      if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
574 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
575 138 storres
      return 0;
576 138 storres
    }
577 138 storres
  }
578 138 storres
  else
579 138 storres
  {
580 138 storres
    if (! pobyso_is_constant_expression(subFun1So))
581 138 storres
    {
582 138 storres
      if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So);
583 138 storres
      if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So);
584 138 storres
      return 0;
585 138 storres
    }
586 138 storres
  }
587 138 storres
  return 1;
588 138 storres
} /* End pobyso_is_polynomial_term. */
589 133 storres
/** @see pobyso.h#pobyso_new_monomial. */
590 133 storres
pobyso_func_exp_t
591 134 storres
pobyso_new_monomial(pobyso_func_exp_t coefficientSo, long degree)
592 133 storres
{
593 134 storres
  sollya_obj_t degreeSo   = NULL;
594 134 storres
  sollya_obj_t varToPowSo = NULL;
595 134 storres
  sollya_obj_t monomialSo = NULL;
596 134 storres
  int initial_verbosity_level = 0;
597 134 storres
598 134 storres
  /* Arguments check. */
599 134 storres
  if (coefficientSo == NULL)
600 133 storres
  {
601 133 storres
    pobyso_error_message("pobyso_parse_string",
602 133 storres
                        "NULL_POINTER_ARGUMENT",
603 133 storres
                        "The expression is a NULL pointer");
604 134 storres
    return NULL;
605 133 storres
  }
606 134 storres
  if (! pobyso_is_constant_expression(coefficientSo))
607 133 storres
  {
608 134 storres
    return NULL;
609 133 storres
  }
610 133 storres
  if (degree < 0)
611 133 storres
  {
612 134 storres
    pobyso_error_message("pobyso_new_monomial",
613 134 storres
                        "NEGATIVE_DEGREE_ARGUMENT",
614 134 storres
                        "The degree is a negative integer");
615 134 storres
    return NULL;
616 133 storres
  }
617 134 storres
  /* If degree == 0, just return a copy of the coefficient. Do not
618 134 storres
   * return the coefficient itself to avoid "double clear" issues. */
619 134 storres
  if (degree == 0)
620 134 storres
  {
621 134 storres
    initial_verbosity_level = pobyso_set_verbosity_off();
622 134 storres
    monomialSo = sollya_lib_copy_obj(coefficientSo);
623 134 storres
    pobyso_set_verbosity_to(initial_verbosity_level);
624 134 storres
  }
625 134 storres
  degreeSo    = sollya_lib_constant_from_int64(degree);
626 134 storres
  varToPowSo  = sollya_lib_build_function_pow(sollya_lib_free_variable(),
627 134 storres
                                              degreeSo);
628 134 storres
  /* Do not use the "build" version to avoid "eating up" the coefficient. */
629 134 storres
  monomialSo = sollya_lib_mul(coefficientSo,varToPowSo);
630 134 storres
  sollya_lib_clear_obj(varToPowSo);
631 134 storres
  /* Do not clear degreeSa: it was "eaten" by sollya-lib_build_function. */
632 134 storres
  return monomialSo;
633 133 storres
} /* End pobyso_new_monomial. */
634 133 storres
635 149 storres
/* @see pobyso.h#pobyso_off */
636 149 storres
pobyso_on_off_t
637 149 storres
pobyso_off()
638 149 storres
{
639 149 storres
  return sollya_lib_off();
640 149 storres
} /* End pobyso_off. */
641 149 storres
642 149 storres
/* @see pobyso.h#pobyso_off */
643 149 storres
pobyso_on_off_t
644 149 storres
pobyso_on()
645 149 storres
{
646 149 storres
  return sollya_lib_on();
647 149 storres
} /* End pobyso_on. */
648 149 storres
649 149 storres
650 149 storres
/* @see pobyso.h#pobyso_parse_string */
651 130 storres
pobyso_func_exp_t
652 26 storres
pobyso_parse_string(const char* expression)
653 26 storres
{
654 128 storres
  int expressionLength, i;
655 127 storres
  char *expressionWithSemiCo;
656 136 storres
  sollya_obj_t expressionSo;
657 128 storres
658 127 storres
  /* Arguments check. */
659 26 storres
  if (expression == NULL)
660 26 storres
  {
661 26 storres
    pobyso_error_message("pobyso_parse_string",
662 26 storres
                        "NULL_POINTER_ARGUMENT",
663 26 storres
                        "The expression is a NULL pointer");
664 134 storres
    return  NULL;
665 26 storres
  }
666 127 storres
  expressionLength = strlen(expression);
667 127 storres
  if (expressionLength == 0)
668 127 storres
  {
669 127 storres
    pobyso_error_message("pobyso_parse_string",
670 127 storres
                        "EMPTY_STRING_ARGUMENT",
671 127 storres
                        "The expression is an empty string");
672 134 storres
    return NULL;
673 127 storres
  }
674 128 storres
  /* Search from the last char of the expression until, whichever happens
675 128 storres
   * first:
676 128 storres
   * a ";" is found;
677 128 storres
   * a char  != ';' is found the the ";" is appended.
678 128 storres
   * If the head of the string is reached before any of these two events happens
679 128 storres
   * return an error.
680 128 storres
   */
681 128 storres
  for (i = expressionLength - 1 ; i >= 0 ; i--)
682 127 storres
  {
683 128 storres
    if (expression[i] == ';') /* Nothing special to do:
684 128 storres
                                 try to parse the string*/
685 127 storres
    {
686 136 storres
      expressionSo = sollya_lib_parse_string(expression);
687 136 storres
      if (sollya_lib_obj_is_error(expressionSo))
688 136 storres
      {
689 136 storres
        sollya_lib_clear_obj(expressionSo);
690 136 storres
        return NULL;
691 136 storres
      }
692 136 storres
      else
693 136 storres
      {
694 136 storres
        return expressionSo;
695 136 storres
      }
696 127 storres
    }
697 128 storres
    else
698 128 storres
    {
699 128 storres
      if (expression[i] == ' ' || expression[i] == '\t') /* A blank,
700 128 storres
                                                           just continue. */
701 128 storres
      {
702 128 storres
        continue;
703 128 storres
      }
704 128 storres
      else /* a character != ';' and from a blank: create the ';'ed string. */
705 128 storres
      {
706 128 storres
        /* Create a new string for the expression, add the ";" and
707 128 storres
         * and call sollya_lib_parse_string. */
708 128 storres
        expressionWithSemiCo = calloc(i + 3, sizeof(char));
709 128 storres
        if (expressionWithSemiCo == NULL)
710 128 storres
        {
711 128 storres
          pobyso_error_message("pobyso_parse_string",
712 128 storres
                                "MEMORY_ALLOCATION_ERROR",
713 128 storres
                                "Could not allocate the expression string");
714 134 storres
          return NULL;
715 128 storres
        }
716 132 storres
        strncpy(expressionWithSemiCo, expression, i+1);
717 128 storres
        expressionWithSemiCo[i + 1] = ';';
718 128 storres
        expressionWithSemiCo[i + 2] = '\0';
719 136 storres
        expressionSo = sollya_lib_parse_string(expressionWithSemiCo);
720 128 storres
        free(expressionWithSemiCo);
721 136 storres
        if (sollya_lib_obj_is_error(expressionSo))
722 136 storres
        {
723 136 storres
          sollya_lib_clear_obj(expressionSo);
724 136 storres
          return NULL;
725 136 storres
        }
726 136 storres
        else
727 136 storres
        {
728 136 storres
          return expressionSo;
729 136 storres
        }
730 128 storres
      } /* End character != ';' and from a blank. */
731 128 storres
    /* Create a new string for the expression, add the ";" and
732 128 storres
     * and call sollya_lib_parse_string. */
733 128 storres
    } /* End else. */
734 128 storres
  } /* End for. */
735 128 storres
  /* We get here, it is because we did not find any char == anything different
736 128 storres
   * from ' ' or '\t': the string is empty.
737 128 storres
   */
738 128 storres
  pobyso_error_message("pobyso_parse_string",
739 128 storres
                       "ONLY_BLANK_ARGUMENT",
740 128 storres
                        "The expression is only made of blanks");
741 134 storres
  return NULL;
742 26 storres
} /* pobyso_parse_string */
743 26 storres
744 149 storres
/** @see pobyso.h#pobyso_range_from_bounds */
745 149 storres
pobyso_range_t
746 149 storres
pobyso_range_from_bounds(mpfr_t lowerBound, mpfr_t upperBound)
747 149 storres
{
748 149 storres
  /* Supferficial check of arguments. */
749 149 storres
  if (mpfr_cmp(lowerBound, upperBound) > 0)
750 149 storres
  {
751 149 storres
      return(NULL);
752 149 storres
  }
753 149 storres
  return(sollya_lib_range_from_bounds(lowerBound, upperBound));
754 149 storres
}
755 149 storres
756 132 storres
pobyso_func_exp_t
757 132 storres
pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function,
758 132 storres
                                      long int degree,
759 132 storres
                                      pobyso_range_t interval,
760 132 storres
                                      pobyso_func_exp_t weight,
761 132 storres
                                      double quality,
762 132 storres
                                      pobyso_range_t bounds)
763 132 storres
{
764 134 storres
  sollya_obj_t  degreeSo = NULL;
765 134 storres
  sollya_obj_t qualitySo = NULL;
766 132 storres
767 134 storres
  degreeSo = sollya_lib_constant_from_int(degree);
768 134 storres
  qualitySo = sollya_lib_constant_from_double(quality);
769 132 storres
770 134 storres
  sollya_lib_clear_obj(degreeSo);
771 134 storres
  sollya_lib_clear_obj(qualitySo);
772 132 storres
  return NULL;
773 132 storres
} /* End pobyso_remez_canonical_monomials_base. */
774 132 storres
775 132 storres
int
776 132 storres
pobyso_set_canonical_on()
777 132 storres
{
778 132 storres
  pobyso_on_off_t currentCanonicalModeSo;
779 132 storres
  pobyso_on_off_t on;
780 132 storres
781 132 storres
  currentCanonicalModeSo = sollya_lib_get_canonical();
782 132 storres
  if (sollya_lib_is_on(currentCanonicalModeSo))
783 132 storres
  {
784 134 storres
    sollya_lib_clear_obj(currentCanonicalModeSo);
785 134 storres
    return POBYSO_ON;
786 132 storres
  }
787 132 storres
  else
788 132 storres
  {
789 134 storres
    on = sollya_lib_on();
790 134 storres
    sollya_lib_set_canonical(on);
791 134 storres
    sollya_lib_clear_obj(on);
792 134 storres
    sollya_lib_clear_obj(currentCanonicalModeSo);
793 134 storres
    return POBYSO_OFF;
794 132 storres
  }
795 132 storres
} /* End pobyso_set_canonical_on. */
796 132 storres
797 134 storres
int
798 128 storres
pobyso_set_verbosity_off()
799 26 storres
{
800 134 storres
  sollya_obj_t verbosityLevelZeroSo;
801 134 storres
  sollya_obj_t currentVerbosityLevelSo = NULL;
802 134 storres
  int currentVerbosityLevel = 0;
803 128 storres
804 134 storres
  currentVerbosityLevelSo = sollya_lib_get_verbosity();
805 134 storres
  sollya_lib_get_constant_as_int(&currentVerbosityLevel,
806 134 storres
                                 currentVerbosityLevelSo);
807 134 storres
  verbosityLevelZeroSo = sollya_lib_constant_from_int(0);
808 134 storres
  sollya_lib_set_verbosity(verbosityLevelZeroSo);
809 134 storres
  sollya_lib_clear_obj(verbosityLevelZeroSo);
810 134 storres
  sollya_lib_clear_obj(currentVerbosityLevelSo);
811 134 storres
  return currentVerbosityLevel;
812 26 storres
} /* End of pobyso_set_verbosity_off. */
813 26 storres
814 134 storres
int
815 134 storres
pobyso_set_verbosity_to(int newVerbosityLevel)
816 26 storres
{
817 134 storres
  int initialVerbosityLevel = 0;
818 134 storres
  sollya_obj_t initialVerbosityLevelSo = NULL;
819 134 storres
  sollya_obj_t newVerbosityLevelSo = NULL;
820 134 storres
821 134 storres
  initialVerbosityLevelSo = sollya_lib_get_verbosity();
822 134 storres
  sollya_lib_get_constant_as_int(&initialVerbosityLevel,
823 134 storres
                                 initialVerbosityLevelSo);
824 134 storres
  sollya_lib_clear_obj(initialVerbosityLevelSo);
825 134 storres
  if (newVerbosityLevel < 0)
826 26 storres
  {
827 26 storres
    pobyso_error_message("pobyso_set_verbosity_to",
828 134 storres
                        "INVALID_VALUE",
829 134 storres
                        "The new verbosity level is a negative number");
830 134 storres
    return initialVerbosityLevel;
831 26 storres
  }
832 134 storres
  newVerbosityLevelSo = sollya_lib_constant_from_int(newVerbosityLevel);
833 134 storres
  sollya_lib_set_verbosity(newVerbosityLevelSo);
834 134 storres
  sollya_lib_clear_obj(newVerbosityLevelSo);
835 134 storres
  return initialVerbosityLevel;
836 26 storres
} /* End of pobyso_set_verbosity_to. */
837 26 storres
838 134 storres
/**
839 134 storres
 * @see pobyso.h#pobyso_subpoly
840 134 storres
 */
841 134 storres
pobyso_func_exp_t
842 136 storres
pobyso_subpoly(pobyso_func_exp_t polynomialSo, long expsNum, long int* expsList)
843 132 storres
{
844 136 storres
  sollya_obj_t  expsListSo    = NULL;
845 136 storres
  sollya_obj_t* expsList_pso  = NULL;
846 136 storres
  sollya_obj_t subpoly        = NULL;
847 136 storres
  int i, j;
848 134 storres
849 134 storres
  /* Arguments check. */
850 134 storres
  if (polynomialSo == NULL) return NULL;
851 134 storres
  if (expsNum < 0) return NULL;
852 134 storres
  if (expsNum == 0) return sollya_lib_copy_obj(polynomialSo);
853 136 storres
  if (expsList == 0) return NULL;
854 136 storres
  /* Create a list of Sollya constants. */
855 136 storres
  expsList_pso = (sollya_obj_t*) malloc(expsNum * sizeof(sollya_obj_t));
856 136 storres
  if (expsList_pso == NULL)
857 132 storres
  {
858 134 storres
    pobyso_error_message("pobyso_subpoly",
859 134 storres
                          "MEMORY_ALLOCATION_ERROR",
860 136 storres
                          "Could not allocate the Sollya exponents list");
861 134 storres
    return NULL;
862 132 storres
  }
863 136 storres
  /* Fill up the list. */
864 134 storres
  for (i = 0 ; i < expsNum ; i++)
865 134 storres
  {
866 136 storres
    /* Abort if an exponent is negative. */
867 136 storres
    if (expsList[i] < 0 )
868 136 storres
    {
869 136 storres
      for (j = 0 ; j < i ; j++)
870 136 storres
      {
871 136 storres
        sollya_lib_clear_obj(expsList_pso[j]);
872 136 storres
      }
873 136 storres
      free(expsList_pso);
874 136 storres
      return NULL;
875 136 storres
    }
876 136 storres
    expsList_pso[i] = sollya_lib_constant_from_int64(expsList[i]);
877 134 storres
  } /* End for */
878 136 storres
  expsListSo = sollya_lib_list(expsList_pso, expsNum);
879 136 storres
  for (i = 0 ; i < expsNum ; i++)
880 134 storres
  {
881 136 storres
    sollya_lib_clear_obj(expsList_pso[i]);
882 136 storres
  }
883 136 storres
  free(expsList_pso);
884 136 storres
  if (expsListSo == NULL)
885 136 storres
  {
886 134 storres
    pobyso_error_message("pobyso_subpoly",
887 134 storres
                          "LIST_CREATIONERROR",
888 134 storres
                          "Could not create the exponents list");
889 134 storres
    return NULL;
890 134 storres
  }
891 134 storres
  subpoly = sollya_lib_subpoly(polynomialSo, expsListSo);
892 136 storres
  pobyso_autoprint(expsListSo);
893 134 storres
  sollya_lib_clear_obj(expsListSo);
894 134 storres
  return subpoly;
895 134 storres
} /* pobyso_subpoly. */
896 132 storres
897 26 storres
/* Attic from the sollya_lib < 4. */
898 26 storres
#if 0
899 26 storres
chain*
900 26 storres
pobyso_create_canonical_monomials_base(const unsigned int degree)
901 26 storres
{
902 26 storres
  int i    = 0;
903 26 storres
  chain *monomials  = NULL;
904 26 storres
  node  *monomial   = NULL;
905 26 storres

906 26 storres
  for(i = degree ; i >= 0  ; i--)
907 26 storres
  {
908 26 storres
     monomial   = makePow(makeVariable(), makeConstantDouble((double)i));
909 26 storres
     monomials  = addElement(monomials, monomial);
910 26 storres
     fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
911 26 storres
  }
912 26 storres
  if (monomials == NULL)
913 26 storres
  {
914 26 storres
    pobyso_error_message("pobyso_create_canonical_monomial_base",
915 26 storres
                        "CHAIN_CREATION_ERROR",
916 26 storres
                        "Could not create the monomials chain");
917 26 storres
    return(NULL);
918 26 storres
  }
919 26 storres
  return(monomials);
920 26 storres
} /* End pobyso_create_canonical_monomials_base. */
921 26 storres

922 26 storres
chain*
923 26 storres
pobyso_create_chain_from_int_array(int* intArray,
924 26 storres
                                    const unsigned int arrayLength)
925 26 storres
{
926 26 storres
  int i = 0;
927 26 storres
  chain *newChain = NULL;
928 26 storres
  int *currentInt;
929 26 storres

930 26 storres
  if (arrayLength == 0) return(NULL);
931 26 storres
  if (intArray == NULL)
932 26 storres
  {
933 26 storres
    pobyso_error_message("pobyso_create_chain_from_int_array",
934 26 storres
                        "NULL_POINTER_ARGUMENT",
935 26 storres
                        "The array is a NULL pointer");
936 26 storres
    return(NULL);
937 26 storres
  }
938 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
939 26 storres
  {
940 26 storres
    currentInt = malloc(sizeof(int));
941 26 storres
    if (currentInt == NULL)
942 26 storres
    {
943 26 storres
      pobyso_error_message("pobyso_create_chain_from_int_array",
944 26 storres
                          "MEMORY_ALLOCATION_ERROR",
945 26 storres
                          "Could not allocate one of the integers");
946 26 storres
      freeChain(newChain, free);
947 26 storres
      return(NULL);
948 26 storres
    }
949 26 storres
    *currentInt = intArray[i];
950 26 storres
    newChain = addElement(newChain, currentInt);
951 26 storres
  }
952 26 storres
  return(newChain);
953 26 storres
} // End pobyso_create_chain_from_int_array. */
954 26 storres

955 26 storres
chain*
956 26 storres
pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
957 26 storres
                                        const unsigned int arrayLength)
958 26 storres
{
959 26 storres
  int i = 0;
960 26 storres
  chain *newChain = NULL;
961 26 storres
  unsigned int *currentInt;
962 26 storres

963 26 storres
  /* Argument checking. */
964 26 storres
  if (arrayLength == 0) return(NULL);
965 26 storres
  if (intArray == NULL)
966 26 storres
  {
967 26 storres
    pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
968 26 storres
                        "NULL_POINTER_ARGUMENT",
969 26 storres
                        "The array is a NULL pointer");
970 26 storres
    return(NULL);
971 26 storres
  }
972 26 storres
  for (i = arrayLength - 1 ; i >= 0 ; i--)
973 26 storres
  {
974 26 storres
    currentInt = malloc(sizeof(unsigned int));
975 26 storres
    if (currentInt == NULL)
976 26 storres
    {
977 26 storres
      pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
978 26 storres
                          "MEMORY_ALLOCATION_ERROR",
979 26 storres
                          "Could not allocate one of the integers");
980 26 storres
      freeChain(newChain, free);
981 26 storres
      return(NULL);
982 26 storres
    }
983 26 storres
    *currentInt = intArray[i];
984 26 storres
    newChain = addElement(newChain, currentInt);
985 26 storres
  }
986 26 storres
  return(newChain);
987 26 storres
} // End pobyso_create_chain_from_unsigned_int_array. */
988 26 storres

989 26 storres
node*
990 26 storres
pobyso_differentiate(node *functionNode)
991 26 storres
{
992 26 storres
  /* Argument checking. */
993 26 storres
  node *differentialNode;
994 26 storres
  if (functionNode == NULL)
995 26 storres
  {
996 26 storres
    pobyso_error_message("pobyso_differentiate",
997 26 storres
                        "NULL_POINTER_ARGUMENT",
998 26 storres
                        "The function to differentiate is a NULL pointer");
999 26 storres
    return(NULL);
1000 26 storres
  }
1001 26 storres
  differentialNode = differentiate(functionNode);
1002 26 storres
  if (differentialNode == NULL)
1003 26 storres
  {
1004 26 storres
    pobyso_error_message("pobyso_differentiate",
1005 26 storres
                        "INTERNAL ERROR",
1006 26 storres
                        "Sollya could not differentiate the function");
1007 26 storres
  }
1008 26 storres
  return(differentialNode);
1009 26 storres
} // End pobyso_differentiate
1010 26 storres

1011 26 storres

1012 26 storres
int
1013 26 storres
pobyso_dirty_infnorm(mpfr_t infNorm,
1014 26 storres
                      node *functionNode,
1015 26 storres
                      mpfr_t lowerBound,
1016 26 storres
                      mpfr_t upperBound,
1017 26 storres
                      mp_prec_t precision)
1018 26 storres
{
1019 26 storres
  int functionCallResult;
1020 26 storres
  /* Arguments checking. */
1021 26 storres
  if (functionNode == NULL)
1022 26 storres
  {
1023 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
1024 26 storres
                        "NULL_POINTER_ARGUMENT",
1025 26 storres
                        "The function to compute is a NULL pointer");
1026 26 storres
    return(1);
1027 26 storres
  }
1028 26 storres
  if (mpfr_cmp(lowerBound, upperBound) > 0)
1029 26 storres
  {
1030 26 storres
    pobyso_error_message("pobyso_dirty_infnorm",
1031 26 storres
                        "INCOHERENT_INPUT_DATA",
1032 26 storres
                        "The lower bond is greater than the upper bound");
1033 26 storres
    return(1);
1034 26 storres
  }
1035 26 storres
  /* Particular cases. */
1036 26 storres
  if (mpfr_cmp(lowerBound, upperBound) == 0)
1037 26 storres
  {
1038 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
1039 26 storres
                                                  functionNode,
1040 26 storres
                                                  lowerBound,
1041 26 storres
                                                  precision);
1042 26 storres
    return(functionCallResult);
1043 26 storres
  }
1044 26 storres
  if (isConstant(functionNode))
1045 26 storres
  {
1046 26 storres
    functionCallResult = pobyso_evaluate_faithful(infNorm,
1047 26 storres
                                                  functionNode,
1048 26 storres
                                                  lowerBound,
1049 26 storres
                                                  precision);
1050 26 storres
    if (!functionCallResult)
1051 26 storres
    {
1052 26 storres
      mpfr_abs(infNorm, infNorm, MPFR_RNDN);
1053 26 storres
    }
1054 26 storres
    return(functionCallResult);
1055 26 storres
  }
1056 26 storres
  uncertifiedInfnorm(infNorm,
1057 26 storres
                      functionNode,
1058 26 storres
                      lowerBound,
1059 26 storres
                      upperBound,
1060 26 storres
                      POBYSO_DEFAULT_POINTS,
1061 26 storres
                      precision);
1062 26 storres
  return(0);
1063 26 storres
} /* End pobyso_dirty_infnorm. */
1064 26 storres

1065 26 storres
int
1066 26 storres
pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
1067 26 storres
                          node *nodeToEvaluate,
1068 26 storres
                          mpfr_t argument,
1069 26 storres
                          mpfr_prec_t precision)
1070 26 storres
{
1071 26 storres
  /* Check input arguments. */
1072 26 storres
  if (nodeToEvaluate == NULL)
1073 26 storres
  {
1074 26 storres
    pobyso_error_message("pobyso_evaluate_faithful",
1075 26 storres
                        "NULL_POINTER_ARGUMENT",
1076 26 storres
                        "nodeToEvaluate is a NULL pointer");
1077 26 storres
    return(1);
1078 26 storres
  }
1079 26 storres
  evaluateFaithful(faithfulEvaluation,
1080 26 storres
                   nodeToEvaluate,
1081 26 storres
                   argument,
1082 26 storres
                   precision);
1083 26 storres
  return(0);
1084 26 storres
} /* End pobyso_evaluate_faithfull. */
1085 26 storres

1086 26 storres
chain*
1087 26 storres
pobyso_find_zeros(node *function,
1088 26 storres
                  mpfr_t *lower_bound,
1089 26 storres
                  mpfr_t *upper_bound)
1090 26 storres
{
1091 26 storres
  mp_prec_t currentPrecision;
1092 26 storres
  mpfr_t currentDiameter;
1093 26 storres
  rangetype bounds;
1094 26 storres

1095 26 storres
  currentPrecision = getToolPrecision();
1096 26 storres
  mpfr_init2(currentDiameter, currentPrecision);
1097 26 storres

1098 26 storres
  bounds.a = lower_bound;
1099 26 storres
  bounds.b = upper_bound;
1100 26 storres

1101 26 storres
  if (bounds.a == NULL || bounds.b == NULL)
1102 26 storres
  {
1103 26 storres
    pobyso_error_message("pobyso_find_zeros",
1104 26 storres
                        "MEMORY_ALLOCATION_ERROR",
1105 26 storres
                        "Could not allocate one of the bounds");
1106 26 storres
    return(NULL);
1107 26 storres
  }
1108 26 storres
  return(findZerosFunction(function,
1109 26 storres
                            bounds,
1110 26 storres
                            currentPrecision,
1111 26 storres
                            currentDiameter));
1112 26 storres
} /* End pobyso_find_zeros. */
1113 26 storres

1114 26 storres
void
1115 26 storres
pobyso_free_chain_of_nodes(chain *theChainOfNodes)
1116 26 storres
{
1117 26 storres
  node *currentNode           = NULL;
1118 26 storres
  chain *currentChainElement  = NULL;
1119 26 storres
  chain *nextChainElement     = NULL;
1120 26 storres

1121 26 storres
  nextChainElement = theChainOfNodes;
1122 26 storres

1123 26 storres
  while(nextChainElement != NULL)
1124 26 storres
  {
1125 26 storres
    currentChainElement = nextChainElement;
1126 26 storres
    currentNode = (node*)(currentChainElement->value);
1127 26 storres
    nextChainElement = nextChainElement->next;
1128 26 storres
    free_memory(currentNode);
1129 26 storres
    free((void*)(currentChainElement));
1130 26 storres
  }
1131 26 storres
} /* End pobyso_free_chain_of_nodes. */
1132 26 storres

1133 26 storres
void
1134 26 storres
pobyso_free_range(rangetype range)
1135 26 storres
{
1136 26 storres

1137 26 storres
  mpfr_clear(*(range.a));
1138 26 storres
  mpfr_clear(*(range.b));
1139 26 storres
  free(range.a);
1140 26 storres
  free(range.b);
1141 26 storres
} /* End pobyso_free_range. */
1142 26 storres

1143 26 storres
node*
1144 26 storres
pobyso_fp_minimax_canonical_monomials_base(node *function,
1145 26 storres
                                      int degree,
1146 26 storres
                                      chain *formats,
1147 26 storres
                                      chain *points,
1148 26 storres
                                      mpfr_t lowerBound,
1149 26 storres
                                      mpfr_t upperBound,
1150 26 storres
                                      int fpFixedArg,
1151 26 storres
                                      int absRel,
1152 26 storres
                                      node *constPart,
1153 26 storres
                                      node *minimax)
1154 26 storres
{
1155 26 storres
  return(NULL);
1156 26 storres
} /* End pobyso_fp_minimax_canonical_monomials_base. */
1157 26 storres

1158 26 storres
node*
1159 26 storres
pobyso_parse_function(char *functionString,
1160 26 storres
                      char *freeVariableNameString)
1161 26 storres
{
1162 26 storres
  if (functionString == NULL || freeVariableNameString == NULL)
1163 26 storres
  {
1164 26 storres
    pobyso_error_message("pobyso_parse_function",
1165 26 storres
                        "NULL_POINTER_ARGUMENT",
1166 26 storres
                        "One of the arguments is a NULL pointer");
1167 26 storres
    return(NULL);
1168 26 storres
  }
1169 26 storres
  return(parseString(functionString));
1170 26 storres

1171 26 storres
} /* End pobyso_parse_function */
1172 26 storres

1173 26 storres
node*
1174 26 storres
pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
1175 26 storres
                                                  unsigned int mode,
1176 26 storres
                                                  mpfr_t lowerBound,
1177 26 storres
                                                  mpfr_t upperBound,
1178 26 storres
                                                  mpfr_t eps)
1179 26 storres
{
1180 26 storres
  node *weight              = NULL;
1181 26 storres
  node *bestApproxPolyNode  = NULL;
1182 26 storres
  node *bestApproxHorner    = NULL;
1183 26 storres
  node *errorNode           = NULL;
1184 26 storres
  rangetype degreeRange;
1185 26 storres
  mpfr_t quality;
1186 26 storres
  mpfr_t currentError;
1187 26 storres
  unsigned int degree;
1188 26 storres

1189 26 storres
  /* Check the parameters. */
1190 26 storres
  if (functionNode == NULL)
1191 26 storres
  {
1192 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
1193 26 storres
                        "NULL_POINTER_ARGUMENT",
1194 26 storres
                        "functionNode is a NULL pointer");
1195 26 storres
    return(NULL);
1196 26 storres
  }
1197 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
1198 26 storres
  {
1199 26 storres
    pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
1200 26 storres
                        "INCOHERENT_INPUT_DATA",
1201 26 storres
                        "the lower_bound >= upper_bound");
1202 26 storres
    return(NULL);
1203 26 storres
  }
1204 26 storres
  /* Set the weight. */
1205 26 storres
  if (mode == POBYSO_ABSOLUTE)
1206 26 storres
  {
1207 26 storres
    /* Set the weight to 1 for the ABSOLUTE_MODE. */
1208 26 storres
    weight = makeConstantDouble(1.0);
1209 26 storres
  }
1210 26 storres
  else
1211 26 storres
  {
1212 26 storres
    if (mode == POBYSO_RELATIVE)
1213 26 storres
    {
1214 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
1215 26 storres
                          "NOT_IMPLEMENTED",
1216 26 storres
                          "the search for relative error is not implemented yet");
1217 26 storres
      return(NULL);
1218 26 storres
    }
1219 26 storres
    else
1220 26 storres
    {
1221 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
1222 26 storres
                          "INCOHERENT_INPUT_DATA",
1223 26 storres
                          "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
1224 26 storres
      return(NULL);
1225 26 storres
    }
1226 26 storres
  }
1227 26 storres
  //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
1228 26 storres
  degreeRange = guessDegree(functionNode,
1229 26 storres
                            weight,
1230 26 storres
                            lowerBound,
1231 26 storres
                            upperBound,
1232 26 storres
                            eps,
1233 26 storres
                            POBYSO_GUESS_DEGREE_BOUND);
1234 26 storres
  degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
1235 26 storres
  //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
1236 26 storres
  fprintf(stderr, "Guessed degree: ");
1237 26 storres
  mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
1238 26 storres
  fprintf(stderr, " - as int: %u\n", degree);
1239 26 storres
  /* Reduce the degree by 1 in the foolish hope it could work. */
1240 26 storres
  if (degree > 0) degree--;
1241 26 storres
  /* Both elements of degreeRange where "inited" within guessDegree. */
1242 26 storres
  mpfr_clear(*(degreeRange.a));
1243 26 storres
  mpfr_clear(*(degreeRange.b));
1244 26 storres
  free(degreeRange.a);
1245 26 storres
  free(degreeRange.b);
1246 26 storres
  /* Mimic the default behavior of interactive Sollya. */
1247 26 storres
  mpfr_init(quality);
1248 26 storres
  mpfr_set_d(quality, 1e-5, MPFR_RNDN);
1249 26 storres
  mpfr_init2(currentError, getToolPrecision());
1250 26 storres
  mpfr_set_inf(currentError,1);
1251 26 storres

1252 26 storres
  /* Try to refine the initial guess: loop with increasing degrees until
1253 26 storres
   * we find a satisfactory one. */
1254 26 storres
  while(mpfr_cmp(currentError, eps) > 0)
1255 26 storres
  {
1256 26 storres
    /* Get rid of the previous polynomial, if any. */
1257 26 storres
    if (bestApproxPolyNode != NULL)
1258 26 storres
    {
1259 26 storres
      free_memory(bestApproxPolyNode);
1260 26 storres
    }
1261 26 storres
    fprintf(stderr, "Degree: %u\n", degree);
1262 26 storres
    fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
1263 26 storres
    /* Try to find a polynomial with the guessed degree. */
1264 26 storres
    bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
1265 26 storres
                                                            weight,
1266 26 storres
                                                            degree,
1267 26 storres
                                                            lowerBound,
1268 26 storres
                                                            upperBound,
1269 26 storres
                                                            quality);
1270 26 storres

1271 26 storres
    if (bestApproxPolyNode == NULL)
1272 26 storres
    {
1273 26 storres
      pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
1274 26 storres
                          "INTERNAL_ERROR",
1275 26 storres
                          "could not compute the bestApproxPolyNode");
1276 26 storres
      mpfr_clear(currentError);
1277 26 storres
      mpfr_clear(quality);
1278 26 storres
      return(NULL);
1279 26 storres
    }
1280 26 storres

1281 26 storres
    setDisplayMode(DISPLAY_MODE_DECIMAL);
1282 26 storres
    fprintTree(stderr, bestApproxPolyNode);
1283 26 storres
    fprintf(stderr, "\n\n");
1284 26 storres

1285 26 storres
    errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
1286 26 storres
    /* Check the error with the computed polynomial. */
1287 26 storres
    uncertifiedInfnorm(currentError,
1288 26 storres
                        errorNode,
1289 26 storres
                        lowerBound,
1290 26 storres
                        upperBound,
1291 26 storres
                        POBYSO_INF_NORM_NUM_POINTS,
1292 26 storres
                        getToolPrecision());
1293 26 storres
    fprintf(stderr, "Inf norm: ");
1294 26 storres
    mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
1295 26 storres
    fprintf(stderr, "\n\n");
1296 26 storres
    /* Free the errorNode but not the bestApproxPolyNode (we need it if
1297 26 storres
     * we exit the loop at the next iteration). */
1298 26 storres
    free_memory(errorNode);
1299 26 storres
    degree++;
1300 26 storres
  }
1301 26 storres
  /* Use an intermediate variable, since horner() creates a new node
1302 26 storres
   * and does not reorder the argument "in place". This allows for the memory
1303 26 storres
   * reclaim of bestApproxHorner.
1304 26 storres
   */
1305 26 storres
  bestApproxHorner = horner(bestApproxPolyNode);
1306 26 storres
  free_memory(bestApproxPolyNode);
1307 26 storres
  mpfr_clear(currentError);
1308 26 storres
  mpfr_clear(quality);
1309 26 storres
  free_memory(weight);
1310 26 storres
  return(bestApproxHorner);
1311 26 storres
} /* End pobyso_remez_approx_canonical_monomials_base_for_error */
1312 26 storres

1313 26 storres
node*
1314 26 storres
pobyso_remez_canonical_monomials_base(node *function,
1315 26 storres
                                     node *weight,
1316 26 storres
                                     unsigned int degree,
1317 26 storres
                                     mpfr_t lowerBound,
1318 26 storres
                                     mpfr_t upperBound,
1319 26 storres
                                     mpfr_t quality)
1320 26 storres
{
1321 26 storres
  node  *bestApproxPoly = NULL;
1322 26 storres
  chain *monomials      = NULL;
1323 26 storres
  chain *curMonomial    = NULL;
1324 26 storres

1325 26 storres
  mpfr_t satisfying_error;
1326 26 storres
  mpfr_t target_error;
1327 26 storres

1328 26 storres
  /* Argument checking */
1329 26 storres
  /* Function tree. */
1330 26 storres
  if (function == NULL)
1331 26 storres
  {
1332 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1333 26 storres
                        "NULL_POINTER_ARGUMENT",
1334 26 storres
                        "the \"function\" argument is a NULL pointer");
1335 26 storres
    return(NULL);
1336 26 storres
  }
1337 26 storres
  if (weight == NULL)
1338 26 storres
  {
1339 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1340 26 storres
                        "NULL_POINTER_ARGUMENT",
1341 26 storres
                        "the \"weight\" argument is a NULL pointer");
1342 26 storres
    return(NULL);
1343 26 storres
  }
1344 26 storres
  /* Check the bounds. */
1345 26 storres
  if (mpfr_cmp(lowerBound, upperBound) >= 0)
1346 26 storres
  {
1347 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1348 26 storres
                        "INCOHERENT_IMPUT_DATA",
1349 26 storres
                        "the lower_bound >= upper_bound");
1350 26 storres
    return(NULL);
1351 26 storres
  }
1352 26 storres
  /* The quality must be a non null positive number. */
1353 26 storres
  if (mpfr_sgn(quality) <= 0)
1354 26 storres
  {
1355 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1356 26 storres
                        "INCOHERENT_INPUT_DATA",
1357 26 storres
                        "the quality <= 0");
1358 26 storres
  }
1359 26 storres
  /* End argument checking. */
1360 26 storres
  /* Create the monomials nodes chains. */
1361 26 storres
  monomials = pobyso_create_canonical_monomials_base(degree);
1362 26 storres
  fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
1363 26 storres
  if (monomials == NULL || (lengthChain(monomials) != degree + 1))
1364 26 storres
  {
1365 26 storres
    pobyso_error_message("pobyso_remez_canonical_monomials_base",
1366 26 storres
                        "CHAIN_CREATION_ERROR",
1367 26 storres
                        "could not create the monomials chain");
1368 26 storres
    return(NULL);
1369 26 storres
  }
1370 26 storres
  curMonomial = monomials;
1371 26 storres

1372 26 storres
  while (curMonomial != NULL)
1373 26 storres
  {
1374 26 storres
    fprintf(stderr, "monomial tree: ");
1375 26 storres
    //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
1376 26 storres
    fprintTree(stderr, (node*)(curMonomial->value));
1377 26 storres
    fprintf(stderr, "\n");
1378 26 storres
    curMonomial = curMonomial->next;
1379 26 storres
  }
1380 26 storres

1381 26 storres
  /* Deal with NULL weight. */
1382 26 storres
  if (weight == NULL)
1383 26 storres
  {
1384 26 storres
    weight = makeConstantDouble(1.0);
1385 26 storres
  }
1386 26 storres
  /* Compute the best polynomial with the required quality.
1387 26 storres
     The behavior is as if satisfying error and target_error had
1388 26 storres
     not been used.*/
1389 26 storres
  mpfr_init(satisfying_error);
1390 26 storres
  mpfr_init(target_error);
1391 26 storres
  mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
1392 26 storres
  mpfr_set_inf(target_error, 1);
1393 26 storres

1394 26 storres

1395 26 storres
  fprintf(stderr, "satisfying_error: ");
1396 26 storres
  mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
1397 26 storres
  fprintf(stderr, ".\n");
1398 26 storres
  fprintf(stderr, "target_error: ");
1399 26 storres
  mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
1400 26 storres
  fprintf(stderr, ".\n");
1401 26 storres

1402 26 storres
  fprintf(stderr,
1403 26 storres
          "current precision: %li\n", getToolPrecision());
1404 26 storres
  /* Call the Sollya function. */
1405 26 storres
  bestApproxPoly = remez(function,
1406 26 storres
                          weight,
1407 26 storres
                          monomials,
1408 26 storres
                          lowerBound,
1409 26 storres
                          upperBound,
1410 26 storres
                          quality,
1411 26 storres
                          satisfying_error,
1412 26 storres
                          target_error,
1413 26 storres
                          getToolPrecision());
1414 26 storres

1415 26 storres
  mpfr_clear(satisfying_error);
1416 26 storres
  mpfr_clear(target_error);
1417 26 storres
  pobyso_free_chain_of_nodes(monomials);
1418 26 storres

1419 26 storres
  return(bestApproxPoly);
1420 26 storres
} /* End pobyso_remez_canonical_monomials_base. */
1421 26 storres

1422 26 storres
#endif
1423 26 storres
1424 26 storres
void
1425 26 storres
pobyso_error_message(char *functionName, char *messageName, char* message)
1426 26 storres
{
1427 26 storres
  fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
1428 26 storres
} /* End pobyso_error_message */