[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [plt-scheme] Divisions



> 
> I predict that the next question from Paulo is how to get rid of the extra
> multiplication. :-)

Indeed... :) I'd like it to be faster...
Anyway, meanwhile I've implemented a version of it:

;Returns a list with three elements: 
;the gcd of a and b, 
;the multiplicative coeficient of (max a b) 
;and the multiplicative coeficient of (min a b)
(define ext-gcd
  (lambda (a b)
    (define ext-gcd-aux
      (lambda (a b q x1 x2 y1 y2)
        (let ((r (remainder a b))
              (id (floor (/ a b))))
          (if (zero? r)
              (list b x2 y2)
              (ext-gcd-aux b
                           r
                           id
                           x2 
                           (- x1 (* id x2))
                           y2 
                           (- y1 (* id y2)))))))
    (ext-gcd-aux (max a b) (min a b) 0 1 0 0 1)))

I get these result for bigger numbers:
> (time (gcd 2341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306234173284723941923061923741973249172394614765128937419237419732491723946147651289374234173284723941923061923741973249172394614765128937412341234 249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394))
cpu time: 0 real time: 2 gc time: 0
2

> (time (ext-gcd 2341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306234173284723941923061923741973249172394614765128937419237419732491723946147651289374234173284723941923061923741973249172394614765128937412341234 249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394))
cpu time: 840 real time: 837 gc time: 360
(2
 -41197154655173832587064954268624743092022135716832520557927007541172097204423073946484596355229020829444480605448311257646728435925159716050005045674243507999429188158810442586737506581272270680938206111665025452731472528937542081002243375535894793162505640858394864255183049357413457699045437450549580895716926260667997483128623671605759595274524313896951925534893931992537469778599
 387172625675389562580840348830899029571167923796088777060487379495035263506624857239058810939356454542673806923367141508267365940642301568498249241331563765015878882823006933367052477949223113000170543414488365791523414499418544773850680864916716438593558573703454787650807285998237649105627749723486941457868225846107392623651642772599817583127041489006709558854091060011515652825716898389697137950671179649098530290340272)

With that simple hack of a = b * q + r that I forgot about... I
get much better results using:
(define ext-gcd-ndiv
  (lambda (a b)
    (define ext-gcd-aux
      (lambda (a b q x1 x2 y1 y2)
        (let* ((id (quotient a b))
               (r (- a (* b id))))
          (if (zero? r)
              (list b x2 y2)
              (ext-gcd-aux b
                           r
                           id
                           x2 
                           (- x1 (* id x2))
                           y2 
                           (- y1 (* id y2)))))))
    (ext-gcd-aux (max a b) (min a b) 0 1 0 0 1)))

> (time (ext-gcd-ndiv 2341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394192306234173284723941923061923741973249172394614765128937419237419732491723946147651289374234173284723941923061923741973249172394614765128937412341234 249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394249172394614765128937423417328472394192306192374197324917239461476512893742341732847239419230619237419732491723946147651289374234173284723941923061923741973249172394614765128937423417328472394))
cpu time: 10 real time: 11 gc time: 0
(2
 -41197154655173832587064954268624743092022135716832520557927007541172097204423073946484596355229020829444480605448311257646728435925159716050005045674243507999429188158810442586737506581272270680938206111665025452731472528937542081002243375535894793162505640858394864255183049357413457699045437450549580895716926260667997483128623671605759595274524313896951925534893931992537469778599
 387172625675389562580840348830899029571167923796088777060487379495035263506624857239058810939356454542673806923367141508267365940642301568498249241331563765015878882823006933367052477949223113000170543414488365791523414499418544773850680864916716438593558573703454787650807285998237649105627749723486941457868225846107392623651642772599817583127041489006709558854091060011515652825716898389697137950671179649098530290340272)

:)

> 
> Since mzscheme uses the GMP library one could use
> 

Hummm, didn't know bout that... :) 
Seems that it is a pretty interesting detail... Always thought
that mzscheme implemented its own bignum library!

>     void mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t n, mpz_t d)
>     Divide n by d, forming a quotient q and/or remainder r.
> 
> which calculates the quotient and the remainder at the same time, to
> implement a quotient-remainder function.
> 
>     (let-values ([(q r) (quotient-remainder 42 5)]))
>        ...)
> 
> The file to hack must be "mzscheme/src/numarith.c".
> 

I'm going to check it. :)
hehehe, it seems fun!
Anyway, do you think that if it is hardcoded, it will get much
faster as the scheme implemened one?
Anyway, let's try... what will we lose anyway?
eheh

> If anyone has the energy (Paulo, Matthew?) one could add the extended gcd
> algorithm at the same time:
> 
>     http://www.swox.com/gmp/manual/Extended-GCD.html#Extended%20GCD
>

I know it... :)
I'm studying number theory and I'm implementing some number
theoretic functions and some crypt algorithms... it's extremely
fun! ehehe
 
> Yours,
> --
> Jens Axel Søgaard
> 

Best regards,

-- 
Paulo J. Matos : pocm(_at_)rnl.ist.utl.pt
Instituto Superior Tecnico - Lisbon    
Software & Computer Engineering - A.I.
 - > http://www.rnl.ist.utl.pt/~pocm 
 ---	
	Yes, God had a deadline...
		So, He wrote it all in Lisp!