ex-2/fancier: use optimised Lambert W

ex-2/fancier.c

@@ -40,6 +40,23 @@
  * the computation of Euler's constant, Jean-Pierre Demailly.
  */
 
+
+/* Approximation of the Lambert W(x) function at large values
+ * where x = 5π/9 10^(2D + 1). The expression of x has been
+ * inlined and optimized to avoid overflowing on large D.
+ *
+ * Note: from x>10 this is an overestimate of W(x), which is
+ * good for the application below.
+ *
+ * source: "On the Lambert W function", Corless et al.
+ */
+double lambertW_large(size_t digits) {
+  double L1 = log(5*M_PI/9) + (2*digits + 1)*log(10);
+  double L2 = log(L1);
+  return L1 - L2  + L2/L1 + (L2 - 2)/(2*L1*L1)
+       + L2*(3*L2*L2*L2 - 22*L2*L2 + 36*L2 - 12)/(12*L1*L1*L1*L1);
+}
+
 /* `mpq_dec_str(z, n)`
  * returns the decimal representation of the
  * rational z as a string, up to n digits.
@@ -525,7 +542,7 @@ int main(int argc, char** argv) {
   // requested decimal digits
   size_t D = atol(argv[1]);
   // Brent-McMillan number N
-  size_t N = gsl_sf_lambert_W0(5*M_PI/9*pow(10, 2*D + 1))/16 + 1;
+  size_t N = lambertW_large(D)/16 + 1;
 
   mpq_t exact, indicator;
   mpq_inits(exact, indicator, NULL);