How to improve pow()

[Bonita Montero]

I've written my own pow( double, int ) and it's about
10% slower than that of MSVCRT. Where's the mistake ?

double powPow( double d, int e )
{
unsigned eAbs = e >= 0 ? e : -e;
double v = 1.0, sq = d;
for( ; eAbs; eAbs >>= 1, sq *= sq )
if( eAbs & 1 )
v *= sq;
return e >= 0 ? v : 1.0 / v;
}

[Lynn McGuire]

You need to unroll the loop for the common cases of e == 1, e == 2, e ==
3, etc.

Lynn

[Bonita Montero]

LOL.

[Bonita Montero]

I've got it ! I've optimized away one jump and the performance is
now the same:

double powPow( double d, int e )
{

auto run = [&]<bool Neg>( bool_constant<Neg> ) -> double
{
unsigned eAbs;
if constexpr( !Neg )
eAbs = e;
else
eAbs = -e;

double v = 1.0, sq = d;
for( ; eAbs; eAbs >>= 1, sq *= sq )
if( eAbs & 1 )
v *= sq;

if constexpr( Neg )
v = 1.0 / v;
return v;
};
if( e >= 0 )
return run( false_type() );
else
return run( true_type() );
}

[Bonita Montero]

Now with loop-unrolling, but not the one Lynn suggested:

double powPow( double d, int e )
{
auto run = [&]<bool Neg>( bool_constant<Neg> ) -> double
{

auto unroll_n = [&]<typename Fn, size_t ... Indices>( Fn fn,
std::index_sequence<Indices ...> )
{
return (fn.template operator ()<Indices>() && ...);
};

unsigned eAbs;
if constexpr( !Neg )
eAbs = e;
else
eAbs = -e;
double v = 1.0, sq = d;

while( unroll_n(
[&]<size_t I>() -> bool
{
if( !eAbs )
return false;

if( eAbs & 1 )
v *= sq;

sq *= sq;
eAbs >>= 1;
return true;
}, std::make_index_sequence<5>() ) );

if constexpr( Neg )
v = 1.0 / v;
return v;
};
if( e >= 0 )
return run( false_type() );
else
return run( true_type() );
}

Now it's 13% faster than MSVCRT.

[Bonita Montero]

This is +110% faster than MSVCRT's implementation:

ouble powPow( double d, int e )
{

constexpr size_t UNROLL_FACTOR = 1;

auto run = [&]<bool Neg>( bool_constant<Neg> ) -> double
{

unsigned eAbs = !Neg ? e : -e;

double v = 1.0, sq = d;

for( ; eAbs; eAbs >>= 1, sq *= sq )

v *= (int)(eAbs & 1) * sq + (int)(eAbs & 1 ^ 1);
v = !Neg ? v : 1.0 / v;

return v;
};
if( e >= 0 )
return run( false_type() );
else
return run( true_type() );
}

There are no conditional moves for floating-point operations,
so I did the decision whether to add the squared value with
that: (int)(eAbs & 1) * sq + (int)(eAbs & 1 ^ 1).

[Alf P. Steinbach]

Maybe you could measure also my old code with slightly different
approach (no `constexpr`), at <url:
https://github.com/progrock-libraries/kickstart/blob/master/source/library/kickstart/main_library/core/ns%E2%96%B8language/operations/intpow.hpp>:


#include <assert.h>

#include <type_traits> // is_floating_point_v

