drmikeando

To me the neat bit isn't that it got the exponential decay right - that's pretty standard, its that it realised there were two different timescales for the decay and got ball-park numbers for them pretty well.

This is the kind of model you would expect from a simple cylindrical model of the coffee cup with some inbuilt heat capacity of its own.

However, those decay coefficients are going to be very dependent of the physical parameters of your coffee cup - in particular the geometry and thermal parameters of the porcelain. There's a lot of assumptions and variability to account for that the models will have to deal with.

While this is a great article, I feels it buries the lede.

For me, the key insight was from the last paragraph of the article:

C++23 introduces "deducing this", which is a way to avoid the performance cost of dynamic dispatch without needing to use tricks like CRRT, by writing:

    class Base {
    public:
      auto foo(this auto&& self) -> int { return 77 + self.bar(); }
    };

    class Derived : public Base {
    public:
      auto bar() -> int { return 88; }
    };

My mind was exploded by this somewhat similar technique https://en.wikipedia.org/wiki/Tanh-sinh_quadrature - it uses a similar transformation of domain, but uses some properties of optimal quadrature for infinite sequences in the complex plane to produce quickly convergent integrals for many simple and pathological cases.

You can view this result as the convolution of the signal with an exponentially decaying sine and cosine.

That is, `y(t') = integral e^kt x(t' - t) dt`, with k complex and negative real part.

If you discretize that using simple integration and t' = i dt, t = j dt you get

    y_i = dt sum_j e^(k j dt) x_{i - j}
    y_{i+1} = dt sum_j e^(k j dt) x_{i+1 - j}
            = (dt e^(k dt) sum_j' e^(k j' dt) x_{i - j'}) + x_i 
            = dt e^(k dt) y_i + x_i

    z_{i+1} = dt e^(k dt) z_i + A x_i

We can write `z_i = e^{w dt i} r_i` where w is the imaginary part of k

   e^{w dt (i+1)} r_{i+1} = dt e^(k dt) e^{w dt i} r_i + A x_i
             r_{i+1} = dt e^((k - w) dt) r_i + e^{-w dt (i+1) } A x_i
                     = (1-alpha) r_i + p_i x_i

The neat thing about recognising this as a convolution integral, is that we can use shaping other than exponential decay - we can implement a box filter using only two states, or a triangular filter (this is a bit trickier and takes more states). While they're tricky to derive, they tend to run really quickly.

IMO the reason the compiler doesn't add special cases for the simplest version is that it doesn't know which of its _many_ special cases to use. If you actually use the unoptimised version of the code like

    void withSwitch(vector<int>& Values, bool v) {
      if (v) {
        multiply1(Values, 2.0);
      } else {
        multiply1(Values, 3.0);
      }
    }

You can see the code here: https://godbolt.org/z/5beeYe77a