Some faster floorDiv calculations

This suggests potentially faster calculations for floorDiv by 86400 and 86400e9,
respectively the number of seconds and nanoseconds in a normal day.

***The usual disclaimers***begin***

This code is published in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty
of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Also, I have no reason to suppose any of this is in any way novel,
because it applies some basic mathematical concepts,
so I’m not making any claims to copyright.

I haven’t found online examples of this actual code, but that doesn’t mean
there aren’t any online examples which are essentially very similar.

Henry S Warren’s book Hacker’s Delight has a chapter on integer division
by constants, and his discussion there is useful for understanding this code,
but the approach here is somewhat different: his code works for all integers
of a specified number of bits, whereas the code here uses faster approximations
for a subset of a range of integers, switching to the slower “accurate” version
for values outside that subset. The idea is we can use faster approximations
for “usual” values, and deal with the possibility of “unusual” values
by using the slower “accurate” version for those “unusual” values.

***The usual disclaimers***end***

floorDiv86400e9 seems much faster for all 64bit signed integers,
and floorDiv86400e9Approximate seems even faster.

floorDiv86400 seems materially faster for
-466_018_172_160 < x < 466_018_172_160
which if x is a number of seconds represents about
-14_760 years <= x <= +14_760 years,
which seems a more than adequate range for most practical date calculations.
So being slightly slower for unusual values of x seems a reasonable price
to pay for being materially faster for usual values of x.

In very rough benchmarking on a low specification netbook
the average nanoseconds per calculation were:
20 using: floorDiv86400e9Approximate
38 using: floorDiv86400, if x is in the fast calculation range
53 using: floorDiv86400e9
92 using: x >= 0 ? x / 86400 : (x + 1) / 86400 – 1
102 using: floorDiv86400, if x is outside the fast calculation range
267 using: x >= 0 ? x / 86400_000_000_000 : (x + 1) / 86400_000_000_000 – 1

The code follows (remember the disclaimers).
Below the code are what are hopefully proofs that the code will work as intended.

    public static long floorDiv86400e9Approximate(long x) {
        // Returns a value which is -1, +0, or +1 from the correct value.
        // It is useful for when we can use fast approximate values
        // in initial calculations before getting the exact values
        // with just one call to floorDiv86400e9.
        // The result is in the 32bit signed integer range
        // but it is returned as a 64bit signed integer.
        return ((x >> 16) * 417) >> 39;
    }

    public static long floorDiv86400e9(long x) {
        // The result is in the 32bit signed integer range
        // but it is returned as a 64bit signed integer.
        long r = x >> 16; // floorDiv by 2**16, highest power of 2 in 86400e9
        // we have used floorDiv by 2**16 so the next line will not overflow
        long n = (r * 417) >> 39;
        r = r - n * 1_318_359_375L;  // 1_318_359_375L * 2**16 == 86_400e9
        if (r < 0) {
            return n - 1;
        } else if (r >= 1_318_359_375L) {
            return n + 1;
        }
        return n;
        // If you need both the floorDiv result and the floorMod result
        // then it may be worth using the following code "inline":
        //   long n = ((x >> 16) * 417) >> 39;
        //   long r = x - n * 86400_000_000_000L;
        //   if (r < 0) {
        //       n = n - 1;
        //       r = r + 86400_000_000_000L;
        //   } else if (r >= 86400_000_000_000L) {
        //       n = n + 1;
        //       r = r - 86400_000_000_000L;
        //   }
    }

    public static long floorDiv86400(long x) {
        if (x >= 466_018_172_160L) {
            return x / 86_400;
        } else if (x <= -466_018_172_160L) {
            return (x + 1) / 86_400 - 1;
        }
        // If here we can use a faster approximate calculation.
        // If x represents seconds, then this approximation works
        // for a range of about -14_760 years to +14_760 years,
        // a more than adequate range for most practical date calculations.
        return ((x >> 7) * 1_628_906_115L + 814_453_057L) >> 40;
    }


