00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012 #include <stdio.h>
00013
00014 #include "CouenneTypes.hpp"
00015 #include "CouennePrecisions.hpp"
00016
00017 namespace Couenne {
00018
00019
00020
00021 CouNumber Q (register int k, CouNumber x) {
00022
00023 register CouNumber xp = x, Q = 1.;
00024
00025 k *= 2;
00026
00027 for (register int i=2; i<=k; i++) {
00028
00029 Q += (CouNumber) i * xp;
00030 xp *= x;
00031 }
00032
00033 return Q;
00034 }
00035
00036
00037
00038
00039
00040
00041
00042 CouNumber rootQ (int k) {
00043
00044 if (k==1) return - 0.5;
00045 else {
00046
00047 register CouNumber
00048 l = - 1.0 + 0.5 / k,
00049 u = - 0.5,
00050
00051 Qm,
00052 midpoint;
00053 do {
00054
00055 midpoint = 0.5 * (l+u);
00056 Qm = Q (k, midpoint);
00057
00058
00059
00060 if (Qm<0) {l = midpoint; Qm = - Qm;}
00061 else {u = midpoint; }
00062
00063 } while (Qm > 1e-15);
00064
00065 return midpoint;
00066 }
00067 }
00068
00069 }
00070
00071 #ifdef DEBUG_ROOTQ
00072 int main () {
00073
00074 register int k;
00075 CouNumber x, q;
00076
00077 for (k=6; --k;) {
00078
00079 printf ("root, %3d -> %.15f\n", 2*k+1, rootQ (k));
00080
00081
00082
00083
00084
00085
00086
00087
00088 }
00089 }
00090
00091 #endif