root / pobysoC-4.0 / src / pobyso.c @ 133
Historique | Voir | Annoter | Télécharger (25,12 ko)
1 |
/** @file pobyso.c
|
---|---|
2 |
* Name & purpose
|
3 |
*
|
4 |
* @author S.T.
|
5 |
* @date 2011-10-12
|
6 |
*
|
7 |
*
|
8 |
*/
|
9 |
|
10 |
/******************************************************************************/
|
11 |
/* Headers, applying the "particular to general" convention.*/
|
12 |
|
13 |
#include "pobyso.h" |
14 |
|
15 |
/* includes of local headers */
|
16 |
|
17 |
/* includes of project headers */
|
18 |
|
19 |
/* includes of system headers */
|
20 |
#include <string.h> |
21 |
#include <stdlib.h> |
22 |
#include <stdio.h> |
23 |
|
24 |
/* Other declarations */
|
25 |
|
26 |
/* Internal prototypes */
|
27 |
void
|
28 |
pobyso_error_message(char *functionName, char *messageName, char* message); |
29 |
/* Types, constants and macros definitions */
|
30 |
|
31 |
/* Global variables */
|
32 |
|
33 |
/* Functions */
|
34 |
|
35 |
/* @see pobyso.h#pobyso_autprint */
|
36 |
void
|
37 |
pobyso_autoprint(sollya_obj_t solObj, ...) |
38 |
{ |
39 |
va_list va; |
40 |
va_start(va, solObj); |
41 |
sollya_lib_v_autoprint(solObj, va); |
42 |
va_end(va); |
43 |
} /* End pobyso_autoprint. */
|
44 |
|
45 |
/* @see pobyso.h#pobyso_is_constant_expression*/
|
46 |
int
|
47 |
pobyso_is_constant_expression(sollya_obj_t obj_to_test) |
48 |
{ |
49 |
mpfr_t dummy; |
50 |
int test;
|
51 |
sollya_obj_t verbosity_level = NULL;
|
52 |
/* Test argument. */
|
53 |
if (obj_to_test == NULL) |
54 |
{ |
55 |
pobyso_error_message("pobyso_parse_string",
|
56 |
"NULL_POINTER_ARGUMENT",
|
57 |
"The expression is a NULL pointer");
|
58 |
return 0; |
59 |
} |
60 |
verbosity_level = pobyso_set_verbosity_off(); |
61 |
|
62 |
if (! sollya_lib_obj_is_function(obj_to_test))
|
63 |
{ |
64 |
pobyso_set_verbosity_to(verbosity_level); |
65 |
sollya_lib_clear_obj(verbosity_level); |
66 |
return 0; |
67 |
} |
68 |
mpfr_init2(dummy,64);
|
69 |
/* Call to previous Sollya function resets verbosity. */
|
70 |
verbosity_level = pobyso_set_verbosity_off(); |
71 |
test = sollya_lib_get_constant(dummy, obj_to_test); |
72 |
mpfr_clear(dummy); |
73 |
pobyso_set_verbosity_to(verbosity_level); |
74 |
sollya_lib_clear_obj(verbosity_level); |
75 |
return test;
|
76 |
} /* End pobyso_is_constant_expression. */
|
77 |
|
78 |
/** @see pobyso.h#pobyso_new_monomial. */
|
79 |
pobyso_func_exp_t |
80 |
pobyso_new_monomial(pobyso_func_exp_t coefficientSa, long degree)
|
81 |
{ |
82 |
sollya_obj_t degreeSa = NULL;
|
83 |
sollya_obj_t varToPowSa = NULL;
|
84 |
sollya_obj_t monomialSa = NULL;
|
85 |
if (coefficientSa == NULL) |
86 |
{ |
87 |
pobyso_error_message("pobyso_parse_string",
|
88 |
"NULL_POINTER_ARGUMENT",
|
89 |
"The expression is a NULL pointer");
|
90 |
return sollya_lib_error();
|
91 |
} |
92 |
if (! pobyso_is_constant_expression(coefficientSa))
|
93 |
{ |
94 |
return sollya_lib_error();
|
95 |
} |
96 |
if (degree < 0) |
97 |
{ |
98 |
return sollya_lib_error();
|
99 |
} |
100 |
degreeSa = sollya_lib_constant_from_int64(degree); |
101 |
varToPowSa = sollya_lib_build_function_pow(sollya_lib_free_variable(), |
102 |
degreeSa); |
103 |
monomialSa = sollya_lib_mul(coefficientSa,varToPowSa); |
104 |
sollya_lib_clear_obj(varToPowSa); |
105 |
return monomialSa;
|
106 |
} /* End pobyso_new_monomial. */
|
107 |
|
108 |
/* @see pobyso.h#pobyso_parse_string*/
|
109 |
pobyso_func_exp_t |
110 |
pobyso_parse_string(const char* expression) |
111 |
{ |
112 |
int expressionLength, i;
|
113 |
char *expressionWithSemiCo;
|
114 |
sollya_obj_t expressionSa; |
115 |
|
116 |
/* Arguments check. */
|
117 |
if (expression == NULL) |
118 |
{ |
119 |
pobyso_error_message("pobyso_parse_string",
|
120 |
"NULL_POINTER_ARGUMENT",
|
121 |
"The expression is a NULL pointer");
|
122 |
return(sollya_lib_error());
|
123 |
} |
124 |
expressionLength = strlen(expression); |
125 |
if (expressionLength == 0) |
126 |
{ |
127 |
pobyso_error_message("pobyso_parse_string",
|
128 |
"EMPTY_STRING_ARGUMENT",
|
129 |
"The expression is an empty string");
|
130 |
return sollya_lib_error();
|
131 |
} |
132 |
/* Search from the last char of the expression until, whichever happens
|
133 |
* first:
|
134 |
* a ";" is found;
|
135 |
* a char != ';' is found the the ";" is appended.
|
136 |
* If the head of the string is reached before any of these two events happens
|
137 |
* return an error.
|
138 |
*/
|
139 |
for (i = expressionLength - 1 ; i >= 0 ; i--) |
140 |
{ |
141 |
if (expression[i] == ';') /* Nothing special to do: |
142 |
try to parse the string*/
|
143 |
{ |
144 |
expressionSa = sollya_lib_parse_string(expression); |
145 |
return expressionSa;
|
146 |
} |
147 |
else
|
148 |
{ |
149 |
if (expression[i] == ' ' || expression[i] == '\t') /* A blank, |
150 |
just continue. */
|
151 |
{ |
152 |
continue;
|
153 |
} |
154 |
else /* a character != ';' and from a blank: create the ';'ed string. */ |
155 |
{ |
156 |
/* Create a new string for the expression, add the ";" and
|
157 |
* and call sollya_lib_parse_string. */
|
158 |
expressionWithSemiCo = calloc(i + 3, sizeof(char)); |
159 |
if (expressionWithSemiCo == NULL) |
160 |
{ |
161 |
pobyso_error_message("pobyso_parse_string",
|
162 |
"MEMORY_ALLOCATION_ERROR",
|
163 |
"Could not allocate the expression string");
|
164 |
return sollya_lib_error();
|
165 |
} |
166 |
strncpy(expressionWithSemiCo, expression, i+1);
|
167 |
expressionWithSemiCo[i + 1] = ';'; |
168 |
expressionWithSemiCo[i + 2] = '\0'; |
169 |
expressionSa = sollya_lib_parse_string(expressionWithSemiCo); |
170 |
free(expressionWithSemiCo); |
171 |
return expressionSa;
|
172 |
} /* End character != ';' and from a blank. */
|
173 |
/* Create a new string for the expression, add the ";" and
|
174 |
* and call sollya_lib_parse_string. */
|
175 |
} /* End else. */
|
176 |
} /* End for. */
|
177 |
/* We get here, it is because we did not find any char == anything different
|
178 |
* from ' ' or '\t': the string is empty.
|
179 |
*/
|
180 |
pobyso_error_message("pobyso_parse_string",
|
181 |
"ONLY_BLANK_ARGUMENT",
|
182 |
"The expression is only made of blanks");
|
183 |
return(sollya_lib_error());
|
184 |
|
185 |
} /* pobyso_parse_string */
|
186 |
|
187 |
pobyso_func_exp_t |
188 |
pobyso_remez_canonical_monomials_base(pobyso_func_exp_t function, |
189 |
long int degree, |
190 |
pobyso_range_t interval, |
191 |
pobyso_func_exp_t weight, |
192 |
double quality,
|
193 |
pobyso_range_t bounds) |
194 |
{ |
195 |
sollya_obj_t degreeSa = NULL;
|
196 |
sollya_obj_t qualitySa = NULL;
|
197 |
|
198 |
degreeSa = sollya_lib_constant_from_int(degree); |
199 |
qualitySa = sollya_lib_constant_from_double(quality); |
200 |
|
201 |
return NULL; |
202 |
} /* End pobyso_remez_canonical_monomials_base. */
|
203 |
|
204 |
int
|
205 |
pobyso_set_canonical_on() |
206 |
{ |
207 |
pobyso_on_off_t currentCanonicalModeSo; |
208 |
int currentCanonicalMode;
|
209 |
pobyso_on_off_t on; |
210 |
|
211 |
on = sollya_lib_on(); |
212 |
currentCanonicalModeSo = sollya_lib_get_canonical(); |
213 |
if (sollya_lib_is_on(currentCanonicalModeSo))
|
214 |
{ |
215 |
currentCanonicalMode = POBYSO_ON; |
216 |
} |
217 |
else
|
218 |
{ |
219 |
currentCanonicalMode = POBYSO_OFF; |
220 |
} |
221 |
sollya_lib_set_canonical(on); |
222 |
sollya_lib_clear_obj(on); |
223 |
sollya_lib_clear_obj(currentCanonicalModeSo); |
224 |
return currentCanonicalMode;
|
225 |
} /* End pobyso_set_canonical_on. */
|
226 |
|
227 |
pobyso_sollya_verbosity_t |
228 |
pobyso_set_verbosity_off() |
229 |
{ |
230 |
sollya_obj_t verbosityLevelZero; |
231 |
sollya_obj_t currentVerbosityLevel = NULL;
|
232 |
|
233 |
currentVerbosityLevel = sollya_lib_get_verbosity(); |
234 |
verbosityLevelZero = sollya_lib_constant_from_int(0);
|
235 |
sollya_lib_set_verbosity(verbosityLevelZero); |
236 |
sollya_lib_clear_obj(verbosityLevelZero); |
237 |
return(currentVerbosityLevel);
|
238 |
} /* End of pobyso_set_verbosity_off. */
|
239 |
|
240 |
pobyso_sollya_verbosity_t |
241 |
pobyso_set_verbosity_to(pobyso_sollya_verbosity_t newVerbosityLevel) |
242 |
{ |
243 |
sollya_obj_t currentVerbosityLevel = NULL;
|
244 |
if (newVerbosityLevel == NULL) |
245 |
{ |
246 |
pobyso_error_message("pobyso_set_verbosity_to",
|
247 |
"NULL_POINTER_ARGUMENT",
|
248 |
"The new verbosity level is a null pointer");
|
249 |
return(sollya_lib_error());
|
250 |
} |
251 |
currentVerbosityLevel = sollya_lib_get_verbosity(); |
252 |
sollya_lib_set_verbosity(newVerbosityLevel); |
253 |
return(currentVerbosityLevel);
|
254 |
} /* End of pobyso_set_verbosity_to. */
|
255 |
|
256 |
pobyso_sollya_verbosity_t |
257 |
pobyso_set_verbosity_to_int(int intVerbosityLevel)
|
258 |
{ |
259 |
sollya_obj_t currentVerbosityLevel = NULL;
|
260 |
sollya_obj_t newVerbosityLevel = NULL;
|
261 |
if (intVerbosityLevel < 0) |
262 |
{ |
263 |
pobyso_error_message("pobyso_set_verbosity_to",
|
264 |
"INVALID_VALUE",
|
265 |
"The new verbosity level is a negative number");
|
266 |
return(sollya_lib_error());
|
267 |
} |
268 |
newVerbosityLevel = sollya_lib_constant_from_int(intVerbosityLevel); |
269 |
currentVerbosityLevel = sollya_lib_get_verbosity(); |
270 |
sollya_lib_set_verbosity(newVerbosityLevel); |
271 |
return(currentVerbosityLevel);
|
272 |
} /* End of pobyso_set_verbosity_to_int. */
|
273 |
|
274 |
/* Attic from the sollya_lib < 4. */
|
275 |
#if 0
|
276 |
chain*
|
277 |
pobyso_create_canonical_monomials_base(const unsigned int degree)
|
278 |
{
|
279 |
int i = 0;
|
280 |
chain *monomials = NULL;
|
281 |
node *monomial = NULL;
|
282 |
|
283 |
for(i = degree ; i >= 0 ; i--)
|
284 |
{
|
285 |
monomial = makePow(makeVariable(), makeConstantDouble((double)i));
|
286 |
monomials = addElement(monomials, monomial);
|
287 |
fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
|
288 |
}
|
289 |
if (monomials == NULL)
|
290 |
{
|
291 |
pobyso_error_message("pobyso_create_canonical_monomial_base",
|
292 |
"CHAIN_CREATION_ERROR",
|
293 |
"Could not create the monomials chain");
|
294 |
return(NULL);
|
295 |
}
|
296 |
return(monomials);
|
297 |
} /* End pobyso_create_canonical_monomials_base. */
|
298 |
|
299 |
chain*
|
300 |
pobyso_create_chain_from_int_array(int* intArray,
|
301 |
const unsigned int arrayLength)
|
302 |
{
|
303 |
int i = 0;
|
304 |
chain *newChain = NULL;
|
305 |
int *currentInt;
|
306 |
|
307 |
if (arrayLength == 0) return(NULL);
|
308 |
if (intArray == NULL)
|
309 |
{
|
310 |
pobyso_error_message("pobyso_create_chain_from_int_array",
|
311 |
"NULL_POINTER_ARGUMENT",
|
312 |
"The array is a NULL pointer");
|
313 |
return(NULL);
|
314 |
}
|
315 |
for (i = arrayLength - 1 ; i >= 0 ; i--)
|
316 |
{
|
317 |
currentInt = malloc(sizeof(int));
|
318 |
if (currentInt == NULL)
|
319 |
{
|
320 |
pobyso_error_message("pobyso_create_chain_from_int_array",
|
321 |
"MEMORY_ALLOCATION_ERROR",
|
322 |
"Could not allocate one of the integers");
|
323 |
freeChain(newChain, free);
|
324 |
return(NULL);
|
325 |
}
|
326 |
*currentInt = intArray[i];
|
327 |
newChain = addElement(newChain, currentInt);
|
328 |
}
|
329 |
return(newChain);
|
330 |
} // End pobyso_create_chain_from_int_array. */
|
331 |
|
332 |
chain*
|
333 |
pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
|
334 |
const unsigned int arrayLength)
|
335 |
{
|
336 |
int i = 0;
|
337 |
chain *newChain = NULL;
|
338 |
unsigned int *currentInt;
|
339 |
|
340 |
/* Argument checking. */
|
341 |
if (arrayLength == 0) return(NULL);
|
342 |
if (intArray == NULL)
|
343 |
{
|
344 |
pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
345 |
"NULL_POINTER_ARGUMENT",
|
346 |
"The array is a NULL pointer");
|
347 |
return(NULL);
|
348 |
}
|
349 |
for (i = arrayLength - 1 ; i >= 0 ; i--)
|
350 |
{
|
351 |
currentInt = malloc(sizeof(unsigned int));
|
352 |
if (currentInt == NULL)
|
353 |
{
|
354 |
pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
355 |
"MEMORY_ALLOCATION_ERROR",
|
356 |
"Could not allocate one of the integers");
|
357 |
freeChain(newChain, free);
|
358 |
return(NULL);
|
359 |
}
|
360 |
*currentInt = intArray[i];
|
361 |
newChain = addElement(newChain, currentInt);
|
362 |
}
|
363 |
return(newChain);
|
364 |
} // End pobyso_create_chain_from_unsigned_int_array. */
|
365 |
|
366 |
node*
|
367 |
pobyso_differentiate(node *functionNode)
|
368 |
{
|
369 |
/* Argument checking. */
|
370 |
node *differentialNode;
|
371 |
if (functionNode == NULL)
|
372 |
{
|
373 |
pobyso_error_message("pobyso_differentiate",
|
374 |
"NULL_POINTER_ARGUMENT",
|
375 |
"The function to differentiate is a NULL pointer");
|
376 |
return(NULL);
|
377 |
}
|
378 |
differentialNode = differentiate(functionNode);
|
379 |
if (differentialNode == NULL)
|
380 |
{
|
381 |
pobyso_error_message("pobyso_differentiate",
|
382 |
"INTERNAL ERROR",
|
383 |
"Sollya could not differentiate the function");
|
384 |
}
|
385 |
return(differentialNode);
|
386 |
} // End pobyso_differentiate
|
387 |
|
388 |
|
389 |
int
|
390 |
pobyso_dirty_infnorm(mpfr_t infNorm,
|
391 |
node *functionNode,
|
392 |
mpfr_t lowerBound,
|
393 |
mpfr_t upperBound,
|
394 |
mp_prec_t precision)
|
395 |
{
|
396 |
int functionCallResult;
|
397 |
/* Arguments checking. */
|
398 |
if (functionNode == NULL)
|
399 |
{
|
400 |
pobyso_error_message("pobyso_dirty_infnorm",
|
401 |
"NULL_POINTER_ARGUMENT",
|
402 |
"The function to compute is a NULL pointer");
|
403 |
return(1);
|
404 |
}
|
405 |
if (mpfr_cmp(lowerBound, upperBound) > 0)
|
406 |
{
|
407 |
pobyso_error_message("pobyso_dirty_infnorm",
|
408 |
"INCOHERENT_INPUT_DATA",
|
409 |
"The lower bond is greater than the upper bound");
|
410 |
return(1);
|
411 |
}
|
412 |
/* Particular cases. */
|
413 |
if (mpfr_cmp(lowerBound, upperBound) == 0)
|
414 |
{
|
415 |
functionCallResult = pobyso_evaluate_faithful(infNorm,
|
416 |
functionNode,
|
417 |
lowerBound,
|
418 |
precision);
|
419 |
return(functionCallResult);
|
420 |
}
|
421 |
if (isConstant(functionNode))
|
422 |
{
|
423 |
functionCallResult = pobyso_evaluate_faithful(infNorm,
|
424 |
functionNode,
|
425 |
lowerBound,
|
426 |
precision);
|
427 |
if (!functionCallResult)
|
428 |
{
|
429 |
mpfr_abs(infNorm, infNorm, MPFR_RNDN);
|
430 |
}
|
431 |
return(functionCallResult);
|
432 |
}
|
433 |
uncertifiedInfnorm(infNorm,
|
434 |
functionNode,
|
435 |
lowerBound,
|
436 |
upperBound,
|
437 |
POBYSO_DEFAULT_POINTS,
|
438 |
precision);
|
439 |
return(0);
|
440 |
} /* End pobyso_dirty_infnorm. */
|
441 |
|
442 |
int
|
443 |
pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
|
444 |
node *nodeToEvaluate,
|
445 |
mpfr_t argument,
|
446 |
mpfr_prec_t precision)
|
447 |
{
|
448 |
/* Check input arguments. */
|
449 |
if (nodeToEvaluate == NULL)
|
450 |
{
|
451 |
pobyso_error_message("pobyso_evaluate_faithful",
|
452 |
"NULL_POINTER_ARGUMENT",
|
453 |
"nodeToEvaluate is a NULL pointer");
|
454 |
return(1);
|
455 |
}
|
456 |
evaluateFaithful(faithfulEvaluation,
|
457 |
nodeToEvaluate,
|
458 |
argument,
|
459 |
precision);
|
460 |
return(0);
|
461 |
} /* End pobyso_evaluate_faithfull. */
|
462 |
|
463 |
chain*
|
464 |
pobyso_find_zeros(node *function,
|
465 |
mpfr_t *lower_bound,
|
466 |
mpfr_t *upper_bound)
|
467 |
{
|
468 |
mp_prec_t currentPrecision;
|
469 |
mpfr_t currentDiameter;
|
470 |
rangetype bounds;
|
471 |
|
472 |
currentPrecision = getToolPrecision();
|
473 |
mpfr_init2(currentDiameter, currentPrecision);
|
474 |
|
475 |
bounds.a = lower_bound;
|
476 |
bounds.b = upper_bound;
|
477 |
|
478 |
if (bounds.a == NULL || bounds.b == NULL)
|
479 |
{
|
480 |
pobyso_error_message("pobyso_find_zeros",
|
481 |
"MEMORY_ALLOCATION_ERROR",
|
482 |
"Could not allocate one of the bounds");
|
483 |
return(NULL);
|
484 |
}
|
485 |
return(findZerosFunction(function,
|
486 |
bounds,
|
487 |
currentPrecision,
|
488 |
currentDiameter));
|
489 |
} /* End pobyso_find_zeros. */
|
490 |
|
491 |
void
|
492 |
pobyso_free_chain_of_nodes(chain *theChainOfNodes)
|
493 |
{
|
494 |
node *currentNode = NULL;
|
495 |
chain *currentChainElement = NULL;
|
496 |
chain *nextChainElement = NULL;
|
497 |
|
498 |
nextChainElement = theChainOfNodes;
|
499 |
|
500 |
while(nextChainElement != NULL)
|
501 |
{
|
502 |
currentChainElement = nextChainElement;
|
503 |
currentNode = (node*)(currentChainElement->value);
|
504 |
nextChainElement = nextChainElement->next;
|
505 |
free_memory(currentNode);
|
506 |
free((void*)(currentChainElement));
|
507 |
}
|
508 |
} /* End pobyso_free_chain_of_nodes. */
|
509 |
|
510 |
void
|
511 |
pobyso_free_range(rangetype range)
|
512 |
{
|
513 |
|
514 |
mpfr_clear(*(range.a));
|
515 |
mpfr_clear(*(range.b));
|
516 |
free(range.a);
|
517 |
free(range.b);
|
518 |
} /* End pobyso_free_range. */
|
519 |
|
520 |
node*
|
521 |
pobyso_fp_minimax_canonical_monomials_base(node *function,
|
522 |
int degree,
|
523 |
chain *formats,
|
524 |
chain *points,
|
525 |
mpfr_t lowerBound,
|
526 |
mpfr_t upperBound,
|
527 |
int fpFixedArg,
|
528 |
int absRel,
|
529 |
node *constPart,
|
530 |
node *minimax)
|
531 |
{
|
532 |
return(NULL);
|
533 |
} /* End pobyso_fp_minimax_canonical_monomials_base. */
|
534 |
|
535 |
node*
|
536 |
pobyso_parse_function(char *functionString,
|
537 |
char *freeVariableNameString)
|
538 |
{
|
539 |
if (functionString == NULL || freeVariableNameString == NULL)
|
540 |
{
|
541 |
pobyso_error_message("pobyso_parse_function",
|
542 |
"NULL_POINTER_ARGUMENT",
|
543 |
"One of the arguments is a NULL pointer");
|
544 |
return(NULL);
|
545 |
}
|
546 |
return(parseString(functionString));
|
547 |
|
548 |
} /* End pobyso_parse_function */
|
549 |
|
550 |
node*
|
551 |
pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
|
552 |
unsigned int mode,
|
553 |
mpfr_t lowerBound,
|
554 |
mpfr_t upperBound,
|
555 |
mpfr_t eps)
|
556 |
{
|
557 |
node *weight = NULL;
|
558 |
node *bestApproxPolyNode = NULL;
|
559 |
node *bestApproxHorner = NULL;
|
560 |
node *errorNode = NULL;
|
561 |
rangetype degreeRange;
|
562 |
mpfr_t quality;
|
563 |
mpfr_t currentError;
|
564 |
unsigned int degree;
|
565 |
|
566 |
/* Check the parameters. */
|
567 |
if (functionNode == NULL)
|
568 |
{
|
569 |
pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
570 |
"NULL_POINTER_ARGUMENT",
|
571 |
"functionNode is a NULL pointer");
|
572 |
return(NULL);
|
573 |
}
|
574 |
if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
575 |
{
|
576 |
pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
577 |
"INCOHERENT_INPUT_DATA",
|
578 |
"the lower_bound >= upper_bound");
|
579 |
return(NULL);
|
580 |
}
|
581 |
/* Set the weight. */
|
582 |
if (mode == POBYSO_ABSOLUTE)
|
583 |
{
|
584 |
/* Set the weight to 1 for the ABSOLUTE_MODE. */
|
585 |
weight = makeConstantDouble(1.0);
|
586 |
}
|
587 |
else
|
588 |
{
|
589 |
if (mode == POBYSO_RELATIVE)
|
590 |
{
|
591 |
pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
592 |
"NOT_IMPLEMENTED",
|
593 |
"the search for relative error is not implemented yet");
|
594 |
return(NULL);
|
595 |
}
|
596 |
else
|
597 |
{
|
598 |
pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
599 |
"INCOHERENT_INPUT_DATA",
|
600 |
"the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
|
601 |
return(NULL);
|
602 |
}
|
603 |
}
|
604 |
//fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
|
605 |
degreeRange = guessDegree(functionNode,
|
606 |
weight,
|
607 |
lowerBound,
|
608 |
upperBound,
|
609 |
eps,
|
610 |
POBYSO_GUESS_DEGREE_BOUND);
|
611 |
degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
|
612 |
//fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
|
613 |
fprintf(stderr, "Guessed degree: ");
|
614 |
mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
|
615 |
fprintf(stderr, " - as int: %u\n", degree);
|
616 |
/* Reduce the degree by 1 in the foolish hope it could work. */
|
617 |
if (degree > 0) degree--;
|
618 |
/* Both elements of degreeRange where "inited" within guessDegree. */
|
619 |
mpfr_clear(*(degreeRange.a));
|
620 |
mpfr_clear(*(degreeRange.b));
|
621 |
free(degreeRange.a);
|
622 |
free(degreeRange.b);
|
623 |
/* Mimic the default behavior of interactive Sollya. */
|
624 |
mpfr_init(quality);
|
625 |
mpfr_set_d(quality, 1e-5, MPFR_RNDN);
|
626 |
mpfr_init2(currentError, getToolPrecision());
|
627 |
mpfr_set_inf(currentError,1);
|
628 |
|
629 |
/* Try to refine the initial guess: loop with increasing degrees until
|
630 |
* we find a satisfactory one. */
|
631 |
while(mpfr_cmp(currentError, eps) > 0)
|
632 |
{
|
633 |
/* Get rid of the previous polynomial, if any. */
|
634 |
if (bestApproxPolyNode != NULL)
|
635 |
{
|
636 |
free_memory(bestApproxPolyNode);
|
637 |
}
|
638 |
fprintf(stderr, "Degree: %u\n", degree);
|
639 |
fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
|
640 |
/* Try to find a polynomial with the guessed degree. */
|
641 |
bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
|
642 |
weight,
|
643 |
degree,
|
644 |
lowerBound,
|
645 |
upperBound,
|
646 |
quality);
|
647 |
|
648 |
if (bestApproxPolyNode == NULL)
|
649 |
{
|
650 |
pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
651 |
"INTERNAL_ERROR",
|
652 |
"could not compute the bestApproxPolyNode");
|
653 |
mpfr_clear(currentError);
|
654 |
mpfr_clear(quality);
|
655 |
return(NULL);
|
656 |
}
|
657 |
|
658 |
setDisplayMode(DISPLAY_MODE_DECIMAL);
|
659 |
fprintTree(stderr, bestApproxPolyNode);
|
660 |
fprintf(stderr, "\n\n");
|
661 |
|
662 |
errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
|
663 |
/* Check the error with the computed polynomial. */
|
664 |
uncertifiedInfnorm(currentError,
|
665 |
errorNode,
|
666 |
lowerBound,
|
667 |
upperBound,
|
668 |
POBYSO_INF_NORM_NUM_POINTS,
|
669 |
getToolPrecision());
|
670 |
fprintf(stderr, "Inf norm: ");
|
671 |
mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
|
672 |
fprintf(stderr, "\n\n");
|
673 |
/* Free the errorNode but not the bestApproxPolyNode (we need it if
|
674 |
* we exit the loop at the next iteration). */
|
675 |
free_memory(errorNode);
|
676 |
degree++;
|
677 |
}
|
678 |
/* Use an intermediate variable, since horner() creates a new node
|
679 |
* and does not reorder the argument "in place". This allows for the memory
|
680 |
* reclaim of bestApproxHorner.
|
681 |
*/
|
682 |
bestApproxHorner = horner(bestApproxPolyNode);
|
683 |
free_memory(bestApproxPolyNode);
|
684 |
mpfr_clear(currentError);
|
685 |
mpfr_clear(quality);
|
686 |
free_memory(weight);
|
687 |
return(bestApproxHorner);
|
688 |
} /* End pobyso_remez_approx_canonical_monomials_base_for_error */
|
689 |
|
690 |
node*
|
691 |
pobyso_remez_canonical_monomials_base(node *function,
|
692 |
node *weight,
|
693 |
unsigned int degree,
|
694 |
mpfr_t lowerBound,
|
695 |
mpfr_t upperBound,
|
696 |
mpfr_t quality)
|
697 |
{
|
698 |
node *bestApproxPoly = NULL;
|
699 |
chain *monomials = NULL;
|
700 |
chain *curMonomial = NULL;
|
701 |
|
702 |
mpfr_t satisfying_error;
|
703 |
mpfr_t target_error;
|
704 |
|
705 |
/* Argument checking */
|
706 |
/* Function tree. */
|
707 |
if (function == NULL)
|
708 |
{
|
709 |
pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
710 |
"NULL_POINTER_ARGUMENT",
|
711 |
"the \"function\" argument is a NULL pointer");
|
712 |
return(NULL);
|
713 |
}
|
714 |
if (weight == NULL)
|
715 |
{
|
716 |
pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
717 |
"NULL_POINTER_ARGUMENT",
|
718 |
"the \"weight\" argument is a NULL pointer");
|
719 |
return(NULL);
|
720 |
}
|
721 |
/* Check the bounds. */
|
722 |
if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
723 |
{
|
724 |
pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
725 |
"INCOHERENT_IMPUT_DATA",
|
726 |
"the lower_bound >= upper_bound");
|
727 |
return(NULL);
|
728 |
}
|
729 |
/* The quality must be a non null positive number. */
|
730 |
if (mpfr_sgn(quality) <= 0)
|
731 |
{
|
732 |
pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
733 |
"INCOHERENT_INPUT_DATA",
|
734 |
"the quality <= 0");
|
735 |
}
|
736 |
/* End argument checking. */
|
737 |
/* Create the monomials nodes chains. */
|
738 |
monomials = pobyso_create_canonical_monomials_base(degree);
|
739 |
fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
|
740 |
if (monomials == NULL || (lengthChain(monomials) != degree + 1))
|
741 |
{
|
742 |
pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
743 |
"CHAIN_CREATION_ERROR",
|
744 |
"could not create the monomials chain");
|
745 |
return(NULL);
|
746 |
}
|
747 |
curMonomial = monomials;
|
748 |
|
749 |
while (curMonomial != NULL)
|
750 |
{
|
751 |
fprintf(stderr, "monomial tree: ");
|
752 |
//mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
|
753 |
fprintTree(stderr, (node*)(curMonomial->value));
|
754 |
fprintf(stderr, "\n");
|
755 |
curMonomial = curMonomial->next;
|
756 |
}
|
757 |
|
758 |
/* Deal with NULL weight. */
|
759 |
if (weight == NULL)
|
760 |
{
|
761 |
weight = makeConstantDouble(1.0);
|
762 |
}
|
763 |
/* Compute the best polynomial with the required quality.
|
764 |
The behavior is as if satisfying error and target_error had
|
765 |
not been used.*/
|
766 |
mpfr_init(satisfying_error);
|
767 |
mpfr_init(target_error);
|
768 |
mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
|
769 |
mpfr_set_inf(target_error, 1);
|
770 |
|
771 |
|
772 |
fprintf(stderr, "satisfying_error: ");
|
773 |
mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
|
774 |
fprintf(stderr, ".\n");
|
775 |
fprintf(stderr, "target_error: ");
|
776 |
mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
|
777 |
fprintf(stderr, ".\n");
|
778 |
|
779 |
fprintf(stderr,
|
780 |
"current precision: %li\n", getToolPrecision());
|
781 |
/* Call the Sollya function. */
|
782 |
bestApproxPoly = remez(function,
|
783 |
weight,
|
784 |
monomials,
|
785 |
lowerBound,
|
786 |
upperBound,
|
787 |
quality,
|
788 |
satisfying_error,
|
789 |
target_error,
|
790 |
getToolPrecision());
|
791 |
|
792 |
mpfr_clear(satisfying_error);
|
793 |
mpfr_clear(target_error);
|
794 |
pobyso_free_chain_of_nodes(monomials);
|
795 |
|
796 |
return(bestApproxPoly);
|
797 |
} /* End pobyso_remez_canonical_monomials_base. */
|
798 |
|
799 |
#endif
|
800 |
|
801 |
void
|
802 |
pobyso_error_message(char *functionName, char *messageName, char* message) |
803 |
{ |
804 |
fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
|
805 |
} /* End pobyso_error_message */
|