root / pobysoC / pobyso.c @ 3
Historique | Voir | Annoter | Télécharger (18,78 ko)
| 1 | 3 | storres | /** @file pobyso.c
|
|---|---|---|---|
| 2 | 3 | storres | * Name & purpose
|
| 3 | 3 | storres | *
|
| 4 | 3 | storres | * @author S.T.
|
| 5 | 3 | storres | * @date 2011-10-12
|
| 6 | 3 | storres | *
|
| 7 | 3 | storres | */
|
| 8 | 3 | storres | /******************************************************************************/
|
| 9 | 3 | storres | /* Headers, applying the "particular to general" convention.*/
|
| 10 | 3 | storres | |
| 11 | 3 | storres | #include "pobyso.h" |
| 12 | 3 | storres | |
| 13 | 3 | storres | /* includes of local headers */
|
| 14 | 3 | storres | |
| 15 | 3 | storres | /* includes of project headers */
|
| 16 | 3 | storres | |
| 17 | 3 | storres | /* includes of system headers */
|
| 18 | 3 | storres | |
| 19 | 3 | storres | /* Other declarations */
|
| 20 | 3 | storres | |
| 21 | 3 | storres | /* Internal prototypes */
|
| 22 | 3 | storres | void
|
| 23 | 3 | storres | pobyso_error_message(char *functionName, char *messageName, char* message); |
| 24 | 3 | storres | /* Types, constants and macros definitions */
|
| 25 | 3 | storres | |
| 26 | 3 | storres | /* Global variables */
|
| 27 | 3 | storres | |
| 28 | 3 | storres | /* Functions */
|
| 29 | 3 | storres | void
|
| 30 | 3 | storres | pobyso_autoprint(sollya_obj_t solObj, ...) |
| 31 | 3 | storres | {
|
| 32 | 3 | storres | va_list va; |
| 33 | 3 | storres | va_start(va, solObj); |
| 34 | 3 | storres | sollya_lib_v_autoprint(solObj, va); |
| 35 | 3 | storres | va_end(va); |
| 36 | 3 | storres | } /* End pobyso_autoprint. */
|
| 37 | 3 | storres | |
| 38 | 3 | storres | sollya_obj_t |
| 39 | 3 | storres | pobyso_parse_string(const char* expression) |
| 40 | 3 | storres | {
|
| 41 | 3 | storres | if (expression == NULL) |
| 42 | 3 | storres | {
|
| 43 | 3 | storres | pobyso_error_message("pobyso_parse_string",
|
| 44 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 45 | 3 | storres | "The expression is a NULL pointer");
|
| 46 | 3 | storres | return(NULL); |
| 47 | 3 | storres | } |
| 48 | 3 | storres | return(sollya_lib_parse_string(expression));
|
| 49 | 3 | storres | } /* pobyso_parse_string */
|
| 50 | 3 | storres | |
| 51 | 3 | storres | void
|
| 52 | 3 | storres | pobyso_set_verbosity_off(sollya_obj_t* currentVerbosityLevel) |
| 53 | 3 | storres | {
|
| 54 | 3 | storres | sollya_obj_t verbosityLevelZero; |
| 55 | 3 | storres | if (currentVerbosityLevel != NULL) |
| 56 | 3 | storres | {
|
| 57 | 3 | storres | *currentVerbosityLevel = sollya_lib_get_verbosity(); |
| 58 | 3 | storres | } |
| 59 | 3 | storres | verbosityLevelZero = sollya_lib_constant_from_int(0);
|
| 60 | 3 | storres | sollya_lib_set_verbosity(verbosityLevelZero); |
| 61 | 3 | storres | sollya_lib_clear_obj(verbosityLevelZero); |
| 62 | 3 | storres | } /* End of pobyso_set_verbosity_off. */
|
| 63 | 3 | storres | |
| 64 | 3 | storres | void
|
| 65 | 3 | storres | pobyso_set_verbosity_to(sollya_obj_t newVerbosityLevel) |
| 66 | 3 | storres | {
|
| 67 | 3 | storres | if (newVerbosityLevel == NULL) |
| 68 | 3 | storres | {
|
| 69 | 3 | storres | pobyso_error_message("pobyso_set_verbosity_to",
|
| 70 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 71 | 3 | storres | "The new verbosity level is a NULL pointer");
|
| 72 | 3 | storres | return;
|
| 73 | 3 | storres | } |
| 74 | 3 | storres | sollya_lib_set_verbosity(newVerbosityLevel); |
| 75 | 3 | storres | } /* End of pobyso_set_verbosity_to. */
|
| 76 | 3 | storres | |
| 77 | 3 | storres | /* Attic from the sollya_lib < 4. */
|
| 78 | 3 | storres | #if 0
|
| 79 | 3 | storres | chain*
|
| 80 | 3 | storres | pobyso_create_canonical_monomials_base(const unsigned int degree)
|
| 81 | 3 | storres | {
|
| 82 | 3 | storres | int i = 0;
|
| 83 | 3 | storres | chain *monomials = NULL;
|
| 84 | 3 | storres | node *monomial = NULL;
|
| 85 | 3 | storres | |
| 86 | 3 | storres | for(i = degree ; i >= 0 ; i--)
|
| 87 | 3 | storres | {
|
| 88 | 3 | storres | monomial = makePow(makeVariable(), makeConstantDouble((double)i));
|
| 89 | 3 | storres | monomials = addElement(monomials, monomial);
|
| 90 | 3 | storres | fprintf(stderr, "pobyso_create_canonical_monomials_base: %u\n", i);
|
| 91 | 3 | storres | }
|
| 92 | 3 | storres | if (monomials == NULL)
|
| 93 | 3 | storres | {
|
| 94 | 3 | storres | pobyso_error_message("pobyso_create_canonical_monomial_base",
|
| 95 | 3 | storres | "CHAIN_CREATION_ERROR",
|
| 96 | 3 | storres | "Could not create the monomials chain");
|
| 97 | 3 | storres | return(NULL);
|
| 98 | 3 | storres | }
|
| 99 | 3 | storres | return(monomials);
|
| 100 | 3 | storres | } /* End pobyso_create_canonical_monomials_base. */
|
| 101 | 3 | storres | |
| 102 | 3 | storres | chain*
|
| 103 | 3 | storres | pobyso_create_chain_from_int_array(int* intArray,
|
| 104 | 3 | storres | const unsigned int arrayLength)
|
| 105 | 3 | storres | {
|
| 106 | 3 | storres | int i = 0;
|
| 107 | 3 | storres | chain *newChain = NULL;
|
| 108 | 3 | storres | int *currentInt;
|
| 109 | 3 | storres | |
| 110 | 3 | storres | if (arrayLength == 0) return(NULL);
|
| 111 | 3 | storres | if (intArray == NULL)
|
| 112 | 3 | storres | {
|
| 113 | 3 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
| 114 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 115 | 3 | storres | "The array is a NULL pointer");
|
| 116 | 3 | storres | return(NULL);
|
| 117 | 3 | storres | }
|
| 118 | 3 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
| 119 | 3 | storres | {
|
| 120 | 3 | storres | currentInt = malloc(sizeof(int));
|
| 121 | 3 | storres | if (currentInt == NULL)
|
| 122 | 3 | storres | {
|
| 123 | 3 | storres | pobyso_error_message("pobyso_create_chain_from_int_array",
|
| 124 | 3 | storres | "MEMORY_ALLOCATION_ERROR",
|
| 125 | 3 | storres | "Could not allocate one of the integers");
|
| 126 | 3 | storres | freeChain(newChain, free);
|
| 127 | 3 | storres | return(NULL);
|
| 128 | 3 | storres | }
|
| 129 | 3 | storres | *currentInt = intArray[i];
|
| 130 | 3 | storres | newChain = addElement(newChain, currentInt);
|
| 131 | 3 | storres | }
|
| 132 | 3 | storres | return(newChain);
|
| 133 | 3 | storres | } // End pobyso_create_chain_from_int_array. */
|
| 134 | 3 | storres | |
| 135 | 3 | storres | chain*
|
| 136 | 3 | storres | pobyso_create_chain_from_unsigned_int_array(unsigned int* intArray,
|
| 137 | 3 | storres | const unsigned int arrayLength)
|
| 138 | 3 | storres | {
|
| 139 | 3 | storres | int i = 0;
|
| 140 | 3 | storres | chain *newChain = NULL;
|
| 141 | 3 | storres | unsigned int *currentInt;
|
| 142 | 3 | storres | |
| 143 | 3 | storres | /* Argument checking. */
|
| 144 | 3 | storres | if (arrayLength == 0) return(NULL);
|
| 145 | 3 | storres | if (intArray == NULL)
|
| 146 | 3 | storres | {
|
| 147 | 3 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
| 148 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 149 | 3 | storres | "The array is a NULL pointer");
|
| 150 | 3 | storres | return(NULL);
|
| 151 | 3 | storres | }
|
| 152 | 3 | storres | for (i = arrayLength - 1 ; i >= 0 ; i--)
|
| 153 | 3 | storres | {
|
| 154 | 3 | storres | currentInt = malloc(sizeof(unsigned int));
|
| 155 | 3 | storres | if (currentInt == NULL)
|
| 156 | 3 | storres | {
|
| 157 | 3 | storres | pobyso_error_message("pobyso_create_chain_from_unsigned_int_array",
|
| 158 | 3 | storres | "MEMORY_ALLOCATION_ERROR",
|
| 159 | 3 | storres | "Could not allocate one of the integers");
|
| 160 | 3 | storres | freeChain(newChain, free);
|
| 161 | 3 | storres | return(NULL);
|
| 162 | 3 | storres | }
|
| 163 | 3 | storres | *currentInt = intArray[i];
|
| 164 | 3 | storres | newChain = addElement(newChain, currentInt);
|
| 165 | 3 | storres | }
|
| 166 | 3 | storres | return(newChain);
|
| 167 | 3 | storres | } // End pobyso_create_chain_from_unsigned_int_array. */
|
| 168 | 3 | storres | |
| 169 | 3 | storres | node*
|
| 170 | 3 | storres | pobyso_differentiate(node *functionNode)
|
| 171 | 3 | storres | {
|
| 172 | 3 | storres | /* Argument checking. */
|
| 173 | 3 | storres | node *differentialNode;
|
| 174 | 3 | storres | if (functionNode == NULL)
|
| 175 | 3 | storres | {
|
| 176 | 3 | storres | pobyso_error_message("pobyso_differentiate",
|
| 177 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 178 | 3 | storres | "The function to differentiate is a NULL pointer");
|
| 179 | 3 | storres | return(NULL);
|
| 180 | 3 | storres | }
|
| 181 | 3 | storres | differentialNode = differentiate(functionNode);
|
| 182 | 3 | storres | if (differentialNode == NULL)
|
| 183 | 3 | storres | {
|
| 184 | 3 | storres | pobyso_error_message("pobyso_differentiate",
|
| 185 | 3 | storres | "INTERNAL ERROR",
|
| 186 | 3 | storres | "Sollya could not differentiate the function");
|
| 187 | 3 | storres | }
|
| 188 | 3 | storres | return(differentialNode);
|
| 189 | 3 | storres | } // End pobyso_differentiate
|
| 190 | 3 | storres | |
| 191 | 3 | storres | |
| 192 | 3 | storres | int
|
| 193 | 3 | storres | pobyso_dirty_infnorm(mpfr_t infNorm,
|
| 194 | 3 | storres | node *functionNode,
|
| 195 | 3 | storres | mpfr_t lowerBound,
|
| 196 | 3 | storres | mpfr_t upperBound,
|
| 197 | 3 | storres | mp_prec_t precision)
|
| 198 | 3 | storres | {
|
| 199 | 3 | storres | int functionCallResult;
|
| 200 | 3 | storres | /* Arguments checking. */
|
| 201 | 3 | storres | if (functionNode == NULL)
|
| 202 | 3 | storres | {
|
| 203 | 3 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
| 204 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 205 | 3 | storres | "The function to compute is a NULL pointer");
|
| 206 | 3 | storres | return(1);
|
| 207 | 3 | storres | }
|
| 208 | 3 | storres | if (mpfr_cmp(lowerBound, upperBound) > 0)
|
| 209 | 3 | storres | {
|
| 210 | 3 | storres | pobyso_error_message("pobyso_dirty_infnorm",
|
| 211 | 3 | storres | "INCOHERENT_INPUT_DATA",
|
| 212 | 3 | storres | "The lower bond is greater than the upper bound");
|
| 213 | 3 | storres | return(1);
|
| 214 | 3 | storres | }
|
| 215 | 3 | storres | /* Particular cases. */
|
| 216 | 3 | storres | if (mpfr_cmp(lowerBound, upperBound) == 0)
|
| 217 | 3 | storres | {
|
| 218 | 3 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
| 219 | 3 | storres | functionNode,
|
| 220 | 3 | storres | lowerBound,
|
| 221 | 3 | storres | precision);
|
| 222 | 3 | storres | return(functionCallResult);
|
| 223 | 3 | storres | }
|
| 224 | 3 | storres | if (isConstant(functionNode))
|
| 225 | 3 | storres | {
|
| 226 | 3 | storres | functionCallResult = pobyso_evaluate_faithful(infNorm,
|
| 227 | 3 | storres | functionNode,
|
| 228 | 3 | storres | lowerBound,
|
| 229 | 3 | storres | precision);
|
| 230 | 3 | storres | if (!functionCallResult)
|
| 231 | 3 | storres | {
|
| 232 | 3 | storres | mpfr_abs(infNorm, infNorm, MPFR_RNDN);
|
| 233 | 3 | storres | }
|
| 234 | 3 | storres | return(functionCallResult);
|
| 235 | 3 | storres | }
|
| 236 | 3 | storres | uncertifiedInfnorm(infNorm,
|
| 237 | 3 | storres | functionNode,
|
| 238 | 3 | storres | lowerBound,
|
| 239 | 3 | storres | upperBound,
|
| 240 | 3 | storres | POBYSO_DEFAULT_POINTS,
|
| 241 | 3 | storres | precision);
|
| 242 | 3 | storres | return(0);
|
| 243 | 3 | storres | } /* End pobyso_dirty_infnorm. */
|
| 244 | 3 | storres | |
| 245 | 3 | storres | int
|
| 246 | 3 | storres | pobyso_evaluate_faithful(mpfr_t faithfulEvaluation,
|
| 247 | 3 | storres | node *nodeToEvaluate,
|
| 248 | 3 | storres | mpfr_t argument,
|
| 249 | 3 | storres | mpfr_prec_t precision)
|
| 250 | 3 | storres | {
|
| 251 | 3 | storres | /* Check input arguments. */
|
| 252 | 3 | storres | if (nodeToEvaluate == NULL)
|
| 253 | 3 | storres | {
|
| 254 | 3 | storres | pobyso_error_message("pobyso_evaluate_faithful",
|
| 255 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 256 | 3 | storres | "nodeToEvaluate is a NULL pointer");
|
| 257 | 3 | storres | return(1);
|
| 258 | 3 | storres | }
|
| 259 | 3 | storres | evaluateFaithful(faithfulEvaluation,
|
| 260 | 3 | storres | nodeToEvaluate,
|
| 261 | 3 | storres | argument,
|
| 262 | 3 | storres | precision);
|
| 263 | 3 | storres | return(0);
|
| 264 | 3 | storres | } /* End pobyso_evaluate_faithfull. */
|
| 265 | 3 | storres | |
| 266 | 3 | storres | chain*
|
| 267 | 3 | storres | pobyso_find_zeros(node *function,
|
| 268 | 3 | storres | mpfr_t *lower_bound,
|
| 269 | 3 | storres | mpfr_t *upper_bound)
|
| 270 | 3 | storres | {
|
| 271 | 3 | storres | mp_prec_t currentPrecision;
|
| 272 | 3 | storres | mpfr_t currentDiameter;
|
| 273 | 3 | storres | rangetype bounds;
|
| 274 | 3 | storres | |
| 275 | 3 | storres | currentPrecision = getToolPrecision();
|
| 276 | 3 | storres | mpfr_init2(currentDiameter, currentPrecision);
|
| 277 | 3 | storres | |
| 278 | 3 | storres | bounds.a = lower_bound;
|
| 279 | 3 | storres | bounds.b = upper_bound;
|
| 280 | 3 | storres | |
| 281 | 3 | storres | if (bounds.a == NULL || bounds.b == NULL)
|
| 282 | 3 | storres | {
|
| 283 | 3 | storres | pobyso_error_message("pobyso_find_zeros",
|
| 284 | 3 | storres | "MEMORY_ALLOCATION_ERROR",
|
| 285 | 3 | storres | "Could not allocate one of the bounds");
|
| 286 | 3 | storres | return(NULL);
|
| 287 | 3 | storres | }
|
| 288 | 3 | storres | return(findZerosFunction(function,
|
| 289 | 3 | storres | bounds,
|
| 290 | 3 | storres | currentPrecision,
|
| 291 | 3 | storres | currentDiameter));
|
| 292 | 3 | storres | } /* End pobyso_find_zeros. */
|
| 293 | 3 | storres | |
| 294 | 3 | storres | void
|
| 295 | 3 | storres | pobyso_free_chain_of_nodes(chain *theChainOfNodes)
|
| 296 | 3 | storres | {
|
| 297 | 3 | storres | node *currentNode = NULL;
|
| 298 | 3 | storres | chain *currentChainElement = NULL;
|
| 299 | 3 | storres | chain *nextChainElement = NULL;
|
| 300 | 3 | storres | |
| 301 | 3 | storres | nextChainElement = theChainOfNodes;
|
| 302 | 3 | storres | |
| 303 | 3 | storres | while(nextChainElement != NULL)
|
| 304 | 3 | storres | {
|
| 305 | 3 | storres | currentChainElement = nextChainElement;
|
| 306 | 3 | storres | currentNode = (node*)(currentChainElement->value);
|
| 307 | 3 | storres | nextChainElement = nextChainElement->next;
|
| 308 | 3 | storres | free_memory(currentNode);
|
| 309 | 3 | storres | free((void*)(currentChainElement));
|
| 310 | 3 | storres | }
|
| 311 | 3 | storres | } /* End pobyso_free_chain_of_nodes. */
|
| 312 | 3 | storres | |
| 313 | 3 | storres | void
|
| 314 | 3 | storres | pobyso_free_range(rangetype range)
|
| 315 | 3 | storres | {
|
| 316 | 3 | storres | |
| 317 | 3 | storres | mpfr_clear(*(range.a));
|
| 318 | 3 | storres | mpfr_clear(*(range.b));
|
| 319 | 3 | storres | free(range.a);
|
| 320 | 3 | storres | free(range.b);
|
| 321 | 3 | storres | } /* End pobyso_free_range. */
|
| 322 | 3 | storres | |
| 323 | 3 | storres | node*
|
| 324 | 3 | storres | pobyso_fp_minimax_canonical_monomials_base(node *function,
|
| 325 | 3 | storres | int degree,
|
| 326 | 3 | storres | chain *formats,
|
| 327 | 3 | storres | chain *points,
|
| 328 | 3 | storres | mpfr_t lowerBound,
|
| 329 | 3 | storres | mpfr_t upperBound,
|
| 330 | 3 | storres | int fpFixedArg,
|
| 331 | 3 | storres | int absRel,
|
| 332 | 3 | storres | node *constPart,
|
| 333 | 3 | storres | node *minimax)
|
| 334 | 3 | storres | {
|
| 335 | 3 | storres | return(NULL);
|
| 336 | 3 | storres | } /* End pobyso_fp_minimax_canonical_monomials_base. */
|
| 337 | 3 | storres | |
| 338 | 3 | storres | node*
|
| 339 | 3 | storres | pobyso_parse_function(char *functionString,
|
| 340 | 3 | storres | char *freeVariableNameString)
|
| 341 | 3 | storres | {
|
| 342 | 3 | storres | if (functionString == NULL || freeVariableNameString == NULL)
|
| 343 | 3 | storres | {
|
| 344 | 3 | storres | pobyso_error_message("pobyso_parse_function",
|
| 345 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 346 | 3 | storres | "One of the arguments is a NULL pointer");
|
| 347 | 3 | storres | return(NULL);
|
| 348 | 3 | storres | }
|
| 349 | 3 | storres | return(parseString(functionString));
|
| 350 | 3 | storres | |
| 351 | 3 | storres | } /* End pobyso_parse_function */
|
| 352 | 3 | storres | |
| 353 | 3 | storres | node*
|
| 354 | 3 | storres | pobyso_remez_approx_canonical_monomials_base_for_error(node *functionNode,
|
| 355 | 3 | storres | unsigned int mode,
|
| 356 | 3 | storres | mpfr_t lowerBound,
|
| 357 | 3 | storres | mpfr_t upperBound,
|
| 358 | 3 | storres | mpfr_t eps)
|
| 359 | 3 | storres | {
|
| 360 | 3 | storres | node *weight = NULL;
|
| 361 | 3 | storres | node *bestApproxPolyNode = NULL;
|
| 362 | 3 | storres | node *bestApproxHorner = NULL;
|
| 363 | 3 | storres | node *errorNode = NULL;
|
| 364 | 3 | storres | rangetype degreeRange;
|
| 365 | 3 | storres | mpfr_t quality;
|
| 366 | 3 | storres | mpfr_t currentError;
|
| 367 | 3 | storres | unsigned int degree;
|
| 368 | 3 | storres | |
| 369 | 3 | storres | /* Check the parameters. */
|
| 370 | 3 | storres | if (functionNode == NULL)
|
| 371 | 3 | storres | {
|
| 372 | 3 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
| 373 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 374 | 3 | storres | "functionNode is a NULL pointer");
|
| 375 | 3 | storres | return(NULL);
|
| 376 | 3 | storres | }
|
| 377 | 3 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
| 378 | 3 | storres | {
|
| 379 | 3 | storres | pobyso_error_message("remezApproxCanonicalMonomialsBaseForError",
|
| 380 | 3 | storres | "INCOHERENT_INPUT_DATA",
|
| 381 | 3 | storres | "the lower_bound >= upper_bound");
|
| 382 | 3 | storres | return(NULL);
|
| 383 | 3 | storres | }
|
| 384 | 3 | storres | /* Set the weight. */
|
| 385 | 3 | storres | if (mode == POBYSO_ABSOLUTE)
|
| 386 | 3 | storres | {
|
| 387 | 3 | storres | /* Set the weight to 1 for the ABSOLUTE_MODE. */
|
| 388 | 3 | storres | weight = makeConstantDouble(1.0);
|
| 389 | 3 | storres | }
|
| 390 | 3 | storres | else
|
| 391 | 3 | storres | {
|
| 392 | 3 | storres | if (mode == POBYSO_RELATIVE)
|
| 393 | 3 | storres | {
|
| 394 | 3 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
| 395 | 3 | storres | "NOT_IMPLEMENTED",
|
| 396 | 3 | storres | "the search for relative error is not implemented yet");
|
| 397 | 3 | storres | return(NULL);
|
| 398 | 3 | storres | }
|
| 399 | 3 | storres | else
|
| 400 | 3 | storres | {
|
| 401 | 3 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
| 402 | 3 | storres | "INCOHERENT_INPUT_DATA",
|
| 403 | 3 | storres | "the mode is node of POBYSO_ABOLUTE or POBYSO_RELATIVE");
|
| 404 | 3 | storres | return(NULL);
|
| 405 | 3 | storres | }
|
| 406 | 3 | storres | }
|
| 407 | 3 | storres | //fprintf(stderr, "\n\n\n******* I'm here! ********\n\n\n");
|
| 408 | 3 | storres | degreeRange = guessDegree(functionNode,
|
| 409 | 3 | storres | weight,
|
| 410 | 3 | storres | lowerBound,
|
| 411 | 3 | storres | upperBound,
|
| 412 | 3 | storres | eps,
|
| 413 | 3 | storres | POBYSO_GUESS_DEGREE_BOUND);
|
| 414 | 3 | storres | degree = mpfr_get_ui(*(degreeRange.a), MPFR_RNDN);
|
| 415 | 3 | storres | //fprintf(stderr, "\n\n\n******* I'm back! ********\n\n\n");
|
| 416 | 3 | storres | fprintf(stderr, "Guessed degree: ");
|
| 417 | 3 | storres | mpfr_out_str(stderr, 10, 17, *(degreeRange.a), MPFR_RNDN);
|
| 418 | 3 | storres | fprintf(stderr, " - as int: %u\n", degree);
|
| 419 | 3 | storres | /* Reduce the degree by 1 in the foolish hope it could work. */
|
| 420 | 3 | storres | if (degree > 0) degree--;
|
| 421 | 3 | storres | /* Both elements of degreeRange where "inited" within guessDegree. */
|
| 422 | 3 | storres | mpfr_clear(*(degreeRange.a));
|
| 423 | 3 | storres | mpfr_clear(*(degreeRange.b));
|
| 424 | 3 | storres | free(degreeRange.a);
|
| 425 | 3 | storres | free(degreeRange.b);
|
| 426 | 3 | storres | /* Mimic the default behavior of interactive Sollya. */
|
| 427 | 3 | storres | mpfr_init(quality);
|
| 428 | 3 | storres | mpfr_set_d(quality, 1e-5, MPFR_RNDN);
|
| 429 | 3 | storres | mpfr_init2(currentError, getToolPrecision());
|
| 430 | 3 | storres | mpfr_set_inf(currentError,1);
|
| 431 | 3 | storres | |
| 432 | 3 | storres | /* Try to refine the initial guess: loop with increasing degrees until
|
| 433 | 3 | storres | * we find a satisfactory one. */
|
| 434 | 3 | storres | while(mpfr_cmp(currentError, eps) > 0)
|
| 435 | 3 | storres | {
|
| 436 | 3 | storres | /* Get rid of the previous polynomial, if any. */
|
| 437 | 3 | storres | if (bestApproxPolyNode != NULL)
|
| 438 | 3 | storres | {
|
| 439 | 3 | storres | free_memory(bestApproxPolyNode);
|
| 440 | 3 | storres | }
|
| 441 | 3 | storres | fprintf(stderr, "Degree: %u\n", degree);
|
| 442 | 3 | storres | fprintf(stderr, "Calling pobyso_remez_canonical_monomials_base...\n");
|
| 443 | 3 | storres | /* Try to find a polynomial with the guessed degree. */
|
| 444 | 3 | storres | bestApproxPolyNode = pobyso_remez_canonical_monomials_base(functionNode,
|
| 445 | 3 | storres | weight,
|
| 446 | 3 | storres | degree,
|
| 447 | 3 | storres | lowerBound,
|
| 448 | 3 | storres | upperBound,
|
| 449 | 3 | storres | quality);
|
| 450 | 3 | storres | |
| 451 | 3 | storres | if (bestApproxPolyNode == NULL)
|
| 452 | 3 | storres | {
|
| 453 | 3 | storres | pobyso_error_message("computeRemezApproxCanonicalMonomialsBaseForError",
|
| 454 | 3 | storres | "INTERNAL_ERROR",
|
| 455 | 3 | storres | "could not compute the bestApproxPolyNode");
|
| 456 | 3 | storres | mpfr_clear(currentError);
|
| 457 | 3 | storres | mpfr_clear(quality);
|
| 458 | 3 | storres | return(NULL);
|
| 459 | 3 | storres | }
|
| 460 | 3 | storres | |
| 461 | 3 | storres | setDisplayMode(DISPLAY_MODE_DECIMAL);
|
| 462 | 3 | storres | fprintTree(stderr, bestApproxPolyNode);
|
| 463 | 3 | storres | fprintf(stderr, "\n\n");
|
| 464 | 3 | storres | |
| 465 | 3 | storres | errorNode = makeSub(copyTree(functionNode), copyTree(bestApproxPolyNode));
|
| 466 | 3 | storres | /* Check the error with the computed polynomial. */
|
| 467 | 3 | storres | uncertifiedInfnorm(currentError,
|
| 468 | 3 | storres | errorNode,
|
| 469 | 3 | storres | lowerBound,
|
| 470 | 3 | storres | upperBound,
|
| 471 | 3 | storres | POBYSO_INF_NORM_NUM_POINTS,
|
| 472 | 3 | storres | getToolPrecision());
|
| 473 | 3 | storres | fprintf(stderr, "Inf norm: ");
|
| 474 | 3 | storres | mpfr_out_str(stderr, 10, 17, currentError, MPFR_RNDN);
|
| 475 | 3 | storres | fprintf(stderr, "\n\n");
|
| 476 | 3 | storres | /* Free the errorNode but not the bestApproxPolyNode (we need it if
|
| 477 | 3 | storres | * we exit the loop at the next iteration). */
|
| 478 | 3 | storres | free_memory(errorNode);
|
| 479 | 3 | storres | degree++;
|
| 480 | 3 | storres | }
|
| 481 | 3 | storres | /* Use an intermediate variable, since horner() creates a new node
|
| 482 | 3 | storres | * and does not reorder the argument "in place". This allows for the memory
|
| 483 | 3 | storres | * reclaim of bestApproxHorner.
|
| 484 | 3 | storres | */
|
| 485 | 3 | storres | bestApproxHorner = horner(bestApproxPolyNode);
|
| 486 | 3 | storres | free_memory(bestApproxPolyNode);
|
| 487 | 3 | storres | mpfr_clear(currentError);
|
| 488 | 3 | storres | mpfr_clear(quality);
|
| 489 | 3 | storres | free_memory(weight);
|
| 490 | 3 | storres | return(bestApproxHorner);
|
| 491 | 3 | storres | } /* End pobyso_remez_approx_canonical_monomials_base_for_error */
|
| 492 | 3 | storres | |
| 493 | 3 | storres | node*
|
| 494 | 3 | storres | pobyso_remez_canonical_monomials_base(node *function,
|
| 495 | 3 | storres | node *weight,
|
| 496 | 3 | storres | unsigned int degree,
|
| 497 | 3 | storres | mpfr_t lowerBound,
|
| 498 | 3 | storres | mpfr_t upperBound,
|
| 499 | 3 | storres | mpfr_t quality)
|
| 500 | 3 | storres | {
|
| 501 | 3 | storres | node *bestApproxPoly = NULL;
|
| 502 | 3 | storres | chain *monomials = NULL;
|
| 503 | 3 | storres | chain *curMonomial = NULL;
|
| 504 | 3 | storres | |
| 505 | 3 | storres | mpfr_t satisfying_error;
|
| 506 | 3 | storres | mpfr_t target_error;
|
| 507 | 3 | storres | |
| 508 | 3 | storres | /* Argument checking */
|
| 509 | 3 | storres | /* Function tree. */
|
| 510 | 3 | storres | if (function == NULL)
|
| 511 | 3 | storres | {
|
| 512 | 3 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
| 513 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 514 | 3 | storres | "the \"function\" argument is a NULL pointer");
|
| 515 | 3 | storres | return(NULL);
|
| 516 | 3 | storres | }
|
| 517 | 3 | storres | if (weight == NULL)
|
| 518 | 3 | storres | {
|
| 519 | 3 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
| 520 | 3 | storres | "NULL_POINTER_ARGUMENT",
|
| 521 | 3 | storres | "the \"weight\" argument is a NULL pointer");
|
| 522 | 3 | storres | return(NULL);
|
| 523 | 3 | storres | }
|
| 524 | 3 | storres | /* Check the bounds. */
|
| 525 | 3 | storres | if (mpfr_cmp(lowerBound, upperBound) >= 0)
|
| 526 | 3 | storres | {
|
| 527 | 3 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
| 528 | 3 | storres | "INCOHERENT_IMPUT_DATA",
|
| 529 | 3 | storres | "the lower_bound >= upper_bound");
|
| 530 | 3 | storres | return(NULL);
|
| 531 | 3 | storres | }
|
| 532 | 3 | storres | /* The quality must be a non null positive number. */
|
| 533 | 3 | storres | if (mpfr_sgn(quality) <= 0)
|
| 534 | 3 | storres | {
|
| 535 | 3 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
| 536 | 3 | storres | "INCOHERENT_INPUT_DATA",
|
| 537 | 3 | storres | "the quality <= 0");
|
| 538 | 3 | storres | }
|
| 539 | 3 | storres | /* End argument checking. */
|
| 540 | 3 | storres | /* Create the monomials nodes chains. */
|
| 541 | 3 | storres | monomials = pobyso_create_canonical_monomials_base(degree);
|
| 542 | 3 | storres | fprintf(stderr, "monomials chain length = %d\n", lengthChain(monomials));
|
| 543 | 3 | storres | if (monomials == NULL || (lengthChain(monomials) != degree + 1))
|
| 544 | 3 | storres | {
|
| 545 | 3 | storres | pobyso_error_message("pobyso_remez_canonical_monomials_base",
|
| 546 | 3 | storres | "CHAIN_CREATION_ERROR",
|
| 547 | 3 | storres | "could not create the monomials chain");
|
| 548 | 3 | storres | return(NULL);
|
| 549 | 3 | storres | }
|
| 550 | 3 | storres | curMonomial = monomials;
|
| 551 | 3 | storres | |
| 552 | 3 | storres | while (curMonomial != NULL)
|
| 553 | 3 | storres | {
|
| 554 | 3 | storres | fprintf(stderr, "monomial tree: ");
|
| 555 | 3 | storres | //mpfr_out_str(stderr, 10, 17, *((mpfr_t*)((node*)(curMonomial->value))->value), MPFR_RNDN);
|
| 556 | 3 | storres | fprintTree(stderr, (node*)(curMonomial->value));
|
| 557 | 3 | storres | fprintf(stderr, "\n");
|
| 558 | 3 | storres | curMonomial = curMonomial->next;
|
| 559 | 3 | storres | }
|
| 560 | 3 | storres | |
| 561 | 3 | storres | /* Deal with NULL weight. */
|
| 562 | 3 | storres | if (weight == NULL)
|
| 563 | 3 | storres | {
|
| 564 | 3 | storres | weight = makeConstantDouble(1.0);
|
| 565 | 3 | storres | }
|
| 566 | 3 | storres | /* Compute the best polynomial with the required quality.
|
| 567 | 3 | storres | The behavior is as if satisfying error and target_error had
|
| 568 | 3 | storres | not been used.*/
|
| 569 | 3 | storres | mpfr_init(satisfying_error);
|
| 570 | 3 | storres | mpfr_init(target_error);
|
| 571 | 3 | storres | mpfr_set_str(satisfying_error, "0", 10, MPFR_RNDN);
|
| 572 | 3 | storres | mpfr_set_inf(target_error, 1);
|
| 573 | 3 | storres | |
| 574 | 3 | storres | |
| 575 | 3 | storres | fprintf(stderr, "satisfying_error: ");
|
| 576 | 3 | storres | mpfr_out_str(stderr, 10, 17, satisfying_error, MPFR_RNDN);
|
| 577 | 3 | storres | fprintf(stderr, ".\n");
|
| 578 | 3 | storres | fprintf(stderr, "target_error: ");
|
| 579 | 3 | storres | mpfr_out_str(stderr, 10, 17, target_error,MPFR_RNDN);
|
| 580 | 3 | storres | fprintf(stderr, ".\n");
|
| 581 | 3 | storres | |
| 582 | 3 | storres | fprintf(stderr,
|
| 583 | 3 | storres | "current precision: %li\n", getToolPrecision());
|
| 584 | 3 | storres | /* Call the Sollya function. */
|
| 585 | 3 | storres | bestApproxPoly = remez(function,
|
| 586 | 3 | storres | weight,
|
| 587 | 3 | storres | monomials,
|
| 588 | 3 | storres | lowerBound,
|
| 589 | 3 | storres | upperBound,
|
| 590 | 3 | storres | quality,
|
| 591 | 3 | storres | satisfying_error,
|
| 592 | 3 | storres | target_error,
|
| 593 | 3 | storres | getToolPrecision());
|
| 594 | 3 | storres | |
| 595 | 3 | storres | mpfr_clear(satisfying_error);
|
| 596 | 3 | storres | mpfr_clear(target_error);
|
| 597 | 3 | storres | pobyso_free_chain_of_nodes(monomials);
|
| 598 | 3 | storres | |
| 599 | 3 | storres | return(bestApproxPoly);
|
| 600 | 3 | storres | } /* End pobyso_remez_canonical_monomials_base. */
|
| 601 | 3 | storres | |
| 602 | 3 | storres | #endif
|
| 603 | 3 | storres | |
| 604 | 3 | storres | void
|
| 605 | 3 | storres | pobyso_error_message(char *functionName, char *messageName, char* message) |
| 606 | 3 | storres | {
|
| 607 | 3 | storres | fprintf(stderr, "?%s: %s.\n%s.\n", functionName, messageName, message);
|
| 608 | 3 | storres | } /* End pobyso_error_message */ |