Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WISH] Better implementation of Num.quo_num #3473

Closed
vicuna opened this issue Feb 18, 2005 · 2 comments
Closed

[WISH] Better implementation of Num.quo_num #3473

vicuna opened this issue Feb 18, 2005 · 2 comments
Assignees

Comments

@vicuna
Copy link

vicuna commented Feb 18, 2005

Original bug ID: 3473
Reporter: administrator
Assigned to: @pierreweis
Status: closed (set by @pierreweis on 2009-04-29T18:41:47Z)
Resolution: fixed
Priority: normal
Severity: feature
Fixed in version: 3.11.0+beta
Category: ~DO NOT USE (was: OCaml general)
Monitored by: @Chris00

Bug description

Hi,

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)----



File attachments

@vicuna
Copy link
Author

vicuna commented May 24, 2008

Comment author: @Chris00

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).

@vicuna
Copy link
Author

vicuna commented Sep 10, 2008

Comment author: @pierreweis

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).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants