From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Delivered-To: caml-list@yquem.inria.fr Received: from concorde.inria.fr (concorde.inria.fr [192.93.2.39]) by yquem.inria.fr (Postfix) with ESMTP id 5D58CBC8B for ; Sun, 13 Feb 2005 20:05:10 +0100 (CET) Received: from pauillac.inria.fr (pauillac.inria.fr [128.93.11.35]) by concorde.inria.fr (8.13.0/8.13.0) with ESMTP id j1DJ5953024757 for ; Sun, 13 Feb 2005 20:05:10 +0100 Received: from nez-perce.inria.fr (nez-perce.inria.fr [192.93.2.78]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id UAA15249 for ; Sun, 13 Feb 2005 20:05:09 +0100 (MET) Received: from swip.net (mailfe06.swip.net [212.247.154.161]) by nez-perce.inria.fr (8.13.0/8.13.0) with ESMTP id j1DJ58Lo011496 for ; Sun, 13 Feb 2005 20:05:09 +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" Subject: Num.quo_num From: Christophe TROESTLER 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 Content-Type: Multipart/Mixed; boundary="--Next_Part(Sun_Feb_13_19_24_10_2005_886)--" Content-Transfer-Encoding: 8bit X-Miltered: at concorde with ID 420FA4E5.001 by Joe's j-chkmail (http://j-chkmail.ensmp.fr)! X-Miltered: at nez-perce with ID 420FA4E5.000 by Joe's j-chkmail (http://j-chkmail.ensmp.fr)! X-Spam: no; 0.00; christophe:01 troestler:01 christophe:01 troestler:01 umh:01 comlab:01 oucl:01 vastly:01 computed:01 printf:01 mult:01 struct:01 rec:01 printf:01 argv:01 X-Attachments: name="pidigits.ml" name="pidigits-num.ml" X-Spam-Checker-Version: SpamAssassin 3.0.2 (2004-11-16) on yquem.inria.fr X-Spam-Status: No, score=0.0 required=5.0 tests=none autolearn=disabled version=3.0.2 X-Spam-Level: ----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-Transfer-Encoding: 8bit Content-Disposition: inline; filename="pidigits.ml" (* 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 ( !$ ) = Big_int.big_int_of_int let ( +$ ) = Big_int.add_big_int let ( *$ ) = Big_int.mult_big_int let ( =$ ) = Big_int.eq_big_int let zero = Big_int.zero_big_int and one = Big_int.unit_big_int and three = !$ 3 and four = !$ 4 and ten = !$ 10 and neg_ten = !$(-10) (* Linear Fractional (aka Möbius) Transformations *) module LFT = struct let floor_ev (q,r,s,t) x = div_big_int (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 (!$ k, !$(2*den), zero, !$ 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_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 = 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)-- Content-Type: Text/Plain; charset=us-ascii Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="pidigits-num.ml" (* 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)----