source: branches/Couenne/Couenne/src/util/rootQ.c @ 498

Last change on this file since 498 was 498, checked in by pbelotti, 13 years ago

fixed serious bug with x(2k+1)

File size: 1.6 KB
Line 
1/*
2 * Name:    rootQ.c
3 * Author:  Pietro Belotti
4 * Purpose: find roots of polynomial Q^k(x) (see Liberti and Pantelides, 2003)
5 *
6 * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL)
7 */
8
9#include <math.h>
10
11#include <CouenneTypes.h>
12#include <CouennePrecisions.h>
13#include <stdio.h>
14
15
16/* compute Q(x)*/
17
18CouNumber Q (register int k, CouNumber x) {
19
20  register CouNumber xp, Q;
21  register int i;
22
23  k *= 2;
24
25  xp = x;
26  Q = 1;
27
28  for (i=2; i<=k; i++) {
29
30    Q += (CouNumber) i * xp;
31    xp *= x;
32  }
33
34  return Q;
35}
36
37
38/*
39 * Find roots of polynomial $Q^k(x) = \sum_{i=1}^{2k} i x^{i-1}$. Used
40 *  in convexification of powers with odd exponent
41 */
42
43CouNumber rootQ (int k) {
44
45  if (k==1) return - 0.5;
46  else {
47
48    register CouNumber l  = - 1.0 + 0.5 / k, 
49                       u  = - 0.5,
50                       Ql = Q (k, l), Qu = Q (k, u), Qm,
51                       midpoint;
52    do {
53
54      midpoint = 0.5 * (l+u); /* (- Ql * u + Qu * l) / (Qu - Ql); */
55      Qm = Q (k, midpoint);
56
57      /*      printf ("[%.4f, %.4f] --> %.4f: %.24f\n", l, u, midpoint, Qm); */
58
59      if (Qm<0) {l = midpoint; Ql = Qm; Qm = - Qm;}
60      else      {u = midpoint; Qu = Qm;}
61
62    } while (Qm > 1e-15);
63
64    return midpoint;
65  }
66}
67
68#ifdef DEBUG_ROOTQ
69int main () {
70
71  register int k;
72  CouNumber x, q;
73
74  for (k=6; --k;) {
75
76    printf ("root, %3d -> %.15f\n", 2*k+1, rootQ (k));
77    /*
78    printf ("k=%3d: ", 2*k+1);
79    for (x = -1.0; x < 0.4; x += 0.1) {
80      //      Q (k, x, &q, NULL, NULL);
81      printf ("[%.2f, %.3f] ", x, rootQ);
82    }
83    printf ("\n");
84    */
85  }
86}
87#endif
Note: See TracBrowser for help on using the repository browser.