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(¤tVerbosityLevel, |
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 */ |