Make XcFun compatible with the Microsoft visual C++ compiler.

1) Change and and or to && and ||
2) Function pointer casts obey different rules
3) Microsoft do not know math
4) Auto arrays not supported??
+ polish. This compiles under windows and the tests pass..

include/xcfun.h

@@ -90,6 +90,14 @@ extern "C" {
   // Get the current value of the setting.
   double xc_get_param(xc_functional fun, int param);
 
+  /* Transform the output of xc_eval to a different mode,
+     for example because your program wants AB mode but you
+     want to take advantage of knowing that you only care
+     about non-polarizing N mode derivatives. */
+  void xc_transform(int order, 
+		    int from_mode, const double *from_data,
+		    int to_mode, double *to_data);
+
 #ifdef __cplusplus
 } // End of extern "C"
 #endif

src/evaluators.cpp

@@ -544,8 +544,8 @@ evaluator xc_evaluator_lookup(int mode, int type, int order)
       eval_setup_mgga<XC_MGGA_MAX_ORDER>();
       eval_setup_mlgga<XC_MLGGA_MAX_ORDER>();
     }
-  assert(mode>=0 and mode < XC_NR_MODES);
-  assert(type>=0 and type < XC_NR_TYPES);
-  assert(order>=0 and order <= XC_MAX_ORDER);
+  assert(mode>=0 && mode < XC_NR_MODES);
+  assert(type>=0 && type < XC_NR_TYPES);
+  assert(order>=0 && order <= XC_MAX_ORDER);
   return eval_tab->tab[mode][type][order];
 }

src/fortran.c

@@ -51,9 +51,10 @@ void FSYM(xceval)(int *fun,  int *order,
 		  double *first_result,
 		  double *second_result)
 {
+  int dpitch, rpitch;
   assert(*fun >= 0 && *fun < MAX_FORTRAN_FUNCTIONALS);
-  int dpitch = second_density - first_density;
-  int rpitch = second_result - first_result;
+  dpitch = second_density - first_density;
+  rpitch = second_result - first_result;
   xc_eval_vec(fortran_functionals[*fun], *order, *nr_points, 
 	      first_density, dpitch,
 	      first_result, rpitch);

src/functional.cpp

@@ -26,7 +26,7 @@ void functional::describe(enum xc_parameters weight_param, int type,
 			  const char *oneliner,
 			  const char *reference)
 {
-  assert(type >= 0 and type < XC_NR_TYPES);
+  assert(type >= 0 && type < XC_NR_TYPES);
   m_type = type;
   m_name = weight_param;
   xcint_set_short_description(weight_param, oneliner);

src/functional.h

@@ -3,7 +3,9 @@
 
 // Everything required to implement a functional
 #include "taylor.h"
+#ifdef XCFUN_CONTRACTIONS
 #include "ctaylor.h"
+#endif
 #define XCFUN_INTERNAL
 #include "xcfun.h"
 #include "config.h"
@@ -68,7 +70,12 @@ class functional
   {
     assert(Nvar<=XC_MAX_NVAR);
     assert(Ndeg<=XC_MAX_ORDER);
+    // See comment at eval about this cast
+#ifdef _MSC_VER
+    ftab[Nvar][Ndeg] = *reinterpret_cast<void **>(&f);
+#else
     ftab[Nvar][Ndeg] = reinterpret_cast<void *>(f);
+#endif
   }
 #ifdef XCFUN__CONTRACTIONS
   template<int Ndeg>
@@ -86,10 +93,22 @@ class functional
     assert(Nvar<=XC_MAX_NVAR);
     assert(Ndeg<=XC_MAX_ORDER);
     assert(ftab[Nvar][Ndeg]);
+#ifdef _MSC_VER
+    /* Here we need to cast the void pointer back to a function
+       pointer of the appropriate type. This is supposed to work
+       under C++0x, may not be supported under all compilers.
+       Using solution from 
+       http://stackoverflow.com/questions/1096341/function-pointers-casting-in-c
+     */
+    taylor<double,Nvar,Ndeg> 
+      (*f)(const densvars< taylor<double,Nvar,Ndeg> > &);
+    *reinterpret_cast<void**>(&f) = ftab[Nvar][Ndeg];
+#else
     taylor<double,Nvar,Ndeg> 
       (*f)(const densvars< taylor<double,Nvar,Ndeg> > &) = 
       reinterpret_cast<taylor<double,Nvar,Ndeg> 
       (*)(const densvars< taylor<double,Nvar,Ndeg> > &)>(ftab[Nvar][Ndeg]);
+#endif
     return f(dv);
   }
 #ifdef XCFUN_CONTRACTIONS

src/parameters.cpp

@@ -31,30 +31,30 @@ static functional *param_functionals[XC_NR_PARAMS] = {0};
 extern "C"
 const char *xc_name(int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);    
   return param_symbols[param];
 }
 extern "C"
 const char *xc_short_description(int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   xcint_assure_setup();
   return param_short[param];
 }
 extern "C"
 const char *xc_long_description(int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   xcint_assure_setup();
   return param_long[param];
 }
 
 double xcint_default(int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   xcint_assure_setup();
   return param_default[param];
@@ -64,36 +64,36 @@ double xcint_default(int param)
 // Short description is a string without newlines!
 void xcint_set_short_description(int param, const char *text)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
     xcint_die("Invalid parameter nr (version mismatch?)",param);
   param_short[param] = text;
 }
 // Long description should end with a final newline, and may have many lines.
 void xcint_set_long_description(int param, const char *text)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   param_long[param] = text;
 }
 
 void xcint_set_default(int param, double value)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   param_default[param] = value;
 }
 
 // Short description is a string without newlines!
 void xcint_set_functional(int param, functional *f)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   param_functionals[param] = f;
 }
 // Long description should end with a final newline, and may have many lines.
 functional *xcint_functional(int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter nr (version mismatch?)",param);
   return param_functionals[param];
 }

