Threaded Mode | Linear Mode

robve · 04-08-2021, 06:17 PM

(04-08-2021 04:27 PM)Albert Chan Wrote:
(04-08-2021 01:32 PM)robve Wrote: Using exp(t) in the inner loop instead of the strength-reduced expt recurrence should be the most accurate. Try it.

I did. It make no difference.

That's really nice. I was stating this more in "theory" that it would be most accurate, but practically makes no difference, which is really nice. On the SHARP Pocket Computers it makes a difference only because variables hold 10 digits whereas expressions are evaluated with 12 digit precision. So using EXP is a tiny bit better than repeated multiplication with a variable, e.g. we get the integral 2.0 exactly instead of 1.999999998.

(04-08-2021 04:27 PM)Albert Chan Wrote: Same conclusion when I compared strength-reduced exp vs supposed more accurate expm1.
I even tried "improved u", when t get closer to 1. Again, no difference.

u = exp(1/t - t) = exp((1-t)*(1+1/t))

Think of this as Newton's method.
Once starting to converge, every iteration doubled accuracy.
In other words, correction term only required half precision. (guess had the other half)

Same with quadrature, required precision keep reducing.
Inaccuracy of strength-reduced expt does not contribute to the sum.

I didn't think about it that way. That is a nice way to confirm.

I would also like to confirm that removing if (r == 0 || r == 1) break; makes no difference with all integrals tested. Logically the loop terminates before these conditions are met. This check is redundant and can be removed as you observed correctly.

One thing I am still curious about is reusing or setting fp and fm to zero. I will look into the decaying factor idea. Perhaps there is also a difference between Inf and NaN, because with Inf the limit is clearly not finite and with NaN we don't necessarily know, but reusing in that case makes sense, I think. But I may be wrong if it makes no difference. I just don't know yet.

- Rob

Albert Chan · 04-08-2021, 06:59 PM

Just to show expt make no difference, this version eliminated t = exp(h), replaced by cheaper sqrt().

Code:

function qthsh3(f, a, b, c, n, eps) -- tanh-sinh quadrature

    c = c or 1                      -- decay factor, 0 to 1

    n = n or 6                      -- 6 is optimal, 7 take longer

    eps = eps or 1e-9               -- relative error

    local s, d = f((b+a)/2), (b-a)/2

    local k, fp, fm, t0, err = 0, 0, 0, 7.38905609893065

    repeat

        local p, t = 0, sqrt(t0)

        local t1 = k==0 and t or t0

        t0 = t                      -- e, e^.5, e^.25, ...

        repeat

            local u = exp(1/t-t)

            local r = 2*u/(1+u)

            local x, y = d*r

            if a+x ~= a then

                y = f(a+x)

                if y-y == 0 then fp=y else fp=c*fp end

            end

            if b-x ~= b then

                y = f(b-x)

                if y-y == 0 then fm=y else fm=c*fm end

            end

            y = (t+1/t)*r/(1+u) * (fp+fm)   -- apply weight

            p, t = p+y, t*t1

        until abs(y) <= eps*abs(p)

        s, err, k = s+p, abs(s-p), k+1

    until err <= 10*eps*abs(s) or k>n

    return d*ldexp(s,1-k), err/abs(s)

end

lua> f = function(x) return 1/sqrt(-log(x)) end

lua> qthsh (f, 0, 1) -- exp version
1.7724538240988943 6.586211614189097e-009
lua> qthsh2(f, 0, 1) -- expm1 version
1.7724538240988927 6.58621067462469e-009
lua> qthsh3(f, 0, 1) -- sqrt version
1.7724538240988927 6.586210486711807e-009

Albert Chan · 04-08-2021, 10:43 PM

(04-08-2021 06:17 PM)robve Wrote: Perhaps there is also a difference between Inf and NaN, because with
Inf the limit is clearly not finite and with NaN we don't necessarily know.

Unless formula is very simple, there is not enough information to claim returning inf implied infinite true result.

Cas> f(u) := atanh(u*u*(3-2*u)) * 6*u*(1-u)
Cas> u = 1 - 10^-9
Cas> f(float(u)) → ∞
Cas> float(f(u)) → 1.23123199576e-07

robve · 04-09-2021, 02:54 AM

(04-08-2021 10:43 PM)Albert Chan Wrote:
(04-08-2021 06:17 PM)robve Wrote: Perhaps there is also a difference between Inf and NaN, because with
Inf the limit is clearly not finite and with NaN we don't necessarily know.

Unless formula is very simple, there is not enough information to claim returning inf implied infinite true result.

Cas> f(u) := atanh(u*u*(3-2*u)) * 6*u*(1-u)
Cas> u = 1 - 10^-9
Cas> f(float(u)) → ∞
Cas> float(f(u)) → 1.23123199576e-07

The function is a black box and we cannot trust ill-conditioned functions. With NaN and Inf we don't know what the cause is. Tracing the points toward an endpoint may be suggestive of the "missing" values by extrapolation (reusing the previous point). However, a small change in the bounds can make a huge difference in the result and cause Inf or NaN.

