Changeset 1236


Ignore:
Timestamp:
Jan 16, 2015 10:36:37 AM (5 years ago)
Author:
tkr
Message:

Switching to use of CoinRational?

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cgl/src/CglGMI/CglGMI.cpp

    r1122 r1236  
    2727#include "CglGMI.hpp"
    2828#include "CoinFinite.hpp"
     29#include "CoinRational.hpp"
    2930
    3031//-------------------------------------------------------------------
     
    890891  long numerator = 0, denominator = 0;
    891892  // 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();
    895897  }
    896898  else{
     
    904906      continue;
    905907    }
    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()));
    909912    }
    910913    else{
     
    929932  return true;
    930933} /* scaleCutIntegral */
    931 
    932 /************************************************************************/
    933 /* arguments:
    934  * val = double precision value that must be converted
    935  * maxdelta = max allowed difference between val and the rational computed
    936  * maxdnom = max allowed denominator
    937  * numerator = the numerator will be stored here if successful
    938  * denominator = the denominator will be stored here if successful
    939  * 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 GPL
    944  * from ZIB. More information can be obtained from the authors.
    945  *
    946  * Copyright (C) 2012 Konrad-Zuse-Zentrum                           
    947  *                    fuer Informationstechnik Berlin
    948  */
    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 table
    968    * multiplied by powers of 10 is tried as denominator
    969    */
    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 rational
    998    * 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 - (g0-1.0)/h0 : val - (g0+1.0)/h0);
    1009  
    1010   while ((fabs(delta0) > maxdelta) && (fabs(delta1) > maxdelta)) {
    1011     if ((b-a) < 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 - (g0-1.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 */
    1062934
    1063935/************************************************************************/
Note: See TracChangeset for help on using the changeset viewer.