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

18 | CouNumber 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 | |
---|

43 | CouNumber 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 |
---|

69 | int 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 |
---|