// Important to not introduce possible future name conflicts with
<math.h>, hence
// the “lx” (short for “language extension”) ns▸
namespace kickstart::language::lx::_definitions {
using std::is_floating_point_v;

namespace impl
{
// Essentially this is Horner's rule adapted to calculating a
power, so that the
// number of floating point multiplications is at worst O(log₂n).
template< class Number_type >
constexpr inline auto intpow_( const Number_type base, const
int exponent )
-> Number_type
{
Number_type result = 1;
Number_type weight = base;
for( int n = exponent; n != 0; weight *= weight ) {
if( n % 2 != 0 ) {
result *= weight;
}
n /= 2;
}
return result;
}
} // namespace impl

template< class Number_type >
constexpr inline auto intpow( const Number_type base, const int
exponent )
-> Number_type
{
// TODO: proper failure handling
if( exponent < 0 and not is_floating_point_v<Number_type> ) {
assert( "For negative powers the base must be floating
point" and false );
}
return (0?0
: exponent > 0? impl::intpow_<Number_type>( base,
exponent )
: exponent == 0? 1
: Number_type(
1.0/impl::intpow_<Number_type>( base, -exponent ) )
);
}


//-----------------------------------------------------------
@exported:
namespace d = _definitions;
namespace exported_names { using
d::intpow;
} // namespace exported names
} // namespace kickstart::language::lx::_definitions

namespace kickstart::language::lx { using namespace
_definitions::exported_names; }


In the above copy+paste I removed an include of a header that defines a
`bool` replacement type called `Truth`, b/c I can't see that it's used.

The main difference from your code is that I discriminate the different
cases in a public function that just calls the internal implementation
with appropriate adjustments of arguments and result. Also, I use IMO
more clear notation like `/= 2` instead of `>>= 1`. Should yield same
machine code though.

I never measured this so it's interesting how it would stack up to your
code, and MSVC's.


- Alf

[Bonita Montero]

I think that it's very rare that someone uses an integer type for
pow( x, intType ). And with a double you have more precision within
the calculcaton itself and you'd just to cast back the result. So
ther's no need to have that generic.

[Juha Nieminen]

test.cc:3:18: error: 'size_t' does not name a type
test.cc:4:34: error: 'bool_constant' has not been declared
test.cc:14:28: error: 'false_type' was not declared in this scope
test.cc:16:28: error: 'true_type' was not declared in this scope

[Scott Lurndal]

Instead of insults, how about you simply disassemble the MS runtime
version of the function and figure out for yourself how it works?

A competent programmer would have no problem doing that.

[Scott Lurndal]

Typical unreadable code from Christoph/Bonnie.

[Bonita Montero]

There's nothing unreadable with that.
It's simply that you're both too stupid to adapt my style.

[Bonita Montero]

That's not necessary any more since I've developed some code
that is more than twice as fast as the MSVC implementation.

[Mut...@dastardlyhq.com]

Its quite clear you love complexity.

[Bonita Montero]

I've written complex code, but this absolutely isn't complex.
You're just overstrained.

[Cholo Lennon]

Don't feed the troll, she has serious mental issues. After reading her
answer I completely blocked her, she's not OK.

[Mut...@dastardlyhq.com]

On Wed, 7 Sep 2022 18:07:43 +0200

I'm not talking about the mathematics, I'm talking about the needless use of
a local lambda function. Why not do it inline? And wtf is the bool_constant
nonsense about with false_type and true_type? You don't need any of that.

[Juha Nieminen]

Mut...@dastardlyhq.com wrote:
> I'm not talking about the mathematics, I'm talking about the needless use of
> a local lambda function. Why not do it inline? And wtf is the bool_constant
> nonsense about with false_type and true_type? You don't need any of that.

Yeah. Normal people do a simpler run<true>() and run<false>().

Less writing, less confusing, simpler.

[Mut...@dastardlyhq.com]

You don't even need that. She's setting the bool based on whether e >= 0
then setting eAbs to e or -e based on the bool. Why not just do "eAbs = e"?
And since e isn't changed you can use it in the "v = " test at the bottom too.

And there's zero reason for a lambda and UNROLL_FACTOR does nothing. Could be
rewritten as follows:

double powPow( double d, int e )
{

double v = 1.0, sq = d;

for(unsigned aAbs=e; eAbs; eAbs >>= 1, sq *= sq )

v *= (int)(eAbs & 1) * sq + (int)(eAbs & 1 ^ 1);

return (e >= 0 ? v : 1.0 / v);
}

[Richard Damon]

Is that going to work for negative e?

if e == -1, then aAbs == UNSIGNED_MAX, and the return value will be
something like 1 / (d^UNSIGNED_MAX) instead of 1/d.

[Scott Lurndal]

Not to mention the unused constant UNROLL_FACTOR.

[Mut...@dastardlyhq.com]

I didn't compile and run it , but something along those lines will do the
exact same job as Bonitas code and probably more efficiently.

[Bonita Montero]

> a local lambda function. ...

That's not needles, that's my style which saves a lot of code.

[Bonita Montero]

There's nothing confusing with my code and the way I did it saves
one non-predictible branch and results in half the code size if
you do that manually.

[Bonita Montero]

Am 08.09.2022 um 15:25 schrieb Scott Lurndal:

> Not to mention the unused constant UNROLL_FACTOR.

That's a line of code I forgot to remove from by before version.

[Bonita Montero]

I don't understand why that is trolling.
But I think people being provoked by that
have like you have mental issues, not me.

[Juha Nieminen]

Mut...@dastardlyhq.com wrote:
>>Is that going to work for negative e?
>>
>>if e == -1, then aAbs == UNSIGNED_MAX, and the return value will be
>>something like 1 / (d^UNSIGNED_MAX) instead of 1/d.
>
> I didn't compile and run it , but something along those lines will do the
> exact same job as Bonitas code and probably more efficiently.

If you need the absolute value of an integer, I don't think there's any
way around having to use a conditional.

Or, better yet, just std::abs(). (The compiler will do a better job than
you at optimizing it. Eg. if there's a CPU instruction that directly
converts a register into its absolute value without conditionals, certainly
the compiler will use that.)

[Öö Tiib]

He just "optimised" sign check out and posted it as "simpler". But unfortunately
it does not matter how simply and quickly we get wrong answers.

[Malcolm McLean]

On Friday, 9 September 2022 at 06:29:23 UTC+1, Juha Nieminen wrote:
> Mut...@dastardlyhq.com wrote:
> >>Is that going to work for negative e?
> >>
> >>if e == -1, then aAbs == UNSIGNED_MAX, and the return value will be
> >>something like 1 / (d^UNSIGNED_MAX) instead of 1/d.
> >
> > I didn't compile and run it , but something along those lines will do the
> > exact same job as Bonitas code and probably more efficiently.
> If you need the absolute value of an integer, I don't think there's any
> way around having to use a conditional.
>

(int) sqrt(N*N);

[Mut...@dastardlyhq.com]

There's nothing clever or simple about writing code that has pointless
syntatic contructs and calls that don't need to be there. If you worked in
a team with other people you'd understand than when you have to debug/modify
someone elses code you don't want to have to wade through a load of ego
tripping show off code first.

[Öö Tiib]

On Friday, 9 September 2022 at 12:16:36 UTC+3, Mut...@dastardlyhq.com wrote:
> On Fri, 9 Sep 2022 00:20:33 -0700 (PDT)
> =?UTF-8?B?w5bDtiBUaWli?= <oot...@hot.ee> wrote:
> >On Friday, 9 September 2022 at 08:29:23 UTC+3, Juha Nieminen wrote:
> >> Mut...@dastardlyhq.com wrote:
> >> >>Is that going to work for negative e?
> >> >>
> >> >>if e == -1, then aAbs == UNSIGNED_MAX, and the return value will be
> >> >>something like 1 / (d^UNSIGNED_MAX) instead of 1/d.
> >> >
> >> > I didn't compile and run it , but something along those lines will do the
> >> > exact same job as Bonitas code and probably more efficiently.
> >> If you need the absolute value of an integer, I don't think there's any
> >> way around having to use a conditional.
> >>
> >> Or, better yet, just std::abs(). (The compiler will do a better job than
> >> you at optimizing it. Eg. if there's a CPU instruction that directly
> >> converts a register into its absolute value without conditionals, certainly
> >> the compiler will use that.)
> >
> >He just "optimised" sign check out and posted it as "simpler". But
> >unfortunately
> >it does not matter how simply and quickly we get wrong answers.
>
> There's nothing clever or simple about writing code that has pointless
> syntatic contructs and calls that don't need to be there.

That is nothing I would dispute but irrelevant about your code snippet. Sign
check is not "syntatic contruct" (sic) but necessity in given algorithm.

> If you worked in
> a team with other people you'd understand than when you have to debug/modify
> someone elses code you don't want to have to wade through a load of ego
> tripping show off code first.

Sign check is also not "ego tripping show off code" but necessity in given
algorithm.

[Mut...@dastardlyhq.com]

On Fri, 9 Sep 2022 02:39:33 -0700 (PDT)

I never said my code was ready to run, I simply removed the pointless lambda
and boolean that didn't need to be there.

[Öö Tiib]

Not true ... your major complaint was about sign check: "She's setting the bool

based on whether e >= 0 then setting eAbs to e or -e based on the bool. Why
not just do "eAbs = e"?"

That sign check is needed for the algorithm to work. Explain to your rubber
duck <https://en.wikipedia.org/wiki/Rubber_duck_debugging> that it is
unneeded and see yourself how you two reach conclusion that it is needed.

[Mut...@dastardlyhq.com]

Clearly you're looking for an argument. I'm not interested, find another foil.

[Öö Tiib]

I am totally uninterested in argument ... I was just pointing out that you were
incorrect. You want to argue because you are incapable to admit that you
were incorrect. But describing that sign check as "syntatic contruct" and
"ego tripping show off code" just digs you deeper. Indeed, better stop.

[Mr Flibble]

On Tue, 6 Sep 2022 19:08:32 +0200
Bonita Montero <Bonita....@gmail.com> wrote:

> I've written my own pow( double, int ) and it's about
> 10% slower than that of MSVCRT. Where's the mistake ?
>
> double powPow( double d, int e )
> {
> unsigned eAbs = e >= 0 ? e : -e;
> double v = 1.0, sq = d;
> for( ; eAbs; eAbs >>= 1, sq *= sq )
> if( eAbs & 1 )
> v *= sq;
> return e >= 0 ? v : 1.0 / v;
> }

You can improve it in two ways:

1) parameter e should be a double not an int
2) use logarithms and exponents