For example, 1/(1+cos(x)^x) integrated on [-pi/2,pi/2] is OK but produces NaN with the inexact bounds [-1.570796327,1.570796327] see the spreadsheet row 64. Earlier I mentioned that fp and fm reuse appears to be better than zero for this function and for two other functions. But there is more noise than signal to use these examples as sufficient evidence to make a determination. I can't tell yet if a decay approach can help to control singularities at an endpoint. It may be moot attempt. If so, there is not much left that remains to do for further improvements of Tanh-Sinh.

- Rob

Albert Chan · 04-10-2021, 12:05 AM

(04-08-2021 01:58 PM)Albert Chan Wrote:
(04-08-2021 11:20 AM)emece67 Wrote: integrals (1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x)), from -inf to +inf) needed one more level.
This particular case worsened dramatically, from 13 correct digits with strength reduction to only 4 with these changes)

I believe there is nothing to worry about.
Integrand suffered from subtraction cancellations.

Amazingly, qthsh (I picked sqrt version, for speed), does not get garbage from this.

lua> f = function(x) n=n+1; return 1-sqrt(2)*cosh(x)/sqrt(cosh(2*x)) end
lua> f01 = function(t) local y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y) end

With default setting of n=6, eps=1e-9, we have:

lua> Q = require 'quad'
lua> n=0; print(Q.qthsh(f01, 0, 1, 1), n) -- Previous Point
-0.6931471805996943 80
lua> n=0; print(Q.qthsh(f01, 0, 1, 0), n) -- Default Zero
-0.6931471805582439 77

Showing running integral sum, for each level k (c=0, n=5, eps=0)

lua> n=0; Q.qthsh(f01, 0, 1, 0, 5, 0)
0 -1.6568543106755638
1 -0.8783982505319032
2 -0.7037071116271125
3 -0.693194301031653
4 -0.6931471818102728
5 -0.6931471805599556
lua> n, -log(2)
99 -0.6931471805599453

Albert Chan · (This post was last modified: 04-10-2021 02:25 PM by Albert Chan.)

(04-08-2021 11:20 AM)emece67 Wrote: integrals (1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x)), from -inf to +inf) needed one more level.
This particular case worsened dramatically, from 13 correct digits with strength reduction to only 4 with these changes)

The reason is not the alogrithm, but how mpmath numbers work.
Because mpmath does not do overflow / underflow, the integrand may generate noise.
(with IEEE double, generated Inf/NaN will allow the algorithm to fix it ...)

With limits of ±∞, even slight noise will make it hard to converge. (*)

>>> from mpmath import *
>>> f = lambda x: 1 - sqrt(2)*cosh(x)/sqrt(cosh(2*x))
>>> for x in range(15,20): print x, f(x)
...
15 -9.37028232783632e-14
16 -1.26565424807268e-14
17 -1.99840144432528e-15
18 -2.22044604925031e-16
19 -2.22044604925031e-16

We can explicitly remove the noise, then integrate.

>>> from double_exponential import *
>>> f2 = lambda x: 0.0 if abs(x) > 19 else f(x)

>>> double_exponential(f2, -inf, inf)
(mpf('-0.69314718055994773'), mpf('3.8011316316755028e-10'), 61, 4, 2)

This is almost as good as the "sharpened" version of f:

>>> def f3(x): y=1/cosh(2*x); return -y/(1 + sqrt(1+y))
...
>>> double_exponential(f3, -inf, inf)
(mpf('-0.69314718055994584'), mpf('3.8011260805603797e-10'), 69, 4, 2)

(*) Tracing the function calls, errors started at point #10:
old version: f(6704.2155162232693) → mpf('0.0')
new version: f(6704.215516223293) → mpf('-2.2204460492503131e-16')

The old version just get lucky ...

>>> f3(6704.2155162232693)
mpf('-6.1999991241406458e-5824')

robve · 04-10-2021, 04:13 PM

(04-10-2021 12:05 AM)Albert Chan Wrote:
(04-08-2021 01:58 PM)Albert Chan Wrote: I believe there is nothing to worry about.
Integrand suffered from subtraction cancellations.

Amazingly, qthsh (I picked sqrt version, for speed), does not get garbage from this.

lua> f = function(x) n=n+1; return 1-sqrt(2)*cosh(x)/sqrt(cosh(2*x)) end
lua> f01 = function(t) local y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y) end

With default setting of n=6, eps=1e-9, we have:

lua> Q = require 'quad'
lua> n=0; print(Q.qthsh(f01, 0, 1, 1), n) -- Previous Point
-0.6931471805996943 80
lua> n=0; print(Q.qthsh(f01, 0, 1, 0), n) -- Default Zero
-0.6931471805582439 77

Amazing. Confirmed with the C version, reusing previous point:

