实现单纯性算法时遇到的精度问题

代码主要是实现线性规划中的单纯性算法的，分别有三种情况：1.有解；2.无解；3.无限解
在一个样例输入中
1 1
26
4362
27
正确输出是167.769230769230770000
我的程序输出为167.000000000000000
检查了一下，没有发现强制类型转换，应该是没有截断的

#include <iostream>
#include <iomanip>
#include <vector>
#include <exception>
#include <cmath>

using std::cin;
using std::cout;
using std::endl;
using std::vector;
using std::exception;


bool check_if_better(const vector<vector<double>>& restriction_matrix,const vector<double>& variable_values, const vector<double>& pleasure_num, double& prev_value)
{
    //this function checks both if it's better and is valid;
    //It also updates the prev_value into "pres_value";
    
    for(int i=0;i<restriction_matrix.size();i++)
    {
        double result = 0;
        for(int j=0;j<restriction_matrix[i].size()-1;j++)
        {
            result += restriction_matrix[i].at(j) * variable_values[j];
        }
        
        if(result > restriction_matrix[i].at(restriction_matrix[i].size() - 1)) // check whether it satisfies the restrictions
        {
            return false;
        }
    }
    
    double pres_value = 0;
    for(int i=0;i<pleasure_num.size();i++)
    {
        pres_value += pleasure_num[i] * variable_values[i];
    }
    
    if(pres_value < prev_value)
    {
        return false;
    }
    else
    {
        prev_value = pres_value;
        return true;
    }
}

void swap_equations(vector< vector<double>>& matrix, int equation_1, int equation_2)
{
    if(equation_1 == equation_2)
    {
        return;
    }
    std::swap(matrix[equation_1], matrix[equation_2]);
}

void scale(vector< vector<double>>& matrix, int equation_index, double scale_num)
{
    if(scale_num == 0 || scale_num == 1)
    {
        return;
    }
    if(scale_num < 0)
    {
        scale_num = scale_num * -1;
    }
    for(int i=0;i<matrix[equation_index].size();i++)
    {
        matrix[equation_index].at(i) /= scale_num;
    }
}

void subtract(vector< vector<double>>& matrix, int equation_1, int equation_2, double muti)
{
    if(muti == 0)
    {
        return;
    }
    for(int i=0;i<matrix[equation_1].size();i++)
    {
        matrix[equation_1].at(i) -= matrix[equation_2].at(i) * muti;
    }
}

void Gaussian_Elimination(vector<vector<double>> matrix)
{
    int pres_variable = 0;
    for(int pres_equation = 0;pres_equation<matrix.size();pres_equation++)
    {
        double pres_number = matrix[pres_equation].at(pres_variable);
        
        int change_index = pres_equation;
        while(pres_number == 0 && change_index < matrix.size())
        {
            change_index++;
            if(change_index < matrix.size())
            {
                pres_number = matrix[pres_equation].at(change_index);
            }
        }
        
        if(change_index >= matrix.size())
        {
            continue;
        }
        
        swap_equations(matrix, change_index, pres_equation);
        scale(matrix, pres_equation, matrix[pres_equation].at(pres_variable));
        
        for(int i=0;i<matrix.size();i++)
        {
            double muti = matrix[i].at(pres_variable); // matrix[pres_equation].at(pres_variable);
            
            if(i == pres_equation)
            {
                continue;
            }
            
            subtract(matrix, i, pres_equation, muti);
        }
        
        pres_variable++;
    }
    
    for(int i=0;i<matrix.size();i++)
    {
        int num_of_zero = 0;
        for(int j=0;j<matrix[i].size()-1;j++)
        {
            if(matrix[i].at(j) == 0)
            {
                num_of_zero++;
            }
        }
        if(num_of_zero == matrix[i].size()-1 && matrix[i].at(matrix[i].size()-1) != 0)
        {
            cout<<"No solution"<<endl;
            exit(0);
        }
    }

    int variable_num = matrix[0].size()-1;
    //int rank = 0;
    
    for(int i=0;i<variable_num;i++)
    {
        int num_of_zero = 0;
        for(int j=0;j<matrix.size();j++)
        {
            if(matrix[j].at(i) == 0)
            {
                num_of_zero++;
            }
        }
        //cout<<num_of_zero<<endl;
        if(num_of_zero == matrix.size())
        {
            cout<<"Infinity"<<endl;
            exit(0);
        }
    }
}

