[Agda] Data.Rational

Sergei Meshveliani mechvel at botik.ru
Wed Mar 12 13:37:35 CET 2014


On Tue, 2014-03-11 at 20:58 +0100, Andreas Abel wrote:

> I have no personal experience, but I listened to a talk by Anders 
> M"ortberg who presented his joint work with Cyril Cohen. He said that
> if you want to be efficient you should *not* try to cancel the common 
> primfactors in fractions each time.  The cost of making things 
> coprime repeatedly seems to be higher than computing with larger 
> integers. 
 

My experience is very different.

1. GHC Rational
---------------
As I recall, I tried this in earlier GHC versions (in 1990-ies)

   sumF n =  1 + 1/2 + 1/3 + ... + 1/n,   with taking large  n

-- something of this sort.
My DoCon is written in Haskell, and runs under GHC.
But its arithmetic for  Fraction Integer  differs from  Ratio Integer
of standard.
It applies the approach 

  (C)     -- compute with cancelling each appearing fraction to coprime.

Another approach is 

  (Delay) -- compute with delaying cancellation to the end

(I expect that any Computer Algebra developer, of course, chooses (C)).


And DoCon was somewhat 100 times faster on  sumF  than  Rational of
GHC.  
As I recall, this was precisely due to the (very strange) (Delay)
approach in old GHC.
I do not know, may be this remains so in the current GHC.

Anyway, GHC is great. And I do not care for this Rational, because DoCon
uses its own algorithms for the  Fraction domain constructor. 


2. Example with rational functions
----------------------------------
This experience (of about 1992) was as follows.
Find the inverse matrix 

    inverse M
    where
    M = [[ x+y , (x + 1)/(y^2 - x) ,  1/(2*x^2 - y) ,  y^2 + 3 ] , 
         ...,
         ...,
         ...
        ]

M is  4 by 4  matrix which elements are  fractions of polynomials  over
Integer coefficients in the variables [x, y]
-- the domain of  Fraction ℤ[x,y].

The three remaining rows in M to be filled with as _small expressions_
as the first one. If we do this at random, most probably M will occur
invertible. 
`inverse' is computed by the Gauss elimination -- by reducing M to
triangular. (This a contrived and simple solution, there is possible a
faster method). 

I tried several examples, and have observed that   
* running (Delay) costs more than 100 times more than (C),
* in (Delay) there appear intermediate denominators having thousands
monomials, like this:  x^40*y^35 + ... + 234*x + y  -- an expression
taking a couple of pages to print. 
* (C) computed in about 30 sec (on about 5 MHz machine),
  and (Delay) could not compute at all, because after 1 hour it was
going to fill all the memory with large polynomials.

(1) and (2) are typical examples from mathematical practice, they have
been taken at random: "take any popular computation method and see". 

May be you can specially choose an example where (Delay) is faster.



3. In our current example
-------------------------

we have  Fraction Integer  instead of  Fraction ℤ[x,y].

I expect the comparison result for   inverse M    will be similar.
Because  ℤ  is in _unary_ representation, and in (Delay), there will
appear denominators of thousands  `suc'  applications  -- similar as
there appeared large polynomials in [x,y].

n-ary  representation for ℤ will make the cost ratio smaller,
but I do not think that (Delay) will occur faster -- see the point (1)
in this letter.


4. CA books and developers
--------------------------

Look into the classic of   "The Art of Computer Programming" by D.Knuth.
Knuth is an expect in this area, he is a mathematician.
He considers many methods for each chosen problem.
But he does not even mention such a strange idea as (Delay) for
fractions.

Also ask, for example, the Axiom CA developers (the FriCAS e-mail list).
I am sure that they consider only the (C) approach to the arithmetic of
Fraction.

All right, I can mistake.
I only present here my experience. 

Anyway, for _this subject_ one needs just to consult 
   mathematicians who do algorithms for Computer Algebra libraries. 

Regards,

------
Sergei


 