qthsh a NAN!!!
qthsh b NAN!!!
qthsh a NAN!!!
qthsh b NAN!!!
qthsh a NAN!!!
k=5 err=1.86609e-09 pts=80
qthsh = -0.693147180599693
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
k=6 err=5.67726e-07 pts=139
wp34s = -0.693147748281705

Using zero as default when NaN occurs:

qthsh a NAN!!!
qthsh b NAN!!!
k=5 err=1.80629e-09 pts=77
qthsh = -0.693147180558242
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
wp34s a NAN!!!
wp34s b NAN!!!
k=6 err=4.70357e-12 pts=133
wp34s = -0.693147180559957

Zero as default improves the result. Note that for the wp34s I had to set eps=1e-15 to prevent early termination. With qthsh with eps=1e-12:

qthsh a NAN!!!
qthsh b NAN!!!
k=6 err=1.10518e-14 pts=155
qthsh = -0.693147180559942

- Rob

robve · (This post was last modified: 04-11-2021 01:07 PM by robve.)

Here is my latest contribution with integrated Tanh-Sinh, Exp-Sinh, Sinh-Sinh rules. The Tanh-Sinh rule has the same Michalski & Mosig abscissas and weights, but reformulated to combine the three quadratures into one routine (see the updated doc for the formulas):

Code:

// Tanh-Sinh bounded [a,b], Exp-Sinh (one inf either side), Sinh-Sinh [-inf,inf]

double quad(double (*f)(double), double a, double b, double n, double eps) {

  const double tol = 10*eps;

  double c = 0, d = 1, s, sign = 1, e, v, h = 2;

  int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2

  if (b < a) { // swap bounds

    v = b;

    b = a;

    a = v;

    sign = -1;

  }

  if (isfinite(a) && isfinite(b)) {

    c = (a+b)/2;

    d = (b-a)/2;

    s = f(c);

  }

  else if (isfinite(a)) {

    mode = 1; // Exp-Sinh

    c = a;

    s = f(a + 1);

  }

  else if (isfinite(b)) {

    mode = 1; // Exp-Sinh

    c = b;

    s = f(b - 1);

    d = -d;

    sign = -sign;

  }

  else {

    mode = 2; // Sinh-Sinh

    s = f(0);

  }

  do {

    double p = 0, q, t, eh;

    h /= 2;

    eh = exp(h);

    t = eh/2;

    if (k > 0)

      eh *= eh;

    do {

      double r, w, x, y;

      q = 0;

      r = w = exp(t-.25/t); // = exp(sinh(j*h))

      if (mode != 1) {      // if Tanh-Sinh or Sinh-Sinh

        w += 1/w;           // = 2*cosh(sinh(j*h))

        r -= 1/r;           // = 2*sinh(sinh(j*h))

        if (mode == 0) {    // if Tanh-Sinh

          r /= w;           // = tanh(sinh(j*h))

          w = 4/(w*w);      // = 1/cosh(sinh(j*h))^2

        }

        else {              // if Sinh-Sinh

          r /= 2;           // = sinh(sinh(j*h))

          w /= 2;           // = cosh(sinh(j*h))

        }

        x = c - d*r;

        if (x > a) {

          y = f(x);

          if (isfinite(y))  // if f(x) is finite, add to local sum

            q += y*w;

        }

      }

      else {                // Exp-Sinh

        x = c + d/r;

        if (x > a) {

          y = f(x);

          if (isfinite(y))  // if f(x) is finite, add to local sum

            q += y/w;

        }

      }

      x = c + d*r;

      if (x < b) {

        y = f(x);

        if (isfinite(y))    // if f(x) is finite, add to local sum

          q += y*w;

      }

      q *= t+.25/t;         // q *= cosh(j*h)

      p += q;

      t *= eh;

    } while (fabs(q) > eps*fabs(p));

    v = 2*s;

    s += p;

    ++k;

  } while (fabs(v-s) > tol*fabs(s) && k <= n);

  e = fabs(v-s)/(fabs(s)+eps);

  return sign*d*s*h; // result with estimated relative error e

}

This was highly optimized while keeping the code compact. Optimizations include strength reductions, branch eliminations, and variable reuse when possible to reduce CPU register pressure.

The partially optimized code for the inner loop is perhaps a bit easier to comprehend:

Code:

    t = eh = exp(h);

    if (k > 0)

      eh *= eh;

    do {

      double ch, u, r, x, w;

      ch = (t+1/t)/2; // = cosh(j*h)

      w = r = u = exp((t-1/t)/2); // = exp(sinh(j*h))

      if (mode != 1) {

        u += 1/u; // = 2*cosh(sinh(j*h))

        r -= 1/r;   // = 2*sinh(sinh(j*h))

        if (mode == 0) {

          r /= u; // = tanh(sinh(j*h))

          w = 4/(u*u); // = 1/cosh(sinh(j*h))

        }

        else {

          r /= 2;   // = sinh(sinh(j*h))

          w = u/2; // = cosh(sinh(j*h))

        }

      }

      x = d*r;

      q = 0;

      if (c+x < b) {

        double y = f(c+x);

        if (isfinite(y))

          q = y*w;

      }

      if (mode == 1)

        x = -d/r;

      if (c-x > a) {

        double y = f(c-x);

        if (isfinite(y))

          q += mode == 1 ? y/w : y*w;

      }

      q *= ch;

      p += q;

      t *= eh;

    } while (fabs(q) > eps*fabs(p));

