🐛 sign(0) was returning 1 instead of 0 with use_dll

The C99 copysign() function was used in the generated C output, but that
function does not correctly handle zero. Replace it by a custom sign()
function.
master
Sébastien Villemot 2023-04-07 15:45:44 +02:00
parent 5a4297088d
commit b1e4884237
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
5 changed files with 36 additions and 34 deletions

View File

@ -890,7 +890,7 @@ DataTree::minLagForSymbol(int symb_id) const
}
void
DataTree::writePowerDeriv(ostream &output) const
DataTree::writeCHelpersDefinition(ostream &output) const
{
if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
output << "/*" << endl
@ -909,13 +909,21 @@ DataTree::writePowerDeriv(ostream &output) const
<< " return dxp;" << endl
<< " }" << endl
<< "}" << endl;
if (isUnaryOpUsed(UnaryOpcode::sign))
output << "double sign(double x)" << endl
<< "{" << endl
<< " return (x > 0) ? 1 : ((x < 0) ? -1 : 0);" << endl
<< "}" << endl;
}
void
DataTree::writePowerDerivHeader(ostream &output) const
DataTree::writeCHelpersDeclaration(ostream &output) const
{
if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
output << "double getPowerDeriv(double x, double p, int k);" << endl;
if (isUnaryOpUsed(UnaryOpcode::sign))
output << "double sign(double x);" << endl;
}
vector<string>

View File

@ -277,10 +277,10 @@ public:
//! Returns the minimum lag (as a negative number) of the given symbol in the whole data tree (and not only in the equations !!)
/*! Returns 0 if the symbol is not used */
int minLagForSymbol(int symb_id) const;
//! Write getPowerDeriv in C (function body)
void writePowerDeriv(ostream &output) const;
//! Write getPowerDeriv in C (prototype)
void writePowerDerivHeader(ostream &output) const;
//! Writes definitions of C function helpers (getPowerDeriv(), sign())
void writeCHelpersDefinition(ostream &output) const;
//! Writes declarations of C function helpers (getPowerDeriv(), sign())
void writeCHelpersDeclaration(ostream &output) const;
//! Thrown when trying to access an unknown variable by deriv_id
class UnknownDerivIDException
{

View File

@ -2950,10 +2950,10 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "abs";
break;
case UnaryOpcode::sign:
if (isCOutput(output_type))
output << "copysign";
else
output << "sign";
/* C does not have a sign() function, and copysign() is not suitable
because it does not handle zero correctly, so we define our own sign()
helper function, see DataTree::writeCHelpersDefinition() */
output << "sign";
break;
case UnaryOpcode::steadyState:
ExprNodeOutputType new_output_type;
@ -3054,8 +3054,6 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
&& arg->precedence(output_type, temporary_terms) < precedence(output_type, temporary_terms)))
{
output << LEFT_PAR(output_type);
if (op_code == UnaryOpcode::sign && isCOutput(output_type))
output << "1.0,";
close_parenthesis = true;
}

View File

@ -973,7 +973,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext,
output << "#include <math.h>" << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< endl;
writePowerDerivHeader(output);
writeCHelpersDeclaration(output);
output << endl
<< prototype_tt << endl
<< "{" << endl
@ -1012,7 +1012,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext,
output << "#include <math.h>" << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< endl;
writePowerDerivHeader(output);
writeCHelpersDeclaration(output);
output << endl
<< prototype_main << endl
<< "{" << endl
@ -1043,8 +1043,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext,
output << "#include " << it.filename() << endl;
output << endl;
// Write function definition if BinaryOpcode::powerDeriv is used
writePowerDeriv(output);
writeCHelpersDefinition(output);
output << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
@ -2645,18 +2644,13 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
NB: The prefix (static/dynamic) is added to the filename (even though its
the same source between static and dynamic) to avoid a race condition when
static and dynamic are compiled in parallel. */
open_file(model_src_dir / (prefix + "power_deriv.h"));
writePowerDerivHeader(output);
output.close();
filesystem::path power_deriv_src {model_src_dir / (prefix + "power_deriv.c")};
open_file(power_deriv_src);
filesystem::path helpers_src {model_src_dir / (prefix + "helpers.c")};
open_file(helpers_src);
output << "#include <math.h>" << endl << endl;
writePowerDeriv(output);
writeCHelpersDefinition(output);
output.close();
auto power_deriv_object {compileMEX(model_src_dir, (prefix + "power_deriv"),
mexext, { power_deriv_src },
matlabroot, false)};
auto helpers_object {compileMEX(model_src_dir, (prefix + "helpers"),
mexext, { helpers_src }, matlabroot, false)};
size_t ttlen {0};
@ -2724,8 +2718,9 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
open_file(source_tt);
output << "#include <math.h>" << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< R"(#include ")" << prefix << R"(power_deriv.h")" << endl
<< endl
<< endl;
writeCHelpersDeclaration(output);
output << endl
<< prototype_tt << endl
<< "{" << endl
<< tt_sparse_output[i].str()
@ -2747,7 +2742,7 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
output << "#include <math.h>" << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< endl;
writePowerDerivHeader(output);
writeCHelpersDeclaration(output);
output << endl
<< prototype_main << endl
<< "{" << endl
@ -2836,7 +2831,7 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
<< "}" << endl;
output.close();
vector<filesystem::path> mex_input_files { power_deriv_object, main_object_file, source_mex };
vector<filesystem::path> mex_input_files { helpers_object, main_object_file, source_mex };
for (int j {0}; j <= i; j++)
mex_input_files.push_back(tt_object_files[j]);
compileMEX(mex_dir, funcname, mexext, mex_input_files, matlabroot);
@ -2864,8 +2859,9 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
open_file(source_mex);
output << "#include <math.h>" << endl
<< R"(#include "mex.h")" << endl
<< R"(#include "../)" << prefix << R"(power_deriv.h")" << endl
<< endl
<< endl;
writeCHelpersDeclaration(output);
output << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " if (nrhs != " << nargin << ")" << endl
@ -2928,7 +2924,7 @@ ModelTree::writeSparseModelCFiles(const string &basename, const string &mexext,
}
output << "}" << endl;
output.close();
compileMEX(block_dir, funcname, mexext, { source_mex, power_deriv_object }, matlabroot);
compileMEX(block_dir, funcname, mexext, { source_mex, helpers_object }, matlabroot);
}
}
}

View File

@ -920,7 +920,7 @@ StaticModel::writeRamseyMultipliersDerivativesCFile(const string &basename, cons
output << "#include <math.h>" << endl << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< endl;
writePowerDeriv(output);
writeCHelpersDefinition(output);
output << endl
<< "void ramsey_multipliers_static_g1(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict g1m_v)" << endl
<< "{" << endl;