Changeset 1236
 Timestamp:
 Jan 16, 2015 10:36:37 AM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cgl/src/CglGMI/CglGMI.cpp
r1122 r1236 27 27 #include "CglGMI.hpp" 28 28 #include "CoinFinite.hpp" 29 #include "CoinRational.hpp" 29 30 30 31 // … … 890 891 long numerator = 0, denominator = 0; 891 892 // Initialize gcd and lcm 892 if (nearestRational(cutRhs, maxdelta, maxdnom, numerator, denominator)) { 893 gcd = labs(numerator); 894 lcm = denominator; 893 CoinRational r = CoinRational(cutRhs, maxdelta, maxdnom); 894 if (r.getNumerator() != 0){ 895 gcd = labs(r.getNumerator()); 896 lcm = r.getDenominator(); 895 897 } 896 898 else{ … … 904 906 continue; 905 907 } 906 if(nearestRational(cutElem[i], maxdelta, maxdnom, numerator, denominator)) { 907 gcd = computeGcd(gcd,labs(numerator)); 908 lcm *= denominator/(computeGcd(lcm,denominator)); 908 CoinRational r = CoinRational(cutElem[i], maxdelta, maxdnom); 909 if (r.getNumerator() != 0){ 910 gcd = computeGcd(gcd, r.getNumerator()); 911 lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator())); 909 912 } 910 913 else{ … … 929 932 return true; 930 933 } /* scaleCutIntegral */ 931 932 /************************************************************************/933 /* arguments:934 * val = double precision value that must be converted935 * maxdelta = max allowed difference between val and the rational computed936 * maxdnom = max allowed denominator937 * numerator = the numerator will be stored here if successful938 * denominator = the denominator will be stored here if successful939 * returns true if successful, false if not.940 *941 * This function is based on SCIPrealToRational() from SCIP, scip@zib.de.942 * The copyright of SCIP and of this function belongs to ZIB.943 * We explicitly obtained the rights to license this function under GPL944 * from ZIB. More information can be obtained from the authors.945 *946 * Copyright (C) 2012 KonradZuseZentrum947 * fuer Informationstechnik Berlin948 */949 bool CglGMI::nearestRational(double val, double maxdelta, long maxdnom,950 long& numerator, long& denominator)951 {952 953 /// Denominators that should be tried for the integral scaling phase.954 /// These values are taken from SCIP.955 static const double simplednoms[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0,956 8.0, 9.0, 11.0, 12.0, 13.0, 14.0,957 15.0, 16.0, 17.0, 18.0, 19.0, 25.0,958 1.0};959 960 double a, b;961 double g0, g1, gx;962 double h0, h1, hx;963 double delta0, delta1;964 double epsilon;965 int i;966 967 /* try the simple denominators first: each value of the simplednoms table968 * multiplied by powers of 10 is tried as denominator969 */970 for (i = 0; simplednoms[i] > 0.0; ++i) {971 double num, dnom;972 double ratval0, ratval1;973 double diff;974 975 /* try powers of 10 (including 10^0) */976 dnom = simplednoms[i];977 while (dnom <= maxdnom) {978 num = floor(val * dnom);979 ratval0 = num/dnom;980 ratval1 = (num+1.0)/dnom;981 diff = fabs(val  ratval0);982 if (diff < maxdelta) {983 numerator = (long)num;984 denominator = (long)dnom;985 return true;986 }987 diff = fabs(val  ratval1);988 if (diff < maxdelta) {989 numerator = (long)(num+1.0);990 denominator = (long)dnom;991 return true;992 }993 dnom *= 10.0;994 }995 }996 997 /* the simple denominators didn't work: calculate rational998 * representation with arbitrary denominator */999 epsilon = maxdelta/2.0;1000 1001 b = val;1002 a = floor(b + epsilon);1003 g0 = a;1004 h0 = 1.0;1005 g1 = 1.0;1006 h1 = 0.0;1007 delta0 = val  g0/h0;1008 delta1 = (delta0 < 0.0 ? val  (g01.0)/h0 : val  (g0+1.0)/h0);1009 1010 while ((fabs(delta0) > maxdelta) && (fabs(delta1) > maxdelta)) {1011 if ((ba) < epsilon  h0 < 0  h1 < 0)1012 return false;1013 1014 b = 1.0 / (b  a);1015 a = floor(b + epsilon);1016 1017 if (a < 0.0)1018 return false;1019 gx = g0;1020 hx = h0;1021 1022 g0 = a * g0 + g1;1023 h0 = a * h0 + h1;1024 1025 g1 = gx;1026 h1 = hx;1027 1028 if (h0 > maxdnom)1029 return false;1030 1031 delta0 = val  g0/h0;1032 delta1 = (delta0 < 0.0 ? val  (g01.0)/h0 : val  (g0+1.0)/h0);1033 }1034 1035 if (fabs(g0) > (LONG_MAX >> 4)  h0 > (LONG_MAX >> 4))1036 return false;1037 1038 if (h0 > 0.5)1039 return false;1040 1041 if (delta0 < maxdelta) {1042 if (fabs(delta1) > maxdelta)1043 return false;1044 numerator = (long)(g0  1.0);1045 denominator = (long)h0;1046 }1047 else if (delta0 > maxdelta) {1048 if (fabs(delta1) > maxdelta)1049 return false;1050 numerator = (long)(g0 + 1.0);1051 denominator = (long)h0;1052 }1053 else{1054 numerator = (long)g0;1055 denominator = (long)h0;1056 }1057 if ((denominator < 1) 1058 (fabs(val  (double)(numerator)/(double)(denominator)) > maxdelta))1059 return false;1060 return true;1061 } /* nearestRational */1062 934 1063 935 /************************************************************************/
Note: See TracChangeset
for help on using the changeset viewer.