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