Add Bessel functions of second kind (yn)

Sync with AMReX-Codes/amrex#4234.

Author: WeiqunZhang

README.md

@@ -20,7 +20,8 @@ in the form of string.  It supports `+`, `-`, `*`, `/`, `**` (power), `^`
 can be computed with `min` and `max`, respectively.  It supports the
 Heaviside step function, `heaviside(x1,x2)` that gives `0`, `x2`, `1`, for
 `x1 < 0`, `x1 = 0` and `x1 > 0`, respectively.  It supports the Bessel
-function of the first kind of order `n` `jn(n,x)`. Complete elliptic
+function of the first kind of order `n` `jn(n,x)`, and the Bessel function
+of the second kind of order ``n`` ``yn(n,x)``. Complete elliptic
 integrals of the first and second kind, `comp_ellint_1(k)` and
 `comp_ellint_2(k)`, are supported.  There is `if(a,b,c)` that gives `b` or
 `c` depending on the value of `a`.  A number of comparison operators are

Src/amrexpr_Parser_Y.H

@@ -90,6 +90,7 @@ enum parser_f2_t {  // Built-in functions with two arguments
     PARSER_OR,
     PARSER_HEAVISIDE,
     PARSER_JN,
+    PARSER_YN,
     PARSER_MIN,
     PARSER_MAX,
     PARSER_FMOD
@@ -115,6 +116,7 @@ std::string_view parser_f2_s[] =
     "or",
     "heaviside",
     "jn",
+    "yn",
     "min",
     "max",
     "fmod"
@@ -399,6 +401,21 @@ T parser_math_jn (int a, T b)
 #endif
 }
 
+template <typename T>
+AMREXPR_GPU_HOST_DEVICE AMREXPR_NO_INLINE
+T parser_math_yn (int a, T b)
+{
+#if defined AMREXPR_USE_SYCL || defined __MINGW32__
+    amrex::ignore_unused(a,b);
+    // neither yn(f) nor std::cyl_bessel_y work yet
+    // https://github.com/oneapi-src/oneAPI-spec/issues/308
+    AMREXPR_ALWAYS_ASSERT_WITH_MESSAGE(false, "parser: yn in SYCL not supported yet");
+    return 0.0;
+#else
+    return yn(a, b);
+#endif
+}
+
 AMREXPR_GPU_HOST_DEVICE AMREXPR_FORCE_INLINE double
 parser_call_f1 (enum parser_f1_t type, double a)
 {
@@ -459,6 +476,8 @@ parser_call_f2 (enum parser_f2_t type, double a, double b)
         return (a < 0.0) ? 0.0 : ((a > 0.0) ? 1.0 : b);
     case PARSER_JN:
         return parser_math_jn<double>(int(a),b);
+    case PARSER_YN:
+        return parser_math_yn<double>(int(a),b);
     case PARSER_MIN:
         return (a < b) ? a : b;
     case PARSER_MAX:

Src/amrexpr_parser.l

@@ -70,6 +70,7 @@ EXP	([Ee][-+]?[0-9]+)
 "pow"         { amrexpr_parserlval.f2 = amrexpr::PARSER_POW;       return F2; }
 "heaviside"   { amrexpr_parserlval.f2 = amrexpr::PARSER_HEAVISIDE; return F2; }
 "jn"          { amrexpr_parserlval.f2 = amrexpr::PARSER_JN;        return F2; }
+"yn"          { amrexpr_parserlval.f2 = amrexpr::PARSER_YN;        return F2; }
 "min"         { amrexpr_parserlval.f2 = amrexpr::PARSER_MIN;       return F2; }
 "max"         { amrexpr_parserlval.f2 = amrexpr::PARSER_MAX;       return F2; }
 "fmod"        { amrexpr_parserlval.f2 = amrexpr::PARSER_FMOD;      return F2; }