Mantis Bug Tracker

View Issue Details Jump to Notes ] Issue History ] Print ]
IDProjectCategoryView StatusDate SubmittedLast Update
0003473OCamlOCaml generalpublic2005-02-18 15:472009-04-29 20:41
Reporteradministrator 
Assigned Toweis 
PrioritynormalSeverityfeatureReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version3.11.0+beta 
Summary0003473: [WISH] Better implementation of Num.quo_num
DescriptionHi,

Attached is a message sent to the list telling that a better
implementation of Num.quo_num is desirable. I forward it to the bug
tracker so that you do not forget about it when you have some time. :)

Best regards,
ChriS

----------------
Return-path: <Christophe.Troestler@umh.ac.be>
Envelope-to: trch@localhost
Received: from poincare
    ([127.0.0.1] helo=localhost ident=trch)
    by poincare with esmtp (Exim 3.36 #1 (Debian))
    id 1D0PD7-0004TP-00
    for <trch@localhost>; Sun, 13 Feb 2005 20:15:13 +0100
Received: from pop.tiscali.be [213.205.33.245]
    by localhost with POP3 (fetchmail-6.2.5)
    for trch@localhost (single-drop); Sun, 13 Feb 2005 20:15:13 +0100 (CET)
Received: from mail-in.tiscali.be (62.235.13.184) by mail-8-bnl.tiscali.it (7.2.052)
        id 41DE61B50032B6BD for del-con@tiscali.be; Sun, 13 Feb 2005 20:05:37 +0100
Received: from [193.190.193.73] (helo=hedwig1.umh.ac.be)
    by mail-in.tiscali.be with esmtp (Tiscali.be http://www.tiscali.be [^])
    id 1D0P3n-0005Bp-9t
    for <del-con@tiscali.be>; Sun, 13 Feb 2005 20:05:37 +0100
Received: from swip.net (mailfe06.swip.net [212.247.154.161])
    by hedwig1.umh.ac.be (8.13.1/8.13.1) with ESMTP id j1DJ8rg0770238
    for <Christophe.Troestler@umh.ac.be>; Sun, 13 Feb 2005 20:08:57 +0100
X-T2-Posting-ID: 2IIuXFmkcTpj3lKEFKW25A==
Received: from [83.176.160.123] (HELO poincare)
  by mailfe06.swip.net (CommuniGate Pro SMTP 4.2.9)
  with ESMTP id 292824364; Sun, 13 Feb 2005 20:05:07 +0100
Received: from poincare
    ([127.0.0.1] helo=localhost ident=trch)
    by poincare with esmtp (Exim 3.36 #1 (Debian))
    id 1D0OPi-00042V-00; Sun, 13 Feb 2005 19:24:10 +0100
Date: Sun, 13 Feb 2005 19:24:10 +0100 (CET)
Message-Id: <20050213.192410.61034584.Christophe.Troestler@umh.ac.be>
To: "O'Caml Mailing List" <caml-list@inria.fr>
From: Christophe TROESTLER <Christophe.Troestler@umh.ac.be>
Organization: Universite de Mons-Hainaut (http://math.umh.ac.be/an/ [^])
X-Spook: AIMSX JUWTF enforcers president [Hello to all my friends and fans
 in domestic surveillance] ANDVT JPL advisors 64 Vauxhall Cross Soviet
X-Blessing: Om Ah Hum Vajra Guru Pema Siddhi Hum
X-Operating-System: GNU/Linux (http://www.linux.org/ [^])
X-Mailer-URL: http://www.mew.org/ [^]
X-Mailer: Mew version 4.2 on Emacs 21.3 / Mule 5.0 (SAKAKI)
Mime-Version: 1.0
X-Scanned-By: MIMEDefang 2.1 (www dot roaringpenguin dot com slash mimedefang)
Subject: Num.quo_num
Content-Type: Multipart/Mixed;
 boundary="--Next_Part(Sun_Feb_13_19_24_10_2005_886)--"
X-Spam-Status: No, hits=1.1 required=5.0 tests=SARE_HEAD_HDR_XT2PID,
    SARE_RECV_SPAM_DOMN81

----Next_Part(Sun_Feb_13_19_24_10_2005_886)--
Content-Type: Text/Plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Hi,

Attached are two implementations of the algorithm for finding the
digits of pi described here:
http://web.comlab.ox.ac.uk/oucl/work/jeremy.gibbons/publications/spigot.pdf [^]
The two implementations differ only by the fact that the first one
uses Big_int while the second one uses Num. However, the efficiency
figures are rather different (here for 100 digits) :

Big_int
  real 0m0.009s
  user 0m0.007s
  sys 0m0.003s

Num
  real 0m0.487s
  user 0m0.484s
  sys 0m0.000s

Upon investigation, it appears that the culpright is the following:

  let quo_num x y = floor_num (div_num x y)

Indeed "Num.div_num" calls "Num.num_of_ratio" wich in turn calls
"Ratio.normalize_ratio" which I guess accounts for the huge time
difference. Now, "Num.normalize_ratio" is actually useless for
"quo_num", since it does "num_of_big_int" anyway. What about writing
a specialized and vastly more efficient implementation for "quo_num"?

Regards,
ChriS

----Next_Part(Sun_Feb_13_19_24_10_2005_886)--
Content-Type: Text/Plain; charset=iso-8859-1
Content-Disposition: inline; filename="pidigits.ml"
Content-Transfer-Encoding: quoted-printable
X-MIME-Autoconverted: from 8bit to quoted-printable by hedwig1.umh.ac.be id j1DJ8rg0770238

(* Pi digits computed with the sreaming algorithm given on pages 4, 6
   & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
   Gibbons, August 2004. *)

open Printf
open Big_int

let ( !$ ) =3D Big_int.big_int_of_int
let ( +$ ) =3D Big_int.add_big_int
let ( *$ ) =3D Big_int.mult_big_int
let ( =3D$ ) =3D Big_int.eq_big_int
let zero =3D Big_int.zero_big_int
and one =3D Big_int.unit_big_int
and three =3D !$ 3
and four =3D !$ 4
and ten =3D !$ 10
and neg_ten =3D !$(-10)

(* Linear Fractional (aka M=F6bius) Transformations *)
module LFT =3D
struct
  let floor_ev (q,r,s,t) x =3D div_big_int (q *$ x +$ r) (s *$ x +$ t)

  let unit =3D (one, zero, zero, one)

  let comp (q,r,s,t) (q',r',s',t') =3D
    (q *$ q' +$ r *$ s', q *$ r' +$ r *$ t',
     s *$ q' +$ t *$ s', s *$ r' +$ t *$ t')
end

let next z =3D LFT.floor_ev z three
and safe z n =3D (n =3D$ LFT.floor_ev z four)
and prod z n =3D LFT.comp (ten, neg_ten *$ n, zero, one) z
and cons z k =3D
  let den =3D 2*k+1 in LFT.comp z (!$ k, !$(2*den), zero, !$ den)

let rec digit k z n row col =3D
  if n > 0 then
    let y =3D next z in
    if safe z y then
      if col =3D 10 then (
    let row =3D row + 10 in
    printf "\t:%i\n%s" row (string_of_big_int y);
    digit k (prod z y) (n-1) row 1
      )
      else (
    print_string(string_of_big_int y);
    digit k (prod z y) (n-1) row (col+1)
      )
    else digit (k+1) (cons z k) n row col
  else
    printf "%*s\t:%i\n" (10 - col) "" (row + col)

let digits n =3D digit 1 LFT.unit n 0 0

let () =3D
  digits(int_of_string Sys.argv.(1))

----Next_Part(Sun_Feb_13_19_24_10_2005_886)--
Content-Type: Text/Plain; charset=us-ascii
Content-Disposition: inline; filename="pidigits-num.ml"
Content-Transfer-Encoding: 7bit

(* Pi digits computed with the sreaming algorithm given on pages 4, 6
   & 7 of "Unbounded Spigot Algorithms for the Digits of Pi", Jeremy
   Gibbons, August 2004. *)

open Printf
open Num

let zero = num_of_int 0
and one = num_of_int 1
and three = num_of_int 3
and four = num_of_int 4
and ten = num_of_int 10
and neg_ten = num_of_int(-10)

(* Linear Fractional Transformation *)
module LFT =
struct
  let floor_ev (q,r,s,t) x = quo_num (q */ x +/ r) (s */ x +/ t)

  let unit = (one, zero, zero, one)

  let comp (q,r,s,t) (q',r',s',t') =
    (q */ q' +/ r */ s', q */ r' +/ r */ t',
     s */ q' +/ t */ s', s */ r' +/ t */ t')
end

let next z = LFT.floor_ev z three
and safe z n = (n =/ LFT.floor_ev z four)
and prod z n = LFT.comp (ten, neg_ten */ n, zero, one) z
and cons z k =
  let den = 2*k+1 in
  LFT.comp z (num_of_int k, num_of_int(2*den), zero, num_of_int den)

let rec digit k z n row col =
  if n > 0 then
    let y = next z in
    if safe z y then
      if col = 10 then (
    let row = row + 10 in
    printf "\t:%i\n%s" row (string_of_num y);
    digit k (prod z y) (n-1) row 1
      )
      else (
    print_string(string_of_num y);
    digit k (prod z y) (n-1) row (col+1)
      )
    else digit (k+1) (cons z k) n row col
  else
    printf "%*s\t:%i\n" (10 - col) "" (row + col)

let digits n = digit 1 LFT.unit n 0 0

let () =
  digits(int_of_string Sys.argv.(1))

----Next_Part(Sun_Feb_13_19_24_10_2005_886)----

----------------
----------------
TagsNo tags attached.
Attached Files? file icon quo_num.ml [^] (2,297 bytes) 2008-05-25 00:33 [Show Content]

- Relationships

-  Notes
(0004510)
Christophe Troestler (reporter)
2008-05-25 00:32
edited on: 2008-05-25 00:35

I have attached a simple change to [quo_num] (refactorizing [div_num]) that makes the code using [Num] basically as fast as the one using [Big_int] (a ten fold improvement for n=1000).

(0004607)
weis (developer)
2008-09-10 18:25

Thank you for the suggestion.

I wrote a better implementation for quo_num and mod_num: now the num and big_int versions of your program run at approximately the same speed.

$ time pi_big_int 10000
3141592653 :10
...
5525637567 :10000

real 0m33.755s

$ time pi_num 10000
3141592653 :10
.....
5525637567 :10000

real 0m33.451s

It will be avaiblable in the next version of the compiler (version >= 3.11).

- Issue History
Date Modified Username Field Change
2005-11-18 10:13 administrator New Issue
2008-05-25 00:32 Christophe Troestler Note Added: 0004510
2008-05-25 00:33 Christophe Troestler File Added: quo_num.ml
2008-05-25 00:35 Christophe Troestler Note Edited: 0004510
2008-09-10 18:25 weis Note Added: 0004607
2008-09-10 18:34 weis Status acknowledged => resolved
2008-09-10 18:34 weis Resolution open => fixed
2008-09-10 18:34 weis Assigned To => weis
2009-04-29 20:41 weis Status resolved => closed
2009-04-29 20:41 weis Fixed in Version => 3.11.0+beta


Copyright © 2000 - 2011 MantisBT Group
Powered by Mantis Bugtracker