From db0dbfdaa43454beb105bf1517cd25b5b678b63b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien.villemot@ens.fr>
Date: Tue, 23 Feb 2010 14:39:49 +0100
Subject: [PATCH] Preprocessor: full support for normcdf() function, in
 evaluation and in C files (closes #84)

---
 ExprNode.cc | 34 ++++++++++++++++++++++------------
 1 file changed, 22 insertions(+), 12 deletions(-)

diff --git a/ExprNode.cc b/ExprNode.cc
index 76457a7b..a063df93 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -2945,8 +2945,7 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
   switch (op_code)
     {
     case oNormcdf:
-      cerr << "NORMCDF: eval not implemented" << endl;
-      exit(EXIT_FAILURE);
+      return (0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
     }
   // Suppress GCC warning
   exit(EXIT_FAILURE);
@@ -3008,9 +3007,6 @@ void
 TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                            const temporary_terms_type &temporary_terms) const
 {
-  // TrinaryOpNode not implemented for C output
-  assert(!IS_C(output_type));
-
   // If current node is a temporary term
   temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
   if (it != temporary_terms.end())
@@ -3022,15 +3018,29 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
   switch (op_code)
     {
     case oNormcdf:
-      output << "normcdf(";
+      if (IS_C(output_type))
+        {
+          // In C, there is no normcdf() primitive, so use erf()
+          output << "(0.5*(1+erf(((";
+          arg1->writeOutput(output, output_type, temporary_terms);
+          output << ")-(";
+          arg2->writeOutput(output, output_type, temporary_terms);
+          output << "))/(";
+          arg3->writeOutput(output, output_type, temporary_terms);
+          output << ")/M_SQRT2)))";
+        }
+      else
+        {
+          output << "normcdf(";
+          arg1->writeOutput(output, output_type, temporary_terms);
+          output << ",";
+          arg2->writeOutput(output, output_type, temporary_terms);
+          output << ",";
+          arg3->writeOutput(output, output_type, temporary_terms);
+          output << ")";
+        }
       break;
     }
-  arg1->writeOutput(output, output_type, temporary_terms);
-  output << ",";
-  arg2->writeOutput(output, output_type, temporary_terms);
-  output << ",";
-  arg3->writeOutput(output, output_type, temporary_terms);
-  output << ")";
 }
 
 void
-- 
GitLab