Proposal: Better power for Rational

[Daniel Fischer]

On Saturday 25 September 2010 17:05:59, Isaac Dupree wrote:
> On 09/24/10 20:34, Daniel Fischer wrote:
> > Proposal: A better implementation of powers for Rational
>
> generally +1
>
> > For well-formed Rationals, the numerator and denominator are known to
> > be coprime, hence all powers of the numerator and denominator are also
> > coprime.
>
> Is it worth putting this stuff in the Data.Ratio code comments to
> explain why what you're doing is valid and useful, or is it already
> well-commented enough in a general sense about why "gcd" is sometimes
> necessary, yet expensive?
>

There are currently not many comments explaining such things there.
I guess adding a comment there why the special implementation for Rationals 
(and not other Ratio t) exists would be a good thing.

> > To avoid superfluous work, I propose a special power function for
> > Rationals and a rewrite rule to replace calls to (^) for Rational
> > bases by the special function. It might also be beneficial to export
> > the specialised function from Data.Ratio to be used if the rule
> > doesn't fire.
>
> Can you do the same for ^^ ? That is, a ratPowCanBeNegative (implement
> in terms of ratPow, or directly, as you please) and a RULE.

Sure, that'd be easy. I'm not sure whether it would make much difference.
The code for (^^) in GHC.Real is

-- | raise a number to an integral power
{-# SPECIALISE (^^) ::
        Rational -> Int -> Rational #-}
(^^)            :: (Fractional a, Integral b) => a -> b -> a
x ^^ n          =  if n >= 0 then x^n else recip (x^(negate n))

The SPECIALISE pragma is odd, btw.

Would a rule for (^^) be likely to fire in cases where the rule for (^) 
wouldn't?

-- testing

Yes, it would.

So some method of handling (^^) for Rational bases is called for.

Question for the experts, what would be more reliable,
do a <- [Int, Integer, Word, ... ]
   [{-# SPECIALISE (^^) :: Rational -> a -> Rational #-}]

or

{-# RULES
"^^/Rational"   (^^) = rationalPower
  #-}
rationalPower :: Integral a => Rational -> a -> Rational
rationalPower r e
    | e < 0     = ratPow (recip r) (-e)
    | otherwise = ratPow r e

?

Would specialising ratPow/rationalPower for common exponent types (Int, 
Integer) give additional benefit?

>
> (better names would be good if these are going to be exported though!

Aye. Unfortunately, I suck at finding good names. Suggestions welcome.

> But I don't think they need to be exported, unless hmm, is removing
> 'gcd' an *asymptotic* speedup for large integers?)
>

Yes. (Well, beaten by Felipe.)

> > Proposed function and rule:
> >
> > ratPow :: Integral a =>  Rational ->  a ->  Rational
> > ratPow _ e
> >
> >      | e<  0     = error "Negative exponent"
> >
> > ratPow _ 0  = 1 :% 1
> > ratPow r 1  = r
> > ratPow (0:%y) _ = 0 :% 1
> > ratPow (x:%1) e = (x^e) :% 1
> > ratPow (x:%y) e = (x^e) :% (y^e)
>
> Wondering why is there an extra case for x:%1 when the x:%y case handles
> that correctly (albeit slower)?

Well, I'm not sure whether that special case should be removed or a special 
case for numerator 1 should be added. It would require extensive 
benchmarking to be sure whether it's an actual improvement.

But perhaps
> Integer-base ^ does not have this
> 1-base optimization (maybe that's just because '1' maybe isn't
> multiplicative identity for general Num, and GHC.Real.^ is written for
> general Num base,

the special case should be added here.
Since x ^ 0 = 1, it seems to be assumed that 1 (or, fromInteger 1) is the 
multiplicative identity.
Then it might also make sense to special case base 0.

> or 1-base isn't common for general exponentiation but
> in Rationals it's common to have a Rational that's a whole number?),

Well, I don't know how common.
Of course, a lot of Rationals are created via fromInteger or fromIntegral, 
so n % 1 is probably overrepresented in code.
But if you calculate a lot of x^e for Rational x - and if you do only a 
handful of exponentiations, small performance differences don't matter 
anyway -, probably only a small fraction of those are whole numbers (or 
their reciprocals).
So we buy a performance gain for some exponentiations with a small extra 
cost of a test (== 1) for all exponentiations.
Whether one effect outweighs the other, and which one it would be, frankly, 
I have no idea.

> and you don't test for 1-base numerator.

Perhaps I should, but where do I get a large representative sample of 
Rationals appearing in real code?

> I think your choice is sensible
> enough overall; would like to hear what you think.
>
> > Like the elimination of `gcd` from recip (#4336), this would yield a
> > great performance boost.
>
> Did you measure the performance, or is it just obvious?  (Either is okay
> with me.)

Both. I've not done many tests, but all showed a 10× - 20× difference (for 
larger numbers).

>
> -Isaac

Cheers,
Daniel