int main()
{
    int restriction_num,variable_num;
    cin>>restriction_num>>variable_num;
    vector< vector<double>> restriction_matrix(restriction_num, vector<double>(variable_num + 1));
    
    //store the restrictions;
    for(int i=0;i<restriction_num;i++)
    {
        for(int j=0;j<variable_num;j++)
        {
            cin>>restriction_matrix[i].at(j);
        }
    }
    
    //store result;
    for(int i=0;i<restriction_num;i++)
    {
        cin>>restriction_matrix[i].at(variable_num);
    }
    
    //store pleasure num;
    vector<double> pleasure_num(variable_num);
    for(int i=0;i<variable_num;i++)
    {
        cin>>pleasure_num[i];
    }
    
    //this step checks whether the solution is infinite;
    if(restriction_num < variable_num)
    {
        cout<<"Infinity"<<endl;
        return 0;
    }
    
    //this step checks whether there is a solution
    Gaussian_Elimination(restriction_matrix);
    
    vector<double> variable_value(variable_num);
    double pres_result = 0;

    for(int pres_variable=0;pres_variable<variable_num;pres_variable++)
    {
        while(true)
        {
            variable_value[pres_variable]++;
            
            if(!check_if_better(restriction_matrix, variable_value, pleasure_num, pres_result))
            {
                variable_value[pres_variable]--; //this means that the new value is invalid or isn't better than the old value, we cancel the adding step;
                break;
            }
        }
    }
    
    //define another matrix to store the changed values and for Gaussian_Elimination
    
    cout<<std::fixed<<std::showpoint<<std::setprecision(15);
    cout<<"Bounded solution"<<endl;
    for(int i=0;i<variable_value.size();i++)
    {
        cout<<variable_value[i]<<' ';
    }
    cout<<endl;
    
    return 0;
}

TechWhizKid参考GPT回答：

#include <iostream>
#include <iomanip>
#include <vector>
#include <exception>
#include <cmath>
 
using std::cin;
using std::cout;
using std::endl;
using std::vector;
using std::exception;
 
 
bool check_if_better(const vector<vector<double>>& restriction_matrix,const vector<double>& variable_values, const vector<double>& pleasure_num, double& prev_value)
{
    //this function checks both if it's better and is valid;
    //It also updates the prev_value into "pres_value";
    
    for(int i=0;i<restriction_matrix.size();i++)
    {
        double result = 0;
        for(int j=0;j<restriction_matrix[i].size()-1;j++)
        {
            result += restriction_matrix[i].at(j) * variable_values[j];
        }
        
        if(result > restriction_matrix[i].at(restriction_matrix[i].size() - 1)) // check whether it satisfies the restrictions
        {
            return false;
        }
    }
    
    double pres_value = 0;
    for(int i=0;i<pleasure_num.size();i++)
    {
        pres_value += pleasure_num[i] * variable_values[i];
    }
    
    if(pres_value < prev_value)
    {
        return false;
    }
    else
    {
        prev_value = pres_value;
        return true;
    }
}
 
void swap_equations(vector< vector<double>>& matrix, int equation_1, int equation_2)
{
    if(equation_1 == equation_2)
    {
        return;
    }
    std::swap(matrix[equation_1], matrix[equation_2]);
}
 
void scale(vector< vector<double>>& matrix, int equation_index, double scale_num)
{
    if(scale_num == 0 || scale_num == 1)
    {
        return;
    }
    if(scale_num < 0)
    {
        scale_num = scale_num * -1;
    }
    for(int i=0;i<matrix[equation_index].size();i++)
    {
        matrix[equation_index].at(i) /= scale_num;
    }
}
 
void subtract(vector< vector<double>>& matrix, int equation_1, int equation_2, double muti)
{
    if(muti == 0)
    {
        return;
    }
    for(int i=0;i<matrix[equation_1].size();i++)
    {
        matrix[equation_1].at(i) -= matrix[equation_2].at(i) * muti;
    }
}
 
