[Agda] summing rationals

[Sergei Meshveliani]

On Sun, 2017-06-25 at 09:46 +0200, Ulf Norell wrote:
> The rationals are not very well-developed in agda-prelude.
> Nevertheless I get unsurprising performance compared to ghc (~2min for
> n = 20,000 on my machine) with the following benchmark:
> [..]
 
What is the machine that runs your test? 
(how much frequency [GHz] has it?).
 
What timing GHC shows on your machine for  n = 200,000  for

     import Data.Ratio ((%))
        
     f :: Integer -> Rational     -- :: Ratio Integer
     f n =
         sum [1 % k | k <- [1 .. n]]

     > ghc --make -O Test
?         


DoCon-A (Agda) has  21 sec  for  n = 2000,  for a  3 GHz  personal
computer.

 

> The reason that DoCon performs so well is that it has a more clever
> implementation of _+_. agda-prelude and GHC implements _+_ as
> 
> 
>   (p :% q) + (p1 :% q1) = (p * q1 + p1 * q) % (q * q1)
>     where a % b = let d = gcd a b in div a d :% div d b
>
> DoCon implements _+_ as
> 
> 
>   (p :% q) + (p1 :% q1) = div n dd :% (q' * q1' * div d dd)
>     where
>       d = gcd q q1
>       q' = div q d
>       q1' = div q1 d
>       n = p * q1' + p1 * q'
>       dd = gcd n d
> 
> 
> I can't claim to understand why this is so much faster, but it is.
> 

As I recall, I just copied there the method from the old book 
"The art of computer programming" by D.Knuth.  
As I understand, the main idea is that for large numbers multiplication
causes the result size growth (number of (macro)digits), while divMod
and gcd do not. 
And the sizes of  (p * q1), (p1 * q), (q * q1)  effect the cost of the
further operations applied to obtain the final fraction.
Hence it is better to cancel p, p1, q, q1 by the corresponding gcd-s as
soon as possible.

But for _random_ numbers a naive method of GHC is, probably, faster on
average than Knuth's. Because for two random numbers m, n, it often
holds  gcd m n = 1.

But in Harmonic, and in many real mathematical tasks, fractions are
often accumulated in the middle of loops, they are "less random". 

--------------
I do not find, so far, why DoCon-A shows a greater cost order than GHC.

May be this is _-fr_ that probably imply _∸_ ?  May be, _∸_ or `pred' is
slow, I do not know.

Thanks for comments,

------
Sergei





> On Sat, Jun 24, 2017 at 9:06 PM, Sergei Meshveliani <mechvel at botik.ru>
> wrote:
>         Dear all,
>         
>         Consider the following performance test for rational
>         arithmetic.
>         Call this test call "Harmonic":
>         
>           f n =  sum [1/k | k <- [1 .. n]]   for n : Nat, 1/k :
>         Rational,
>         
>           sum = foldr _+_ 0    where _+_ = fraction addition
>         
>           test n =  let fn = f n
>                     in  fn - fn
>         
>           Read from file   n = 100, 1000, 2000, 100.000 ...,
>         
>           print  (test n)  and check time.
>         
>         (fn - fn)  is applied here in order to avoid printing large
>         values.
>         
>         
>           --------------------------- Haskell code -------
>           import Data.Ratio ((%))
>         
>           f :: Integer -> Rational     -- :: Ratio Integer
>           f n =
>               sum [1 % k | k <- [1 .. n]]
>           -------------------------------------------------
>         
>         
>         I compare the three libraries:
>         
>           (GHC) GHC  (Glasgow Haskell 7.10.2),
>         
>           (D)   DoCon-2.12  written in Haskell, running under GHC,
>         
>           (DA)  DoCon-A 2.00-pre   written in Agda,
>                                    running under Agda 2.5.2 -like.
>         
>         (D) and (DA) keep fractions with numerator and denominator
>         mutually
>         coprime, and cancel fractions most eagerly by gcd.
>         
>         (GHC)  does not cancel by gcd so eagerly.
>         
>         The performance comparison on Harmonic is:
>         
>                         (D) >> (GHC) > (DA).
>         
>         (D) takes  n = 200.000  easily, while GHC "hangs" on 20.000
>             -- for a 3 GHz machine.
>         
>         (GHC) is considerably faster (DA)
>               (DA hangs for long on 4000).
>         
>         The method in (DA) is the same as (D).
>         gcd is almost as fast as GHC.gcd, this has been tested.
>         Division (with remainder) is via  dimMod  of Agda-prelude, it
>         is fast.
>         
>         And I wonder why (DA) is so slow.
>         (DA) carries all the needed proofs (for coprimality, and
>         such).
>         But I hope, this does not effect performance any essentially.
>         
>         Can anybody provide the code for testing Harmonic for
>         Rational numbers
>         in  agda-prelude ?
>         
>         The (DA) test for (Fraction Integer) is like this:
>         
>         ------------------------------------------------------------------------
>         ...
>         open Fraction.F0 Int1.integralRing-with∣? Int1.gcdRing using
>                                          (Fraction; fraction;
>         fraction-showInst)
>         open Fraction.F1 Int1.integralRing-with∣? Int1.gcdRing using
>         (frRingoid)
>         open Ringoid frRingoid using (_-_)
>         open Fraction
>         
>         inv-n :  ℕ → Fraction              -- \n → 1/(1+n)
>         inv-n n =
>                 fraction (+ 1) (+ (suc n)) +suc≢+0
>         
>         f : ℕ → Fraction    -- sum [1/n | n <- [1 .. n+1]]
>         f n =
>             sum frRingoid (map inv-n (natSegment 0 n))
>         
>         main : IO Unit
>         main = (readFiniteFile "data.txt") >>= putStrLn ∘ toCostring ∘
>         g
>           where
>           g : String -> String
>           g str =
>             -- n  is first read from file  data.txt
>         
>             concatStr
>             ("n = "                 ∷  showNat n                 ∷
>         ",   " ∷
>              "num (f n - f n) =  "  ∷  showInt (num (fn - fn))   ∷ "\n
>         \n"
>              ∷ []
>             )
>             where
>             n  = stringToNat $ filterStr decimalDigit? str
>             fn = f n
>         -------------------------------------------------------------------------
>         
>         
>         So far, I do not find the performance bottleneck in (DA).
>         
>         Thanks,
>         
>         ------
>         Sergei
>         
>         
>         
>         _______________________________________________
>         Agda mailing list
>         Agda at lists.chalmers.se
>         https://lists.chalmers.se/mailman/listinfo/agda
> 
>