> On 11.03.2014 17:57, Li Nuo wrote:
> > For the efficient problem, what you mentioned is right, but I think
> > maybe I didn't make it more precise how this approach works. In fact,
> > what we have done is to have both definitions, but for the operations
> > for the normal forms (Ra), we define them by lifting from the ones for
> > the non-reduced forms(Ra0) via using embedding functions from R -> R0
> > and reducing functions from Ra0 -> Ra. In practise it only helps you
> > defining functions and proving properties, and you can use the lifted
> > operations which will reduce the fractions whenever you do one operation
> > on Ra (because every lifted function just computes a non-reduced
> > fraction and reduce it).
> >
> > I hope this explains what I meant before.....
> >
> > Cheers,
> > Nuo
> >
> >
> >
> > 2014-03-11 15:29 GMT+00:00 Sergei Meshveliani <mechvel at botik.ru
> > <mailto:mechvel at botik.ru>>:
> >
> >     On Tue, 2014-03-11 at 14:28 +0100, Andreas Abel wrote:
> >      > Such a development would also be welcome for the standard library...
> >      >
> >      > On 11.03.2014 11:43, Li Nuo wrote:
> >      > > Hi,
> >      > >
> >      > > I think you can use an alternative approach to make the definition
> >      > > simpler. You could define the setoid version (namely fractions
> >     which
> >      > > doesn't have to be reduced) and use what we call "definable
> >     quotient
> >      > > structure" to relate the setoid definition with the normal form
> >      > > definition (which you can find here
> >      > > http://www.cs.nott.ac.uk/~txa/publ/defquotients.pdf). Generally
> >      > > speaking, it contains a set of operations and properties like
> >      > > normalisation functions (in this case, the fraction-reducing
> >     functions)
> >      > > and representative selection functions.
> >      > >
> >      > > The advantage of it is to lift functions in a more systematic way,
> >      > > benefiting from "quotients" without "quotient types". Here in
> >     this case,
> >      > > it reduces the complexity of checking coprime property all
> >     through which
> >      > > is unnecessary when doing operations on fractions(It is only
> >     needed if
> >      > > we want to check whether they are normal form fractions i.e.
> >     reduced
> >      > > fractions).
> >      > > [..]
> >
> >
> >      >> Here in this case,
> >      >> it reduces the complexity of checking coprime property all through
> >      >> which is unnecessary when doing operations on fractions
> >
> >
> >     I do not precisely follow the subject. But I would note for any
> >     occasion:
> >
> >     if one needs to compute efficiently with fractions, then a good
> >     tactic is to cancel each fraction to coprime as earlier as possible.
> >
> >     If one does not follow this rule, then computation will remain
> >     theoretically correct, but for most practicable examples it will become
> >     somewhat hundreds of times more expensive (the size of intermediate
> >     denominators in a loop will explode).
> >
> >     This feature does not depend on the choice a programming language or
> >     programming technique.
> >
> >     As to the below example of            ℚ- : ℚ → ℚ,
> >     its main part is a matter of proving
> >
> >       -- (1)
> >       lemma :  {m d : ℕ}  →  Coprime ∣ - (+ m) ∣ d  →  Coprime ∣ (+ m) ∣ d
> >
> >     And I think, it is by applying  PropositionalEquality.subst  for
> >     (\x -> Coprime x d)
> >
> >     with using the identity   ∣ - x ∣ ≡ ∣ x ∣   for  x : ℤ
> >
> >     (for absolute value). This latter identity can be taken from
> >     Data.Integer.
> >
> >     Probably, the lemmata like (1) need to be in Standard library.
> >
> >     Regards,
> >
> >     ------
> >     Sergei
> >
> >
> >      > >
> >      > >
> >      > > 2014-03-11 4:38 GMT+00:00 Amr Sabry <sabry at soic.indiana.edu
> >     <mailto:sabry at soic.indiana.edu>
> >      > > <mailto:sabry at soic.indiana.edu <mailto:sabry at soic.indiana.edu>>>:
> >      > >
> >      > >     Hi,
> >      > >
> >      > >     I need some decent Rational library and since the current
> >     Data.Rational
> >      > >     looks like it needs quite a bit of extension, I started
> >     writing a few
> >      > >     small functions. I am getting stuck very quickly though.
> >     Here is a small
> >      > >     example that should work:
> >      > >
> >      > >     -- unary negation
> >      > >     ℚ- : ℚ → ℚ
> >      > >     ℚ- p with ℚ.numerator p | ℚ.denominator-1 p | toWitness
> >     (ℚ.isCoprime p)
> >      > >     ... | -[1+ n ] | d | c = (+ ℕ.suc n ÷ ℕ.suc d) {fromWitness
> >     {!!}}
> >      > >     ... | + 0 | d | _ = p
> >      > >     ... | + ℕ.suc n | d | c = (-[1+ n ] ÷ ℕ.suc d) {fromWitness
> >     {!!}}
> >      > >
> >      > >     In the first line, c has type:
> >      > >
> >      > >        c  : {i : ℕ} → Σ (i ∣ ℕ.suc n) (λ x → i ∣ ℕ.suc d) → i ≡ 1
> >      > >
> >      > >     The goal is:
> >      > >
> >      > >        ?0 : Coprime (ℕ.suc n) (ℕ.suc d)
> >      > >
> >      > >     where Coprime is defined in Data.Nat.Coprimality as:
> >      > >
> >      > >        Coprime m n = ∀ {i} → (i ∣ m) × (i ∣ n) → (i ≡ 1)
> >      > >
> >      > >     As far as I can see, the type of the goal expands directly
> >     to the type
> >      > >     of c and in fact, "giving" c to the hole works but then the
> >     file doesn't
> >      > >     reload after that?
> >      > >
> >      > >     Any suggestions? --Amr
> >      > >
> >      > >     _______________________________________________
> >      > >     Agda mailing list
> >      > > Agda at lists.chalmers.se <mailto:Agda at lists.chalmers.se>
> >     <mailto:Agda at lists.chalmers.se <mailto:Agda at lists.chalmers.se>>
> >      > > https://lists.chalmers.se/mailman/listinfo/agda
> >      > >
> >      > >
> >      > >
> >      > >
> >      > > _______________________________________________
> >      > > Agda mailing list
> >      > > Agda at lists.chalmers.se <mailto:Agda at lists.chalmers.se>
> >      > > https://lists.chalmers.se/mailman/listinfo/agda
> >      > >
> >      >
> >      >
> >
> >
> >     _______________________________________________
> >     Agda mailing list
> >     Agda at lists.chalmers.se <mailto:Agda at lists.chalmers.se>
> >     https://lists.chalmers.se/mailman/listinfo/agda
> >
> >
> 
> 




More information about the Agda mailing list