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