root / pobysoC-4.0 / src / pobyso.c
Historique | Voir | Annoter | Télécharger (44,6 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 | 287 | storres | /* @see pobyso.h#pobyso_clear_obj */
|
43 | 287 | storres | void
|
44 | 287 | storres | pobyso_clear_obj(sollya_obj_t objSo) |
45 | 287 | storres | { |
46 | 287 | storres | if (objSo == NULL) {return;} |
47 | 287 | storres | sollya_lib_clear_obj(objSo); |
48 | 287 | storres | } /* End pobyso_clear_obj. */
|
49 | 287 | storres | |
50 | 147 | storres | /* @see pobyso.h#pobyso_dirty_find_zeros */
|
51 | 147 | storres | mpfr_t* |
52 | 147 | storres | pobyso_dirty_find_zeros_bounds(pobyso_func_exp_t funcExpSo, |
53 | 147 | storres | mpfr_t lowerBound, |
54 | 147 | storres | mpfr_t upperBound, |
55 | 147 | storres | int* zerosCount)
|
56 | 147 | storres | { |
57 | 147 | storres | pobyso_range_t rangeSo; |
58 | 147 | storres | sollya_obj_t zerosListSo = NULL;
|
59 | 147 | storres | sollya_obj_t* zerosArraySo = NULL;
|
60 | 147 | storres | mpfr_t* zerosArrayMp = NULL;
|
61 | 147 | storres | pobyso_precision_t prec; |
62 | 147 | storres | |
63 | 147 | storres | int endEll = 0; |
64 | 147 | storres | int i,j;
|
65 | 147 | storres | |
66 | 147 | storres | /* Arguments check. */
|
67 | 147 | storres | *zerosCount = -1;
|
68 | 147 | storres | if (funcExpSo == NULL ) |
69 | 147 | storres | { |
70 | 147 | storres | pobyso_error_message("pobyso_dirty_find_zeros",
|
71 | 147 | storres | "NULL_POINTER_ARGUMENT",
|
72 | 147 | storres | "The funcExpSo argument is a NULL pointer");
|
73 | 147 | storres | return NULL; |
74 | 147 | storres | } |
75 | 147 | storres | if (mpfr_cmp(lowerBound, upperBound) > 0) |
76 | 147 | storres | { |
77 | 147 | storres | pobyso_error_message("pobyso_dirty_find_zeros",
|
78 | 147 | storres | "INVALID_INTERVAL_ARGUMENT",
|
79 | 147 | storres | "The lower bound is larger than the upper bound");
|
80 | 147 | storres | return NULL; |
81 | 147 | storres | } |
82 | 147 | storres | /* Make a range out of the bounds. */
|
83 | 147 | storres | rangeSo = sollya_lib_range_from_bounds(lowerBound, upperBound); |
84 | 147 | storres | if (rangeSo == NULL) |
85 | 147 | storres | { |
86 | 147 | storres | return NULL; |
87 | 147 | storres | } |
88 | 147 | storres | zerosListSo = sollya_lib_dirtyfindzeros(funcExpSo, rangeSo); |
89 | 147 | storres | if (zerosListSo == NULL) |
90 | 147 | storres | { |
91 | 147 | storres | return NULL; |
92 | 147 | storres | } |
93 | 147 | storres | sollya_lib_clear_obj(rangeSo); |
94 | 147 | storres | /* Transform the Sollya list into an MPFR list. */
|
95 | 147 | storres | if (! sollya_lib_get_list_elements(&zerosArraySo,
|
96 | 147 | storres | zerosCount, |
97 | 147 | storres | &endEll, |
98 | 147 | storres | zerosListSo)) |
99 | 147 | storres | { |
100 | 147 | storres | sollya_lib_clear_obj(zerosListSo); |
101 | 147 | storres | *zerosCount = -1;
|
102 | 147 | storres | return NULL; |
103 | 147 | storres | } |
104 | 147 | storres | sollya_lib_clear_obj(zerosListSo); |
105 | 147 | storres | zerosArrayMp = (mpfr_t*) malloc(*zerosCount * sizeof(mpfr_t));
|
106 | 147 | storres | if (zerosArrayMp == NULL) |
107 | 147 | storres | { |
108 | 147 | storres | pobyso_error_message("pobyso_dirty_find_zeros",
|
109 | 147 | storres | "MEMORY_ALLOCATION_ERROR",
|
110 | 147 | storres | "Could not allocate zeroes array");
|
111 | 147 | storres | *zerosCount = -1;
|
112 | 147 | storres | return NULL; |
113 | 147 | storres | } |
114 | 147 | storres | for (i = 0 ; i < *zerosCount ; i++) |
115 | 147 | storres | { |
116 | 147 | storres | if (! sollya_lib_get_prec_of_constant(&prec, zerosArraySo[i]))
|
117 | 147 | storres | { |
118 | 147 | storres | /* Clean up the already allocated MPFRs. */
|
119 | 147 | storres | for (j = 0 ; j < i ; j++) mpfr_clear(zerosArrayMp[j]); |
120 | 147 | storres | /* Clean up the zerosArrayMp array itself. */
|
121 | 147 | storres | free(zerosArrayMp); |
122 | 147 | storres | /* Clean up what is left in the zerosArraySo. */
|
123 | 147 | storres | for (j = i ; j < *zerosCount ; j++) sollya_lib_clear_obj(zerosArraySo[j]);
|
124 | 147 | storres | /* Clean up the zerosArraySo array itself. */
|
125 | 147 | storres | sollya_lib_free(zerosArraySo); |
126 | 147 | storres | *zerosCount = -1;
|
127 | 147 | storres | return NULL; |
128 | 147 | storres | } |
129 | 147 | storres | mpfr_init2(zerosArrayMp[i], prec); |
130 | 147 | storres | if (! sollya_lib_get_constant(zerosArrayMp[i], zerosArraySo[i]))
|
131 | 147 | storres | { |
132 | 147 | storres | /* Clean up the already allocated MPFRs. */
|
133 | 147 | storres | for (j = 0 ; j <= i ; j++) mpfr_clear(zerosArrayMp[j]); |
134 | 147 | storres | /* Clean up the zerosArrayMp array itself. */
|
135 | 147 | storres | free(zerosArrayMp); |
136 | 147 | storres | /* Clean up what is left in the zerosArraySo. */
|
137 | 147 | storres | for (j = i ; j < *zerosCount ; j++) sollya_lib_clear_obj(zerosArraySo[j]);
|
138 | 147 | storres | /* Clean up the zerosArraySo array itself. */
|
139 | 147 | storres | sollya_lib_free(zerosArraySo); |
140 | 147 | storres | *zerosCount = -1;
|
141 | 147 | storres | return NULL; |
142 | 147 | storres | } |
143 | 147 | storres | sollya_lib_clear_obj(zerosArraySo[i]); |
144 | 147 | storres | } /* End for i. */
|
145 | 147 | storres | sollya_lib_free(zerosArraySo); |
146 | 147 | storres | return zerosArrayMp;
|
147 | 147 | storres | } /* End pobyso_dirty_find_zeros. */
|
148 | 147 | storres | |
149 | 140 | storres | /* @see pobyso.h#pobyso_evaluate_constant */
|
150 | 139 | storres | int
|
151 | 139 | storres | pobyso_evaluate_constant(pobyso_func_exp_t functionSo, |
152 | 139 | storres | mpfr_t argumentMp, |
153 | 139 | storres | mpfr_t evaluationMp) |
154 | 139 | storres | { |
155 | 139 | storres | sollya_obj_t argumentSo = NULL;
|
156 | 139 | storres | sollya_obj_t evaluationSo = NULL;
|
157 | 140 | storres | mpfr_t evaluationTmpMp1; |
158 | 140 | storres | mpfr_t evaluationTmpMp2; |
159 | 140 | storres | mpfr_prec_t evalPrec = 0;
|
160 | 139 | storres | |
161 | 139 | storres | /* Test argument. */
|
162 | 139 | storres | if (functionSo == NULL || argumentMp == NULL || evaluationMp == NULL) |
163 | 139 | storres | { |
164 | 139 | storres | pobyso_error_message("pobyso_evaluate_constant",
|
165 | 139 | storres | "NULL_POINTER_ARGUMENT",
|
166 | 139 | storres | "One of the arguments is a NULL pointer");
|
167 | 139 | storres | return 1; |
168 | 139 | storres | } |
169 | 139 | storres | if (! sollya_lib_obj_is_function(functionSo))
|
170 | 139 | storres | { |
171 | 139 | storres | pobyso_error_message("pobyso_evaluate_constant",
|
172 | 139 | storres | "INVALID_TYPE_ARGUMENT",
|
173 | 139 | storres | "The functionSo argument is not a functional expression");
|
174 | 139 | storres | return 1; |
175 | 139 | storres | } |
176 | 140 | storres | /* Function evaluation and checks. */
|
177 | 139 | storres | argumentSo = sollya_lib_constant(argumentMp); |
178 | 139 | storres | evaluationSo = sollya_lib_evaluate(functionSo, argumentSo); |
179 | 147 | storres | /* Not needed any more. */
|
180 | 147 | storres | //pobyso_autoprint(evaluationSo);
|
181 | 139 | storres | sollya_lib_clear_obj(argumentSo); |
182 | 140 | storres | /* The range case: we return the mean of the bounds. The result
|
183 | 140 | storres | * is not faithfully rounded. */
|
184 | 140 | storres | if (sollya_lib_obj_is_range(evaluationSo))
|
185 | 140 | storres | { |
186 | 147 | storres | //pobyso_autoprint(evaluationSo);
|
187 | 140 | storres | if (sollya_lib_get_prec_of_range(&evalPrec, evaluationSo))
|
188 | 140 | storres | { |
189 | 140 | storres | mpfr_init2(evaluationTmpMp1, evalPrec); |
190 | 140 | storres | mpfr_init2(evaluationTmpMp2, evalPrec); |
191 | 140 | storres | if (sollya_lib_get_bounds_from_range(evaluationTmpMp1,
|
192 | 140 | storres | evaluationTmpMp2, |
193 | 140 | storres | evaluationSo)) |
194 | 140 | storres | { |
195 | 140 | storres | /* Compute the mean of the bounds. */
|
196 | 140 | storres | mpfr_clear(evaluationMp); |
197 | 140 | storres | mpfr_init2(evaluationMp, evalPrec); |
198 | 140 | storres | mpfr_add(evaluationMp,evaluationTmpMp1,evaluationTmpMp2,MPFR_RNDN); |
199 | 140 | storres | mpfr_div_2ui(evaluationMp, evaluationMp, 1, MPFR_RNDN);
|
200 | 140 | storres | mpfr_clear(evaluationTmpMp1); |
201 | 140 | storres | mpfr_clear(evaluationTmpMp2); |
202 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
203 | 147 | storres | /* It may happen, in this case, when the bounds are -Infty and
|
204 | 147 | storres | * +Infty, that the average is NaN. */
|
205 | 147 | storres | if (mpfr_nan_p(evaluationMp))
|
206 | 147 | storres | { |
207 | 147 | storres | return POBYSO_NAN;
|
208 | 147 | storres | } |
209 | 147 | storres | else
|
210 | 147 | storres | { |
211 | 147 | storres | return POBYSO_UNFAITHFUL;
|
212 | 147 | storres | } |
213 | 140 | storres | } |
214 | 140 | storres | else /* Could not get the values of the bounds. */ |
215 | 140 | storres | { |
216 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
217 | 140 | storres | return 1; |
218 | 140 | storres | } |
219 | 140 | storres | } |
220 | 140 | storres | else /* Could not get the precision of the range. */ |
221 | 140 | storres | { |
222 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
223 | 140 | storres | return 1; |
224 | 140 | storres | } |
225 | 140 | storres | } /* End the evaluation is a range. */
|
226 | 140 | storres | /* From now on, we assume that the evaluation is constant. */
|
227 | 147 | storres | if (sollya_lib_get_prec_of_constant(&evalPrec, evaluationSo))
|
228 | 140 | storres | { |
229 | 140 | storres | mpfr_init2(evaluationTmpMp1, evalPrec); |
230 | 140 | storres | if (sollya_lib_get_constant(evaluationTmpMp1, evaluationSo))
|
231 | 140 | storres | { |
232 | 140 | storres | /* Deal with the NaN case. */
|
233 | 140 | storres | if (mpfr_nan_p(evaluationTmpMp1))
|
234 | 140 | storres | { |
235 | 140 | storres | mpfr_clear(evaluationTmpMp1); |
236 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
237 | 140 | storres | return POBYSO_NAN;
|
238 | 140 | storres | } |
239 | 140 | storres | else /* The evaluation is not NaN. */ |
240 | 140 | storres | { |
241 | 140 | storres | mpfr_clear(evaluationMp); |
242 | 140 | storres | mpfr_init2(evaluationMp, evalPrec); |
243 | 140 | storres | mpfr_set(evaluationMp, evaluationTmpMp1, MPFR_RNDN); |
244 | 140 | storres | mpfr_clear(evaluationTmpMp1); |
245 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
246 | 140 | storres | return 0; |
247 | 140 | storres | } |
248 | 140 | storres | } |
249 | 140 | storres | else /* Could not recover the value of the evaluation. */ |
250 | 140 | storres | { |
251 | 140 | storres | mpfr_clear(evaluationTmpMp1); |
252 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
253 | 140 | storres | return 1; |
254 | 140 | storres | } |
255 | 140 | storres | } |
256 | 140 | storres | else /* Could not get the precision of the evaluation. */ |
257 | 140 | storres | { |
258 | 140 | storres | sollya_lib_clear_obj(evaluationSo); |
259 | 140 | storres | return 0; |
260 | 140 | storres | } |
261 | 139 | storres | } |
262 | 139 | storres | /* End pobyso_evaluate_constant. */
|
263 | 139 | storres | |
264 | 134 | storres | /* @see pobyso.h#pobyso_get_verbosity */
|
265 | 134 | storres | int
|
266 | 134 | storres | pobyso_get_verbosity() |
267 | 134 | storres | { |
268 | 134 | storres | sollya_obj_t verbositySo = NULL;
|
269 | 134 | storres | int verbosity = 0; |
270 | 134 | storres | |
271 | 134 | storres | verbositySo = sollya_lib_get_verbosity(); |
272 | 134 | storres | sollya_lib_get_constant_as_int(&verbosity,verbositySo); |
273 | 134 | storres | sollya_lib_clear_obj(verbositySo); |
274 | 134 | storres | return verbosity;
|
275 | 134 | storres | } /* End pobyso_get_verbosity. */
|
276 | 134 | storres | |
277 | 137 | storres | /** @see pobyso.h#pobyso_is_constant_expression
|
278 | 137 | storres | * Strategy: rely on sollya_lib_get_constant. It return 1, when the
|
279 | 137 | storres | * expression is constant.
|
280 | 137 | storres | */
|
281 | 133 | storres | int
|
282 | 133 | storres | pobyso_is_constant_expression(sollya_obj_t obj_to_test) |
283 | 133 | storres | { |
284 | 133 | storres | mpfr_t dummy; |
285 | 133 | storres | int test;
|
286 | 134 | storres | int initial_verbosity_level = 0; |
287 | 134 | storres | |
288 | 133 | storres | /* Test argument. */
|
289 | 133 | storres | if (obj_to_test == NULL) |
290 | 133 | storres | { |
291 | 137 | storres | pobyso_error_message("pobyso_is_constant_expression",
|
292 | 133 | storres | "NULL_POINTER_ARGUMENT",
|
293 | 133 | storres | "The expression is a NULL pointer");
|
294 | 133 | storres | return 0; |
295 | 133 | storres | } |
296 | 134 | storres | initial_verbosity_level = pobyso_set_verbosity_off(); |
297 | 139 | storres | /* In Sollya, constant are functional expressions. */
|
298 | 133 | storres | if (! sollya_lib_obj_is_function(obj_to_test))
|
299 | 133 | storres | { |
300 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
301 | 133 | storres | return 0; |
302 | 133 | storres | } |
303 | 133 | storres | mpfr_init2(dummy,64);
|
304 | 133 | storres | /* Call to previous Sollya function resets verbosity. */
|
305 | 137 | storres | /* Todo: change verbosity suppression strategy with a message call back. */
|
306 | 134 | storres | pobyso_set_verbosity_off(); |
307 | 139 | storres | /* Try to convert the would be constant into an MPFR number. */
|
308 | 139 | storres | /* If OK, we indeed have a constant. If the conversion fails, return 0. */
|
309 | 133 | storres | test = sollya_lib_get_constant(dummy, obj_to_test); |
310 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
311 | 138 | storres | if (test)
|
312 | 137 | storres | { |
313 | 147 | storres | if(mpfr_number_p(dummy) || mpfr_inf_p(dummy))
|
314 | 138 | storres | { |
315 | 138 | storres | mpfr_clear(dummy); |
316 | 138 | storres | return test;
|
317 | 138 | storres | } |
318 | 147 | storres | else /* We do not consider NaNs as constants. */ |
319 | 138 | storres | { |
320 | 138 | storres | mpfr_clear(dummy); |
321 | 138 | storres | return 0; |
322 | 138 | storres | } |
323 | 137 | storres | } |
324 | 138 | storres | else
|
325 | 137 | storres | { |
326 | 139 | storres | mpfr_clear(dummy); |
327 | 137 | storres | return 0; |
328 | 137 | storres | } |
329 | 138 | storres | } /* End pobyso_is_constant_expression. */
|
330 | 137 | storres | |
331 | 137 | storres | /** @see pobyso.h#pobyso_is_monomial. */
|
332 | 137 | storres | int
|
333 | 137 | storres | pobyso_is_int(pobyso_func_exp_t exprSo) |
334 | 137 | storres | { |
335 | 137 | storres | mpfr_t float1M; |
336 | 137 | storres | mpfr_t float2M; |
337 | 137 | storres | mpfr_t tempFloat1M; |
338 | 137 | storres | mpfr_t tempFloat2M; |
339 | 137 | storres | mpfr_prec_t prec; |
340 | 137 | storres | int64_t asInt; |
341 | 137 | storres | sollya_obj_t newConstantSo = NULL;
|
342 | 137 | storres | /* Arguments check. */
|
343 | 137 | storres | if (exprSo == NULL) |
344 | 137 | storres | { |
345 | 137 | storres | pobyso_error_message("pobyso_is_free_var_posze_int_power",
|
346 | 137 | storres | "NULL_POINTER_ARGUMENT",
|
347 | 137 | storres | "The expression is a NULL pointer");
|
348 | 137 | storres | return 0; |
349 | 137 | storres | } |
350 | 137 | storres | //fprintf(stdout, "Not NULL.\n"); pobyso_autoprint(exprSo);
|
351 | 137 | storres | if (! pobyso_is_constant_expression(exprSo)) return 0; |
352 | 137 | storres | if (! sollya_lib_get_constant_as_int64(&asInt, exprSo)) return 0; |
353 | 137 | storres | if (asInt == INT64_MIN || asInt == INT64_MAX) return 0; |
354 | 137 | storres | /* Some constant integer expression can't have their precision computed
|
355 | 137 | storres | * (e.g. cos(pi). */
|
356 | 137 | storres | if (! sollya_lib_get_prec_of_constant(&prec, exprSo))
|
357 | 137 | storres | { |
358 | 137 | storres | mpfr_init2(tempFloat1M, 165);
|
359 | 137 | storres | mpfr_init2(tempFloat2M, 165);
|
360 | 137 | storres | mpfr_abs(tempFloat1M, tempFloat1M, MPFR_RNDN); |
361 | 137 | storres | mpfr_log2(tempFloat2M, tempFloat1M, MPFR_RNDU); |
362 | 137 | storres | mpfr_rint_ceil(tempFloat1M, tempFloat2M, MPFR_RNDU); |
363 | 137 | storres | prec = mpfr_get_si(tempFloat1M, MPFR_RNDN) + 10;
|
364 | 137 | storres | if (prec < 1024) prec = 1024; |
365 | 137 | storres | mpfr_clear(tempFloat1M); |
366 | 137 | storres | mpfr_clear(tempFloat2M); |
367 | 137 | storres | mpfr_init2(float1M, prec); |
368 | 137 | storres | if (!sollya_lib_get_constant(float1M, exprSo))
|
369 | 137 | storres | { |
370 | 137 | storres | mpfr_clear(float1M); |
371 | 137 | storres | return 0; |
372 | 137 | storres | } |
373 | 137 | storres | } |
374 | 137 | storres | else /* Precision could be given. */ |
375 | 137 | storres | { |
376 | 137 | storres | mpfr_init2(float1M, prec); |
377 | 137 | storres | if (! sollya_lib_get_constant(float1M, exprSo))
|
378 | 137 | storres | { |
379 | 137 | storres | mpfr_clear(float1M); |
380 | 137 | storres | return 0; |
381 | 137 | storres | } |
382 | 137 | storres | } |
383 | 137 | storres | if (mpfr_nan_p(float1M) || mpfr_inf_p(float1M))
|
384 | 137 | storres | { |
385 | 137 | storres | mpfr_clear(float1M); |
386 | 137 | storres | return 0; |
387 | 137 | storres | } |
388 | 137 | storres | if ((newConstantSo = sollya_lib_constant_from_int64(asInt)) == NULL) |
389 | 137 | storres | { |
390 | 137 | storres | mpfr_clear(float1M); |
391 | 137 | storres | return 0; |
392 | 137 | storres | } |
393 | 137 | storres | sollya_lib_get_prec_of_constant(&prec, newConstantSo); |
394 | 137 | storres | mpfr_init2(float2M, prec); |
395 | 137 | storres | sollya_lib_get_constant(float2M, newConstantSo); |
396 | 137 | storres | if (mpfr_cmp(float1M, float2M) == 0) |
397 | 137 | storres | { |
398 | 137 | storres | mpfr_clear(float1M); |
399 | 137 | storres | mpfr_clear(float2M); |
400 | 137 | storres | sollya_lib_clear_obj(newConstantSo); |
401 | 137 | storres | return 1; |
402 | 137 | storres | } |
403 | 137 | storres | else
|
404 | 137 | storres | { |
405 | 137 | storres | mpfr_clear(float1M); |
406 | 137 | storres | mpfr_clear(float2M); |
407 | 137 | storres | sollya_lib_clear_obj(newConstantSo); |
408 | 137 | storres | return 0; |
409 | 137 | storres | } |
410 | 137 | storres | } /* End pobyso_is_int. */
|
411 | 137 | storres | |
412 | 137 | storres | /** @see pobyso.h#pobyso_is_monomial.
|
413 | 137 | storres | * Strategy: check that the object is a functional expression and
|
414 | 137 | storres | * if so check that it is made of cte * free_var ^ some_power where :
|
415 | 137 | storres | * - cte is a constant expression (a per pobyso_is_constant_experession;
|
416 | 137 | storres | * - some_power is a positive or null power. t*/
|
417 | 137 | storres | int
|
418 | 137 | storres | pobyso_is_monomial(sollya_obj_t objSo) |
419 | 137 | storres | { |
420 | 138 | storres | int arity;
|
421 | 138 | storres | sollya_obj_t subFun1So = NULL;
|
422 | 138 | storres | sollya_obj_t subFun2So = NULL;
|
423 | 138 | storres | sollya_obj_t subFun3So = NULL;
|
424 | 138 | storres | sollya_base_function_t head = 0;
|
425 | 138 | storres | long int exponent = 0; |
426 | 138 | storres | long int exprIntValue = 0; |
427 | 138 | storres | |
428 | 138 | storres | /* Arguments check. */
|
429 | 138 | storres | if (objSo == NULL) |
430 | 138 | storres | { |
431 | 138 | storres | pobyso_error_message("pobyso_is_monomial",
|
432 | 138 | storres | "NULL_POINTER_ARGUMENT",
|
433 | 138 | storres | "The expression is a NULL pointer");
|
434 | 138 | storres | return 0; |
435 | 138 | storres | } |
436 | 138 | storres | /* The object must be a function. */
|
437 | 138 | storres | if (! sollya_lib_obj_is_function(objSo)) return 0; |
438 | 138 | storres | /* Check if it is the 1 constant. */
|
439 | 138 | storres | if (pobyso_is_int(objSo))
|
440 | 138 | storres | { |
441 | 138 | storres | if (! sollya_lib_get_constant_as_int64(&exprIntValue, objSo))
|
442 | 138 | storres | { |
443 | 138 | storres | return 0; |
444 | 138 | storres | } |
445 | 138 | storres | else
|
446 | 138 | storres | { |
447 | 138 | storres | if (exprIntValue == 1) return 1; |
448 | 138 | storres | else return 0; |
449 | 138 | storres | } |
450 | 138 | storres | } |
451 | 138 | storres | if (! sollya_lib_decompose_function(objSo,
|
452 | 138 | storres | &head, |
453 | 138 | storres | &arity, |
454 | 138 | storres | &subFun1So, |
455 | 138 | storres | &subFun2So, |
456 | 138 | storres | NULL)) return 0; |
457 | 138 | storres | if (arity > 2) |
458 | 138 | storres | { |
459 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
460 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
461 | 138 | storres | return 0; |
462 | 138 | storres | } |
463 | 138 | storres | /* Arity == 1 must be free variable by itself. */
|
464 | 138 | storres | if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE) |
465 | 138 | storres | { |
466 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
467 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
468 | 138 | storres | return 1; |
469 | 138 | storres | } |
470 | 138 | storres | else
|
471 | 138 | storres | { |
472 | 138 | storres | /* Another expression. Must be: free variable ^ poze integer. */
|
473 | 138 | storres | if (arity == 2) |
474 | 138 | storres | { |
475 | 138 | storres | if (head != SOLLYA_BASE_FUNC_POW)
|
476 | 138 | storres | { |
477 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
478 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
479 | 138 | storres | return 0; |
480 | 138 | storres | } |
481 | 138 | storres | if (! pobyso_is_int(subFun2So))
|
482 | 138 | storres | { |
483 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
484 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
485 | 138 | storres | return 0; |
486 | 138 | storres | } |
487 | 138 | storres | if (! sollya_lib_get_constant_as_int64(&exponent, subFun2So))
|
488 | 138 | storres | { |
489 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
490 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
491 | 138 | storres | return 0; |
492 | 138 | storres | } |
493 | 138 | storres | if (exponent < 0) |
494 | 138 | storres | { |
495 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
496 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
497 | 138 | storres | return 0; |
498 | 138 | storres | } |
499 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
500 | 138 | storres | /* Check that the first subfunction is the free variable. */
|
501 | 138 | storres | if (! sollya_lib_decompose_function(subFun1So,
|
502 | 138 | storres | &head, |
503 | 138 | storres | &arity, |
504 | 138 | storres | &subFun2So, |
505 | 138 | storres | &subFun3So, |
506 | 138 | storres | NULL))
|
507 | 138 | storres | { |
508 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
509 | 138 | storres | return 0; |
510 | 138 | storres | } |
511 | 138 | storres | if (arity == 1 && head == SOLLYA_BASE_FUNC_FREE_VARIABLE) |
512 | 138 | storres | { |
513 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
514 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
515 | 138 | storres | if (subFun3So != NULL) sollya_lib_clear_obj(subFun3So); |
516 | 138 | storres | return 1; |
517 | 138 | storres | } |
518 | 138 | storres | else
|
519 | 138 | storres | { |
520 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
521 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
522 | 138 | storres | return 0; |
523 | 138 | storres | } |
524 | 138 | storres | } /* End if arity == 2. */
|
525 | 138 | storres | } /* End else if arity == 1. */
|
526 | 137 | storres | return 0; |
527 | 137 | storres | } /* End pobyso_is_monomial. */
|
528 | 137 | storres | |
529 | 138 | storres | /** @see pobyso.h#pobyso_is_polynomial_term.
|
530 | 138 | storres | * Strategy: check that the object is a functional expression and
|
531 | 138 | storres | * if so check that it is made of cte * monomial.
|
532 | 138 | storres | */
|
533 | 138 | storres | int
|
534 | 138 | storres | pobyso_is_polynomial_term(sollya_obj_t objSo) |
535 | 138 | storres | { |
536 | 138 | storres | int arity;
|
537 | 138 | storres | sollya_obj_t subFun1So = NULL;
|
538 | 138 | storres | sollya_obj_t subFun2So = NULL;
|
539 | 138 | storres | sollya_base_function_t head = 0;
|
540 | 138 | storres | |
541 | 138 | storres | /* Arguments check. */
|
542 | 138 | storres | if (objSo == NULL) |
543 | 138 | storres | { |
544 | 138 | storres | pobyso_error_message("pobyso_is_polynomial_term",
|
545 | 138 | storres | "NULL_POINTER_ARGUMENT",
|
546 | 138 | storres | "The expression is a NULL pointer");
|
547 | 138 | storres | return 0; |
548 | 138 | storres | } |
549 | 138 | storres | /* The object must be a function. */
|
550 | 138 | storres | if (! sollya_lib_obj_is_function(objSo)) return 0; |
551 | 138 | storres | /* A constant expression is degree 0 polynomial term. */
|
552 | 138 | storres | if (pobyso_is_constant_expression(objSo)) return 1; |
553 | 138 | storres | /* A monomial is a polynomial term. */
|
554 | 138 | storres | if (pobyso_is_monomial(objSo)) return 1; |
555 | 138 | storres | /* Decompose the functional expression and study the elements. */
|
556 | 138 | storres | if (! sollya_lib_decompose_function(objSo,
|
557 | 138 | storres | &head, |
558 | 138 | storres | &arity, |
559 | 138 | storres | &subFun1So, |
560 | 138 | storres | &subFun2So, |
561 | 138 | storres | NULL)) return 0; |
562 | 138 | storres | /* Monomial case has been dealt with abobe. */
|
563 | 138 | storres | if (arity != 2) |
564 | 138 | storres | { |
565 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
566 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
567 | 138 | storres | return 0; |
568 | 138 | storres | } |
569 | 138 | storres | /* The expression must be: cte * monomial or monomial * cte. */
|
570 | 138 | storres | if (head != SOLLYA_BASE_FUNC_MUL)
|
571 | 138 | storres | { |
572 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
573 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
574 | 138 | storres | return 0; |
575 | 138 | storres | } |
576 | 138 | storres | if (! pobyso_is_monomial(subFun2So))
|
577 | 138 | storres | { |
578 | 138 | storres | if (! pobyso_is_constant_expression(subFun2So) ||
|
579 | 138 | storres | ! pobyso_is_monomial(subFun1So)) |
580 | 138 | storres | { |
581 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
582 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
583 | 138 | storres | return 0; |
584 | 138 | storres | } |
585 | 138 | storres | } |
586 | 138 | storres | else
|
587 | 138 | storres | { |
588 | 138 | storres | if (! pobyso_is_constant_expression(subFun1So))
|
589 | 138 | storres | { |
590 | 138 | storres | if (subFun1So != NULL) sollya_lib_clear_obj(subFun1So); |
591 | 138 | storres | if (subFun2So != NULL) sollya_lib_clear_obj(subFun2So); |
592 | 138 | storres | return 0; |
593 | 138 | storres | } |
594 | 138 | storres | } |
595 | 138 | storres | return 1; |
596 | 138 | storres | } /* End pobyso_is_polynomial_term. */
|
597 | 133 | storres | /** @see pobyso.h#pobyso_new_monomial. */
|
598 | 133 | storres | pobyso_func_exp_t |
599 | 134 | storres | pobyso_new_monomial(pobyso_func_exp_t coefficientSo, long degree)
|
600 | 133 | storres | { |
601 | 134 | storres | sollya_obj_t degreeSo = NULL;
|
602 | 134 | storres | sollya_obj_t varToPowSo = NULL;
|
603 | 134 | storres | sollya_obj_t monomialSo = NULL;
|
604 | 134 | storres | int initial_verbosity_level = 0; |
605 | 134 | storres | |
606 | 134 | storres | /* Arguments check. */
|
607 | 134 | storres | if (coefficientSo == NULL) |
608 | 133 | storres | { |
609 | 133 | storres | pobyso_error_message("pobyso_parse_string",
|
610 | 133 | storres | "NULL_POINTER_ARGUMENT",
|
611 | 133 | storres | "The expression is a NULL pointer");
|
612 | 134 | storres | return NULL; |
613 | 133 | storres | } |
614 | 134 | storres | if (! pobyso_is_constant_expression(coefficientSo))
|
615 | 133 | storres | { |
616 | 134 | storres | return NULL; |
617 | 133 | storres | } |
618 | 133 | storres | if (degree < 0) |
619 | 133 | storres | { |
620 | 134 | storres | pobyso_error_message("pobyso_new_monomial",
|
621 | 134 | storres | "NEGATIVE_DEGREE_ARGUMENT",
|
622 | 134 | storres | "The degree is a negative integer");
|
623 | 134 | storres | return NULL; |
624 | 133 | storres | } |
625 | 134 | storres | /* If degree == 0, just return a copy of the coefficient. Do not
|
626 | 134 | storres | * return the coefficient itself to avoid "double clear" issues. */
|
627 | 134 | storres | if (degree == 0) |
628 | 134 | storres | { |
629 | 134 | storres | initial_verbosity_level = pobyso_set_verbosity_off(); |
630 | 134 | storres | monomialSo = sollya_lib_copy_obj(coefficientSo); |
631 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
632 | 134 | storres | } |
633 | 134 | storres | degreeSo = sollya_lib_constant_from_int64(degree); |
634 | 134 | storres | varToPowSo = sollya_lib_build_function_pow(sollya_lib_free_variable(), |
635 | 134 | storres | degreeSo); |
636 | 134 | storres | /* Do not use the "build" version to avoid "eating up" the coefficient. */
|
637 | 134 | storres | monomialSo = sollya_lib_mul(coefficientSo,varToPowSo); |
638 | 134 | storres | sollya_lib_clear_obj(varToPowSo); |
639 | 134 | storres | /* Do not clear degreeSa: it was "eaten" by sollya-lib_build_function. */
|
640 | 134 | storres | return monomialSo;
|
641 | 133 | storres | } /* End pobyso_new_monomial. */
|
642 | 133 | storres | |
643 | 149 | storres | /* @see pobyso.h#pobyso_off */
|
644 | 149 | storres | pobyso_on_off_t |
645 | 149 | storres | pobyso_off() |
646 | 149 | storres | { |
647 | 149 | storres | return sollya_lib_off();
|
648 | 149 | storres | } /* End pobyso_off. */
|
649 | 149 | storres | |
650 | 149 | storres | /* @see pobyso.h#pobyso_off */
|
651 | 149 | storres | pobyso_on_off_t |
652 | 149 | storres | pobyso_on() |
653 | 149 | storres | { |
654 | 149 | storres | return sollya_lib_on();
|
655 | 149 | storres | } /* End pobyso_on. */
|
656 | 149 | storres | |
657 | 149 | storres | |
658 | 149 | storres | /* @see pobyso.h#pobyso_parse_string */
|
659 | 130 | storres | pobyso_func_exp_t |
660 | 26 | storres | pobyso_parse_string(const char* expression) |
661 | 26 | storres | { |
662 | 128 | storres | int expressionLength, i;
|
663 | 127 | storres | char *expressionWithSemiCo;
|
664 | 136 | storres | sollya_obj_t expressionSo; |
665 | 128 | storres | |
666 | 127 | storres | /* Arguments check. */
|
667 | 26 | storres | if (expression == NULL) |
668 | 26 | storres | { |
669 | 26 | storres | pobyso_error_message("pobyso_parse_string",
|
670 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
671 | 26 | storres | "The expression is a NULL pointer");
|
672 | 134 | storres | return NULL; |
673 | 26 | storres | } |
674 | 127 | storres | expressionLength = strlen(expression); |
675 | 127 | storres | if (expressionLength == 0) |
676 | 127 | storres | { |
677 | 127 | storres | pobyso_error_message("pobyso_parse_string",
|
678 | 127 | storres | "EMPTY_STRING_ARGUMENT",
|
679 | 127 | storres | "The expression is an empty string");
|
680 | 134 | storres | return NULL; |
681 | 127 | storres | } |
682 | 128 | storres | /* Search from the last char of the expression until, whichever happens
|
683 | 128 | storres | * first:
|
684 | 128 | storres | * a ";" is found;
|
685 | 128 | storres | * a char != ';' is found the the ";" is appended.
|
686 | 128 | storres | * If the head of the string is reached before any of these two events happens
|
687 | 128 | storres | * return an error.
|
688 | 128 | storres | */
|
689 | 128 | storres | for (i = expressionLength - 1 ; i >= 0 ; i--) |
690 | 127 | storres | { |
691 | 128 | storres | if (expression[i] == ';') /* Nothing special to do: |
692 | 128 | storres | try to parse the string*/
|
693 | 127 | storres | { |
694 | 136 | storres | expressionSo = sollya_lib_parse_string(expression); |
695 | 136 | storres | if (sollya_lib_obj_is_error(expressionSo))
|
696 | 136 | storres | { |
697 | 136 | storres | sollya_lib_clear_obj(expressionSo); |
698 | 136 | storres | return NULL; |
699 | 136 | storres | } |
700 | 136 | storres | else
|
701 | 136 | storres | { |
702 | 136 | storres | return expressionSo;
|
703 | 136 | storres | } |
704 | 127 | storres | } |
705 | 128 | storres | else
|
706 | 128 | storres | { |
707 | 128 | storres | if (expression[i] == ' ' || expression[i] == '\t') /* A blank, |
708 | 128 | storres | just continue. */
|
709 | 128 | storres | { |
710 | 128 | storres | continue;
|
711 | 128 | storres | } |
712 | 128 | storres | else /* a character != ';' and from a blank: create the ';'ed string. */ |
713 | 128 | storres | { |
714 | 128 | storres | /* Create a new string for the expression, add the ";" and
|
715 | 128 | storres | * and call sollya_lib_parse_string. */
|
716 | 128 | storres | expressionWithSemiCo = calloc(i + 3, sizeof(char)); |
717 | 128 | storres | if (expressionWithSemiCo == NULL) |
718 | 128 | storres | { |
719 | 128 | storres | pobyso_error_message("pobyso_parse_string",
|
720 | 128 | storres | "MEMORY_ALLOCATION_ERROR",
|
721 | 128 | storres | "Could not allocate the expression string");
|
722 | 134 | storres | return NULL; |
723 | 128 | storres | } |
724 | 132 | storres | strncpy(expressionWithSemiCo, expression, i+1);
|
725 | 128 | storres | expressionWithSemiCo[i + 1] = ';'; |
726 | 128 | storres | expressionWithSemiCo[i + 2] = '\0'; |
727 | 136 | storres | expressionSo = sollya_lib_parse_string(expressionWithSemiCo); |
728 | 128 | storres | free(expressionWithSemiCo); |
729 | 136 | storres | if (sollya_lib_obj_is_error(expressionSo))
|
730 | 136 | storres | { |
731 | 136 | storres | sollya_lib_clear_obj(expressionSo); |
732 | 136 | storres | return NULL; |
733 | 136 | storres | } |
734 | 136 | storres | else
|
735 | 136 | storres | { |
736 | 136 | storres | return expressionSo;
|
737 | 136 | storres | } |
738 | 128 | storres | } /* End character != ';' and from a blank. */
|
739 | 128 | storres | /* Create a new string for the expression, add the ";" and
|
740 | 128 | storres | * and call sollya_lib_parse_string. */
|
741 | 128 | storres | } /* End else. */
|
742 | 128 | storres | } /* End for. */
|
743 | 128 | storres | /* We get here, it is because we did not find any char == anything different
|
744 | 128 | storres | * from ' ' or '\t': the string is empty.
|
745 | 128 | storres | */
|
746 | 128 | storres | pobyso_error_message("pobyso_parse_string",
|
747 | 128 | storres | "ONLY_BLANK_ARGUMENT",
|
748 | 128 | storres | "The expression is only made of blanks");
|
749 | 134 | storres | return NULL; |
750 | 26 | storres | } /* pobyso_parse_string */
|
751 | 26 | storres | |
752 | 150 | storres | pobyso_constant_t |
753 | 150 | storres | pobyso_quiet() |
754 | 150 | storres | { |
755 | 150 | storres | pobyso_constant_t verbositySo = NULL;
|
756 | 150 | storres | pobyso_constant_t zeroSo = NULL;
|
757 | 150 | storres | |
758 | 150 | storres | verbositySo = sollya_lib_get_verbosity(); |
759 | 150 | storres | zeroSo = sollya_lib_constant_from_int64(0);
|
760 | 150 | storres | if (zeroSo != NULL) |
761 | 150 | storres | { |
762 | 150 | storres | sollya_lib_set_verbosity(zeroSo); |
763 | 150 | storres | sollya_lib_clear_obj(zeroSo); |
764 | 150 | storres | return verbositySo;
|
765 | 150 | storres | } |
766 | 150 | storres | else
|
767 | 150 | storres | { |
768 | 150 | storres | sollya_lib_clear_obj(verbositySo); |
769 | 150 | storres | return NULL; |
770 | 150 | storres | } |
771 | 150 | storres | |
772 | 150 | storres | } /* End pobyso_quiet. */
|
773 | 150 | storres | |
774 | 149 | storres | /** @see pobyso.h#pobyso_range_from_bounds */
|
775 | 149 | storres | pobyso_range_t |
776 | 149 | storres | pobyso_range_from_bounds(mpfr_t lowerBound, mpfr_t upperBound) |
777 | 149 | storres | { |
778 | 149 | storres | /* Supferficial check of arguments. */
|
779 | 149 | storres | if (mpfr_cmp(lowerBound, upperBound) > 0) |
780 | 149 | storres | { |
781 | 149 | storres | return(NULL); |
782 | 149 | storres | } |
783 | 149 | storres | return(sollya_lib_range_from_bounds(lowerBound, upperBound));
|
784 | 149 | storres | } |
785 | 149 | storres | |
786 | 132 | storres | pobyso_func_exp_t |
787 | 132 | storres | pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function, |
788 | 132 | storres | long int degree, |
789 | 132 | storres | pobyso_range_t interval, |
790 | 132 | storres | pobyso_func_exp_t weight, |
791 | 132 | storres | double quality,
|
792 | 132 | storres | pobyso_range_t bounds) |
793 | 132 | storres | { |
794 | 134 | storres | sollya_obj_t degreeSo = NULL;
|
795 | 134 | storres | sollya_obj_t qualitySo = NULL;
|
796 | 132 | storres | |
797 | 134 | storres | degreeSo = sollya_lib_constant_from_int(degree); |
798 | 134 | storres | qualitySo = sollya_lib_constant_from_double(quality); |
799 | 132 | storres | |
800 | 134 | storres | sollya_lib_clear_obj(degreeSo); |
801 | 134 | storres | sollya_lib_clear_obj(qualitySo); |
802 | 132 | storres | return NULL; |
803 | 132 | storres | } /* End pobyso_remez_canonical_monomials_base. */
|
804 | 132 | storres | |
805 | 132 | storres | int
|
806 | 132 | storres | pobyso_set_canonical_on() |
807 | 132 | storres | { |
808 | 132 | storres | pobyso_on_off_t currentCanonicalModeSo; |
809 | 132 | storres | pobyso_on_off_t on; |
810 | 132 | storres | |
811 | 132 | storres | currentCanonicalModeSo = sollya_lib_get_canonical(); |
812 | 132 | storres | if (sollya_lib_is_on(currentCanonicalModeSo))
|
813 | 132 | storres | { |
814 | 134 | storres | sollya_lib_clear_obj(currentCanonicalModeSo); |
815 | 134 | storres | return POBYSO_ON;
|
816 | 132 | storres | } |
817 | 132 | storres | else
|
818 | 132 | storres | { |
819 | 134 | storres | on = sollya_lib_on(); |
820 | 134 | storres | sollya_lib_set_canonical(on); |
821 | 134 | storres | sollya_lib_clear_obj(on); |
822 | 134 | storres | sollya_lib_clear_obj(currentCanonicalModeSo); |
823 | 134 | storres | return POBYSO_OFF;
|
824 | 132 | storres | } |
825 | 132 | storres | } /* End pobyso_set_canonical_on. */
|
826 | 132 | storres | |
827 | 134 | storres | int
|
828 | 128 | storres | pobyso_set_verbosity_off() |
829 | 26 | storres | { |
830 | 134 | storres | sollya_obj_t verbosityLevelZeroSo; |
831 | 134 | storres | sollya_obj_t currentVerbosityLevelSo = NULL;
|
832 | 134 | storres | int currentVerbosityLevel = 0; |
833 | 128 | storres | |
834 | 134 | storres | currentVerbosityLevelSo = sollya_lib_get_verbosity(); |
835 | 134 | storres | sollya_lib_get_constant_as_int(¤tVerbosityLevel, |
836 | 134 | storres | currentVerbosityLevelSo); |
837 | 134 | storres | verbosityLevelZeroSo = sollya_lib_constant_from_int(0);
|
838 | 134 | storres | sollya_lib_set_verbosity(verbosityLevelZeroSo); |
839 | 134 | storres | sollya_lib_clear_obj(verbosityLevelZeroSo); |
840 | 134 | storres | sollya_lib_clear_obj(currentVerbosityLevelSo); |
841 | 134 | storres | return currentVerbosityLevel;
|
842 | 26 | storres | } /* End of pobyso_set_verbosity_off. */
|
843 | 26 | storres | |
844 | 134 | storres | int
|
845 | 134 | storres | pobyso_set_verbosity_to(int newVerbosityLevel)
|
846 | 26 | storres | { |
847 | 134 | storres | int initialVerbosityLevel = 0; |
848 | 134 | storres | sollya_obj_t initialVerbosityLevelSo = NULL;
|
849 | 134 | storres | sollya_obj_t newVerbosityLevelSo = NULL;
|
850 | 134 | storres | |
851 | 134 | storres | initialVerbosityLevelSo = sollya_lib_get_verbosity(); |
852 | 134 | storres | sollya_lib_get_constant_as_int(&initialVerbosityLevel, |
853 | 134 | storres | initialVerbosityLevelSo); |
854 | 134 | storres | sollya_lib_clear_obj(initialVerbosityLevelSo); |
855 | 134 | storres | if (newVerbosityLevel < 0) |
856 | 26 | storres | { |
857 | 26 | storres | pobyso_error_message("pobyso_set_verbosity_to",
|
858 | 134 | storres | "INVALID_VALUE",
|
859 | 134 | storres | "The new verbosity level is a negative number");
|
860 | 134 | storres | return initialVerbosityLevel;
|
861 | 26 | storres | } |
862 | 134 | storres | newVerbosityLevelSo = sollya_lib_constant_from_int(newVerbosityLevel); |
863 | 134 | storres | sollya_lib_set_verbosity(newVerbosityLevelSo); |
864 | 134 | storres | sollya_lib_clear_obj(newVerbosityLevelSo); |
865 | 134 | storres | return initialVerbosityLevel;
|
866 | 26 | storres | } /* End of pobyso_set_verbosity_to. */
|
867 | 26 | storres | |
868 | 134 | storres | /**
|
869 | 134 | storres | * @see pobyso.h#pobyso_subpoly
|
870 | 134 | storres | */
|
871 | 134 | storres | pobyso_func_exp_t |
872 | 136 | storres | pobyso_subpoly(pobyso_func_exp_t polynomialSo, long expsNum, long int* expsList) |
873 | 132 | storres | { |
874 | 136 | storres | sollya_obj_t expsListSo = NULL;
|
875 | 136 | storres | sollya_obj_t* expsList_pso = NULL;
|
876 | 136 | storres | sollya_obj_t subpoly = NULL;
|
877 | 136 | storres | int i, j;
|
878 | 134 | storres | |
879 | 134 | storres | /* Arguments check. */
|
880 | 134 | storres | if (polynomialSo == NULL) return NULL; |
881 | 134 | storres | if (expsNum < 0) return NULL; |
882 | 134 | storres | if (expsNum == 0) return sollya_lib_copy_obj(polynomialSo); |
883 | 136 | storres | if (expsList == 0) return NULL; |
884 | 136 | storres | /* Create a list of Sollya constants. */
|
885 | 136 | storres | expsList_pso = (sollya_obj_t*) malloc(expsNum * sizeof(sollya_obj_t));
|
886 | 136 | storres | if (expsList_pso == NULL) |
887 | 132 | storres | { |
888 | 134 | storres | pobyso_error_message("pobyso_subpoly",
|
889 | 134 | storres | "MEMORY_ALLOCATION_ERROR",
|
890 | 136 | storres | "Could not allocate the Sollya exponents list");
|
891 | 134 | storres | return NULL; |
892 | 132 | storres | } |
893 | 136 | storres | /* Fill up the list. */
|
894 | 134 | storres | for (i = 0 ; i < expsNum ; i++) |
895 | 134 | storres | { |
896 | 136 | storres | /* Abort if an exponent is negative. */
|
897 | 136 | storres | if (expsList[i] < 0 ) |
898 | 136 | storres | { |
899 | 136 | storres | for (j = 0 ; j < i ; j++) |
900 | 136 | storres | { |
901 | 136 | storres | sollya_lib_clear_obj(expsList_pso[j]); |
902 | 136 | storres | } |
903 | 136 | storres | free(expsList_pso); |
904 | 136 | storres | return NULL; |
905 | 136 | storres | } |
906 | 136 | storres | expsList_pso[i] = sollya_lib_constant_from_int64(expsList[i]); |
907 | 134 | storres | } /* End for */
|
908 | 136 | storres | expsListSo = sollya_lib_list(expsList_pso, expsNum); |
909 | 136 | storres | for (i = 0 ; i < expsNum ; i++) |
910 | 134 | storres | { |
911 | 136 | storres | sollya_lib_clear_obj(expsList_pso[i]); |
912 | 136 | storres | } |
913 | 136 | storres | free(expsList_pso); |
914 | 136 | storres | if (expsListSo == NULL) |
915 | 136 | storres | { |
916 | 134 | storres | pobyso_error_message("pobyso_subpoly",
|
917 | 134 | storres | "LIST_CREATIONERROR",
|
918 | 134 | storres | "Could not create the exponents list");
|
919 | 134 | storres | return NULL; |
920 | 134 | storres | } |
921 | 134 | storres | subpoly = sollya_lib_subpoly(polynomialSo, expsListSo); |
922 | 136 | storres | pobyso_autoprint(expsListSo); |
923 | 134 | storres | sollya_lib_clear_obj(expsListSo); |
924 | 134 | storres | return subpoly;
|
925 | 134 | storres | } /* pobyso_subpoly. */
|
926 | 132 | storres | |
927 | 26 | storres | /* Attic from the sollya_lib < 4. */
|
928 | 26 | storres | #if 0
|
929 | 26 | storres | chain*
|
930 | 26 | storres | pobyso_create_canonical_monomials_base(const unsigned int degree)
|
931 | 26 | storres | {
|
932 | 26 | storres | int i = 0;
|
933 | 26 | storres | chain *monomials = NULL;
|
934 | 26 | storres | node *monomial = NULL;
|
935 | 26 | storres | |
936 | 26 | storres | for(i = degree ; i >= 0 ; i--)
|
937 | 26 | storres | {
|
938 | 26 | storres | monomial = makePow(makeVariable(), makeConstantDouble((double)i));
|
939 | 26 | storres | monomials = addElement(monomials, monomial);
|
940 | 26 | storres | fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
|
941 | 26 | storres | }
|
942 | 26 | storres | if (monomials == NULL)
|
943 | 26 | storres | {
|
944 | 26 | storres | pobyso_error_message("pobyso_create_canonical_monomial_base",
|
945 | 26 | storres | "CHAIN_CREATION_ERROR",
|
946 | 26 | storres | "Could not create the monomials chain");
|
947 | 26 | storres | return(NULL);
|
948 | 26 | storres | }
|
949 | 26 | storres | return(monomials);
|
950 | 26 | storres | } /* End pobyso_create_canonical_monomials_base. */
|
951 | 26 | storres | |
952 | 26 | storres | chain*
|
953 | 26 | storres | pobyso_create_chain_from_int_array(int* intArray,
|
954 | 26 | storres | const unsigned int arrayLength)
|
955 | 26 | storres | {
|
956 | 26 | storres | int i = 0;
|
957 | 26 | storres | chain *newChain = NULL;
|
958 | 26 | storres | int *currentInt;
|
959 | 26 | storres | |
960 | 26 | storres | if (arrayLength == 0) return(NULL);
|
961 | 26 | storres | if (intArray == NULL)
|
962 | 26 | storres | {
|
963 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
964 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
965 | 26 | storres | "The array is a NULL pointer");
|
966 | 26 | storres | return(NULL);
|
967 | 26 | storres | }
|
968 | 26 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
969 | 26 | storres | {
|
970 | 26 | storres | currentInt = malloc(sizeof(int));
|
971 | 26 | storres | if (currentInt == NULL)
|
972 | 26 | storres | {
|
973 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
974 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
975 | 26 | storres | "Could not allocate one of the integers");
|
976 | 26 | storres | freeChain(newChain, free);
|
977 | 26 | storres | return(NULL);
|
978 | 26 | storres | }
|
979 | 26 | storres | *currentInt = intArray[i];
|
980 | 26 | storres | newChain = addElement(newChain, currentInt);
|
981 | 26 | storres | }
|
982 | 26 | storres | return(newChain);
|
983 | 26 | storres | } // End pobyso_create_chain_from_int_array. */
|
984 | 26 | storres | |
985 | 26 | storres | chain*
|
986 | 26 | storres | pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
|
987 | 26 | storres | const unsigned int arrayLength)
|
988 | 26 | storres | {
|
989 | 26 | storres | int i = 0;
|
990 | 26 | storres | chain *newChain = NULL;
|
991 | 26 | storres | unsigned int *currentInt;
|
992 | 26 | storres | |
993 | 26 | storres | /* Argument checking. */
|
994 | 26 | storres | if (arrayLength == 0) return(NULL);
|
995 | 26 | storres | if (intArray == NULL)
|
996 | 26 | storres | {
|
997 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
998 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
999 | 26 | storres | "The array is a NULL pointer");
|
1000 | 26 | storres | return(NULL);
|
1001 | 26 | storres | }
|
1002 | 26 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
1003 | 26 | storres | {
|
1004 | 26 | storres | currentInt = malloc(sizeof(unsigned int));
|
1005 | 26 | storres | if (currentInt == NULL)
|
1006 | 26 | storres | {
|
1007 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
1008 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
1009 | 26 | storres | "Could not allocate one of the integers");
|
1010 | 26 | storres | freeChain(newChain, free);
|
1011 | 26 | storres | return(NULL);
|
1012 | 26 | storres | }
|
1013 | 26 | storres | *currentInt = intArray[i];
|
1014 | 26 | storres | newChain = addElement(newChain, currentInt);
|
1015 | 26 | storres | }
|
1016 | 26 | storres | return(newChain);
|
1017 | 26 | storres | } // End pobyso_create_chain_from_unsigned_int_array. */
|
1018 | 26 | storres | |
1019 | 26 | storres | node*
|
1020 | 26 | storres | pobyso_differentiate(node *functionNode)
|
1021 | 26 | storres | {
|
1022 | 26 | storres | /* Argument checking. */
|
1023 | 26 | storres | node *differentialNode;
|
1024 | 26 | storres | if (functionNode == NULL)
|
1025 | 26 | storres | {
|
1026 | 26 | storres | pobyso_error_message("pobyso_differentiate",
|
1027 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1028 | 26 | storres | "The function to differentiate is a NULL pointer");
|
1029 | 26 | storres | return(NULL);
|
1030 | 26 | storres | }
|
1031 | 26 | storres | differentialNode = differentiate(functionNode);
|
1032 | 26 | storres | if (differentialNode == NULL)
|
1033 | 26 | storres | {
|
1034 | 26 | storres | pobyso_error_message("pobyso_differentiate",
|
1035 | 26 | storres | "INTERNAL ERROR",
|
1036 | 26 | storres | "Sollya could not differentiate the function");
|
1037 | 26 | storres | }
|
1038 | 26 | storres | return(differentialNode);
|
1039 | 26 | storres | } // End pobyso_differentiate
|
1040 | 26 | storres | |
1041 | 26 | storres | |
1042 | 26 | storres | int
|
1043 | 26 | storres | pobyso_dirty_infnorm(mpfr_t infNorm,
|
1044 | 26 | storres | node *functionNode,
|
1045 | 26 | storres | mpfr_t lowerBound,
|
1046 | 26 | storres | mpfr_t upperBound,
|
1047 | 26 | storres | mp_prec_t precision)
|
1048 | 26 | storres | {
|
1049 | 26 | storres | int functionCallResult;
|
1050 | 26 | storres | /* Arguments checking. */
|
1051 | 26 | storres | if (functionNode == NULL)
|
1052 | 26 | storres | {
|
1053 | 26 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
1054 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1055 | 26 | storres | "The function to compute is a NULL pointer");
|
1056 | 26 | storres | return(1);
|
1057 | 26 | storres | }
|
1058 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) > 0)
|
1059 | 26 | storres | {
|
1060 | 26 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
1061 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
1062 | 26 | storres | "The lower bond is greater than the upper bound");
|
1063 | 26 | storres | return(1);
|
1064 | 26 | storres | }
|
1065 | 26 | storres | /* Particular cases. */
|
1066 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) == 0)
|
1067 | 26 | storres | {
|
1068 | 26 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
1069 | 26 | storres | functionNode,
|
1070 | 26 | storres | lowerBound,
|
1071 | 26 | storres | precision);
|
1072 | 26 | storres | return(functionCallResult);
|
1073 | 26 | storres | }
|
1074 | 26 | storres | if (isConstant(functionNode))
|
1075 | 26 | storres | {
|
1076 | 26 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
1077 | 26 | storres | functionNode,
|
1078 | 26 | storres | lowerBound,
|
1079 | 26 | storres | precision);
|
1080 | 26 | storres | if (!functionCallResult)
|
1081 | 26 | storres | {
|
1082 | 26 | storres | mpfr_abs(infNorm, infNorm, MPFR_RNDN);
|
1083 | 26 | storres | }
|
1084 | 26 | storres | return(functionCallResult);
|
1085 | 26 | storres | }
|
1086 | 26 | storres | uncertifiedInfnorm(infNorm,
|
1087 | 26 | storres | functionNode,
|
1088 | 26 | storres | lowerBound,
|
1089 | 26 | storres | upperBound,
|
1090 | 26 | storres | POBYSO_DEFAULT_POINTS,
|
1091 | 26 | storres | precision);
|
1092 | 26 | storres | return(0);
|
1093 | 26 | storres | } /* End pobyso_dirty_infnorm. */
|
1094 | 26 | storres | |
1095 | 26 | storres | int
|
1096 | 26 | storres | pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
|
1097 | 26 | storres | node *nodeToEvaluate,
|
1098 | 26 | storres | mpfr_t argument,
|
1099 | 26 | storres | mpfr_prec_t precision)
|
1100 | 26 | storres | {
|
1101 | 26 | storres | /* Check input arguments. */
|
1102 | 26 | storres | if (nodeToEvaluate == NULL)
|
1103 | 26 | storres | {
|
1104 | 26 | storres | pobyso_error_message("pobyso_evaluate_faithful",
|
1105 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1106 | 26 | storres | "nodeToEvaluate is a NULL pointer");
|
1107 | 26 | storres | return(1);
|
1108 | 26 | storres | }
|
1109 | 26 | storres | evaluateFaithful(faithfulEvaluation,
|
1110 | 26 | storres | nodeToEvaluate,
|
1111 | 26 | storres | argument,
|
1112 | 26 | storres | precision);
|
1113 | 26 | storres | return(0);
|
1114 | 26 | storres | } /* End pobyso_evaluate_faithfull. */
|
1115 | 26 | storres | |
1116 | 26 | storres | chain*
|
1117 | 26 | storres | pobyso_find_zeros(node *function,
|
1118 | 26 | storres | mpfr_t *lower_bound,
|
1119 | 26 | storres | mpfr_t *upper_bound)
|
1120 | 26 | storres | {
|
1121 | 26 | storres | mp_prec_t currentPrecision;
|
1122 | 26 | storres | mpfr_t currentDiameter;
|
1123 | 26 | storres | rangetype bounds;
|
1124 | 26 | storres | |
1125 | 26 | storres | currentPrecision = getToolPrecision();
|
1126 | 26 | storres | mpfr_init2(currentDiameter, currentPrecision);
|
1127 | 26 | storres | |
1128 | 26 | storres | bounds.a = lower_bound;
|
1129 | 26 | storres | bounds.b = upper_bound;
|
1130 | 26 | storres | |
1131 | 26 | storres | if (bounds.a == NULL || bounds.b == NULL)
|
1132 | 26 | storres | {
|
1133 | 26 | storres | pobyso_error_message("pobyso_find_zeros",
|
1134 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
1135 | 26 | storres | "Could not allocate one of the bounds");
|
1136 | 26 | storres | return(NULL);
|
1137 | 26 | storres | }
|
1138 | 26 | storres | return(findZerosFunction(function,
|
1139 | 26 | storres | bounds,
|
1140 | 26 | storres | currentPrecision,
|
1141 | 26 | storres | currentDiameter));
|
1142 | 26 | storres | } /* End pobyso_find_zeros. */
|
1143 | 26 | storres | |
1144 | 26 | storres | void
|
1145 | 26 | storres | pobyso_free_chain_of_nodes(chain *theChainOfNodes)
|
1146 | 26 | storres | {
|
1147 | 26 | storres | node *currentNode = NULL;
|
1148 | 26 | storres | chain *currentChainElement = NULL;
|
1149 | 26 | storres | chain *nextChainElement = NULL;
|
1150 | 26 | storres | |
1151 | 26 | storres | nextChainElement = theChainOfNodes;
|
1152 | 26 | storres | |
1153 | 26 | storres | while(nextChainElement != NULL)
|
1154 | 26 | storres | {
|
1155 | 26 | storres | currentChainElement = nextChainElement;
|
1156 | 26 | storres | currentNode = (node*)(currentChainElement->value);
|
1157 | 26 | storres | nextChainElement = nextChainElement->next;
|
1158 | 26 | storres | free_memory(currentNode);
|
1159 | 26 | storres | free((void*)(currentChainElement));
|
1160 | 26 | storres | }
|
1161 | 26 | storres | } /* End pobyso_free_chain_of_nodes. */
|
1162 | 26 | storres | |
1163 | 26 | storres | void
|
1164 | 26 | storres | pobyso_free_range(rangetype range)
|
1165 | 26 | storres | {
|
1166 | 26 | storres | |
1167 | 26 | storres | mpfr_clear(*(range.a));
|
1168 | 26 | storres | mpfr_clear(*(range.b));
|
1169 | 26 | storres | free(range.a);
|
1170 | 26 | storres | free(range.b);
|
1171 | 26 | storres | } /* End pobyso_free_range. */
|
1172 | 26 | storres | |
1173 | 26 | storres | node*
|
1174 | 26 | storres | pobyso_fp_minimax_canonical_monomials_base(node *function,
|
1175 | 26 | storres | int degree,
|
1176 | 26 | storres | chain *formats,
|
1177 | 26 | storres | chain *points,
|
1178 | 26 | storres | mpfr_t lowerBound,
|
1179 | 26 | storres | mpfr_t upperBound,
|
1180 | 26 | storres | int fpFixedArg,
|
1181 | 26 | storres | int absRel,
|
1182 | 26 | storres | node *constPart,
|
1183 | 26 | storres | node *minimax)
|
1184 | 26 | storres | {
|
1185 | 26 | storres | return(NULL);
|
1186 | 26 | storres | } /* End pobyso_fp_minimax_canonical_monomials_base. */
|
1187 | 26 | storres | |
1188 | 26 | storres | node*
|
1189 | 26 | storres | pobyso_parse_function(char *functionString,
|
1190 | 26 | storres | char *freeVariableNameString)
|
1191 | 26 | storres | {
|
1192 | 26 | storres | if (functionString == NULL || freeVariableNameString == NULL)
|
1193 | 26 | storres | {
|
1194 | 26 | storres | pobyso_error_message("pobyso_parse_function",
|
1195 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1196 | 26 | storres | "One of the arguments is a NULL pointer");
|
1197 | 26 | storres | return(NULL);
|
1198 | 26 | storres | }
|
1199 | 26 | storres | return(parseString(functionString));
|
1200 | 26 | storres | |
1201 | 26 | storres | } /* End pobyso_parse_function */
|
1202 | 26 | storres | |
1203 | 26 | storres | node*
|
1204 | 26 | storres | pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
|
1205 | 26 | storres | unsigned int mode,
|
1206 | 26 | storres | mpfr_t lowerBound,
|
1207 | 26 | storres | mpfr_t upperBound,
|
1208 | 26 | storres | mpfr_t eps)
|
1209 | 26 | storres | {
|
1210 | 26 | storres | node *weight = NULL;
|
1211 | 26 | storres | node *bestApproxPolyNode = NULL;
|
1212 | 26 | storres | node *bestApproxHorner = NULL;
|
1213 | 26 | storres | node *errorNode = NULL;
|
1214 | 26 | storres | rangetype degreeRange;
|
1215 | 26 | storres | mpfr_t quality;
|
1216 | 26 | storres | mpfr_t currentError;
|
1217 | 26 | storres | unsigned int degree;
|
1218 | 26 | storres | |
1219 | 26 | storres | /* Check the parameters. */
|
1220 | 26 | storres | if (functionNode == NULL)
|
1221 | 26 | storres | {
|
1222 | 26 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
1223 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1224 | 26 | storres | "functionNode is a NULL pointer");
|
1225 | 26 | storres | return(NULL);
|
1226 | 26 | storres | }
|
1227 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
1228 | 26 | storres | {
|
1229 | 26 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
1230 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
1231 | 26 | storres | "the lower_bound >= upper_bound");
|
1232 | 26 | storres | return(NULL);
|
1233 | 26 | storres | }
|
1234 | 26 | storres | /* Set the weight. */
|
1235 | 26 | storres | if (mode == POBYSO_ABSOLUTE)
|
1236 | 26 | storres | {
|
1237 | 26 | storres | /* Set the weight to 1 for the ABSOLUTE_MODE. */
|
1238 | 26 | storres | weight = makeConstantDouble(1.0);
|
1239 | 26 | storres | }
|
1240 | 26 | storres | else
|
1241 | 26 | storres | {
|
1242 | 26 | storres | if (mode == POBYSO_RELATIVE)
|
1243 | 26 | storres | {
|
1244 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
1245 | 26 | storres | "NOT_IMPLEMENTED",
|
1246 | 26 | storres | "the search for relative error is not implemented yet");
|
1247 | 26 | storres | return(NULL);
|
1248 | 26 | storres | }
|
1249 | 26 | storres | else
|
1250 | 26 | storres | {
|
1251 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
1252 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
1253 | 26 | storres | "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
|
1254 | 26 | storres | return(NULL);
|
1255 | 26 | storres | }
|
1256 | 26 | storres | }
|
1257 | 26 | storres | //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
|
1258 | 26 | storres | degreeRange = guessDegree(functionNode,
|
1259 | 26 | storres | weight,
|
1260 | 26 | storres | lowerBound,
|
1261 | 26 | storres | upperBound,
|
1262 | 26 | storres | eps,
|
1263 | 26 | storres | POBYSO_GUESS_DEGREE_BOUND);
|
1264 | 26 | storres | degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
|
1265 | 26 | storres | //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
|
1266 | 26 | storres | fprintf(stderr, "Guessed degree: ");
|
1267 | 26 | storres | mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
|
1268 | 26 | storres | fprintf(stderr, " - as int: %u\n", degree);
|
1269 | 26 | storres | /* Reduce the degree by 1 in the foolish hope it could work. */
|
1270 | 26 | storres | if (degree > 0) degree--;
|
1271 | 26 | storres | /* Both elements of degreeRange where "inited" within guessDegree. */
|
1272 | 26 | storres | mpfr_clear(*(degreeRange.a));
|
1273 | 26 | storres | mpfr_clear(*(degreeRange.b));
|
1274 | 26 | storres | free(degreeRange.a);
|
1275 | 26 | storres | free(degreeRange.b);
|
1276 | 26 | storres | /* Mimic the default behavior of interactive Sollya. */
|
1277 | 26 | storres | mpfr_init(quality);
|
1278 | 26 | storres | mpfr_set_d(quality, 1e-5, MPFR_RNDN);
|
1279 | 26 | storres | mpfr_init2(currentError, getToolPrecision());
|
1280 | 26 | storres | mpfr_set_inf(currentError,1);
|
1281 | 26 | storres | |
1282 | 26 | storres | /* Try to refine the initial guess: loop with increasing degrees until
|
1283 | 26 | storres | * we find a satisfactory one. */
|
1284 | 26 | storres | while(mpfr_cmp(currentError, eps) > 0)
|
1285 | 26 | storres | {
|
1286 | 26 | storres | /* Get rid of the previous polynomial, if any. */
|
1287 | 26 | storres | if (bestApproxPolyNode != NULL)
|
1288 | 26 | storres | {
|
1289 | 26 | storres | free_memory(bestApproxPolyNode);
|
1290 | 26 | storres | }
|
1291 | 26 | storres | fprintf(stderr, "Degree: %u\n", degree);
|
1292 | 26 | storres | fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
|
1293 | 26 | storres | /* Try to find a polynomial with the guessed degree. */
|
1294 | 26 | storres | bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
|
1295 | 26 | storres | weight,
|
1296 | 26 | storres | degree,
|
1297 | 26 | storres | lowerBound,
|
1298 | 26 | storres | upperBound,
|
1299 | 26 | storres | quality);
|
1300 | 26 | storres | |
1301 | 26 | storres | if (bestApproxPolyNode == NULL)
|
1302 | 26 | storres | {
|
1303 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
1304 | 26 | storres | "INTERNAL_ERROR",
|
1305 | 26 | storres | "could not compute the bestApproxPolyNode");
|
1306 | 26 | storres | mpfr_clear(currentError);
|
1307 | 26 | storres | mpfr_clear(quality);
|
1308 | 26 | storres | return(NULL);
|
1309 | 26 | storres | }
|
1310 | 26 | storres | |
1311 | 26 | storres | setDisplayMode(DISPLAY_MODE_DECIMAL);
|
1312 | 26 | storres | fprintTree(stderr, bestApproxPolyNode);
|
1313 | 26 | storres | fprintf(stderr, "\n\n");
|
1314 | 26 | storres | |
1315 | 26 | storres | errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
|
1316 | 26 | storres | /* Check the error with the computed polynomial. */
|
1317 | 26 | storres | uncertifiedInfnorm(currentError,
|
1318 | 26 | storres | errorNode,
|
1319 | 26 | storres | lowerBound,
|
1320 | 26 | storres | upperBound,
|
1321 | 26 | storres | POBYSO_INF_NORM_NUM_POINTS,
|
1322 | 26 | storres | getToolPrecision());
|
1323 | 26 | storres | fprintf(stderr, "Inf norm: ");
|
1324 | 26 | storres | mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
|
1325 | 26 | storres | fprintf(stderr, "\n\n");
|
1326 | 26 | storres | /* Free the errorNode but not the bestApproxPolyNode (we need it if
|
1327 | 26 | storres | * we exit the loop at the next iteration). */
|
1328 | 26 | storres | free_memory(errorNode);
|
1329 | 26 | storres | degree++;
|
1330 | 26 | storres | }
|
1331 | 26 | storres | /* Use an intermediate variable, since horner() creates a new node
|
1332 | 26 | storres | * and does not reorder the argument "in place". This allows for the memory
|
1333 | 26 | storres | * reclaim of bestApproxHorner.
|
1334 | 26 | storres | */
|
1335 | 26 | storres | bestApproxHorner = horner(bestApproxPolyNode);
|
1336 | 26 | storres | free_memory(bestApproxPolyNode);
|
1337 | 26 | storres | mpfr_clear(currentError);
|
1338 | 26 | storres | mpfr_clear(quality);
|
1339 | 26 | storres | free_memory(weight);
|
1340 | 26 | storres | return(bestApproxHorner);
|
1341 | 26 | storres | } /* End pobyso_remez_approx_canonical_monomials_base_for_error */
|
1342 | 26 | storres | |
1343 | 26 | storres | node*
|
1344 | 26 | storres | pobyso_remez_canonical_monomials_base(node *function,
|
1345 | 26 | storres | node *weight,
|
1346 | 26 | storres | unsigned int degree,
|
1347 | 26 | storres | mpfr_t lowerBound,
|
1348 | 26 | storres | mpfr_t upperBound,
|
1349 | 26 | storres | mpfr_t quality)
|
1350 | 26 | storres | {
|
1351 | 26 | storres | node *bestApproxPoly = NULL;
|
1352 | 26 | storres | chain *monomials = NULL;
|
1353 | 26 | storres | chain *curMonomial = NULL;
|
1354 | 26 | storres | |
1355 | 26 | storres | mpfr_t satisfying_error;
|
1356 | 26 | storres | mpfr_t target_error;
|
1357 | 26 | storres | |
1358 | 26 | storres | /* Argument checking */
|
1359 | 26 | storres | /* Function tree. */
|
1360 | 26 | storres | if (function == NULL)
|
1361 | 26 | storres | {
|
1362 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
1363 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1364 | 26 | storres | "the \"function\" argument is a NULL pointer");
|
1365 | 26 | storres | return(NULL);
|
1366 | 26 | storres | }
|
1367 | 26 | storres | if (weight == NULL)
|
1368 | 26 | storres | {
|
1369 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
1370 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
1371 | 26 | storres | "the \"weight\" argument is a NULL pointer");
|
1372 | 26 | storres | return(NULL);
|
1373 | 26 | storres | }
|
1374 | 26 | storres | /* Check the bounds. */
|
1375 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
1376 | 26 | storres | {
|
1377 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
1378 | 26 | storres | "INCOHERENT_IMPUT_DATA",
|
1379 | 26 | storres | "the lower_bound >= upper_bound");
|
1380 | 26 | storres | return(NULL);
|
1381 | 26 | storres | }
|
1382 | 26 | storres | /* The quality must be a non null positive number. */
|
1383 | 26 | storres | if (mpfr_sgn(quality) <= 0)
|
1384 | 26 | storres | {
|
1385 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
1386 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
1387 | 26 | storres | "the quality <= 0");
|
1388 | 26 | storres | }
|
1389 | 26 | storres | /* End argument checking. */
|
1390 | 26 | storres | /* Create the monomials nodes chains. */
|
1391 | 26 | storres | monomials = pobyso_create_canonical_monomials_base(degree);
|
1392 | 26 | storres | fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
|
1393 | 26 | storres | if (monomials == NULL || (lengthChain(monomials) != degree + 1))
|
1394 | 26 | storres | {
|
1395 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
1396 | 26 | storres | "CHAIN_CREATION_ERROR",
|
1397 | 26 | storres | "could not create the monomials chain");
|
1398 | 26 | storres | return(NULL);
|
1399 | 26 | storres | }
|
1400 | 26 | storres | curMonomial = monomials;
|
1401 | 26 | storres | |
1402 | 26 | storres | while (curMonomial != NULL)
|
1403 | 26 | storres | {
|
1404 | 26 | storres | fprintf(stderr, "monomial tree: ");
|
1405 | 26 | storres | //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
|
1406 | 26 | storres | fprintTree(stderr, (node*)(curMonomial->value));
|
1407 | 26 | storres | fprintf(stderr, "\n");
|
1408 | 26 | storres | curMonomial = curMonomial->next;
|
1409 | 26 | storres | }
|
1410 | 26 | storres | |
1411 | 26 | storres | /* Deal with NULL weight. */
|
1412 | 26 | storres | if (weight == NULL)
|
1413 | 26 | storres | {
|
1414 | 26 | storres | weight = makeConstantDouble(1.0);
|
1415 | 26 | storres | }
|
1416 | 26 | storres | /* Compute the best polynomial with the required quality.
|
1417 | 26 | storres | The behavior is as if satisfying error and target_error had
|
1418 | 26 | storres | not been used.*/
|
1419 | 26 | storres | mpfr_init(satisfying_error);
|
1420 | 26 | storres | mpfr_init(target_error);
|
1421 | 26 | storres | mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
|
1422 | 26 | storres | mpfr_set_inf(target_error, 1);
|
1423 | 26 | storres | |
1424 | 26 | storres | |
1425 | 26 | storres | fprintf(stderr, "satisfying_error: ");
|
1426 | 26 | storres | mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
|
1427 | 26 | storres | fprintf(stderr, ".\n");
|
1428 | 26 | storres | fprintf(stderr, "target_error: ");
|
1429 | 26 | storres | mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
|
1430 | 26 | storres | fprintf(stderr, ".\n");
|
1431 | 26 | storres | |
1432 | 26 | storres | fprintf(stderr,
|
1433 | 26 | storres | "current precision: %li\n", getToolPrecision());
|
1434 | 26 | storres | /* Call the Sollya function. */
|
1435 | 26 | storres | bestApproxPoly = remez(function,
|
1436 | 26 | storres | weight,
|
1437 | 26 | storres | monomials,
|
1438 | 26 | storres | lowerBound,
|
1439 | 26 | storres | upperBound,
|
1440 | 26 | storres | quality,
|
1441 | 26 | storres | satisfying_error,
|
1442 | 26 | storres | target_error,
|
1443 | 26 | storres | getToolPrecision());
|
1444 | 26 | storres | |
1445 | 26 | storres | mpfr_clear(satisfying_error);
|
1446 | 26 | storres | mpfr_clear(target_error);
|
1447 | 26 | storres | pobyso_free_chain_of_nodes(monomials);
|
1448 | 26 | storres | |
1449 | 26 | storres | return(bestApproxPoly);
|
1450 | 26 | storres | } /* End pobyso_remez_canonical_monomials_base. */
|
1451 | 26 | storres | |
1452 | 26 | storres | #endif
|
1453 | 26 | storres | |
1454 | 26 | storres | void
|
1455 | 26 | storres | pobyso_error_message(char *functionName, char *messageName, char* message) |
1456 | 26 | storres | { |
1457 | 26 | storres | fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
|
1458 | 26 | storres | } /* End pobyso_error_message */ |