[<prev] [<thread-prev] [day] [month] [year] [list]
Message-ID: <a0492425-fdb6-a15b-3034-5d25c72d8b76@esi.com.au>
Date: Sun, 8 Feb 2026 13:40:41 +1100 (AEDT)
From: Damian McGuckin <damianm@....com.au>
To: Szabolcs Nagy <nsz@...t70.net>
cc: MUSL <musl@...ts.openwall.com>
Subject: Re: Slightly faster 'scalbn(x,i) = x*2**i'
Subsequent to further testing on newer CPUs using more practical values of
the integer argument and a revision to avoid a potential integer overflow
bug that I had introduced, I can safely say that my code runs in much the
same time as MUSL's own code within timing tolerances.
Mine and MUSL's both compile with GCC to the same number of instructions
with GCC using -O3 and -mavx2 with the GCC 11 I am using.
So, no need to change the status quo.
Apologies for the noise - Damian
Below, I provide the code from the other project for those interested.
#include <stdint.h>
#include <assert.h>
#include <stdio.h>
#ifdef LINT
#define __builtin_expect(x, j) (x)
#define UFT 1.0e-38
#else
#define UFT 0x1.0p-1022
#endif
#define ABS fabs
#define unused /*@...sed@*/
#define likely(x) __builtin_expect((x),1)
#define unlikely(x) __builtin_expect((x),0)
#define asuint32(__f) ((union{float _f; uint32_t _i;}){__f})._i
#define asuint64(__f) ((union{double _f; uint64_t _i;}){__f})._i
#define asfloat(__u) ((union{uint32_t _i; float _f;}){__u})._f
#define asdouble(__u) ((union{uint64_t _i; double _f;}){__u})._f
/*
* Compute x*2^n wth near minimalist multiplications
*/
double myldexp(double x, int n)
{
static const uint64_t w = 64, p = 53;
static const uint64_t einf = (((uint64_t) 1) << (w - p)) - 1;
static const int b = (int) (einf / 2);
if (unlikely(n > b))
{
/*
* x*2^n = x * 2^*(b-n%1)*(2^(n/2))^2
* where n is limited to b<n<2*b+p+1
*/
static const double inflate[] = { 1 / UFT, 2 / UFT };
static const int N = (int) (einf + p);
const double y = x * inflate[n & 1];
const uint64_t et = (uint64_t) ((((n < N ? n : N) + (b + 1))) / 2);
const double t = asdouble(et << (p - 1));
return (y * t) * t; /* suboptimal if cost of multiplication is high */
}
else if (unlikely(n < 1 - b))
{
/*
* x*2^n = (x/2^(1-b-n))*2^(1-b) = (x/R)*2^(1-b), R = 2^m, m>0
*
* At this point, 1 <= R. If |x| <= 2^(-p), the products (x/R)*2^(1-b)
* and x*2^(1-b) both underflow to zero for any R, i.e. the former
* product can be computed as the latter. As a result, the computation
* of x/R only really needs to be done when x > 2^(-p). It can be shown
* that (x/R) = (x/r) where r = min(2^m,2^(e-1)) = 2^k for all practical
* purposes because once m>e-1, both underflow to zero. The use of a
* (normal) r rather than (potentially subnormal) R has the advantage
* that x/r = x/2^k can be computed rapidly by subtracting the biased
* exponent fields of x and r in-place, or because r has only zero bits
* in the negative and significand bit field, the encodings of x and r.
*/
const uint64_t _x = asuint64(x);
const uint64_t ex = (_x >> (p - 1)) & einf;
if (likely(b - p < ex && ex < einf)) /* only do this is necessary */
{
const uintt64_t e = ex - 1;
const uint64_t m = (uint64_t) ((1 - b) - n); /* never overflows */
const uint64_t _r = (m < e ? m : e) << (p - 1); /* always valid */
x = asdouble(_x - _r); /* replace x by x / r */
}
return x * UFT;
}
return x * asdouble(((uint64_t) (n + b)) << (p - 1));
}
#ifdef TEST
#include <stdio.h>
#include <math.h>
#include <assert.h>
int main(int argc, char *argv[])
{
double t = 0x1.9192d74f53405p-6;
union { double r; unsigned long _r; } _t = { t };
int e = -1027;
double z;
double big = 0x1.0p+1000;
double tiny = 0x1.0p-1020;
double res = myldexp (t, e);
/* double res = (t * 0x1p-1020) * 0x1p-7; */
double ref = t * 0x1p-1027;
printf ("res=%la\nref=%la\n", res, ref);
printf ("0x%016lx\n", _t._r);
_t.r = ref;
printf ("0x%016lx\n", _t._r);
/* test for same bug seen in Apple libm on ARM */
assert(res == ref);
assert(ldexp(big, -2060) == myldexp(big, -2060));
assert(__builtin_ldexp(tiny, 2000) == (z = myldexp(tiny, 2000)));
printf ("%1a\n", z);
return(0);
}
#endif
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