Proofs:

In the following "/" means Rational division, not integer division.

Suppose we want
[EQ.1] floorDiv(x + R, B) == floor((x + R) / B),
where x is an integer, B is a positive rational number,
and R is any rational number such that: Rmin <= R < Rmax

(If we want floorDiv(x, B) where B is a positive integer,
 then Rmin = 0 and Rmax = 1. But, for example, to convert a
 Chronological Julian Day Number (CJDN) to a Julian calendar year,
 then xxxxx xxxxx xxxxx xxxxx xxxxx.)

Suppose we can find positive integers C and G such that 
(C / G) is a good approximation to B.
(For fast calculations we choose C to be a power of 2,
 but we do not make that assumption in this explanation.)

Set E = ((C / G) / B) - 1,
so that (1 + E) == (C / G) / B, 
where E is small, and may be positive or negative or zero.

floorDiv(x + R, B) == floor((x + R) / B)
                   == floor(((x + R) / B) * (C/G) * (G/C))
                   == floor((x + R) * G * ((C/G)/B) / C)
                   == floor((x * G + R*G) * (1 + E) / C)
[EQ.2]             == floor((x * G + x * G*E + R*G*(1 + E)) / C)
Note that this last expression will work for all the values of R that
floorDiv(x + R, B) will work for, so we haven't really changed anything:
it is not yet an appproximation.

Now suppose that instead of: floor((x * G + x * G*E + R*G*(1 + E)) / C)
we use an approximation:     floor((x * G + K) / C)
where K is a constant chosen for a particular application.

We ask the question: given that constant K, for what values of x does:
[EQ.3] floor((x * G + K) / C) == floor((x * G + x * G*E + R*G*(1 + E)) / C)
                              == floor((x + R) / B) == floorDiv(x + R, B)

If we can find an R with Rmin <= R < Rmax such that:
[EQ.4] x * G + K == x * G + x * G*E + R*G*(1 + E)
then the approximation works for x.
(Note that the approximation may work for an x even if we cannot find such an R,
 so that the actual bounds for a particular case may be wider than the bounds
 estimated here. Such wider bounds can be found by testing.)

So to find if the approximation works for a particular x, 
we need to find an R with Rmin <= R < Rmax such that:
[EQ.5] K == x * G*E + R*G*(1 + E)

But using Jacobi's maxim "always invert" we instead ask the question:
Given K and given R with Rmin <= R < Rmax what value of x satisfies [EQ.5]?
That is, what value of x satisfies:
[EQ.6] x * G*E == K - R*G*(1 + E)
That is, what value of x satisfies:
[EQ.7] x == K/(G*E) - R*(1 + E)/E

By varying R from Rmin to Rmax we can generate all values of x 
that satisfy [EQ.7].

In other words, floor((x * G + K) / C) works for these values of x:
(1) If E > 0 it works for: K/(G*E) - Rmax*(1+E)/E <  X <= K/(G*E) - Rmin*(1+E)/E
(2) If E < 0 it works for: K/(G*E) - Rmin*(1+E)/E <= X <  K/(G*E) - Rmax*(1+E)/E

Similar ranges apply if we change "<=" to "<", etc.

TODO:

* Add proofs of ranges for where the approximation definitely doesn't work,
and at least one example of where the actual range for which x works
exceeds the bounds calculated above.

* Add proofs for floorDiv86400e9 and floorDiv86400e9Approximate.

Advertisements

About Colin Bartlett

I'm interested in arts, mathematics, science. Suliram is a partial conflation of the names of three good actors: Ira Aldridge, Anna May Wong, and another. My intention is to use a personal experience of arts to make some points, but without being too "me me me" about it. And to follow Strunk's Elements of Style. Except that I won't always "be definite": I prefer Niels Bohr's precept that you shouldn't write clearer than you think.
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s