Statistiques
| Révision :

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 */