src/specmath.h

@@ -100,6 +100,7 @@ static taylor<T,Nvar,Ndeg> ufunc(const taylor<T,Nvar,Ndeg> &x, const double &a)
   return res;
 }
 
+#ifdef XCFUN_CONTRACTIONS
 template<class T, int Ndeg>
 static ctaylor<T,Ndeg> ufunc(const ctaylor<T,Ndeg> &x, const double &a)
 {
@@ -113,5 +114,6 @@ static ctaylor<T,Ndeg> ufunc(const ctaylor<T,Ndeg> &x, const double &a)
   x.compose(res,tmp1.c);
   return res;
 }
+#endif
 
 #endif

src/taylor/micromath.h

@@ -0,0 +1,91 @@
+/*
+  Math functions for the truly unbelivably bad Visual C++,
+  which does not define many common math functions but instead
+  give a table in the documentation for how you can implement
+  them yourself. -ulf
+ */
+
+inline double cbrt(double x)
+{
+  return pow(x,1.0/3.0);
+}
+
+inline double asinh(double x)
+{
+  return log(x+ sqrt(x*x+1));
+}
+
+/* VC++ does not have erf() (sigh), so here is a TEMPORARY
+   solution with a very unclear copyright. From
+   http://www.digitalmars.com/archives/cplusplus/3634.html
+   Ulf
+*/
+
+/***************************
+*   erf.cpp
+*   author:  Steve Strand
+*   written: 29-Jan-04
+***************************/
+
+static const double erf_rel_error= 1E-12;        //calculate 12 significant figures
+//you can adjust rel_error to trade off between accuracy and speed
+//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
+
+
+inline double erfc(double x);
+inline double erf(double x)
+//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
+//       = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
+//       = 1-erfc(x)
+{
+    static const double two_sqrtpi=  1.128379167095512574;        // 2/sqrt(pi)
+    if (fabs(x) > 2.2) {
+        return 1.0 - erfc(x);        //use continued fraction when fabs(x) > 2.2
+    }
+    double sum= x, term= x, xsqr= x*x;
+    int j= 1;
+    do {
+        term*= xsqr/j;
+        sum-= term/(2*j+1);
+        ++j;
+        term*= xsqr/j;
+        sum+= term/(2*j+1);
+        ++j;
+    } while (fabs(term)/sum > erf_rel_error);
+    return two_sqrtpi*sum;
+}
+
+inline double erfc(double x)
+//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
+//        = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
+//        = 1-erf(x)
+//expression inside [] is a continued fraction so '+' means add to denominator
+//only
+{
+    static const double one_sqrtpi=  0.564189583547756287;        // 1/sqrt(pi)
+    if (fabs(x) < 2.2) {
+        return 1.0 - erf(x);        //use series when fabs(x) < 2.2
+    }
+    if (x<0) {               //continued fraction only valid for x>0
+        return 2.0 - erfc(-x);
+    }
+    double a=1, b=x;                //last two convergent numerators
+    double c=x, d=x*x+0.5;          //last two convergent denominators
+    double q1, q2= b/d;             //last two convergents (a/c and b/d)
+    double n= 1.0, t;
+    do {
+        t= a*n+b*x;
+        a= b;
+        b= t;
+        t= c*n+d*x;
+        c= d;
+        d= t;
+        n+= 0.5;
+        q1= q2;
+        q2= b/d;
+      } while (fabs(q1-q2)/q2 > erf_rel_error);
+    return one_sqrtpi*exp(-x*x)*q2;
+}
+
+
+

