root / pobysoC-4.0 / src / pobyso.c @ 137
Historique | Voir | Annoter | Télécharger (31,75 ko)
1 | 26 | storres | /** @file pobyso.c
|
---|---|---|---|
2 | 137 | storres | * Integration of Sollya to C programs
|
3 | 26 | storres | *
|
4 | 26 | storres | * @author S.T.
|
5 | 26 | storres | * @date 2011-10-12
|
6 | 26 | storres | *
|
7 | 137 | storres | * @todo write pobyso_is_monomial function <br>
|
8 | 137 | storres | * write pobyso_is_free_var_int_poson_power function
|
9 | 26 | storres | */
|
10 | 26 | storres | /******************************************************************************/
|
11 | 26 | storres | /* Headers, applying the "particular to general" convention.*/
|
12 | 26 | storres | |
13 | 26 | storres | #include "pobyso.h" |
14 | 26 | storres | |
15 | 26 | storres | /* includes of local headers */
|
16 | 26 | storres | |
17 | 26 | storres | /* includes of project headers */
|
18 | 26 | storres | |
19 | 26 | storres | /* includes of system headers */
|
20 | 128 | storres | #include <string.h> |
21 | 128 | storres | #include <stdlib.h> |
22 | 128 | storres | #include <stdio.h> |
23 | 26 | storres | |
24 | 26 | storres | /* Other declarations */
|
25 | 26 | storres | |
26 | 26 | storres | /* Internal prototypes */
|
27 | 26 | storres | void
|
28 | 26 | storres | pobyso_error_message(char *functionName, char *messageName, char* message); |
29 | 26 | storres | /* Types, constants and macros definitions */
|
30 | 26 | storres | |
31 | 26 | storres | /* Global variables */
|
32 | 26 | storres | |
33 | 26 | storres | /* Functions */
|
34 | 127 | storres | |
35 | 134 | storres | /* @see pobyso.h#pobyso_autoprint */
|
36 | 26 | storres | void
|
37 | 137 | storres | pobyso_autoprint(sollya_obj_t objSo) |
38 | 26 | storres | { |
39 | 137 | storres | sollya_lib_autoprint(objSo, NULL);
|
40 | 26 | storres | } /* End pobyso_autoprint. */
|
41 | 26 | storres | |
42 | 134 | storres | /* @see pobyso.h#pobyso_get_verbosity */
|
43 | 134 | storres | int
|
44 | 134 | storres | pobyso_get_verbosity() |
45 | 134 | storres | { |
46 | 134 | storres | sollya_obj_t verbositySo = NULL;
|
47 | 134 | storres | int verbosity = 0; |
48 | 134 | storres | |
49 | 134 | storres | verbositySo = sollya_lib_get_verbosity(); |
50 | 134 | storres | sollya_lib_get_constant_as_int(&verbosity,verbositySo); |
51 | 134 | storres | sollya_lib_clear_obj(verbositySo); |
52 | 134 | storres | return verbosity;
|
53 | 134 | storres | } /* End pobyso_get_verbosity. */
|
54 | 134 | storres | |
55 | 137 | storres | /** @see pobyso.h#pobyso_is_constant_expression
|
56 | 137 | storres | * Strategy: rely on sollya_lib_get_constant. It return 1, when the
|
57 | 137 | storres | * expression is constant.
|
58 | 137 | storres | */
|
59 | 133 | storres | int
|
60 | 133 | storres | pobyso_is_constant_expression(sollya_obj_t obj_to_test) |
61 | 133 | storres | { |
62 | 133 | storres | mpfr_t dummy; |
63 | 133 | storres | int test;
|
64 | 134 | storres | int initial_verbosity_level = 0; |
65 | 134 | storres | |
66 | 133 | storres | /* Test argument. */
|
67 | 133 | storres | if (obj_to_test == NULL) |
68 | 133 | storres | { |
69 | 137 | storres | pobyso_error_message("pobyso_is_constant_expression",
|
70 | 133 | storres | "NULL_POINTER_ARGUMENT",
|
71 | 133 | storres | "The expression is a NULL pointer");
|
72 | 133 | storres | return 0; |
73 | 133 | storres | } |
74 | 134 | storres | initial_verbosity_level = pobyso_set_verbosity_off(); |
75 | 128 | storres | |
76 | 133 | storres | if (! sollya_lib_obj_is_function(obj_to_test))
|
77 | 133 | storres | { |
78 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
79 | 133 | storres | return 0; |
80 | 133 | storres | } |
81 | 133 | storres | mpfr_init2(dummy,64);
|
82 | 133 | storres | /* Call to previous Sollya function resets verbosity. */
|
83 | 137 | storres | /* Todo: change verbosity suppression strategy with a message call back. */
|
84 | 134 | storres | pobyso_set_verbosity_off(); |
85 | 133 | storres | test = sollya_lib_get_constant(dummy, obj_to_test); |
86 | 133 | storres | mpfr_clear(dummy); |
87 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
88 | 133 | storres | return test;
|
89 | 133 | storres | } /* End pobyso_is_constant_expression. */
|
90 | 133 | storres | |
91 | 137 | storres | /** @see pobyso.h#pobyso_is_free_var_int_poson_power. */
|
92 | 137 | storres | int
|
93 | 137 | storres | pobyso_is_free_var_posze_int_power(sollya_obj_t objToTestSo) |
94 | 137 | storres | { |
95 | 137 | storres | int arity;
|
96 | 137 | storres | sollya_base_function_t headTypeSo; |
97 | 137 | storres | sollya_obj_t sub1So = NULL;
|
98 | 137 | storres | sollya_obj_t sub2So = NULL;
|
99 | 137 | storres | |
100 | 137 | storres | /* Argument check. */
|
101 | 137 | storres | if (objToTestSo == NULL) |
102 | 137 | storres | { |
103 | 137 | storres | pobyso_error_message("pobyso_is_free_var_int_poson_power",
|
104 | 137 | storres | "NULL_POINTER_ARGUMENT",
|
105 | 137 | storres | "The expression is a NULL pointer");
|
106 | 137 | storres | return 0; |
107 | 137 | storres | } |
108 | 137 | storres | if (! sollya_lib_obj_is_function(objToTestSo))
|
109 | 137 | storres | { |
110 | 137 | storres | return 0; |
111 | 137 | storres | } |
112 | 137 | storres | if (!sollya_lib_get_head_function(&headTypeSo, objToTestSo)) return 0; |
113 | 137 | storres | if (! sollya_lib_get_subfunctions(objToTestSo,
|
114 | 137 | storres | &arity, |
115 | 137 | storres | &sub1So, |
116 | 137 | storres | &sub2So, |
117 | 137 | storres | NULL)) return 0; |
118 | 137 | storres | /* Power function: arity == 2. */
|
119 | 137 | storres | if (arity != 2) return 0; |
120 | 137 | storres | if (! pobyso_is_constant_expression(sub1So)) return 0; |
121 | 137 | storres | return 1; |
122 | 137 | storres | } /* End pobyso_is_free_var_posze_int_power. */
|
123 | 137 | storres | |
124 | 137 | storres | |
125 | 137 | storres | /** @see pobyso.h#pobyso_is_monomial. */
|
126 | 137 | storres | int
|
127 | 137 | storres | pobyso_is_int(pobyso_func_exp_t exprSo) |
128 | 137 | storres | { |
129 | 137 | storres | mpfr_t float1M; |
130 | 137 | storres | mpfr_t float2M; |
131 | 137 | storres | mpfr_t tempFloat1M; |
132 | 137 | storres | mpfr_t tempFloat2M; |
133 | 137 | storres | mpfr_prec_t prec; |
134 | 137 | storres | int64_t asInt; |
135 | 137 | storres | sollya_obj_t newConstantSo = NULL;
|
136 | 137 | storres | /* Arguments check. */
|
137 | 137 | storres | if (exprSo == NULL) |
138 | 137 | storres | { |
139 | 137 | storres | pobyso_error_message("pobyso_is_free_var_posze_int_power",
|
140 | 137 | storres | "NULL_POINTER_ARGUMENT",
|
141 | 137 | storres | "The expression is a NULL pointer");
|
142 | 137 | storres | return 0; |
143 | 137 | storres | } |
144 | 137 | storres | //fprintf(stdout, "Not NULL.\n"); pobyso_autoprint(exprSo);
|
145 | 137 | storres | if (! pobyso_is_constant_expression(exprSo)) return 0; |
146 | 137 | storres | if (! sollya_lib_get_constant_as_int64(&asInt, exprSo)) return 0; |
147 | 137 | storres | if (asInt == INT64_MIN || asInt == INT64_MAX) return 0; |
148 | 137 | storres | /* Some constant integer expression can't have their precision computed
|
149 | 137 | storres | * (e.g. cos(pi). */
|
150 | 137 | storres | if (! sollya_lib_get_prec_of_constant(&prec, exprSo))
|
151 | 137 | storres | { |
152 | 137 | storres | mpfr_init2(tempFloat1M, 165);
|
153 | 137 | storres | mpfr_init2(tempFloat2M, 165);
|
154 | 137 | storres | mpfr_abs(tempFloat1M, tempFloat1M, MPFR_RNDN); |
155 | 137 | storres | mpfr_log2(tempFloat2M, tempFloat1M, MPFR_RNDU); |
156 | 137 | storres | mpfr_rint_ceil(tempFloat1M, tempFloat2M, MPFR_RNDU); |
157 | 137 | storres | prec = mpfr_get_si(tempFloat1M, MPFR_RNDN) + 10;
|
158 | 137 | storres | if (prec < 1024) prec = 1024; |
159 | 137 | storres | mpfr_clear(tempFloat1M); |
160 | 137 | storres | mpfr_clear(tempFloat2M); |
161 | 137 | storres | mpfr_init2(float1M, prec); |
162 | 137 | storres | if (!sollya_lib_get_constant(float1M, exprSo))
|
163 | 137 | storres | { |
164 | 137 | storres | mpfr_clear(float1M); |
165 | 137 | storres | return 0; |
166 | 137 | storres | } |
167 | 137 | storres | } |
168 | 137 | storres | else /* Precision could be given. */ |
169 | 137 | storres | { |
170 | 137 | storres | mpfr_init2(float1M, prec); |
171 | 137 | storres | if (! sollya_lib_get_constant(float1M, exprSo))
|
172 | 137 | storres | { |
173 | 137 | storres | mpfr_clear(float1M); |
174 | 137 | storres | return 0; |
175 | 137 | storres | } |
176 | 137 | storres | } |
177 | 137 | storres | if (mpfr_nan_p(float1M) || mpfr_inf_p(float1M))
|
178 | 137 | storres | { |
179 | 137 | storres | mpfr_clear(float1M); |
180 | 137 | storres | return 0; |
181 | 137 | storres | } |
182 | 137 | storres | if ((newConstantSo = sollya_lib_constant_from_int64(asInt)) == NULL) |
183 | 137 | storres | { |
184 | 137 | storres | mpfr_clear(float1M); |
185 | 137 | storres | return 0; |
186 | 137 | storres | } |
187 | 137 | storres | sollya_lib_get_prec_of_constant(&prec, newConstantSo); |
188 | 137 | storres | mpfr_init2(float2M, prec); |
189 | 137 | storres | sollya_lib_get_constant(float2M, newConstantSo); |
190 | 137 | storres | if (mpfr_cmp(float1M, float2M) == 0) |
191 | 137 | storres | { |
192 | 137 | storres | mpfr_clear(float1M); |
193 | 137 | storres | mpfr_clear(float2M); |
194 | 137 | storres | sollya_lib_clear_obj(newConstantSo); |
195 | 137 | storres | return 1; |
196 | 137 | storres | } |
197 | 137 | storres | else
|
198 | 137 | storres | { |
199 | 137 | storres | pobyso_autoprint(exprSo); |
200 | 137 | storres | pobyso_autoprint(newConstantSo); |
201 | 137 | storres | mpfr_clear(float1M); |
202 | 137 | storres | mpfr_clear(float2M); |
203 | 137 | storres | sollya_lib_clear_obj(newConstantSo); |
204 | 137 | storres | return 0; |
205 | 137 | storres | } |
206 | 137 | storres | } /* End pobyso_is_int. */
|
207 | 137 | storres | |
208 | 137 | storres | /** @see pobyso.h#pobyso_is_monomial.
|
209 | 137 | storres | * Strategy: check that the object is a functional expression and
|
210 | 137 | storres | * if so check that it is made of cte * free_var ^ some_power where :
|
211 | 137 | storres | * - cte is a constant expression (a per pobyso_is_constant_experession;
|
212 | 137 | storres | * - some_power is a positive or null power. t*/
|
213 | 137 | storres | int
|
214 | 137 | storres | pobyso_is_monomial(sollya_obj_t objSo) |
215 | 137 | storres | { |
216 | 137 | storres | return 0; |
217 | 137 | storres | } /* End pobyso_is_monomial. */
|
218 | 137 | storres | |
219 | 133 | storres | /** @see pobyso.h#pobyso_new_monomial. */
|
220 | 133 | storres | pobyso_func_exp_t |
221 | 134 | storres | pobyso_new_monomial(pobyso_func_exp_t coefficientSo, long degree)
|
222 | 133 | storres | { |
223 | 134 | storres | sollya_obj_t degreeSo = NULL;
|
224 | 134 | storres | sollya_obj_t varToPowSo = NULL;
|
225 | 134 | storres | sollya_obj_t monomialSo = NULL;
|
226 | 134 | storres | int initial_verbosity_level = 0; |
227 | 134 | storres | |
228 | 134 | storres | /* Arguments check. */
|
229 | 134 | storres | if (coefficientSo == NULL) |
230 | 133 | storres | { |
231 | 133 | storres | pobyso_error_message("pobyso_parse_string",
|
232 | 133 | storres | "NULL_POINTER_ARGUMENT",
|
233 | 133 | storres | "The expression is a NULL pointer");
|
234 | 134 | storres | return NULL; |
235 | 133 | storres | } |
236 | 134 | storres | if (! pobyso_is_constant_expression(coefficientSo))
|
237 | 133 | storres | { |
238 | 134 | storres | return NULL; |
239 | 133 | storres | } |
240 | 133 | storres | if (degree < 0) |
241 | 133 | storres | { |
242 | 134 | storres | pobyso_error_message("pobyso_new_monomial",
|
243 | 134 | storres | "NEGATIVE_DEGREE_ARGUMENT",
|
244 | 134 | storres | "The degree is a negative integer");
|
245 | 134 | storres | return NULL; |
246 | 133 | storres | } |
247 | 134 | storres | /* If degree == 0, just return a copy of the coefficient. Do not
|
248 | 134 | storres | * return the coefficient itself to avoid "double clear" issues. */
|
249 | 134 | storres | if (degree == 0) |
250 | 134 | storres | { |
251 | 134 | storres | initial_verbosity_level = pobyso_set_verbosity_off(); |
252 | 134 | storres | monomialSo = sollya_lib_copy_obj(coefficientSo); |
253 | 134 | storres | pobyso_set_verbosity_to(initial_verbosity_level); |
254 | 134 | storres | } |
255 | 134 | storres | degreeSo = sollya_lib_constant_from_int64(degree); |
256 | 134 | storres | varToPowSo = sollya_lib_build_function_pow(sollya_lib_free_variable(), |
257 | 134 | storres | degreeSo); |
258 | 134 | storres | /* Do not use the "build" version to avoid "eating up" the coefficient. */
|
259 | 134 | storres | monomialSo = sollya_lib_mul(coefficientSo,varToPowSo); |
260 | 134 | storres | sollya_lib_clear_obj(varToPowSo); |
261 | 134 | storres | /* Do not clear degreeSa: it was "eaten" by sollya-lib_build_function. */
|
262 | 134 | storres | return monomialSo;
|
263 | 133 | storres | } /* End pobyso_new_monomial. */
|
264 | 133 | storres | |
265 | 127 | storres | /* @see pobyso.h#pobyso_parse_string*/
|
266 | 130 | storres | pobyso_func_exp_t |
267 | 26 | storres | pobyso_parse_string(const char* expression) |
268 | 26 | storres | { |
269 | 128 | storres | int expressionLength, i;
|
270 | 127 | storres | char *expressionWithSemiCo;
|
271 | 136 | storres | sollya_obj_t expressionSo; |
272 | 128 | storres | |
273 | 127 | storres | /* Arguments check. */
|
274 | 26 | storres | if (expression == NULL) |
275 | 26 | storres | { |
276 | 26 | storres | pobyso_error_message("pobyso_parse_string",
|
277 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
278 | 26 | storres | "The expression is a NULL pointer");
|
279 | 134 | storres | return NULL; |
280 | 26 | storres | } |
281 | 127 | storres | expressionLength = strlen(expression); |
282 | 127 | storres | if (expressionLength == 0) |
283 | 127 | storres | { |
284 | 127 | storres | pobyso_error_message("pobyso_parse_string",
|
285 | 127 | storres | "EMPTY_STRING_ARGUMENT",
|
286 | 127 | storres | "The expression is an empty string");
|
287 | 134 | storres | return NULL; |
288 | 127 | storres | } |
289 | 128 | storres | /* Search from the last char of the expression until, whichever happens
|
290 | 128 | storres | * first:
|
291 | 128 | storres | * a ";" is found;
|
292 | 128 | storres | * a char != ';' is found the the ";" is appended.
|
293 | 128 | storres | * If the head of the string is reached before any of these two events happens
|
294 | 128 | storres | * return an error.
|
295 | 128 | storres | */
|
296 | 128 | storres | for (i = expressionLength - 1 ; i >= 0 ; i--) |
297 | 127 | storres | { |
298 | 128 | storres | if (expression[i] == ';') /* Nothing special to do: |
299 | 128 | storres | try to parse the string*/
|
300 | 127 | storres | { |
301 | 136 | storres | expressionSo = sollya_lib_parse_string(expression); |
302 | 136 | storres | if (sollya_lib_obj_is_error(expressionSo))
|
303 | 136 | storres | { |
304 | 136 | storres | sollya_lib_clear_obj(expressionSo); |
305 | 136 | storres | return NULL; |
306 | 136 | storres | } |
307 | 136 | storres | else
|
308 | 136 | storres | { |
309 | 136 | storres | return expressionSo;
|
310 | 136 | storres | } |
311 | 127 | storres | } |
312 | 128 | storres | else
|
313 | 128 | storres | { |
314 | 128 | storres | if (expression[i] == ' ' || expression[i] == '\t') /* A blank, |
315 | 128 | storres | just continue. */
|
316 | 128 | storres | { |
317 | 128 | storres | continue;
|
318 | 128 | storres | } |
319 | 128 | storres | else /* a character != ';' and from a blank: create the ';'ed string. */ |
320 | 128 | storres | { |
321 | 128 | storres | /* Create a new string for the expression, add the ";" and
|
322 | 128 | storres | * and call sollya_lib_parse_string. */
|
323 | 128 | storres | expressionWithSemiCo = calloc(i + 3, sizeof(char)); |
324 | 128 | storres | if (expressionWithSemiCo == NULL) |
325 | 128 | storres | { |
326 | 128 | storres | pobyso_error_message("pobyso_parse_string",
|
327 | 128 | storres | "MEMORY_ALLOCATION_ERROR",
|
328 | 128 | storres | "Could not allocate the expression string");
|
329 | 134 | storres | return NULL; |
330 | 128 | storres | } |
331 | 132 | storres | strncpy(expressionWithSemiCo, expression, i+1);
|
332 | 128 | storres | expressionWithSemiCo[i + 1] = ';'; |
333 | 128 | storres | expressionWithSemiCo[i + 2] = '\0'; |
334 | 136 | storres | expressionSo = sollya_lib_parse_string(expressionWithSemiCo); |
335 | 128 | storres | free(expressionWithSemiCo); |
336 | 136 | storres | if (sollya_lib_obj_is_error(expressionSo))
|
337 | 136 | storres | { |
338 | 136 | storres | sollya_lib_clear_obj(expressionSo); |
339 | 136 | storres | return NULL; |
340 | 136 | storres | } |
341 | 136 | storres | else
|
342 | 136 | storres | { |
343 | 136 | storres | return expressionSo;
|
344 | 136 | storres | } |
345 | 128 | storres | } /* End character != ';' and from a blank. */
|
346 | 128 | storres | /* Create a new string for the expression, add the ";" and
|
347 | 128 | storres | * and call sollya_lib_parse_string. */
|
348 | 128 | storres | } /* End else. */
|
349 | 128 | storres | } /* End for. */
|
350 | 128 | storres | /* We get here, it is because we did not find any char == anything different
|
351 | 128 | storres | * from ' ' or '\t': the string is empty.
|
352 | 128 | storres | */
|
353 | 128 | storres | pobyso_error_message("pobyso_parse_string",
|
354 | 128 | storres | "ONLY_BLANK_ARGUMENT",
|
355 | 128 | storres | "The expression is only made of blanks");
|
356 | 134 | storres | return NULL; |
357 | 26 | storres | } /* pobyso_parse_string */
|
358 | 26 | storres | |
359 | 132 | storres | pobyso_func_exp_t |
360 | 132 | storres | pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function, |
361 | 132 | storres | long int degree, |
362 | 132 | storres | pobyso_range_t interval, |
363 | 132 | storres | pobyso_func_exp_t weight, |
364 | 132 | storres | double quality,
|
365 | 132 | storres | pobyso_range_t bounds) |
366 | 132 | storres | { |
367 | 134 | storres | sollya_obj_t degreeSo = NULL;
|
368 | 134 | storres | sollya_obj_t qualitySo = NULL;
|
369 | 132 | storres | |
370 | 134 | storres | degreeSo = sollya_lib_constant_from_int(degree); |
371 | 134 | storres | qualitySo = sollya_lib_constant_from_double(quality); |
372 | 132 | storres | |
373 | 134 | storres | sollya_lib_clear_obj(degreeSo); |
374 | 134 | storres | sollya_lib_clear_obj(qualitySo); |
375 | 132 | storres | return NULL; |
376 | 132 | storres | } /* End pobyso_remez_canonical_monomials_base. */
|
377 | 132 | storres | |
378 | 132 | storres | int
|
379 | 132 | storres | pobyso_set_canonical_on() |
380 | 132 | storres | { |
381 | 132 | storres | pobyso_on_off_t currentCanonicalModeSo; |
382 | 132 | storres | pobyso_on_off_t on; |
383 | 132 | storres | |
384 | 132 | storres | currentCanonicalModeSo = sollya_lib_get_canonical(); |
385 | 132 | storres | if (sollya_lib_is_on(currentCanonicalModeSo))
|
386 | 132 | storres | { |
387 | 134 | storres | sollya_lib_clear_obj(currentCanonicalModeSo); |
388 | 134 | storres | return POBYSO_ON;
|
389 | 132 | storres | } |
390 | 132 | storres | else
|
391 | 132 | storres | { |
392 | 134 | storres | on = sollya_lib_on(); |
393 | 134 | storres | sollya_lib_set_canonical(on); |
394 | 134 | storres | sollya_lib_clear_obj(on); |
395 | 134 | storres | sollya_lib_clear_obj(currentCanonicalModeSo); |
396 | 134 | storres | return POBYSO_OFF;
|
397 | 132 | storres | } |
398 | 132 | storres | } /* End pobyso_set_canonical_on. */
|
399 | 132 | storres | |
400 | 134 | storres | int
|
401 | 128 | storres | pobyso_set_verbosity_off() |
402 | 26 | storres | { |
403 | 134 | storres | sollya_obj_t verbosityLevelZeroSo; |
404 | 134 | storres | sollya_obj_t currentVerbosityLevelSo = NULL;
|
405 | 134 | storres | int currentVerbosityLevel = 0; |
406 | 128 | storres | |
407 | 134 | storres | currentVerbosityLevelSo = sollya_lib_get_verbosity(); |
408 | 134 | storres | sollya_lib_get_constant_as_int(¤tVerbosityLevel, |
409 | 134 | storres | currentVerbosityLevelSo); |
410 | 134 | storres | verbosityLevelZeroSo = sollya_lib_constant_from_int(0);
|
411 | 134 | storres | sollya_lib_set_verbosity(verbosityLevelZeroSo); |
412 | 134 | storres | sollya_lib_clear_obj(verbosityLevelZeroSo); |
413 | 134 | storres | sollya_lib_clear_obj(currentVerbosityLevelSo); |
414 | 134 | storres | return currentVerbosityLevel;
|
415 | 26 | storres | } /* End of pobyso_set_verbosity_off. */
|
416 | 26 | storres | |
417 | 134 | storres | int
|
418 | 134 | storres | pobyso_set_verbosity_to(int newVerbosityLevel)
|
419 | 26 | storres | { |
420 | 134 | storres | int initialVerbosityLevel = 0; |
421 | 134 | storres | sollya_obj_t initialVerbosityLevelSo = NULL;
|
422 | 134 | storres | sollya_obj_t newVerbosityLevelSo = NULL;
|
423 | 134 | storres | |
424 | 134 | storres | initialVerbosityLevelSo = sollya_lib_get_verbosity(); |
425 | 134 | storres | sollya_lib_get_constant_as_int(&initialVerbosityLevel, |
426 | 134 | storres | initialVerbosityLevelSo); |
427 | 134 | storres | sollya_lib_clear_obj(initialVerbosityLevelSo); |
428 | 134 | storres | if (newVerbosityLevel < 0) |
429 | 26 | storres | { |
430 | 26 | storres | pobyso_error_message("pobyso_set_verbosity_to",
|
431 | 134 | storres | "INVALID_VALUE",
|
432 | 134 | storres | "The new verbosity level is a negative number");
|
433 | 134 | storres | return initialVerbosityLevel;
|
434 | 26 | storres | } |
435 | 134 | storres | newVerbosityLevelSo = sollya_lib_constant_from_int(newVerbosityLevel); |
436 | 134 | storres | sollya_lib_set_verbosity(newVerbosityLevelSo); |
437 | 134 | storres | sollya_lib_clear_obj(newVerbosityLevelSo); |
438 | 134 | storres | return initialVerbosityLevel;
|
439 | 26 | storres | } /* End of pobyso_set_verbosity_to. */
|
440 | 26 | storres | |
441 | 134 | storres | /**
|
442 | 134 | storres | * @see pobyso.h#pobyso_subpoly
|
443 | 134 | storres | */
|
444 | 134 | storres | pobyso_func_exp_t |
445 | 136 | storres | pobyso_subpoly(pobyso_func_exp_t polynomialSo, long expsNum, long int* expsList) |
446 | 132 | storres | { |
447 | 136 | storres | sollya_obj_t expsListSo = NULL;
|
448 | 136 | storres | sollya_obj_t* expsList_pso = NULL;
|
449 | 136 | storres | sollya_obj_t subpoly = NULL;
|
450 | 136 | storres | int i, j;
|
451 | 134 | storres | |
452 | 134 | storres | /* Arguments check. */
|
453 | 134 | storres | if (polynomialSo == NULL) return NULL; |
454 | 134 | storres | if (expsNum < 0) return NULL; |
455 | 134 | storres | if (expsNum == 0) return sollya_lib_copy_obj(polynomialSo); |
456 | 136 | storres | if (expsList == 0) return NULL; |
457 | 136 | storres | /* Create a list of Sollya constants. */
|
458 | 136 | storres | expsList_pso = (sollya_obj_t*) malloc(expsNum * sizeof(sollya_obj_t));
|
459 | 136 | storres | if (expsList_pso == NULL) |
460 | 132 | storres | { |
461 | 134 | storres | pobyso_error_message("pobyso_subpoly",
|
462 | 134 | storres | "MEMORY_ALLOCATION_ERROR",
|
463 | 136 | storres | "Could not allocate the Sollya exponents list");
|
464 | 134 | storres | return NULL; |
465 | 132 | storres | } |
466 | 136 | storres | /* Fill up the list. */
|
467 | 134 | storres | for (i = 0 ; i < expsNum ; i++) |
468 | 134 | storres | { |
469 | 136 | storres | /* Abort if an exponent is negative. */
|
470 | 136 | storres | if (expsList[i] < 0 ) |
471 | 136 | storres | { |
472 | 136 | storres | for (j = 0 ; j < i ; j++) |
473 | 136 | storres | { |
474 | 136 | storres | sollya_lib_clear_obj(expsList_pso[j]); |
475 | 136 | storres | } |
476 | 136 | storres | free(expsList_pso); |
477 | 136 | storres | return NULL; |
478 | 136 | storres | } |
479 | 136 | storres | expsList_pso[i] = sollya_lib_constant_from_int64(expsList[i]); |
480 | 134 | storres | } /* End for */
|
481 | 136 | storres | expsListSo = sollya_lib_list(expsList_pso, expsNum); |
482 | 136 | storres | for (i = 0 ; i < expsNum ; i++) |
483 | 134 | storres | { |
484 | 136 | storres | sollya_lib_clear_obj(expsList_pso[i]); |
485 | 136 | storres | } |
486 | 136 | storres | free(expsList_pso); |
487 | 136 | storres | if (expsListSo == NULL) |
488 | 136 | storres | { |
489 | 134 | storres | pobyso_error_message("pobyso_subpoly",
|
490 | 134 | storres | "LIST_CREATIONERROR",
|
491 | 134 | storres | "Could not create the exponents list");
|
492 | 134 | storres | return NULL; |
493 | 134 | storres | } |
494 | 134 | storres | subpoly = sollya_lib_subpoly(polynomialSo, expsListSo); |
495 | 136 | storres | pobyso_autoprint(expsListSo); |
496 | 134 | storres | sollya_lib_clear_obj(expsListSo); |
497 | 134 | storres | return subpoly;
|
498 | 134 | storres | } /* pobyso_subpoly. */
|
499 | 132 | storres | |
500 | 26 | storres | /* Attic from the sollya_lib < 4. */
|
501 | 26 | storres | #if 0
|
502 | 26 | storres | chain*
|
503 | 26 | storres | pobyso_create_canonical_monomials_base(const unsigned int degree)
|
504 | 26 | storres | {
|
505 | 26 | storres | int i = 0;
|
506 | 26 | storres | chain *monomials = NULL;
|
507 | 26 | storres | node *monomial = NULL;
|
508 | 26 | storres | |
509 | 26 | storres | for(i = degree ; i >= 0 ; i--)
|
510 | 26 | storres | {
|
511 | 26 | storres | monomial = makePow(makeVariable(), makeConstantDouble((double)i));
|
512 | 26 | storres | monomials = addElement(monomials, monomial);
|
513 | 26 | storres | fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
|
514 | 26 | storres | }
|
515 | 26 | storres | if (monomials == NULL)
|
516 | 26 | storres | {
|
517 | 26 | storres | pobyso_error_message("pobyso_create_canonical_monomial_base",
|
518 | 26 | storres | "CHAIN_CREATION_ERROR",
|
519 | 26 | storres | "Could not create the monomials chain");
|
520 | 26 | storres | return(NULL);
|
521 | 26 | storres | }
|
522 | 26 | storres | return(monomials);
|
523 | 26 | storres | } /* End pobyso_create_canonical_monomials_base. */
|
524 | 26 | storres | |
525 | 26 | storres | chain*
|
526 | 26 | storres | pobyso_create_chain_from_int_array(int* intArray,
|
527 | 26 | storres | const unsigned int arrayLength)
|
528 | 26 | storres | {
|
529 | 26 | storres | int i = 0;
|
530 | 26 | storres | chain *newChain = NULL;
|
531 | 26 | storres | int *currentInt;
|
532 | 26 | storres | |
533 | 26 | storres | if (arrayLength == 0) return(NULL);
|
534 | 26 | storres | if (intArray == NULL)
|
535 | 26 | storres | {
|
536 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
537 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
538 | 26 | storres | "The array is a NULL pointer");
|
539 | 26 | storres | return(NULL);
|
540 | 26 | storres | }
|
541 | 26 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
542 | 26 | storres | {
|
543 | 26 | storres | currentInt = malloc(sizeof(int));
|
544 | 26 | storres | if (currentInt == NULL)
|
545 | 26 | storres | {
|
546 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
547 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
548 | 26 | storres | "Could not allocate one of the integers");
|
549 | 26 | storres | freeChain(newChain, free);
|
550 | 26 | storres | return(NULL);
|
551 | 26 | storres | }
|
552 | 26 | storres | *currentInt = intArray[i];
|
553 | 26 | storres | newChain = addElement(newChain, currentInt);
|
554 | 26 | storres | }
|
555 | 26 | storres | return(newChain);
|
556 | 26 | storres | } // End pobyso_create_chain_from_int_array. */
|
557 | 26 | storres | |
558 | 26 | storres | chain*
|
559 | 26 | storres | pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
|
560 | 26 | storres | const unsigned int arrayLength)
|
561 | 26 | storres | {
|
562 | 26 | storres | int i = 0;
|
563 | 26 | storres | chain *newChain = NULL;
|
564 | 26 | storres | unsigned int *currentInt;
|
565 | 26 | storres | |
566 | 26 | storres | /* Argument checking. */
|
567 | 26 | storres | if (arrayLength == 0) return(NULL);
|
568 | 26 | storres | if (intArray == NULL)
|
569 | 26 | storres | {
|
570 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
571 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
572 | 26 | storres | "The array is a NULL pointer");
|
573 | 26 | storres | return(NULL);
|
574 | 26 | storres | }
|
575 | 26 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
576 | 26 | storres | {
|
577 | 26 | storres | currentInt = malloc(sizeof(unsigned int));
|
578 | 26 | storres | if (currentInt == NULL)
|
579 | 26 | storres | {
|
580 | 26 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
581 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
582 | 26 | storres | "Could not allocate one of the integers");
|
583 | 26 | storres | freeChain(newChain, free);
|
584 | 26 | storres | return(NULL);
|
585 | 26 | storres | }
|
586 | 26 | storres | *currentInt = intArray[i];
|
587 | 26 | storres | newChain = addElement(newChain, currentInt);
|
588 | 26 | storres | }
|
589 | 26 | storres | return(newChain);
|
590 | 26 | storres | } // End pobyso_create_chain_from_unsigned_int_array. */
|
591 | 26 | storres | |
592 | 26 | storres | node*
|
593 | 26 | storres | pobyso_differentiate(node *functionNode)
|
594 | 26 | storres | {
|
595 | 26 | storres | /* Argument checking. */
|
596 | 26 | storres | node *differentialNode;
|
597 | 26 | storres | if (functionNode == NULL)
|
598 | 26 | storres | {
|
599 | 26 | storres | pobyso_error_message("pobyso_differentiate",
|
600 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
601 | 26 | storres | "The function to differentiate is a NULL pointer");
|
602 | 26 | storres | return(NULL);
|
603 | 26 | storres | }
|
604 | 26 | storres | differentialNode = differentiate(functionNode);
|
605 | 26 | storres | if (differentialNode == NULL)
|
606 | 26 | storres | {
|
607 | 26 | storres | pobyso_error_message("pobyso_differentiate",
|
608 | 26 | storres | "INTERNAL ERROR",
|
609 | 26 | storres | "Sollya could not differentiate the function");
|
610 | 26 | storres | }
|
611 | 26 | storres | return(differentialNode);
|
612 | 26 | storres | } // End pobyso_differentiate
|
613 | 26 | storres | |
614 | 26 | storres | |
615 | 26 | storres | int
|
616 | 26 | storres | pobyso_dirty_infnorm(mpfr_t infNorm,
|
617 | 26 | storres | node *functionNode,
|
618 | 26 | storres | mpfr_t lowerBound,
|
619 | 26 | storres | mpfr_t upperBound,
|
620 | 26 | storres | mp_prec_t precision)
|
621 | 26 | storres | {
|
622 | 26 | storres | int functionCallResult;
|
623 | 26 | storres | /* Arguments checking. */
|
624 | 26 | storres | if (functionNode == NULL)
|
625 | 26 | storres | {
|
626 | 26 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
627 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
628 | 26 | storres | "The function to compute is a NULL pointer");
|
629 | 26 | storres | return(1);
|
630 | 26 | storres | }
|
631 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) > 0)
|
632 | 26 | storres | {
|
633 | 26 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
634 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
635 | 26 | storres | "The lower bond is greater than the upper bound");
|
636 | 26 | storres | return(1);
|
637 | 26 | storres | }
|
638 | 26 | storres | /* Particular cases. */
|
639 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) == 0)
|
640 | 26 | storres | {
|
641 | 26 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
642 | 26 | storres | functionNode,
|
643 | 26 | storres | lowerBound,
|
644 | 26 | storres | precision);
|
645 | 26 | storres | return(functionCallResult);
|
646 | 26 | storres | }
|
647 | 26 | storres | if (isConstant(functionNode))
|
648 | 26 | storres | {
|
649 | 26 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
650 | 26 | storres | functionNode,
|
651 | 26 | storres | lowerBound,
|
652 | 26 | storres | precision);
|
653 | 26 | storres | if (!functionCallResult)
|
654 | 26 | storres | {
|
655 | 26 | storres | mpfr_abs(infNorm, infNorm, MPFR_RNDN);
|
656 | 26 | storres | }
|
657 | 26 | storres | return(functionCallResult);
|
658 | 26 | storres | }
|
659 | 26 | storres | uncertifiedInfnorm(infNorm,
|
660 | 26 | storres | functionNode,
|
661 | 26 | storres | lowerBound,
|
662 | 26 | storres | upperBound,
|
663 | 26 | storres | POBYSO_DEFAULT_POINTS,
|
664 | 26 | storres | precision);
|
665 | 26 | storres | return(0);
|
666 | 26 | storres | } /* End pobyso_dirty_infnorm. */
|
667 | 26 | storres | |
668 | 26 | storres | int
|
669 | 26 | storres | pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
|
670 | 26 | storres | node *nodeToEvaluate,
|
671 | 26 | storres | mpfr_t argument,
|
672 | 26 | storres | mpfr_prec_t precision)
|
673 | 26 | storres | {
|
674 | 26 | storres | /* Check input arguments. */
|
675 | 26 | storres | if (nodeToEvaluate == NULL)
|
676 | 26 | storres | {
|
677 | 26 | storres | pobyso_error_message("pobyso_evaluate_faithful",
|
678 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
679 | 26 | storres | "nodeToEvaluate is a NULL pointer");
|
680 | 26 | storres | return(1);
|
681 | 26 | storres | }
|
682 | 26 | storres | evaluateFaithful(faithfulEvaluation,
|
683 | 26 | storres | nodeToEvaluate,
|
684 | 26 | storres | argument,
|
685 | 26 | storres | precision);
|
686 | 26 | storres | return(0);
|
687 | 26 | storres | } /* End pobyso_evaluate_faithfull. */
|
688 | 26 | storres | |
689 | 26 | storres | chain*
|
690 | 26 | storres | pobyso_find_zeros(node *function,
|
691 | 26 | storres | mpfr_t *lower_bound,
|
692 | 26 | storres | mpfr_t *upper_bound)
|
693 | 26 | storres | {
|
694 | 26 | storres | mp_prec_t currentPrecision;
|
695 | 26 | storres | mpfr_t currentDiameter;
|
696 | 26 | storres | rangetype bounds;
|
697 | 26 | storres | |
698 | 26 | storres | currentPrecision = getToolPrecision();
|
699 | 26 | storres | mpfr_init2(currentDiameter, currentPrecision);
|
700 | 26 | storres | |
701 | 26 | storres | bounds.a = lower_bound;
|
702 | 26 | storres | bounds.b = upper_bound;
|
703 | 26 | storres | |
704 | 26 | storres | if (bounds.a == NULL || bounds.b == NULL)
|
705 | 26 | storres | {
|
706 | 26 | storres | pobyso_error_message("pobyso_find_zeros",
|
707 | 26 | storres | "MEMORY_ALLOCATION_ERROR",
|
708 | 26 | storres | "Could not allocate one of the bounds");
|
709 | 26 | storres | return(NULL);
|
710 | 26 | storres | }
|
711 | 26 | storres | return(findZerosFunction(function,
|
712 | 26 | storres | bounds,
|
713 | 26 | storres | currentPrecision,
|
714 | 26 | storres | currentDiameter));
|
715 | 26 | storres | } /* End pobyso_find_zeros. */
|
716 | 26 | storres | |
717 | 26 | storres | void
|
718 | 26 | storres | pobyso_free_chain_of_nodes(chain *theChainOfNodes)
|
719 | 26 | storres | {
|
720 | 26 | storres | node *currentNode = NULL;
|
721 | 26 | storres | chain *currentChainElement = NULL;
|
722 | 26 | storres | chain *nextChainElement = NULL;
|
723 | 26 | storres | |
724 | 26 | storres | nextChainElement = theChainOfNodes;
|
725 | 26 | storres | |
726 | 26 | storres | while(nextChainElement != NULL)
|
727 | 26 | storres | {
|
728 | 26 | storres | currentChainElement = nextChainElement;
|
729 | 26 | storres | currentNode = (node*)(currentChainElement->value);
|
730 | 26 | storres | nextChainElement = nextChainElement->next;
|
731 | 26 | storres | free_memory(currentNode);
|
732 | 26 | storres | free((void*)(currentChainElement));
|
733 | 26 | storres | }
|
734 | 26 | storres | } /* End pobyso_free_chain_of_nodes. */
|
735 | 26 | storres | |
736 | 26 | storres | void
|
737 | 26 | storres | pobyso_free_range(rangetype range)
|
738 | 26 | storres | {
|
739 | 26 | storres | |
740 | 26 | storres | mpfr_clear(*(range.a));
|
741 | 26 | storres | mpfr_clear(*(range.b));
|
742 | 26 | storres | free(range.a);
|
743 | 26 | storres | free(range.b);
|
744 | 26 | storres | } /* End pobyso_free_range. */
|
745 | 26 | storres | |
746 | 26 | storres | node*
|
747 | 26 | storres | pobyso_fp_minimax_canonical_monomials_base(node *function,
|
748 | 26 | storres | int degree,
|
749 | 26 | storres | chain *formats,
|
750 | 26 | storres | chain *points,
|
751 | 26 | storres | mpfr_t lowerBound,
|
752 | 26 | storres | mpfr_t upperBound,
|
753 | 26 | storres | int fpFixedArg,
|
754 | 26 | storres | int absRel,
|
755 | 26 | storres | node *constPart,
|
756 | 26 | storres | node *minimax)
|
757 | 26 | storres | {
|
758 | 26 | storres | return(NULL);
|
759 | 26 | storres | } /* End pobyso_fp_minimax_canonical_monomials_base. */
|
760 | 26 | storres | |
761 | 26 | storres | node*
|
762 | 26 | storres | pobyso_parse_function(char *functionString,
|
763 | 26 | storres | char *freeVariableNameString)
|
764 | 26 | storres | {
|
765 | 26 | storres | if (functionString == NULL || freeVariableNameString == NULL)
|
766 | 26 | storres | {
|
767 | 26 | storres | pobyso_error_message("pobyso_parse_function",
|
768 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
769 | 26 | storres | "One of the arguments is a NULL pointer");
|
770 | 26 | storres | return(NULL);
|
771 | 26 | storres | }
|
772 | 26 | storres | return(parseString(functionString));
|
773 | 26 | storres | |
774 | 26 | storres | } /* End pobyso_parse_function */
|
775 | 26 | storres | |
776 | 26 | storres | node*
|
777 | 26 | storres | pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
|
778 | 26 | storres | unsigned int mode,
|
779 | 26 | storres | mpfr_t lowerBound,
|
780 | 26 | storres | mpfr_t upperBound,
|
781 | 26 | storres | mpfr_t eps)
|
782 | 26 | storres | {
|
783 | 26 | storres | node *weight = NULL;
|
784 | 26 | storres | node *bestApproxPolyNode = NULL;
|
785 | 26 | storres | node *bestApproxHorner = NULL;
|
786 | 26 | storres | node *errorNode = NULL;
|
787 | 26 | storres | rangetype degreeRange;
|
788 | 26 | storres | mpfr_t quality;
|
789 | 26 | storres | mpfr_t currentError;
|
790 | 26 | storres | unsigned int degree;
|
791 | 26 | storres | |
792 | 26 | storres | /* Check the parameters. */
|
793 | 26 | storres | if (functionNode == NULL)
|
794 | 26 | storres | {
|
795 | 26 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
796 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
797 | 26 | storres | "functionNode is a NULL pointer");
|
798 | 26 | storres | return(NULL);
|
799 | 26 | storres | }
|
800 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
801 | 26 | storres | {
|
802 | 26 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
803 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
804 | 26 | storres | "the lower_bound >= upper_bound");
|
805 | 26 | storres | return(NULL);
|
806 | 26 | storres | }
|
807 | 26 | storres | /* Set the weight. */
|
808 | 26 | storres | if (mode == POBYSO_ABSOLUTE)
|
809 | 26 | storres | {
|
810 | 26 | storres | /* Set the weight to 1 for the ABSOLUTE_MODE. */
|
811 | 26 | storres | weight = makeConstantDouble(1.0);
|
812 | 26 | storres | }
|
813 | 26 | storres | else
|
814 | 26 | storres | {
|
815 | 26 | storres | if (mode == POBYSO_RELATIVE)
|
816 | 26 | storres | {
|
817 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
818 | 26 | storres | "NOT_IMPLEMENTED",
|
819 | 26 | storres | "the search for relative error is not implemented yet");
|
820 | 26 | storres | return(NULL);
|
821 | 26 | storres | }
|
822 | 26 | storres | else
|
823 | 26 | storres | {
|
824 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
825 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
826 | 26 | storres | "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
|
827 | 26 | storres | return(NULL);
|
828 | 26 | storres | }
|
829 | 26 | storres | }
|
830 | 26 | storres | //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
|
831 | 26 | storres | degreeRange = guessDegree(functionNode,
|
832 | 26 | storres | weight,
|
833 | 26 | storres | lowerBound,
|
834 | 26 | storres | upperBound,
|
835 | 26 | storres | eps,
|
836 | 26 | storres | POBYSO_GUESS_DEGREE_BOUND);
|
837 | 26 | storres | degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
|
838 | 26 | storres | //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
|
839 | 26 | storres | fprintf(stderr, "Guessed degree: ");
|
840 | 26 | storres | mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
|
841 | 26 | storres | fprintf(stderr, " - as int: %u\n", degree);
|
842 | 26 | storres | /* Reduce the degree by 1 in the foolish hope it could work. */
|
843 | 26 | storres | if (degree > 0) degree--;
|
844 | 26 | storres | /* Both elements of degreeRange where "inited" within guessDegree. */
|
845 | 26 | storres | mpfr_clear(*(degreeRange.a));
|
846 | 26 | storres | mpfr_clear(*(degreeRange.b));
|
847 | 26 | storres | free(degreeRange.a);
|
848 | 26 | storres | free(degreeRange.b);
|
849 | 26 | storres | /* Mimic the default behavior of interactive Sollya. */
|
850 | 26 | storres | mpfr_init(quality);
|
851 | 26 | storres | mpfr_set_d(quality, 1e-5, MPFR_RNDN);
|
852 | 26 | storres | mpfr_init2(currentError, getToolPrecision());
|
853 | 26 | storres | mpfr_set_inf(currentError,1);
|
854 | 26 | storres | |
855 | 26 | storres | /* Try to refine the initial guess: loop with increasing degrees until
|
856 | 26 | storres | * we find a satisfactory one. */
|
857 | 26 | storres | while(mpfr_cmp(currentError, eps) > 0)
|
858 | 26 | storres | {
|
859 | 26 | storres | /* Get rid of the previous polynomial, if any. */
|
860 | 26 | storres | if (bestApproxPolyNode != NULL)
|
861 | 26 | storres | {
|
862 | 26 | storres | free_memory(bestApproxPolyNode);
|
863 | 26 | storres | }
|
864 | 26 | storres | fprintf(stderr, "Degree: %u\n", degree);
|
865 | 26 | storres | fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
|
866 | 26 | storres | /* Try to find a polynomial with the guessed degree. */
|
867 | 26 | storres | bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
|
868 | 26 | storres | weight,
|
869 | 26 | storres | degree,
|
870 | 26 | storres | lowerBound,
|
871 | 26 | storres | upperBound,
|
872 | 26 | storres | quality);
|
873 | 26 | storres | |
874 | 26 | storres | if (bestApproxPolyNode == NULL)
|
875 | 26 | storres | {
|
876 | 26 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
877 | 26 | storres | "INTERNAL_ERROR",
|
878 | 26 | storres | "could not compute the bestApproxPolyNode");
|
879 | 26 | storres | mpfr_clear(currentError);
|
880 | 26 | storres | mpfr_clear(quality);
|
881 | 26 | storres | return(NULL);
|
882 | 26 | storres | }
|
883 | 26 | storres | |
884 | 26 | storres | setDisplayMode(DISPLAY_MODE_DECIMAL);
|
885 | 26 | storres | fprintTree(stderr, bestApproxPolyNode);
|
886 | 26 | storres | fprintf(stderr, "\n\n");
|
887 | 26 | storres | |
888 | 26 | storres | errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
|
889 | 26 | storres | /* Check the error with the computed polynomial. */
|
890 | 26 | storres | uncertifiedInfnorm(currentError,
|
891 | 26 | storres | errorNode,
|
892 | 26 | storres | lowerBound,
|
893 | 26 | storres | upperBound,
|
894 | 26 | storres | POBYSO_INF_NORM_NUM_POINTS,
|
895 | 26 | storres | getToolPrecision());
|
896 | 26 | storres | fprintf(stderr, "Inf norm: ");
|
897 | 26 | storres | mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
|
898 | 26 | storres | fprintf(stderr, "\n\n");
|
899 | 26 | storres | /* Free the errorNode but not the bestApproxPolyNode (we need it if
|
900 | 26 | storres | * we exit the loop at the next iteration). */
|
901 | 26 | storres | free_memory(errorNode);
|
902 | 26 | storres | degree++;
|
903 | 26 | storres | }
|
904 | 26 | storres | /* Use an intermediate variable, since horner() creates a new node
|
905 | 26 | storres | * and does not reorder the argument "in place". This allows for the memory
|
906 | 26 | storres | * reclaim of bestApproxHorner.
|
907 | 26 | storres | */
|
908 | 26 | storres | bestApproxHorner = horner(bestApproxPolyNode);
|
909 | 26 | storres | free_memory(bestApproxPolyNode);
|
910 | 26 | storres | mpfr_clear(currentError);
|
911 | 26 | storres | mpfr_clear(quality);
|
912 | 26 | storres | free_memory(weight);
|
913 | 26 | storres | return(bestApproxHorner);
|
914 | 26 | storres | } /* End pobyso_remez_approx_canonical_monomials_base_for_error */
|
915 | 26 | storres | |
916 | 26 | storres | node*
|
917 | 26 | storres | pobyso_remez_canonical_monomials_base(node *function,
|
918 | 26 | storres | node *weight,
|
919 | 26 | storres | unsigned int degree,
|
920 | 26 | storres | mpfr_t lowerBound,
|
921 | 26 | storres | mpfr_t upperBound,
|
922 | 26 | storres | mpfr_t quality)
|
923 | 26 | storres | {
|
924 | 26 | storres | node *bestApproxPoly = NULL;
|
925 | 26 | storres | chain *monomials = NULL;
|
926 | 26 | storres | chain *curMonomial = NULL;
|
927 | 26 | storres | |
928 | 26 | storres | mpfr_t satisfying_error;
|
929 | 26 | storres | mpfr_t target_error;
|
930 | 26 | storres | |
931 | 26 | storres | /* Argument checking */
|
932 | 26 | storres | /* Function tree. */
|
933 | 26 | storres | if (function == NULL)
|
934 | 26 | storres | {
|
935 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
936 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
937 | 26 | storres | "the \"function\" argument is a NULL pointer");
|
938 | 26 | storres | return(NULL);
|
939 | 26 | storres | }
|
940 | 26 | storres | if (weight == NULL)
|
941 | 26 | storres | {
|
942 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
943 | 26 | storres | "NULL_POINTER_ARGUMENT",
|
944 | 26 | storres | "the \"weight\" argument is a NULL pointer");
|
945 | 26 | storres | return(NULL);
|
946 | 26 | storres | }
|
947 | 26 | storres | /* Check the bounds. */
|
948 | 26 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
949 | 26 | storres | {
|
950 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
951 | 26 | storres | "INCOHERENT_IMPUT_DATA",
|
952 | 26 | storres | "the lower_bound >= upper_bound");
|
953 | 26 | storres | return(NULL);
|
954 | 26 | storres | }
|
955 | 26 | storres | /* The quality must be a non null positive number. */
|
956 | 26 | storres | if (mpfr_sgn(quality) <= 0)
|
957 | 26 | storres | {
|
958 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
959 | 26 | storres | "INCOHERENT_INPUT_DATA",
|
960 | 26 | storres | "the quality <= 0");
|
961 | 26 | storres | }
|
962 | 26 | storres | /* End argument checking. */
|
963 | 26 | storres | /* Create the monomials nodes chains. */
|
964 | 26 | storres | monomials = pobyso_create_canonical_monomials_base(degree);
|
965 | 26 | storres | fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
|
966 | 26 | storres | if (monomials == NULL || (lengthChain(monomials) != degree + 1))
|
967 | 26 | storres | {
|
968 | 26 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
969 | 26 | storres | "CHAIN_CREATION_ERROR",
|
970 | 26 | storres | "could not create the monomials chain");
|
971 | 26 | storres | return(NULL);
|
972 | 26 | storres | }
|
973 | 26 | storres | curMonomial = monomials;
|
974 | 26 | storres | |
975 | 26 | storres | while (curMonomial != NULL)
|
976 | 26 | storres | {
|
977 | 26 | storres | fprintf(stderr, "monomial tree: ");
|
978 | 26 | storres | //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
|
979 | 26 | storres | fprintTree(stderr, (node*)(curMonomial->value));
|
980 | 26 | storres | fprintf(stderr, "\n");
|
981 | 26 | storres | curMonomial = curMonomial->next;
|
982 | 26 | storres | }
|
983 | 26 | storres | |
984 | 26 | storres | /* Deal with NULL weight. */
|
985 | 26 | storres | if (weight == NULL)
|
986 | 26 | storres | {
|
987 | 26 | storres | weight = makeConstantDouble(1.0);
|
988 | 26 | storres | }
|
989 | 26 | storres | /* Compute the best polynomial with the required quality.
|
990 | 26 | storres | The behavior is as if satisfying error and target_error had
|
991 | 26 | storres | not been used.*/
|
992 | 26 | storres | mpfr_init(satisfying_error);
|
993 | 26 | storres | mpfr_init(target_error);
|
994 | 26 | storres | mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
|
995 | 26 | storres | mpfr_set_inf(target_error, 1);
|
996 | 26 | storres | |
997 | 26 | storres | |
998 | 26 | storres | fprintf(stderr, "satisfying_error: ");
|
999 | 26 | storres | mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
|
1000 | 26 | storres | fprintf(stderr, ".\n");
|
1001 | 26 | storres | fprintf(stderr, "target_error: ");
|
1002 | 26 | storres | mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
|
1003 | 26 | storres | fprintf(stderr, ".\n");
|
1004 | 26 | storres | |
1005 | 26 | storres | fprintf(stderr,
|
1006 | 26 | storres | "current precision: %li\n", getToolPrecision());
|
1007 | 26 | storres | /* Call the Sollya function. */
|
1008 | 26 | storres | bestApproxPoly = remez(function,
|
1009 | 26 | storres | weight,
|
1010 | 26 | storres | monomials,
|
1011 | 26 | storres | lowerBound,
|
1012 | 26 | storres | upperBound,
|
1013 | 26 | storres | quality,
|
1014 | 26 | storres | satisfying_error,
|
1015 | 26 | storres | target_error,
|
1016 | 26 | storres | getToolPrecision());
|
1017 | 26 | storres | |
1018 | 26 | storres | mpfr_clear(satisfying_error);
|
1019 | 26 | storres | mpfr_clear(target_error);
|
1020 | 26 | storres | pobyso_free_chain_of_nodes(monomials);
|
1021 | 26 | storres | |
1022 | 26 | storres | return(bestApproxPoly);
|
1023 | 26 | storres | } /* End pobyso_remez_canonical_monomials_base. */
|
1024 | 26 | storres | |
1025 | 26 | storres | #endif
|
1026 | 26 | storres | |
1027 | 26 | storres | void
|
1028 | 26 | storres | pobyso_error_message(char *functionName, char *messageName, char* message) |
1029 | 26 | storres | { |
1030 | 26 | storres | fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
|
1031 | 26 | storres | } /* End pobyso_error_message */ |