void Gaussian_Elimination(vector<vector<double>> matrix)
{
    int pres_variable = 0;
    for(int pres_equation = 0;pres_equation<matrix.size();pres_equation++)
    {
        double pres_number = matrix[pres_equation].at(pres_variable);
        
        int change_index = pres_equation;
        while(pres_number == 0 && change_index < matrix.size())
        {
            change_index++;
            if(change_index < matrix.size())
            {
                pres_number = matrix[pres_equation].at(change_index);
            }
        }
        
        if(change_index >= matrix.size())
        {
            continue;
        }
        
        swap_equations(matrix, change_index, pres_equation);
        scale(matrix, pres_equation, matrix[pres_equation].at(pres_variable));
        
        for(int i=0;i<matrix.size();i++)
        {
            double muti = matrix[i].at(pres_variable); // matrix[pres_equation].at(pres_variable);
            
            if(i == pres_equation)
            {
                continue;
            }
            
            subtract(matrix, i, pres_equation, muti);
        }
        
        pres_variable++;
    }
    
    for(int i=0;i<matrix.size();i++)
    {
        int num_of_zero = 0;
        for(int j=0;j<matrix[i].size()-1;j++)
        {
            if(matrix[i].at(j) == 0)
            {
                num_of_zero++;
            }
        }
        if(num_of_zero == matrix[i].size()-1 && matrix[i].at(matrix[i].size()-1) != 0)
        {
            cout<<"No solution"<<endl;
            exit(0);
        }
    }
 
    int variable_num = matrix[0].size()-1;
    //int rank = 0;
    
    for(int i=0;i<variable_num;i++)
    {
        int num_of_zero = 0;
        for(int j=0;j<matrix.size();j++)
        {
            if(matrix[j].at(i) == 0)
            {
                num_of_zero++;
            }
        }
        //cout<<num_of_zero<<endl;
        if(num_of_zero == matrix.size())
        {
            cout<<"Infinity"<<endl;
            exit(0);
        }
    }
}
 
int main()
{
    int restriction_num,variable_num;
    cin>>restriction_num>>variable_num;
    vector< vector<double>> restriction_matrix(restriction_num, vector<double>(variable_num + 1));
    
    //store the restrictions;
    for(int i=0;i<restriction_num;i++)
    {
        for(int j=0;j<variable_num;j++)
        {
            cin>>restriction_matrix[i].at(j);
        }
    }
    
    //store result;
    for(int i=0;i<restriction_num;i++)
    {
        cin>>restriction_matrix[i].at(variable_num);
    }
    
    //store pleasure num;
    vector<double> pleasure_num(variable_num);
    for(int i=0;i<variable_num;i++)
    {
        cin>>pleasure_num[i];
    }
    
    //this step checks whether the solution is infinite;
    if(restriction_num < variable_num)
    {
        cout<<"Infinity"<<endl;
        return 0;
    }
    
    //this step checks whether there is a solution
    Gaussian_Elimination(restriction_matrix);
    
    vector<double> variable_value(variable_num);
    double pres_result = 0;
 
for(int pres_variable=0;pres_variable<variable_num;pres_variable++)
{
    double step = 1.0;
    while(step > 1e-9)
    {
        while(true)
        {
            variable_value[pres_variable] += step;
            
            if(!check_if_better(restriction_matrix, variable_value, pleasure_num, pres_result))
            {
                variable_value[pres_variable] -= step; //this means that the new value is invalid or isn't better than the old value, we cancel the adding step;
                break;
            }
        }
        step /= 10.0;
    }
}

    
    //define another matrix to store the changed values and for Gaussian_Elimination
    
    cout<<std::fixed<<std::showpoint<<std::setprecision(15);
    cout<<"Bounded solution"<<endl;
    for(int i=0;i<variable_value.size();i++)
    {
        cout<<variable_value[i]<<' ';
    }
    cout<<endl;
    
    return 0;
}
 

基于new bing部分指引作答：
从您提供的代码中，我发现可能导致精度问题的一处地方。在check_if_better函数中，您使用double类型的变量pres_value来计算目标函数的值，然后与prev_value比较。然而，在计算机中，浮点数的表示是有限的，因此可能存在舍入误差。这可能会导致在比较pres_value和prev_value时出现微小的差异。

为了解决这个问题，您可以使用一个小的误差范围来比较pres_value和prev_value，而不是直接进行相等性比较。例如，您可以检查它们之间的绝对差异是否小于某个很小的阈值。这将允许您在接受非常接近的值时避免精度问题。

以下是修改后的check_if_better函数的示例代码：

bool check_if_better(const vector<vector<double>>& restriction_matrix, const vector<double>& variable_values, const vector<double>& pleasure_num, double& prev_value)
{
    for (int i = 0; i < restriction_matrix.size(); i++)
    {
        double result = 0;
        for (int j = 0; j < restriction_matrix[i].size() - 1; j++)
        {
            result += restriction_matrix[i].at(j) * variable_values[j];
        }

        if (result > restriction_matrix[i].at(restriction_matrix[i].size() - 1))
        {
            return false;
        }
    }

    double pres_value = 0;
    for (int i = 0; i < pleasure_num.size(); i++)
    {
        pres_value += pleasure_num[i] * variable_values[i];
    }

    // 比较 pres_value 和 prev_value 的绝对差异是否小于一个阈值
    const double epsilon = 1e-9;
    if (pres_value - prev_value < epsilon)
    {
        return false;
    }
    else if (pres_value > prev_value)
    {
        prev_value = pres_value;
        return true;
    }
    else
    {
        return false;
    }
}