/Flibble

[David Brown]

(((x >> 31) << 1) + 1) * x

That assumes 32-bit two's complement integers and common implementation
dependent behaviour of bit shifts. I haven't checked corner cases such
as INT_MIN.

I've no idea if it is faster or not. And of course, a smart enough
compiler can turn it into a conditional if it thinks that is faster.

[Mut...@dastardlyhq.com]

On Fri, 9 Sep 2022 04:05:46 -0700 (PDT)

Get together with Bonita and write something. I'm sure we could all do with
a good laugh.

[Mut...@dastardlyhq.com]

Possibly. I wasn't commenting on the maths which I'm sure is correct. It was
all the surrounding unnecessary syntax I object to mainly because I regularly
see that sort of showboating in code and its irritating when you have to
debug it. Either the author is showing off or they can't think clearly, there
is no other option.

[Alf P. Steinbach]

The point of Bonita's code is to compute an exact integer power when
that can be represented as a `double`, and anyway to do it in a real
fastish way, with O(log n) multiplications where n is the exponent.

`double` exponent value and computation via logarithm and exponentiation
would yield a generally inexact result, and I believe also in a slower
way, even when logarithm and exponentiation are hardware operations.

That belief is strengthened by Bonita's measurements that show distinct
improvement, more than twice as fast, compared to Microsoft's
implementation.