I have verified the optimized routine with several examples, but more testing will be necessary to determine its numerical properties.

- Rob

Edit: minor change to fix mirrored bound checks.

Albert Chan · (This post was last modified: 04-12-2021 01:08 PM by Albert Chan.)

I thought this might be interesting to others, so I am responding PM as a post

robve Wrote:I now see what the difference is between qthsh and quad, it's quite simple: qthsh uses f(a+x) while quad uses f(c-x). When a=0, the check a+x>a always succeeds for nonzero positive x, but c-x>a fails due to cancellation with c=(a+b)/2.

This means that qthsh approaches the zero point aggressively, something we already saw before. It also means that qthsh is not "symmetric" in the bounds, something we also saw before. This could be desirable actually, when functions exhibit interesting behavior close to x=0.

Yes, it is more accurate interpolate from the closest edge.
https://www.hpmuseum.org/forum/thread-16...#pid140084

Even if interpolate closed to center !
(c=(a+b)/2 might have rounding errors; a, b are user inputs, thus can considered exact)

That's why qthsh lua implementation removed a ≤ b limitation.

And, accuracy affect both edges, not just the "zero" side.
Example, (code snippet from previous post)

Code:

      else {                // Exp-Sinh

        x = c + d/r;

        if (x > a) {

          y = f(x);

          if (isfinite(y))  // if f(x) is finite, add to local sum

            q += y/w;

        }

      }

Note that x is interpolated from the finite edge (c), so is very accurate.
With r = exp(sinh(j*h)) > 1, this x is approaching the finite edge, so it is better to test x ≠ c.

If x=c, we throwaway the point. It does not meant we assume f(c)=0.
Instead, we are assuming f(c) * weight = 0, since weight at the edge is miniscule.

This is what "default 0.0" meant.

(04-10-2021 01:45 PM)Albert Chan Wrote: With limits of ±∞, even slight noise will make it hard to converge. (*)

To reduce noise, we should quit summing, if we actually hit the finite edge.
So this is the proposed patch:

Code:

      else {                // Exp-Sinh

        x = c + d/r;

        if (x == c) break;  // x hit finite edge

        y = f(x);

        if (isfinite(y)) q += y/w;

      }

Albert Chan · 04-12-2021, 01:50 PM

(04-12-2021 12:20 PM)Albert Chan Wrote: And, accuracy affect both edges, not just the "zero" side.

An example, with both edge goes to +∞

>>> Q = require 'quad'
>>> f = function(x) return 1/sqrt(1-x*x) end

lua> Q.qthsh(Q.count(f), -1, 1), Q.n
3.141592644894014 59
lua> Q.quad(Q.count(f), -1, 1), Q.n
3.1415926047697873 115

This is not due to default 0.0 vs previous point.
The function has sharp edges. With x in (-1,1), f(x) never returned inf.
With IEEE double, this is the highest it goes:

lua> f(1 - 0x1p-53)
67108864

There is no zero edge here, and midpoint of (-1,1) is exactly 0.
Interpolate from the edges, qthsh use more accurate points.

robve · 04-12-2021, 03:06 PM

(04-12-2021 12:20 PM)Albert Chan Wrote:
robve Wrote:I now see what the difference is between qthsh and quad, it's quite simple: qthsh uses f(a+x) while quad uses f(c-x). When a=0, the check a+x>a always succeeds for nonzero positive x, but c-x>a fails due to cancellation with c=(a+b)/2.

This means that qthsh approaches the zero point aggressively, something we already saw before. It also means that qthsh is not "symmetric" in the bounds, something we also saw before. This could be desirable actually, when functions exhibit interesting behavior close to x=0.

Yes, it is more accurate interpolate from the closest edge.
https://www.hpmuseum.org/forum/thread-16...#pid140084

Let me add to my PM some context: a notable difference are the function evaluation guards at both endpoints in the different implementations which are a+d*r>a in qthsh versus c-d*r and b-d*r<b in qthsh versus c+d*r<b. If the finite endpoints a or b are zero or close to zero, then the corresponding endpoint may be more closely approached by qthsh with additional abscissas. This is important. More accurate results can be achieved with qthsh when the integrand has a significant area close to \( x=0 \), for example \( \int_0^1 1/\sqrt{x}\, dx \) gives:
qthsh points=63 relative error=1e-13
quad (old version using x=c+d*r similar to WP-34S) points=115 relative errror=6e-9