通过引入一个小的误差范围，您可以更好地处理由于浮点数精度问题而引起的输出差异。希望这可以解决您遇到的问题！

根据你的提供的代码，我认为缺少了一个重要的步骤，即 将输入的值存储为double型，因此存在一定的截断误差。检查输出精度时，您可以从4位小数扩大到10位小数等，但实际上它只会使输出精度更高，而不会根本解决截断误差的问题。
建议读取输入时直接使用double类型，以避免将整数值隐式地转换为double类型。例如，您可以将存储值的vector的类型从int改为double。
由于C ++中浮点数具有有限的精度，因此在大于或等于14位小数的情况下，double类型具有保真度。因此，您应该为cout语句设置至少14位小数，而不是增加位数来解决当前精度问题。

来自GPT的回答：

在你的代码中，问题可能出在使用std::setprecision(15)来设置输出精度时。尽管你将精度设置为15位小数，但由于浮点数的内部表示方式，可能会导致结果出现不准确的情况。

为了解决这个问题，你可以考虑使用字符串来存储结果，并在输出时以字符串的形式进行格式化。这样可以确保输出结果的精度和准确性。以下是对你的代码进行修改的示例：

cpp
Copy code
#include <sstream> // 添加头文件

// ...

std::stringstream ss; // 添加 stringstream 对象

// ...

// 修改输出部分的代码
ss << "Bounded solution" << endl;
for (int i = 0; i < variable_value.size(); i++) {
    ss << std::fixed << std::setprecision(15) << variable_value[i] << ' ';
}
ss << endl;

cout << ss.str(); // 使用字符串输出结果

// ...

通过使用字符串流（stringstream），你可以将结果存储为字符串，并在需要时以字符串形式输出，从而确保精度和准确性。

另外，请确保在编译代码时，使用合适的编译选项来启用浮点数的高精度计算，例如使用 -std=c++11 或更高的标准。

请注意，浮点数计算存在舍入误差和精度限制，因此在处理浮点数时，始终要注意这些问题。

不知道你这个问题是否已经解决, 如果还没有解决的话:

• 这有个类似的问题, 你可以参考下: https://ask.csdn.net/questions/7774173
• 这篇博客也不错, 你可以看下代码随想录算法训练营第5天 ：242.有效的字母异位词，349.两个数组的交集，1.两数之和
• 除此之外, 这篇博客: 数据结构经典算法中的 2.连续自然数1~1000，其中少了一个数字，怎样用常数时间复杂度找出那个数字丢失。 部分也许能够解决你的问题, 你可以仔细阅读以下内容或者直接跳转源博客中阅读:
  
  • 可以直接对所有待查询数字求和sum，然后用1~1000的的数字之和total减去之前求和sum得到的结果即为丢失的数字；

  /**
  	 * 1.递归实现1-1000的累加 从最大的开始往下加
  	 */
  	public static int sum(int num) {
  		if (num == 1) {
  			return 1;
  		} else
  			return num + sum(num - 1);
  	}
  
  	/**
  	 * 2.for循环实现1-1000的累加 从最小的开始往上加
  	 */
  	public static int sum2(int num) {
  		int sum = 0;
  		for (int i = 1; i <= num; i++) {
  			sum += i;
  		}
  		return sum;
  	}
  
  	/**
  	 * 3.利用高斯定理实现1-1000的累加 不往上不往下，首尾成对计算和
  	 */
  	public static int sum3(int num) {
  		/*
  		 * 判断成对的依据By对2求余
  		 */
  		int flag = num % 2;
  		// 如果成对的话，直接利用首尾相加*对数
  		if (flag == 0) {
  			// 如果成对的话
  			int group = num / 2;
  			return group * (1 + num);
  		}
  		// 否则的话 凑够偶数对 减去 那个凑个的数
  		return sum3(num + 1) - (num + 1);
  	}

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^

你的代码使用了高斯消元法来解线性规划问题，并通过遍历所有可能的变量值来找到最优解。你的问题可能出现在浮点数的精度上。你可以尝试在输出结果之前使用std::fixed、std::showpoint和std::setprecision来设置输出结果的小数精度，例如：cout<<std::fixed<<std::showpoint<<std::setprecision(15);。这样可以确保输出结果的精度。