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
00018
00019
00020 CouNumber Q (register int k, CouNumber x) {
00021
00022 register CouNumber xp, Q;
00023
00024 k *= 2;
00025
00026 xp = x;
00027 Q = 1;
00028
00029 for (register int i=2; i<=k; i++) {
00030
00031 Q += (CouNumber) i * xp;
00032 xp *= x;
00033 }
00034
00035 return Q;
00036 }
00037
00038
00039
00040
00041
00042
00043
00044 CouNumber rootQ (int k) {
00045
00046 if (k==1) return - 0.5;
00047 else {
00048
00049 register CouNumber l = - 1.0 + 0.5 / k,
00050 u = - 0.5,
00051 Ql = Q (k, l), Qu = Q (k, u), 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; Ql = Qm; Qm = - Qm;}
00061 else {u = midpoint; Qu = Qm;}
00062
00063 } while (Qm > 1e-15);
00064
00065 return midpoint;
00066 }
00067 }
00068
00069 #ifdef DEBUG_ROOTQ
00070 int main () {
00071
00072 register int k;
00073 CouNumber x, q;
00074
00075 for (k=6; --k;) {
00076
00077 printf ("root, %3d -> %.15f\n", 2*k+1, rootQ (k));
00078
00079
00080
00081
00082
00083
00084
00085
00086 }
00087 }
00088 #endif