Likewise, WP-34S "pure" Tanh-Sinh (C version with IEEE 754) with eps=1e-15 has a high relative error of 5e-7 after 25 points "convergence" at k=2. With eps=1e-20 WP-34S has a 2e-10 relative error after 217 points and k=5 convergence.

qthsh converged quickly at level k=3 whereas quad using x=c+d*r converged at level k=4 with a larger error, requiring one more level in an attempt to produce a result within eps=1e-9. This is illustrated by the following chart (click to enlarge):

This shows the point distributions of qthsh and quad for this integrand, using a log scale to emphasize the points approaching the zero endpoint a=0.

(04-12-2021 12:20 PM)Albert Chan Wrote: Even if interpolate closed to center !
(c=(a+b)/2 might have rounding errors; a, b are user inputs, thus can considered exact)

And, accuracy affect both edges, not just the "zero" side.
Example, (code snippet from previous post)

Yes, c can cause rounding errors and loss of precision in a+b. Whenever we add or subtract there can be a (significant) loss. In the end, IEEE 754 double precision does not help. When c is used we should also compare to c (at least in this case). To do so in the code, the Exp-Sinh and Sinh-Sinh quadratures should be split from the Tanh-Sinh computation, which I've already done in yesterday's updated document (that will continue to evolve.)

(04-12-2021 12:20 PM)Albert Chan Wrote: To reduce noise, we should quit summing, if we actually hit the finite edge.

Good point! It makes a lot more sense to stop Exp-Sinh at this point when x hits c for x=c+d/r and avoid using a like qthsh. The weight gets very small when approaching the finite endpoint.

BTW. one trick to get a and b numerically corrected to prevent noise with c:
c=(a+b)/2
d=(b-a)/2
a=c-d
b=c+d
This helps, but it is "iffy" to use for this algorithm, because it (slightly) changes the endpoints! In other contexts these tricks can be important, for example in numerical differentiation e.g. h=1e-3, x=a-h then h is NOT the difference between a and x exactly, but the adjusted h=a-x is.

Note that the Exp-Sinh and Sinh-Sinh quadratures described in the document use a Michalski & Mosig simplification to drop \( \pi/2 \). These rules are transformations by substitution anyway. Having a few more points in the area between the endpoints looks advantageous, judging from the 818 integrals tested in the spreadsheet compared to a "pure" Tanh-Sinh rule. Other parameterizations are possible to control the point distributions, e.g. Michalski & Mosig in their paper start with h=1.5 instead of 1.

An update of the document and code is forthcoming.

- Rob

robve · 04-12-2021, 03:58 PM

(04-12-2021 12:20 PM)Albert Chan Wrote: To reduce noise, we should quit summing, if we actually hit the finite edge.
So this is the proposed patch:

Code:

