simon-git: spigot (master): Simon Tatham
Commits to Tartarus hosted VCS
tartarus-commits at lists.tartarus.org
Thu Aug 8 18:30:19 BST 2019
TL;DR:
1ef963d Add a debug wrapper for pi^2/6.
41ae7a5 Add a debug-only constant 'Stealth1'.
8367d3e MonotoneHelper: new built-in slope-adjustment system.
561484d Fix monotonicity bug in Hg, using the new slope system.
e5f80e5 Migrate Bessel and trig integrals to new slope-adjuster.
34f87c2 Speed up the cosine integral by pre-squaring the input.
Repository: https://git.tartarus.org/simon/spigot.git
On the web: https://git.tartarus.org/?p=simon/spigot.git
Branch updated: master
Committer: Simon Tatham <anakin at pobox.com>
Date: 2019-08-08 18:30:19
commit 1ef963dce44fc41998ad011ebf4597c559e555b0
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=1ef963dce44fc41998ad011ebf4597c559e555b0;hp=ecae2f7d61afec442dbd5d0581fe34bc7dba26af
Author: Simon Tatham <anakin at pobox.com>
Date: Wed Aug 7 19:20:12 2019 +0100
Add a debug wrapper for pi^2/6.
I happened to notice in passing that although I've got a special-
purpose class for computing this constant, it's not exposed in the
expression language _even_ in debug-functions mode. Like all the other
things for internal use within spigot, it should be.
This is the first case where I've needed a debug-only _constant_
function, so I also had to add the ConstantDebugFnWrapper class in
funcs.h (but that was trivial).
consts.cpp | 4 ++++
funcs.h | 5 +++++
2 files changed, 9 insertions(+)
commit 41ae7a5c23588ee3767d49e54abe0275ed58ad3b
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=41ae7a5c23588ee3767d49e54abe0275ed58ad3b;hp=1ef963dce44fc41998ad011ebf4597c559e555b0
Author: Simon Tatham <anakin at pobox.com>
Date: Wed Aug 7 19:20:20 2019 +0100
Add a debug-only constant 'Stealth1'.
This has the value 1, but not obviously: it doesn't announce itself as
a rational, and it generates a non-terminating stream of matrices that
never quite zero in on the ultimate answer. So it's useful for
debugging code paths that are skipped in the easy rational case,
without sacrificing the convenience of having nice easy actual input
values: evaluating a function like sin at (say) (Stealth1 / 5) should
give you the same answer as at 1/5, but forces MonotoneHelper to get
involved.
Of course you could do the same kind of thing already by adding some
expression like (pi-pi), or multiplying by (pi/pi), or some such. But
this approach wastes less CPU on computations unrelated to what you
were actually trying to debug.
(Perhaps in future, if it ever becomes necessary, I could change this
so that it takes parameters indicating how fast to converge to 1, or
different ways in which it could converge, e.g. flipping the interval
sense every time, or converging via continued-fraction-style iterated
reciprocation...)
consts.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 50 insertions(+)
commit 8367d3eb7b235dc2ed9cc3816a434ee50f8b7865
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=8367d3eb7b235dc2ed9cc3816a434ee50f8b7865;hp=41ae7a5c23588ee3767d49e54abe0275ed58ad3b
Author: Simon Tatham <anakin at pobox.com>
Date: Wed Aug 7 07:27:11 2019 +0100
MonotoneHelper: new built-in slope-adjustment system.
This reimplements in a convenient central location the strategy I've
used in several parts of the code when I want to use MonotoneHelper
but the function I'm trying to compute isn't monotonic: find a real
number s such that g(x) = f(x) + sx _is_ monotonic, use MonotoneHelper
with a constructor that computes the adjusted function g, and subtract
kx back off the output of spigot_monotone.
This new approach uses basically the same technique, but pushes it
right inside MonotoneHelper. To activate it, you fill in an optional
extra method in the MonotoneConstructor, which is called the first
time MH generates an interval it wants to compute f at the bounds of,
and should return a suitable s.
MonotoneHelper takes care of all the rest itself: your constructor
just generates a spigot representation of the function you actually
want, with no adjustment needed, and the output of spigot_monotone
will be the target function too, with no postprocessing needed.
(Internally, MonotoneHelper doesn't even do the _literal_ add-and-
subtract strategy: instead, it uses your 'slope adjustment' to get
rigorous bounds on f on each interval [a,b], by saying: I know f(a),
and I know f(b), and I know that on that interval f can slope
downwards by at most s, and that's enough to work out an interval in
which f(x) must lie for all a<=x<=b.)
So, why is this better?
Firstly, it simplifies the call sites, by separating the adjustment
from the description of the function to be adjusted. Now you only have
to construct a spigot representation of the target function itself,
without having to take the slope issue into account _at the same
time_, which makes those representations easier to understand and
contain fewer special cases.
Secondly, it should be faster, because it doesn't require an extra
evaluation in postprocessing of whatever complicated spigot gave rise
to the input value x: all the adjustment _and_ un-adjustment is done
using by-products of the rational approximations computed once within
MonotoneHelper, and nothing needs to be recomputed separately outside
it.
Thirdly, because this compensation is only actually needed when you're
using MonotoneHelper to evaluate your function at an irrational, it
can be completely sidestepped when the spigot_monotone entry point
notices that your input value is rational: then it can just construct
an instance of your end-user class without having to bother with any
of this at all. This is especially important if computing s is
_expensive_, which it's never been yet, but it's about to be, because ...
... fourthly and most importantly, this new system means that if the
target function has _unbounded_ derivative, then you can still compute
it using the slope-based adjustment technique. Because where the old
system had to give an s that worked for the whole domain, now you can
wait to deliver s until you know what bounded interval you need to
make the function monotonic _on_. So any target function which is
continuously differentiable can now be handled by this system, because
on any closed bounded interval, the derivative will also be bounded,
so _some_ s will work, even if you have to think about it carefully
after you find out your interval of the domain.
This commit just introduces the new system in MonotoneHelper, and
nothing actually starts using it yet. There's a default implementation
of the new MonotoneConstructor method, which does no adjustment, so
for the moment, the existing constructor classes carry on working the
way they always have. In following commits I'll actually start using
the new mechanism.
funcs.h | 30 ++++++++++
monotone.cpp | 177 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
2 files changed, 203 insertions(+), 4 deletions(-)
commit 561484def8bc86858a4277c0622f0fe19d40926f
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=561484def8bc86858a4277c0622f0fe19d40926f;hp=8367d3eb7b235dc2ed9cc3816a434ee50f8b7865
Author: Simon Tatham <anakin at pobox.com>
Date: Wed Aug 7 21:59:20 2019 +0100
Fix monotonicity bug in Hg, using the new slope system.
When I turned the generalised hypergeometric function from a debug one
into a public one recently, I jumped the gun a bit: I had completely
forgotten that it's not safe to pass it directly to MonotoneHelper,
because it's not reliably monotonic. (For example, cos(sqrt(x)) is
expressible in terms of an instance of Hg, and that's definitely
wiggly in the positive domain.)
Unlike all the previous non-monotonic functions I've contrived to use
MonotoneHelper with regardless, this one doesn't have a _fixed_
maximum slope, so you can't solve the problem by just adding x, or 2x,
or whatever, in the series and subtracting it again afterwards. You
need to compute the exact compensation value based on the particular
parameters to that specific hypergeometric series - and not only that,
but the derivative of a given instance of Hg can easily be unbounded
(proof: exp is an instance of Hg!) so you _also_ need to base it on
a specific interval.
In other words, this is a problem that needs the MonotoneHelper slope
compensation system that I've just added! So here is a bug fix, which
implements HgConstructor::extra_slope_for_interval in a way that takes
all of that into account.
My job was made easier by the fact that the canonical representation
for a hypergeometric power series is surely _designed_ to be easy to
differentiate: I think that must be why the usual definition has a
mandatory n! in the denominator of each term, because it means that
the x^n/n! in each term turns into x^{n-1}/(n-1)! when it's
differentiated, so all you have to do is adjust the explicit
parameters by one step each and multiply by a constant.
To get a good lower bound, I also added a feature to HgRational that
allows it to only consider power-series terms of a specified sign,
because including only terms of one sign makes the series monotonic in
the input value (at least provided you stay on the same side of zero),
so I can get a bound over the whole interval by only evaluating at the
largest positive and largest negative values.
(Incidentally, one way that this whole procedure can go wrong is that
the attempt to bound f' fails because the series for f' doesn't
converge at one end of the interval. That arises if MonotoneHelper's
first approximating interval to x is too large, so that one end of it
is completely beyond the valid domain of the function. So you return
failure to MonotoneHelper, which tells it to try again with a better
approximating interval. And that's done by catching the domain_error
coming back from the attempt to evaluate f' - which justifies the
precaution I took back in commit c6a66492f of breaking up spigot_error
into subclasses so that I could catch them independently. I _knew_
that would turn out to be useful, and now it has!)
hypergeom.cpp | 129 ++++++++++++++++++++++++++++++++++++++++++++++++++++++----
test.sh | 3 ++
2 files changed, 125 insertions(+), 7 deletions(-)
commit e5f80e584741f0d536d5c6eff007697cb89132bd
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=e5f80e584741f0d536d5c6eff007697cb89132bd;hp=561484def8bc86858a4277c0622f0fe19d40926f
Author: Simon Tatham <anakin at pobox.com>
Date: Wed Aug 7 19:22:43 2019 +0100
Migrate Bessel and trig integrals to new slope-adjuster.
This removes the previous workarounds in which the MonotoneConstructor
instantiates a class that computes f(x)+x, and the function that calls
spigot_monotone then has to subtract x->clone() from the output. Now
the compensation for non-monotonicity in this whole function family is
done using the new shiny system, so all the series representations in
TrigIntSeries look the way you obviously wanted them to look, and the
extra calls to spigot_mobius in BesselJConstructor::construct are gone.
hypergeom.cpp | 46 +++++++++++++++---------------
trigint.cpp | 89 +++++++++++++++++++++++++----------------------------------
2 files changed, 60 insertions(+), 75 deletions(-)
commit 34f87c27cd20e509e7ba0f145d933c443faaf1d3
web diff https://git.tartarus.org/?p=simon/spigot.git;a=commitdiff;h=34f87c27cd20e509e7ba0f145d933c443faaf1d3;hp=e5f80e584741f0d536d5c6eff007697cb89132bd
Author: Simon Tatham <anakin at pobox.com>
Date: Thu Aug 8 18:23:53 2019 +0100
Speed up the cosine integral by pre-squaring the input.
In the previous commit, one non-monotonic function that I _didn't_
migrate to the new slope compensation system was the ordinary
trigonometric cos(). I tried it, and it did give the right answers,
but it turned out that the existing strategy - having the series class
compute cos(sqrt(x)), and feeding it the output of spigot_square()
applied to the user's original input value - was noticeably faster.
I'm not sure of exactly why, but my guess is that computing rational
approximations directly to x^2 is more numerically efficient, in that
the rationals you get that way have more accuracy per digit than the
ones you get by taking an approximation to x itself and squaring it.
That suggests that pre-squaring the input is a trick that will
generally improve performance whenever you can get away with it, i.e.
whenever the target function is both even and computed by a power
series (in which case the series will depend only on x^2 and never
need x itself).
Another function with that property is Cin(), used as a primitive for
the cosine integral Ci. And indeed, testing shows a considerable
speedup in that function as well if I pre-square the input and adjust
the power-series computation to expect that. So this is a pure
optimisation change that improves Cin and its dependent Ci.
trigint.cpp | 21 ++++++++++++++++++---
1 file changed, 18 insertions(+), 3 deletions(-)
More information about the tartarus-commits
mailing list