- Alf

[Bonita Montero]

There's nothing wrong with my syntax - it's just a differnt style.
I used the templated lambda f.e. to prevent to have duplicate code
which I'd have to write to reduce two non-predictible branches to one.
If you've adapted this style the code becomes more readable than the
manual optimizations.

[Bonita Montero]

The implementations are likely to be nearly the same except from my
latest optimitation to prevent the repeated branch. This is because
the performance of my first implementation nearly matches the per-
formance of the Visual C++ runtime - this is unlikely to happen by
accident.

[Bonita Montero]

There's nothing to laugh with my coding style. It's simply somewhat
different because I use a lot of functional programming, mainly to
reduce to size of the code.

[Juha Nieminen]

David Brown <david...@hesbynett.no> wrote:
> (((x >> 31) << 1) + 1) * x

That's clever.

> I've no idea if it is faster or not.

I suppose it depends on how many times the line of code is executed and
how often the CPU mispredicts the branch in the conditional version.

[David Brown]

More obvious source code that has a conditional is usually compiled as a
conditional move, rather than a branch:

x < 0 ? -x : x

I don't think conditional moves are predicted in the same way. But the
efficiency is all going to depend on pipelining, superscaling,
vectorising (if possible), and so on. Modern x86 processors can do
register renaming and should allow code like mine to pass mix well with
other code, but I have no idea how well condition codes / flags and
conditional moves can be mixed with the flow of instructions.