else { // Exp-Sinh x = c + d/r; if (x == c) break; // x hit finite edge y = f(x); if (isfinite(y)) q += y/w; }

And remove the guard x<b for Exp-Sinh and Sinh-Sinh:

Code:

        x = c + d*r;

        y = f(x);

        if (isfinite(y))       // if f(x) is finite, add to local sum

          q += y*w

robve · (This post was last modified: 04-15-2021 05:44 PM by robve.)

I ran 1084 test integrals on the current quad implementation: 818 definite integrals, 208 one-sided inf integrals, and 58 -inf-inf integrals (latest zip file with the results). This shows that:
- reusing previous points is better than zero default function values (58 versus 5) when hitting the endpoints
- default zero function values for inf/nan singularities is better (8 versus 4)

EDIT adding these additional related observations here instead of in a new post:

- A decay factor of .5 or .8 (i.e. fp *= .8 and fm *= .8 when a singularity occurs) produce 8 better and 4 worse qthsh results of the 818 integrals tested. This is the same as using zero by default, which affect the same 8 and 4 integrals. This might just be noise.
- With an initial h=1.5 in qthsh's first iteration instead of 1, 399 integrals are worse, 391 better with respect to the absolute error with the exact integral value.
- With an initial h=0.5 in qthsh's first iteration, 45 are worse, 55 are better with respect to the absolute error with the exact integral value. However, the average number of points to converge went up from 105 to 125.
- Using the pi/2 factor to reintroduce the "pure" Tanh-Sinh rule in qthshs more often produced less accurate results (417) than better results (363) like the WP-34S implementation, taking the same number of points to converge.

The qthsh and quad implementations and the doc will be updated accordingly.

A few more tests are in the pipeline. Almost finished Smile

- Rob

Albert Chan · 04-14-2021, 03:17 PM

Lua implementation of quad(), using qthsh code to handle finite intervals.
Note: Q.quad is just a shorthand for Q["quad"], a function stored in table Q

Code:

function Q.quad(f, a, b, d, n, eps)

    d = d or 1

    n = n or 6

    eps = eps or 1e-9

    local k, c, t0, s, e = 0, 0, 7.38905609893065

    local mode = (a-a==0 and 0 or 1) + (b-b==0 and 0 or 2)

    if mode == 0 then               -- tanh-sinh

        c = d                       -- decay factor, 0 to 1

        d = (b-a)/2

        s = f((a+b)/2)

    elseif mode == 3 then           -- sinh-sinh

        s = f(c)

    else                            -- exp-sinh

        c = (mode==2) and a or b    -- c = finite edge

        if (a>b) == (mode==2) then d = -d end

        s = f(c+d)

    end

    repeat

        local p, t = 0, sqrt(t0)

        local t1 = k==0 and t or t0

        t0 = t                      -- t0 = e, e^.5, e^.25, ...

        if mode == 0 then           -- tanh-sinh

            local fp, fm = 0, 0

            repeat

                local u = exp(1/t-t)

                local r = 2*u/(1+u)

                local x, y = d*r

                if a+x ~= a then

                    y = f(a+x)

                    if y-y == 0 then fp=y else fp=fp*c end

                end

                if b-x ~= b then

                    y = f(b-x)

                    if y-y == 0 then fm=y else fm=fm*c end

                end

                y = (t+1/t)*r/(1+u) * (fp+fm)   -- apply weight

                p, t = p+y, t*t1

            until abs(y) <= eps*abs(p)

        else

            t = t * 0.5

            repeat

                local r = exp(t-0.25/t) -- = exp(sinh(j*h))

                local w, q, y = r, 0

                if mode == 3 then       -- sinh-sinh

                    r = 0.5*r - 0.5/r   -- = sinh(sinh(j*h))

                    w = 0.5*w + 0.5/w   -- = cosh(sinh(j*h))

                    y = f(c - d*r)

                    if y-y==0 then q = y*w end

                else                    -- exp-sinh

                    y = c + d/r         -- y approaching c

                    if y == c then break end

                    y = f(y)

                    if y-y==0 then q = y/w end

                end

                y = f(c + d*r)

                if y-y==0 then q = q + y*w end

                q = q * (t+0.25/t)      -- q *= cosh(j*h)

                p, t = p+q, t*t1

            until abs(q) <= eps*abs(p)

        end

        s, e, k = s+p, abs(s-p), k+1

    until e <= 10*eps*abs(s) or k>n

    if mode==1 or mode==3 and a>b then d=-d end

    return d*ldexp(s,1-k), e/abs(s)

end

---

I noticed exp-sinh is actually doing 2 integrals, simultaneously.

\(\int_a^∞ = \int_a^{a+d} + \int_{a+d}^∞ \)

Both integrals are using same number of points, but the last one is much harder.
If possible, we should pick d so that first one do most of the work.

This version allowed to input d, instead of hard-coded d=1.
If we have a clue about the shape of curve, it helps. Example, Standard Normal Distribution.

lua> Q = require 'quad'
lua> k = 1/sqrt(2*pi)
lua> Inf = math.huge
lua> pdf = function(x) return k*exp(-x*x/2) end

lua> Q.quad(Q.count(pdf), 1, Inf), Q.n
0.15865525392847651 129
lua> Q.quad(Q.count(pdf), 1, Inf, 4), Q.n
0.15865525393145705 69

Above is now as efficient as integrating closed interval of [0, 1]

lua> 1/2 - Q.quad(Q.count(pdf), 0, 1), Q.n
0.15865525393145719 58

Same idea for \(\int_{-∞}^∞\), we wanted \(\int_{-d}^d\) to do the most work.

lua> Q.quad(Q.count(pdf), -Inf, Inf), Q.n
0.999999999999999 119
lua> Q.quad(Q.count(pdf), -Inf, Inf, 4), Q.n
0.9999999999999996 43

Had we optimized specifically for even functions, above only need (n+1)/2 points !

Albert Chan · 04-15-2021, 05:14 PM

(04-12-2021 03:06 PM)robve Wrote:
(04-12-2021 12:20 PM)Albert Chan Wrote: To reduce noise, we should quit summing, if we actually hit the finite edge.

Good point! It makes a lot more sense to stop Exp-Sinh at this point when x hits c for x=c+d/r and avoid using a like qthsh. The weight gets very small when approaching the finite endpoint.

More importantly, infinite "edge" has weight that gets very very big.

For exp-sinh, with limits of [c, d*inf], q before scaled up by cosh(j*h):

q = f(c + d/r)/r + f(c + d*r)*r

When r is huge, second term might generate magnified noise, messing up the sum.

Example, ∫(x^-0.99, x=0..1)

\(\displaystyle \int_0^1 x^{-0.99}\, dx
= \int_0^∞ (e^{-y})^{-0.99} (e^{-y}\, dy)
= \int_0^∞ e^{-0.01y} \, dy
\)

If we apply transformation without last simplifed step, integrand generate much noise.

>>> Q = require 'quad'
>>> f = function(x) return x^-0.99 end
>>> f1 = function(y) y=exp(-y); return f(y)*y end

lua> Q.quad(Q.count(f), 0, 1), Q.n
99.92043566319589 654
lua> Q.quad(Q.count(f1), 0, huge), Q.n
99.93721615545337 459

Now, try simplified form, avoiding the noise

lua> f2 = function(y) return exp(-0.01*y) end
lua> Q.quad(Q.count(f2), 0, huge), Q.n
99.9999999878134 237

Nice, but we can do better ! (see previous post)

exp(0) = 1
exp(-0.01y) = 0.001 → -y/100 = ln(10^-3) = -3*ln(10) ≈ -3*2.3 ≈ -7

\(\int_0^∞ = \int_0^{700} + \int_{700}^∞ :\)

lua> Q.quad(Q.count(f2), 0, huge, 700), Q.n
100 71

robve · 04-15-2021, 07:02 PM

(04-15-2021 05:14 PM)Albert Chan Wrote: Nice, but we can do better ! (see previous post)

exp(0) = 1
exp(-0.01y) = 0.001 → -y/100 = ln(10^-3) = -3*ln(10) ≈ -3*2.3 ≈ -7

\(\int_0^∞ = \int_0^{700} + \int_{700}^∞ :\)

lua> Q.quad(Q.count(f2), 0, huge, 700), Q.n
100 71

Picking a smart "central point" 'd' in Exp-Sinh (the splitting point between the two quadratures) is advantageous.

A higher accuracy can also be accomplished by manually splitting up the integral into two parts by picking the splitting point smartly:

\(\int_0^∞ = \int_0^{700} + \int_{700}^{701} + \int_{701}^∞ \)

In this case we can also get a more accurate result, while keeping the default Exp-Sinh "central point" 700+1 (the default is d=1 versus d=700 in your example):

\( \int_0^{700} e^{-0.01y}\,dy \):
k=3 err=4.82092e-11 pts=58
quad = 99.9088118034444
\( \int_{700}^\infty e^{-0.01y}\,dy \):
k=5 err=8.41427e-09 pts=235
quad = 0.0911881965443389
Combined \( \int_0^\infty e^{-0.01y}\,dy \) = 99.9088118034444 + 0.0911881965443389 = 99.999999999988739

Thus we get three more digits (12 digits) compared to Exp-Sinh with the default d=1 splitting point:
\( \int_0^\infty e^{-0.01y}\,dy \):
k=5 err=8.41426e-09 pts=237
quad = 99.9999999878134

On the other hand, the splitting point d=700 gives the exact result as you point out, which is very nice indeed:

k=3 err=9.08301e-12 pts=71
quad = 100

The absolute error estimates (100*9e-12 ~ 9e-10) is in the same order (100*5e-11+0.09*8e-9 ~ 5e-9). So using this as a "black box" we won't be able to tell the difference, unless we can get the relative error further down, which seems really cheap to do as we only needed 71 points to get about 11 digits.

It makes sense to specify 'd' as an optional parameter or at least say something about it in the document.

Overall, these methods are interesting. But as always one needs to be smart about the type of integrand and consider splitting up integrals into parts. The 1084 integral tests to check these methods do none of that, obviously.

- Rob

robve · 04-15-2021, 09:29 PM

I should add that the choice of d=1 for Exp-Sinh and c=0 for Sinh-Sinh bothered me as well, because the Exp-Sinh and Sinh-Sinh integrals are split in two parts. I have not yet looked into articles for suggestions for these values, although for Sinh-Sinh I would expect that the best c is where f has (some) symmetry e.g. f(c+x) = f(c-x).

In the meantime, I've mostly worked on testing different parameterizations e.g. h=.5 and h=1.5 initially and compared the "pure" Tanh-Sinh with the pi/2 factors etc.

(04-15-2021 05:14 PM)Albert Chan Wrote: Nice, but we can do better ! (see previous post)

exp(0) = 1
exp(-0.01y) = 0.001 → -y/100 = ln(10^-3) = -3*ln(10) ≈ -3*2.3 ≈ -7

The derivation to select d=700 for the integrand exp(-0.01*y) is clear, except for the initial requirement that is not obvious from the integrand and interval:
exp(-0.01y) = 0.001
In this case I can empirically verify this by a rough gradient search for the value of d that produce the lowest errors, but is there a general algebraic formulation using the inverse function to determine where d might be optimal? Any literature on that?

- Rob

Albert Chan · (This post was last modified: 04-16-2021 01:15 PM by Albert Chan.)

(04-15-2021 09:29 PM)robve Wrote: The derivation to select d=700 for the integrand exp(-0.01*y) is clear, except for the
initial requirement that is not obvious from the integrand and interval: exp(-0.01y) = 0.001

My aim is to pick d such that \(\int_0^d\) covered bulk of the area, less work for \(\int_d^∞\)
Based on experience with plots, if value is below 1/1000 of peak, it seems to touch the x-axis.

Actual number is not that critical. Hitting exact 100 is just lucky.

lua> Q = require 'quad'
lua> f2 = function(x) return exp(-0.01*x) end
lua> for d=500,1500,100 do print(d, Q.quad(Q.count(f2), 0, huge, d), Q.n) end

500 99.99999999999997 71
600 100 71
700 100 71
800 99.99999999999997 71
900 99.99999999999999 71
1000 99.99999999999997 71
1100 99.99999999999997 71
1200 99.99999999999997 71
1300 99.99999999999999 71
1400 99.99999999999996 71
1500 99.99999999999997 71

Note that we could accomplish the same by compressing the function, by factor d

lua> x, d = 20, 512
lua> f2(x) * 100 -- exact result, for \(\int_x^∞\)
81.87307530779819

lua> f3 = function(d) return function(x) return f2(x*d)*d end end
lua> Q.quad(Q.count(f3(d)), x/d, huge/d), Q.n
81.87307530779817 69

lua> Q.quad(Q.count(f2), x, huge, d), Q.n
81.87307530779817 69

Both expressions are equivalent (because I picked d = power-of-2)
But setting d is cost free, and less messy.

Albert Chan · 04-17-2021, 12:03 AM

(04-15-2021 07:02 PM)robve Wrote: \( \int_0^{700} e^{-0.01y}\,dy \):
k=3 err=4.82092e-11 pts=58
quad = 99.9088118034444
\( \int_{700}^\infty e^{-0.01y}\,dy \):
k=5 err=8.41427e-09 pts=235
quad = 0.0911881965443389

We can build exp-sinh transformed integrand, and see the curve, in XCas.

XCas> f(x) := exp(-0.01*x)
XCas> y := d * exp(sinh(t))
XCas> g := unapply(f(y) * diff(y,t), (d,t))
XCas> g2(d) := [g(d,-t), g(d,t)]
XCas> int(g2(1), t=0..5) → [0.995016625083, 99.0049833749]
XCas> int(g2(700), t=0..5) → [99.9088118034, 0.0911881965555]

XCas> plot(g2(700), t=0..3, color=[blue, red])

The plot showed a nice bell-shaped curve, easy to get area, simply by summing squares.
see https://www.hpmuseum.org/forum/thread-13...#pid127590

\(\displaystyle\begin{align}
\int_0^∞ f(x)\,dx &= \int_{-∞}^∞ g(t)\,dt
= \int_{-∞}^∞ {g(t) + g(-t) \over 2} \,dt \\
&= \int_0^∞ g(-t)\,dt + \int_0^∞ g(t)\,dt \\
&= \int_0^d f(x)\,dx\;\,+ \int_d^∞ f(x)\,dx
\end{align} \)

Note we take advantage of even functions, and fold the integrals.
see Integration trick: clever use of even and odd parts - John D. Cook

Albert Chan · 04-17-2021, 01:14 PM

Revisiting original quad(), where we interpolate from the center, we may be able to improve it.
double_exponential.py is doing the same way, so the patch might be useful.

This was the snippet from before, started with r = w = exp(t-.25/t); // = exp(sinh(j*h))
For tanh-sinh, r ← (r-1/r) / (r+1/r), which may suffer rounding errors > 1/2 ulp

(04-11-2021 03:05 AM)robve Wrote:

Code:

w += 1/w; // = 2*cosh(sinh(j*h)) r -= 1/r; // = 2*sinh(sinh(j*h)) if (mode == 0) { // if Tanh-Sinh r /= w; // = tanh(sinh(j*h)) w = 4/(w*w); // = 1/cosh(sinh(j*h))^2 } else { // if Sinh-Sinh r /= 2; // = sinh(sinh(j*h)) w /= 2; // = cosh(sinh(j*h)) }

Instead, we calculate from the edge: r ← (r-1/r) / (r+1/r) = 1 - 2/(1+r*r)
This is the proposed patch, where updated r should be much more accurate, at the edges.

Code:

    if (mode == 0) {        // if Tanh-Sinh

        r = 2/(1+r*r);

        w = w * r;

        r = 1 - r;          // = tanh(sinh(j*h))

        w = w * w;          // = 1/cosh(sinh(j*h))^2

    }

    else {                  // if Sinh-Sinh

        r = (r-1/r) / 2;    // = sinh(sinh(j*h))

        w = (w+1/w) / 2;    // = cosh(sinh(j*h))

    }

(04-12-2021 01:50 PM)Albert Chan Wrote: lua> Q = require 'quad'
lua> f = function(x) return 1/sqrt(1-x*x) end

lua> Q.qthsh(Q.count(f), -1, 1), Q.n
3.141592644894014 59
lua> Q.quad(Q.count(f), -1, 1), Q.n
3.1415926047697873 115

This is still not as good as interpolate from the edge, but it improved a lot.
Redoing the above example, with the proposed patch.

lua> O = require 'quad0' -- original quad(), interpolate from center, but more accurate r
lua> O.quad(O.count(f), -1, 1), O.n
3.1415926448906006 59