src/taylor/taylor.h

@@ -185,14 +185,14 @@ class taylor : public polymul::polynomial<T, Nvar, Ndeg>
     assert(Nres >= Ndeg);
     taylor<T, Nvar,Ndeg> tmp = *this;
     tmp[0] = 0;
-    polymul::taylorcompose0(res,tmp,coeff.c);
+    polymul::taylorcompose0<T,Nvar,Nres,Nres>(res,tmp,coeff.c);
   }
   // Like compose, but this->c[0] _must_ be equal to 0. Slightly faster.
   template<int Nres>
   void compose0(taylor<T, Nvar,Nres>& res, const taylor<T,1,Nres> &coeff) const
   {
     assert(Nres >= Ndeg);
-    polymul::taylorcompose0(res,*this,coeff.c);
+    polymul::taylorcompose0<T,Nvar,Nres,Nres>(res,*this,coeff.c);
   }
   template<int Nres>
   taylor<T, Nvar,Nres> compose(const taylor<T,1,Nres> &coeff) const

src/taylor/taylor_math.h

@@ -5,6 +5,10 @@
 #define M_PI 3.14159265358979323846
 #endif
 
+#ifdef _MSC_VER
+#include "micromath.h"
+#endif
+
 //Taylor series of 1/(a+x)
 template<class T,int N>
 static void inv_taylor(taylor<T,1,N>& t, const T &a)

src/xcfun.cpp

@@ -49,8 +49,8 @@ int xcint_input_length(int mode, int type)
 {
   static int tab[XC_NR_MODES][XC_NR_TYPES] = 
     {{1,2,3,4},{1,2,3,4},{2,5,7,9},{2,5,7,9}};
-  assert(mode>=0 and mode <= XC_NR_MODES);
-  assert(type>=0 and type <= XC_NR_TYPES);
+  assert(mode>=0 && mode <= XC_NR_MODES);
+  assert(type>=0 && type <= XC_NR_TYPES);
   return tab[mode][type];
 }
 
@@ -82,7 +82,7 @@ int xc_functional_data::get_type(void) const
 }
 void xc_functional_data::set_mode(int mode)
 {
-  if (!(mode>=0 and mode < XC_NR_MODES))
+  if (!(mode>=0 && mode < XC_NR_MODES))
     xcint_die("Invalid mode to xc_functional::set_mode()",mode);
   this->mode = mode;
   find_max_order();
@@ -91,7 +91,7 @@ void xc_functional_data::set_mode(int mode)
 void xc_functional_data::find_max_order(void)
 {
   max_order = -1;
-  while (max_order < XC_MAX_ORDER and xc_evaluator_lookup(mode,type,max_order+1))
+  while (max_order < XC_MAX_ORDER && xc_evaluator_lookup(mode,type,max_order+1))
     max_order++;
 }
 
@@ -208,10 +208,10 @@ void xc_set_mode(xc_functional fun, int mode)
 
 void xc_set_param(xc_functional fun, int param, double value)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter in xc_set_param()",param);
   fun->parameters[param] = value;
-  if (xc_is_functional(param) and value != 0)
+  if (xc_is_functional(param) && value != 0)
     {
       if (xcint_functional(param)->m_type > fun->type)
 	fun->type = xcint_functional(param)->m_type;
@@ -221,7 +221,7 @@ void xc_set_param(xc_functional fun, int param, double value)
 
 double xc_get_param(xc_functional fun, int param)
 {
-  if (param < 0 or param >= XC_NR_PARAMS)
+  if (param < 0 || param >= XC_NR_PARAMS)
       xcint_die("Invalid parameter in xc_get_param()",param);
   return fun->parameters[param];
 }
@@ -241,7 +241,7 @@ static int run_tests(functional *fun)
   xc_set_mode(xf,fun->test_mode);
   xc_set_param(xf,fun->m_name,1.0);
   int n = xc_output_length(xf, fun->test_order);
-  double out[n];
+  double *out = (double *)malloc(n*sizeof*out);
   double *reference = fun->test_output;
   xc_eval(xf,fun->test_order,fun->test_input,out);
 
@@ -271,6 +271,7 @@ static int run_tests(functional *fun)
       printf("%s ok\n",xc_name(fun->m_name));
     }
   xc_free_functional(xf);
+  free(out);
   return nerr;
 }
 
@@ -286,7 +287,7 @@ int xcfun_test(void)
 	{
 	  nr_run++;
 	  int res = run_tests(f);
-	  if (res != 0 and res != -1)
+	  if (res != 0 && res != -1)
 	    nr_failed++;
 	}
     }