[Mr Flibble]

On Fri, 9 Sep 2022 18:27:48 +0200
"Alf P. Steinbach" <alf.p.s...@gmail.com> wrote:

> On 9 Sep 2022 15:56, Mr Flibble wrote:
> > On Tue, 6 Sep 2022 19:08:32 +0200
> > Bonita Montero <Bonita....@gmail.com> wrote:
> >
> >> I've written my own pow( double, int ) and it's about
> >> 10% slower than that of MSVCRT. Where's the mistake ?
> >>
> >> double powPow( double d, int e )
> >> {
> >> unsigned eAbs = e >= 0 ? e : -e;
> >> double v = 1.0, sq = d;
> >> for( ; eAbs; eAbs >>= 1, sq *= sq )
> >> if( eAbs & 1 )
> >> v *= sq;
> >> return e >= 0 ? v : 1.0 / v;
> >> }
> >
> > You can improve it in two ways:
> >
> > 1) parameter e should be a double not an int
> > 2) use logarithms and exponents
>
> The point of Bonita's code is to compute an exact integer power when
> that can be represented as a `double`, and anyway to do it in a real
> fastish way, with O(log n) multiplications where n is the exponent.

std::pow taking an integer parameter was dropped in C++11: ask yourself
why.

>
> `double` exponent value and computation via logarithm and
> exponentiation would yield a generally inexact result, and I believe
> also in a slower way, even when logarithm and exponentiation are
> hardware operations.

Multiplication and logarithm/exponentiation are both equally inexact
when we are dealing with IEEE 754 floating point.

[Bonita Montero]

Am 09.09.2022 um 19:39 schrieb Mr Flibble:

> std::pow taking an integer parameter was dropped in C++11: ask yourself
> why.

That's because the performance lagely is dependent on span of the
highest and lowest set bit of the exponent. This doesn't make much
difference if you have a integer or a floating point exponent for
the same value.

[Bonita Montero]

Am 09.09.2022 um 19:17 schrieb David Brown:
> On 09/09/2022 18:54, Juha Nieminen wrote:
>> David Brown <david...@hesbynett.no> wrote:
>>> (((x >> 31) << 1) + 1) * x
>>
>> That's clever.
>>
>>> I've no idea if it is faster or not.
>>
>> I suppose it depends on how many times the line of code is executed and
>> how often the CPU mispredicts the branch in the conditional version.
>
> More obvious source code that has a conditional is usually compiled as a
> conditional move, rather than a branch:
>
>     x < 0 ? -x : x
>

> I don't think conditional moves are predicted in the same way. ...

Conditional moves don't affect the flow in the pipeline and
there are no mis-predictions.

[Mut...@dastardlyhq.com]

On Fri, 9 Sep 2022 18:28:52 +0200

It really doesn't.

[Bonita Montero]

Then take my code, write the lambda expanded twice and think about
if that is more readable.

[Alf P. Steinbach]

On 9 Sep 2022 19:39, Mr Flibble wrote:
>
> Multiplication and logarithm/exponentiation are both equally inexact
> when we are dealing with IEEE 754 floating point.

No no no.

- Alf

[Mut...@dastardlyhq.com]

On Sat, 10 Sep 2022 11:53:31 +0200

face <- palm.

Never mind.

[Mr Flibble]

On Sat, 10 Sep 2022 11:59:38 +0200

"Alf P. Steinbach" <alf.p.s...@gmail.com> wrote:

Yes.

/Flibble

[Mr Flibble]

On Sat, 10 Sep 2022 11:59:38 +0200

"Alf P. Steinbach" <alf.p.s...@gmail.com> wrote:

Yes.

#include <iostream>
#include <cmath>

int main()
{
double const n = 0.123456789;
std::cout << (n * n * n) << std::endl;
std::cout << std::pow(n, 3.0) << std::endl;
std::cout << std::exp(3.0 * std::log(n)) << std::endl;
}

