00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012
00013 #include "CouenneTypes.hpp"
00014 #include "CouenneExprPow.hpp"
00015 #include "CouenneExprConst.hpp"
00016 #include "CouenneExprClone.hpp"
00017 #include "CouenneExprMax.hpp"
00018 #include "CouenneExprMin.hpp"
00019 #include "CouenneExprOpp.hpp"
00020 #include "CouennePrecisions.hpp"
00021 #include "CouenneProblem.hpp"
00022
00023 #include "CoinHelperFunctions.hpp"
00024 #include "CoinFinite.hpp"
00025
00026 using namespace Couenne;
00027
00028
00029
00030
00031 void exprPow::getBounds (expression *&lb, expression *&ub) {
00032
00033
00034
00035
00036 assert (arglist_ [0] -> Type () != CONST);
00037
00038
00039
00040
00041
00042 expression *lbbase, *ubbase;
00043 arglist_ [0] -> getBounds (lbbase, ubbase);
00044
00045
00046
00047 if (arglist_ [1] -> Type () == CONST) {
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 CouNumber expon = arglist_ [1] -> Value ();
00060 int rndexp;
00061
00062 bool isInt = fabs (expon - (rndexp = COUENNE_round (expon))) < COUENNE_EPS,
00063 isInvInt = !isInt &&
00064 ((fabs (expon) > COUENNE_EPS) &&
00065 (fabs (1/expon - (rndexp = COUENNE_round (1/expon))) < COUENNE_EPS));
00066
00067 if ((isInt || isInvInt) && (rndexp % 2) && (rndexp > 0)) {
00068
00069
00070
00071
00072 lb = new exprPow (lbbase, new exprConst (expon));
00073 ub = new exprPow (ubbase, new exprConst (expon));
00074 }
00075 else {
00076
00077
00078
00079 expression **all = new expression * [6];
00080
00081 all [0] = new exprOpp (lbbase);
00082 all [2] = new exprConst (0.);
00083 all [4] = ubbase;
00084
00085 if (expon > 0)
00086 all [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00087 else all [1] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00088
00089
00090
00091 if (expon > COUENNE_EPS) all [3] = new exprConst (0.);
00092 else if (isInt || isInvInt) {
00093 if (rndexp % 2)
00094 all [3] = new exprConst (-COUENNE_INFINITY);
00095 else all [3] = new exprMin (new exprClone (all [1]),
00096 new exprPow (new exprClone (lbbase),
00097 new exprConst (expon)));
00098 }
00099 else all [3] = new exprClone (all [1]);
00100
00101
00102
00103 if (expon > COUENNE_EPS) {
00104 if (isInt && !(rndexp % 2))
00105 all [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00106 else all [5] = new exprConst (0.);
00107 }
00108 else {
00109 if (isInt || isInvInt) {
00110 if (rndexp % 2)
00111 all [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00112 else all [5] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00113 }
00114 else all [5] = new exprConst (0.);
00115 }
00116
00117 lb = new exprMin (all, 6);
00118
00119
00120
00121 if (expon > 0) {
00122
00123
00124
00125
00126 ub = new exprMax (new exprPow (new exprClone (lbbase), new exprConst (expon)),
00127 new exprPow (new exprClone (ubbase), new exprConst (expon)));
00128
00129 } else {
00130
00131 expression **alu = new expression * [6];
00132
00133 alu [0] = new exprClone (all [0]);
00134 alu [2] = new exprConst (0.);
00135 alu [4] = new exprClone (ubbase);
00136
00137
00138
00139
00140
00141
00142
00143 alu [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00144
00145
00146
00147
00148 alu [3] = new exprConst (COUENNE_INFINITY);
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 if (isInt || isInvInt)
00164 alu [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00165 else alu [5] = new exprConst (COUENNE_INFINITY);
00166
00167
00168 ub = new exprMin (alu, 6);
00169 }
00170 }
00171 }
00172 else
00173 printf ("exprPow::getBounds(): Warning, exponent not constant\n");
00174
00175
00176
00177
00178
00179
00180
00181 }
00182
00183
00184
00185 void exprPow::getBounds (CouNumber &lb, CouNumber &ub) {
00186
00187 CouNumber lba, uba, k = (*(arglist_ [1])) ();
00188 arglist_ [0] -> getBounds (lba, uba);
00189 int intk;
00190
00191 bool
00192 isInt = fabs (k - (double) (intk = COUENNE_round (k))) < COUENNE_EPS,
00193 isInvInt = !isInt && fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS;
00194
00195 if (!isInt && (!isInvInt || !(intk % 2))) {
00196
00197
00198
00199 if (lba < 0.) lba = 0.;
00200 if (uba < 0.) uba = 0.;
00201 }
00202
00203 if (isInt && !(intk % 2) && (k > 0)) {
00204
00205 if (uba < 0) {
00206 lb = safe_pow (-uba, k);
00207 ub = safe_pow (-lba, k);
00208 } else if (lba > 0) {
00209 lb = safe_pow (lba, k);
00210 ub = safe_pow (uba, k);
00211 } else {
00212 lb = 0;
00213 ub = safe_pow (CoinMax (-lba, uba), k);
00214 }
00215
00216 } else if (k > 0) {
00217
00218 lb = safe_pow (lba, k);
00219 ub = safe_pow (uba, k);
00220
00221 } else if (isInt && !(intk % 2)) {
00222
00223 if (uba < 0) {
00224 lb = safe_pow (-lba, k);
00225 ub = safe_pow (-uba, k);
00226 } else if (lba > 0) {
00227 lb = safe_pow (uba, k);
00228 ub = safe_pow (lba, k);
00229 } else {
00230 lb = safe_pow (CoinMax (-lba, uba), k);
00231 ub = COUENNE_INFINITY;
00232 }
00233
00234 } else {
00235
00236 if (uba < 0) {
00237 lb = safe_pow (uba, k);
00238 ub = safe_pow (lba, k);
00239 } else if (lba > 0) {
00240 lb = safe_pow (uba, k);
00241 ub = safe_pow (lba, k);
00242 } else {
00243 lb = -COIN_DBL_MAX;
00244 ub = COIN_DBL_MAX;
00245 }
00246 }
00247 }