@R.M. Yes, in fact. All your points are valid and I agree with every single one.

@J.M. Accidentally, from List @@ NDSolveValue[{y'[t] == y[t], y[0] == 1}, y, {t, 0, 1}, InterpolationOrder -> All, Method -> "ExplicitRungeKutta"] and guessing & testing what things meant. Solving ODEs in terms of Chebyshev has been around for decades. No convenient way of constructing such an I.F., AFAIK, but it should be possible to write a function to do it. Guessing about internals always seems a bit risky, tho.

sol = NDSolveValue[{y'[t] == y[t], y[0] == 1}, y, {t, 0, 1},
   InterpolationOrder -> All, Method -> "ExplicitRungeKutta"];
{{grid}, coeff, {{idcs, types}}} = (List @@ sol)[[3 ;;]];

y0 = Piecewise[#, Indeterminate] &@Quiet[
   Map[
    {chebeval[coeff[[Last@#1]]][Rescale[t, grid[[#1]], {-1, 1}]],
      grid[[First@#1]] <= t <= grid[[Last@#1]]} &,
    Rest@idcs],
   CompiledFunction::cfsa]

@MichaelE2 Hmm, interesting. Chebyshev is not what is usually used for dense output of RK methods, but it's good that they're branching out of Hermite.

@J.M. The default LSODA switches between Hermite and "LocalSeries", when Interpolation -> All is specified.

@MichaelE2 Missed that, I guess. Old Mathematica was using LSODA under the hood and interpolated everything with Hermite.

chebeval[c0 : {__Real}] :=
  With[{c1 = Reverse[c0], deg = Length@c0 - 1},
   Compile[{x},
    Block[{c = c1, s = 0., s1 = 0., s2 = 0.},
     Do[s = c[[i]] + 2 x*s1 - s2;
      s2 = s1; s1 = s, {i, deg}];
     Last@c + x*s1 - s2], RuntimeAttributes -> {Listable},
    RuntimeOptions -> {"EvaluateSymbolically" -> False}]];

@MichaelE2 I gathered it was the Clenshaw routine you wrote. ;)

@J.M. Yes, but I added the RuntimeOptions -> {"EvaluateSymbolically" -> False}, which is important here.

Probably should've written a generic one, with the coefficients as a argument, but I was just testing...

@MichaelE2 In any case, this was a nice find. :) Hopefully, exploiting the built-in machinery to make your own chebfun is not too far off.

@J.M. Thanks! :D

@J.M. I believe NDSolve does create a Hermite interpolation when the default InterpolationOrder is used. Setting it to All takes longer, probably because of computing the series.

@MichaelE2 Yes, probably.

    can some one spot what I am doing wrong? Why does this work:
    Expectation[x + y, Distributed[{x, y}, DirichletDistribution[{1, 4, 5}]]]
    gives a number. But this does not work:
    Expectation[x + y, Distributed[{x, y}, UniformDistribution[{-1, 1}]]]
? it returns unevaluated

I have to write it like this to work:

Expectation[x + y, {Distributed[x, NormalDistribution[]], Distributed[y, NormalDistribution[]]}]

I meant like this:
Expectation[
 x*y, {Distributed[x, UniformDistribution[{-1, 1}]],
  Distributed[y, UniformDistribution[{-1, 1}]]}]

@Nasser I'm no expert but to my knowledge there are distributions that are defined for vectors such as the the Dirichlet distribution and the multinomial distribution. Your first way of writing it is suitable for these distributions, it means not that each variable is distributed like that but that the vector in its entirety is distributed like that. Then there are other distributions like the beta distribution and the uniform distribution which are defined for numbers, not vectors.

@yode At the moment it's really a series of hacks, kludges and workarounds, and would need serious cleanup. I may consider cleaning it up, but that can take several days.

Since there's no native Voronoi computation on spheres on Mma, it implements stochastic Lloyd's algorithm to spread cells nicely, and then computes Voronoi cell adjacency graph using symbolic regions (essentially numerical semialgebraic sets) and Resolve@Exists.

This adjacency graph is then rewritten into a Boolean equation system for which an instance that solves a four-coloring problem is solved, with additional constraint of every color having exactly 15 cells.

Because region discretization doesn't work for this kind of shells well enough yet in Mma (Method -> "Semialgebraic" is buggy), next the cell corners are computed somewhat similarly to the adjacency graph, and sorted to follow same handedness. These polygons are split to triangle strips around the Voronoi cell centres.

After this, triangles are subdivided couple of times and normalized to the sphere surface, and surface normals are corrected for visual appearance. Colors of the cells are chosen to have uniform lightness and saturation (well, one might debate about that, but that's how it should be) and to be rotated 45 from LCH primaries in a way, they land nicely on both sRGB and CMYK gamuts.

That's it. Not exactly trivial, at least a "screenful" of code (which is currently quite messed up).

Oh, and there's an extra step to reposition almost-the-same points (numerical errors arising from operating in machine precision) together in order to make the model "airtight" for 3D print modelling.

If one wouldn't need to work around certain bugs and Mma would have certain domain-specific functionalities it lacks now, the implementation would possibly be only half a dozen lines long... but it isn't.

About this (closed and deleted) question:

mathematica.stackexchange.com/questions/104223

There's actually a paper explaining how a CAS might deal with this. I'm inclined to believe Mathematica is using some variant of it under the hood, and at least a short pointer to that paper would be an admissible answer.

I've undeleted it, but I'm loathe to unilaterally use my mod powers to reopen it. If people are interested, please vote to reopen.