0.00188168
0.00188168
0.00188168

/Flibble

[Mr Flibble]

On Sat, 10 Sep 2022 11:59:38 +0200

"Alf P. Steinbach" <alf.p.s...@gmail.com> wrote:

Yes.

#include <iostream>
#include <cmath>

#include <limits>
#include <iomanip>

int main()
{
double const n = 0.123456789;
std::cout <<

std::setprecision(std::numeric_limits<double>::digits10 + 1);

std::cout << (n * n * n) << std::endl;
std::cout << std::pow(n, 3.0) << std::endl;
std::cout << std::exp(3.0 * std::log(n)) << std::endl;
}

0.001881676371789155
0.001881676371789155
0.001881676371789155

/Flibble

[Bonita Montero]

Pow uses a row of multiplications of d ^ (2 ^ n) for each bit in the
exponent. With that you have only a logarithmic number of
multipli-cations. When you multiply N times in a row that might be less
exact.

#include <iostream>

int main()
{
double v = 1.5, r = v;
for( int i = 0; i != 30; ++i, r *= 1.5 )
std::cout << hexfloat << pow( 1.5, i ) - r << std::endl;
}

[Alf P. Steinbach]

You're a good programmer but you're lousy in maths :).

Consider:


----------------------------------------------------------------------
#include <float.h> // DBL_DIG
#include <math.h> // exp, log
#include <stdio.h> // printf

auto main() -> int
{
const double x = 3.0;
printf( "Multiplication yields 3^2 = %.*f.\n", DBL_DIG + 1, x*x );
printf( "exp/log yields 3^2 = %.*f.\n", DBL_DIG + 1, exp(
2*log( x ) ) );
}
----------------------------------------------------------------------


Output with Visual C++ 2020:


Multiplication yields 3^2 = 9.0000000000000000.
exp/log yields 3^2 = 9.0000000000000018.


- Alf

[Mr Flibble]

On Sat, 10 Sep 2022 17:33:40 +0200

I may or may not be lousy at maths however you are completely
missing the point: we use double when we want to deal with real
numbers, if we want to deal with just integers then we use integer
types. IN GENERAL multiplication is just as inexect as exp/log when
using IEEE 754 floating point for the set of real numbers representable
by the floating point type.

Not that it materially effects anything, I note that x*x is likely
evaluated at compile time by the optimiser rather than at runtime.

/Flibble

[Mr Flibble]

On Sat, 10 Sep 2022 17:33:40 +0200

Consider this:

#include <float.h> // DBL_DIG
#include <math.h> // exp, log
#include <stdio.h> // printf

auto main() -> int
{
const double x = 3.0;

printf( "Multiplication yields 3^3 = %.*f.\n", DBL_DIG + 1, x*x*x
); printf( "exp/log yields 3^3 = %.*f.\n", DBL_DIG + 1, exp(
3*log( x ) ) );
}

Output:

Multiplication yields 3^3 = 27.0000000000000000.
exp/log yields 3^3 = 27.0000000000000000.

Feel free to retract your assertion that I am lousy at maths.

/Flibble

[Alf P. Steinbach]

Well let's just say that math is not what you're absolutely best at.

#include <float.h> // DBL_DIG
#include <math.h> // exp, log
#include <stdio.h> // printf

auto main() -> int
{
const double x = 3.0;

printf( "Multiplication yields 3^4 = %.*f.\n", DBL_DIG + 1, x*x*x*x );
printf( "exp/log yields 3^4 = %.*f.\n", DBL_DIG + 1, exp(
4*log( x ) ) );
}


Output:


Multiplication yields 3^4 = 81.0000000000000000.
exp/log yields 3^4 = 81.0000000000000284.


The tiny little difference can be arbitrarily amplified by e.g. rounding
up this result, or rounding down one that errs the opposite way. And
that can matter significantly. Often it will not matter, but it can.

And although Bonita has stated that it was not her intention to use her
code for integer base it is eminently usable for that, and that was the
main use case for my `constexpr` templated version `intpow` years ago.

I thought it was nice to be able to write e.g.

constexpr int k = intpow( 2, 10 );

Arguably a missing part of the standard library.


- Alf