[Thomas.Fischbacher@Physik.Uni-Muenchen.DE: Bug#224417: Problem with the old Num implementation of ocaml] #8423

vicuna · 2003-12-18T22:22:15Z

Original bug ID: 1989
Reporter: administrator
Status: closed
Resolution: fixed
Priority: normal
Severity: minor
Category: ~DO NOT USE (was: OCaml general)

Bug description

Salut.

Je vous transmet ce bug report concernant l'ancienne implementation de
bignum. Je ne sais pas si cela vous interesse encore, mais comme
l'ancienne bignum est encore dans la version releaser de ocaml, ...

Amicalement,

Sven Luther

----- Forwarded message from Thomas Fischbacher Thomas.Fischbacher@Physik.Uni-Muenchen.DE -----

Envelope-to: luther@localhost
Delivery-date: Thu, 18 Dec 2003 22:24:08 +0100
X-Sieve: Server Sieve 2.2
Subject: Bug#224417: Problem with the old Num implementation of ocaml
Reply-To: Thomas Fischbacher Thomas.Fischbacher@Physik.Uni-Muenchen.DE,
224417@bugs.debian.org
Resent-From: Thomas Fischbacher Thomas.Fischbacher@Physik.Uni-Muenchen.DE
Resent-To: debian-bugs-dist@lists.debian.org
Resent-Cc: Sven Luther luther@debian.org
Resent-Date: Thu, 18 Dec 2003 20:33:04 UTC
Resent-Message-ID: handler.224417.B.107177910414660@bugs.debian.org
X-Debian-PR-Message: report 224417
X-Debian-PR-Package: ocaml
X-Debian-PR-Keywords:
From: Thomas Fischbacher Thomas.Fischbacher@Physik.Uni-Muenchen.DE
To: submit@bugs.debian.org
Cc: Marc Schaefer mschaefer@uni-duisburg.de,
Klaus Aehlig aehlig@rz.mathematik.uni-muenchen.de,
Wolfram Krause krause@th.physik.uni-frankfurt.de,
steffen.grunewald@aei.mpg.de
X-BOFH: Daemons did it
X-Scanned-By: MIMEDefang 2.37
X-Spam-Checker-Version: SpamAssassin
2.60-master.debian.org_2003_11_25-bugs.debian.org_2003_12_15
(1.212-2003-09-23-exp) on master.debian.org
X-Spam-Status: No, hits=-3.4 required=4.0 tests=HAS_PACKAGE,LARGE_HEX
autolearn=no
version=2.60-master.debian.org_2003_11_25-bugs.debian.org_2003_12_15
X-Spam-Level:
Resent-Sender: Debian BTS debbugs@master.debian.org
X-Spam: no; 0.00; uni-muenchen:01 readable:01 unreadable:01 mime-aware:01 docserver:01 3.06:01 severity:01 grave:01 string:01 int:01 -type:01 integer:01 interpreter:01 unix:01 module:01
X-Attachments: name="ym-debug.ml" name="ym-debug.ml"

Package: ocaml
Version: 3.06-21
Severity: grave

Please note that this bug concerns the "old" Num implementation that was
removed in a recent update. I nevertheless use(d) it as I am right
now doing research for a new string theory paper and wanted to
consistently use one and the same somewhat environment throughout the
process. Nevertheless, I suppose this affects all versions having the old num
implementation (next thing I'll do is to upgrade my entire system to the
new implementation and check the behaviour of that one for some large test
calculations).

The attachment will reproduce a condition where doing the same calculation
(which is purely functional to the outside and does not use global state
or side effects, but involves Num) many times in sequence will give
different results.

As you see, about half of this file is data; I did not manage to find a
smaller test case creating that condition. Although I consider this highly
unlikely, this piece of code may involve an "int"-type related integer
overflow. (Nevertheless, even so, subsequent identical calculations should
not give different results...)

On my system(*), the output of a num_ocaml interpreter built as described
in the Num section of the documentation, but also including the Unix
module is:

0 ==> -12747938802280433395210374808335212773156996019062242486412495785350259890027861091589626807758181401556588623046875/764762496031801254357991532266868085634832685275244705613839817948600630474284529321916381431208558769420593384116260078281979735715302794541235545562947747053568
1 ==> -24799801902139948285742077521010897923261300584495629346339410063333039539280530883995185222932342824043982837509521484375/4156317913216311164989084414493848291493655898235025573988259880155438209099372441966936855604394341138155398826718804773271628998452732579028454051972542103552
2 ==> -39148920061803210956691061036397438426365134774540146675772774556810648122275561412271743926625375084180283661376953125/2294287488095403763073974596800604256904498055825734116841519453845801891422853587965749144293625676308261780152348780234845939207145908383623706636688843241160704
3 ==> -39148920061803210956691061036397438426365134774540146675772774556810648122275561412271743926625375084180283661376953125/2294287488095403763073974596800604256904498055825734116841519453845801891422853587965749144293625676308261780152348780234845939207145908383623706636688843241160704
4 ==> -24799801902139948285742077521010897923261300584495629346339410063333039539280530883995185222932342824043982837509521484375/4156317913216311164989084414493848291493655898235025573988259880155438209099372441966936855604394341138155398826718804773271628998452732579028454051972542103552
5 ==> -24799801902139948285742077521010897923261300584495629346339410063333039539280530883995185222932342824043982837509521484375/4156317913216311164989084414493848291493655898235025573988259880155438209099372441966936855604394341138155398826718804773271628998452732579028454051972542103552
6 ==> -39148920061803210956691061036397438426365134774540146675772774556810648122275561412271743926625375084180283661376953125/2294287488095403763073974596800604256904498055825734116841519453845801891422853587965749144293625676308261780152348780234845939207145908383623706636688843241160704
7 ==> -39148920061803210956691061036397438426365134774540146675772774556810648122275561412271743926625375084180283661376953125/2294287488095403763073974596800604256904498055825734116841519453845801891422853587965749144293625676308261780152348780234845939207145908383623706636688843241160704
8 ==> -39148920061803210956691061036397438426365134774540146675772774556810648122275561412271743926625375084180283661376953125/2294287488095403763073974596800604256904498055825734116841519453845801891422853587965749144293625676308261780152348780234845939207145908383623706636688843241160704
9 ==> 0

: unit = ()

The correct result is zero.

(*) uname -a

Linux djinn 2.4.23 #6 Thu Dec 4 18:13:28 CET 2003 i686 GNU/Linux

I would like to know whether this can be confirmed.

--
regards, tf@cip.physik.uni-muenchen.de (o_
Dr. Thomas Fischbacher - http://www.cip.physik.uni-muenchen.de/~tf //
(lambda (n) ((lambda (p q r) (p p q r)) (lambda (g x y) V_/_
(if (= x 0) y (g g (- x 1) (* x y)))) n 1)) (Debian GNU)
open Num;;

let int1=Int 1;; (* To avoid number consing later on... *)
let int0=Int 0;;
let int2=Int 2;;

type lpower =
{mutable power: int;
mutable vars: int; (* bit vector )
mutable coeffs: int array ( always make denominators rational *)
};;

let power_of_2 n =
let rec trav p ppow =
if ppow>n
then -1
else if ppow<n
then trav (1+p) (ppow+ppow)
else p
in trav 0 1;;

(*
We need an order on lpowers...
*)

let v_cmp cmp a b =
let len_a = Array.length a in
let len_b = Array.length b in
let rec trav pos =
if pos == len_a
then
if len_a == len_b then 0 else -1
else
if pos == len_b
then 1
else
let xa = Array.get a pos in
let xb = Array.get b pos in
let cmp_xa_xb = cmp xa xb in
if cmp_xa_xb != 0 then cmp_xa_xb else trav (1+pos)
in trav 0;;

(*
Example:

v_cmp (fun x y -> if x < y then -1 else if x > y then 1 else 0) [|1;2;3;4;8|] [|1;2;3;4|];;
*)

let lpower_cmp also_cmp_powers lp_a lp_b =
if lp_a.vars < lp_b.vars then -1
else if lp_a.vars > lp_b.vars then 1
else let cmp =
v_cmp (fun x y -> if x < y
then -1
else if x > y then 1
else 0)
lp_a.coeffs lp_b.coeffs in
if cmp != 0 then cmp
else
if also_cmp_powers
then if lp_a.power < lp_b.power then -1
else if lp_a.power > lp_b.power then 1
else 0
else 0;;

(*
When combining lpower vectors,
we first have to count the number of different entries...
*)

let lpower_nr_different_noncancellating v_lp_a v_lp_b =
let len_a = Array.length v_lp_a in
let len_b = Array.length v_lp_b in
let rec trav n_so_far pos_a pos_b =
if pos_a == len_a then n_so_far + len_b-pos_b
else if pos_b == len_b then n_so_far + len_a-pos_a
else
let entry_a = Array.get v_lp_a pos_a in
let entry_b = Array.get v_lp_b pos_b in
let cmp_here = lpower_cmp false entry_a entry_b in
if cmp_here == 0 then
if entry_a.power + entry_b.power == 0
then trav n_so_far (1+pos_a) (1+pos_b)
else trav (1+n_so_far) (1+pos_a) (1+pos_b)
else if cmp_here == -1 then trav (1+n_so_far) (1+pos_a) pos_b
else trav (1+n_so_far) pos_a (1+pos_b)
in trav 0 0 0;;

(*
We will make good use of lpower substitution not only
for doing calculations, but also for constructing
the initial terms.
*)

(*

First of all, we need integer bitmask functions...

*)

let integer_length n =
let rec trav rest_n bit_now =
if rest_n == 0
then bit_now
else
let contrib = 1 lsl bit_now in
if 0 == contrib land rest_n
then trav rest_n (1+bit_now)
else trav (rest_n-contrib) (1+bit_now)
in trav n 0;;

let nr_bits_in n =
let rec trav rest_n bit_now nr =
if rest_n == 0
then nr
else
let contrib = 1 lsl bit_now in
if 0 == contrib land rest_n
then trav rest_n (1+bit_now) nr
else trav (rest_n-contrib) (1+bit_now) (1+nr)
in trav n 0 0;;

(* Given a number and a value 2^bit, tell me the how manyth set bit
in the number it is, starting to count at 0... *)

let nr_bits_set_before num bitpow =
let rec walk nr val_now =
if val_now == bitpow
then nr
else walk (if (0==val_now land num) then nr else nr+1)
(val_now+val_now)
in walk 0 1;;

(* ... and some other helpers... *)

let gcd a b =
let rec work p q =
let mpq = p mod q in
if mpq==0 then q else work q mpq
in
let aa = abs a in
let ab = abs b in
if aa > ab then work aa ab else work ab aa;;

let v_gcd v_int =
let len = Array.length v_int in
if len == 0
then 1
else if len == 1 then abs (Array.get v_int 0)
else let rec trav pos gcd_now =
if gcd_now==1 then 1
else if pos==len then gcd_now
else trav (1+pos) (gcd gcd_now (Array.get v_int pos))
in trav 1 (Array.get v_int 0);;

let lcm a b = abs(a*b/(gcd a b));;

(*

Here, quite some interesting things do happen...

Note that the given var MUST occur in the substitution rule...

We desructively change the extra factors in ref_numer and ref_denom
in order to be able to keep our coeffs integral and stay fast...

Note that one has to convince oneself that, by looking at the given
structures of the terms showing up, we do not leave integer range
by applying #r substitutions of linear functions with coeffs max. 2
(and what one can get from that in consecutive subs) as well as
powers up to 2^r...

Now, how big can we make r before violating these constraints...?

XXX Note: need another variant, where we do not cons the result,
but instead just eval it!

*)

let get_coeff lp coeff_bitpow =
let pos = nr_bits_set_before lp.vars coeff_bitpow in
Array.get lp.coeffs pos;;

let lpower_subs ref_xfactor var_bitpow to_subs lp =
if (var_bitpow land lp.vars) == 0 then lp (* nothing to substitute )
else
let subs_other_vars = var_bitpow lxor to_subs.vars in
let lp_other_vars = var_bitpow lxor lp.vars in
let total_vars = subs_other_vars lor lp_other_vars in
let pos_var_in_lp = nr_bits_set_before lp.vars var_bitpow in
let pos_var_in_subs = nr_bits_set_before to_subs.vars var_bitpow in
let var_coeff_in_lp = Array.get lp.coeffs pos_var_in_lp in
let var_coeff_in_subs = Array.get to_subs.coeffs pos_var_in_subs in
let var_coeff = lcm var_coeff_in_lp var_coeff_in_subs in
let factor_subs = -var_coeff/var_coeff_in_subs in
let factor_lp = var_coeff/var_coeff_in_lp in
let nr_total_vars = nr_bits_in total_vars in ( May still be reduced by additional cancellations! )
let result_coeffs = Array.make nr_total_vars 0 in
let rec fill_result bitpow result_pos true_result_vars =
if bitpow > total_vars
then ( finished -- but may have to take sub-range! *)
let final_coeffs =
if result_pos==nr_total_vars
then result_coeffs
else Array.sub result_coeffs 0 result_pos in
{power=lp.power;vars=true_result_vars;coeffs=final_coeffs}
else if 0 != (lp_other_vars land bitpow)
then if 0 != (subs_other_vars land bitpow)
then let val_to_set = (get_coeff lp bitpow)*factor_lp
+ (get_coeff to_subs bitpow)factor_subs in
if val_to_set==0
then fill_result (bitpow+bitpow) result_pos true_result_vars
else (Array.set result_coeffs result_pos val_to_set;
fill_result (bitpow+bitpow) (1+result_pos) (true_result_vars+bitpow))
else ( contrib from lp, but not from subs... *)
(Array.set result_coeffs result_pos
((get_coeff lp bitpow)factor_lp);
fill_result (bitpow+bitpow) (1+result_pos) (true_result_vars+bitpow))
else if 0 != (subs_other_vars land bitpow)
( Contrib from subs, but not from lp *)
then (Array.set result_coeffs result_pos
((get_coeff to_subs bitpow)factor_subs);
fill_result (bitpow+bitpow) (1+result_pos) (true_result_vars+bitpow))
else ( both contribs zero... )
fill_result (bitpow+bitpow) result_pos true_result_vars in
let result_lp = fill_result 1 0 0 in
if Array.length result_lp.coeffs == 0
( if we substituted X by X, we get out an empty factor...
For term normalization purposes, we make 0^0=1.
)
then
(result_lp.power <- 0;
result_lp.vars <- 0;
result_lp)
else
let coeffs_gcd0 = v_gcd result_lp.coeffs in
let coeffs_gcd = if (Array.get result_lp.coeffs 0) < 0 then (-coeffs_gcd0) else coeffs_gcd0 in
( Printf.printf "coeffs_gcd0=%d coeffs_gcd=%d\n" coeffs_gcd0 coeffs_gcd; DDD *)
(
if coeffs_gcd != 1
then Array.iteri (fun pos coeff -> Array.set result_lp.coeffs pos (coeff/coeffs_gcd)) result_lp.coeffs
else ();
ref_xfactor:= !ref_xfactor */ (power_num ((Int factor_lp)//(Int coeffs_gcd)) (Int (-lp.power)));
result_lp
);;

let v_lpower_subs_uniq numcoeff_plus_v_lp var_bitpow to_subs =
let (nc,v_lp) = numcoeff_plus_v_lp in
let rnc = ref nc in
let len_v_lp = Array.length v_lp in
let v_lp_subs = Array.map (lpower_subs rnc var_bitpow to_subs) v_lp in
(Array.sort (lpower_cmp true) v_lp_subs;
let rec trav_uniq pos_src pos_dst =
if pos_src == len_v_lp
then if pos_dst == pos_src
then v_lp_subs
else Array.sub v_lp_subs 0 pos_dst
else (* still work to be done )
if pos_dst==0
then ( no terms on stack )
(Array.set v_lp_subs pos_dst (Array.get v_lp_subs pos_src);
trav_uniq (1+pos_src) (1+pos_dst))
else
let lp1 = Array.get v_lp_subs (pos_dst-1) in
let lp2 = Array.get v_lp_subs pos_src in
let cmp = lpower_cmp false lp1 lp2 in
if cmp==0 ( combinable terms )
then let total_power = lp1.power + lp2.power in
if total_power == 0
then trav_uniq (1+pos_src) (pos_dst-1) ( Term on stack was eaten )
else (Array.set v_lp_subs (pos_dst-1)
{power=total_power;
vars=lp1.vars;
coeffs=lp1.coeffs}
( This conses more often than strictly necessary...
However, I do not want to worry about when I already
have copied a value and may do destructive
modification...
);
trav_uniq (1+pos_src) pos_dst)
else ( cmp!=0, which means that we have to transfer next term *)
(Array.set v_lp_subs pos_dst lp2;
trav_uniq (1+pos_src) (1+pos_dst)) in
let v_lp_subs_uniq = trav_uniq 0 0 in
(!rnc,v_lp_subs_uniq));;

let eval_nc_v_lp nc_v_lp v_Ex =
let (nc,v_lp)=nc_v_lp in
let nr_factors=Array.length v_lp in
let bitmask = (1 lsl (Array.length v_Ex))-1 in
(* quick hack: use a v_lp to retrieve coeff data... *)
let eval_lp={power=0;vars=bitmask;coeffs=v_Ex} in
let rec trav so_far n =
if n == nr_factors then so_far else
let factor_n = Array.get v_lp n in
let coeffs_n = factor_n.coeffs in
let nr_coeffs=Array.length coeffs_n in
let rec trav_vars rest_vars bitpow sum_now =
if rest_vars==0
then sum_now
else if (bitpow land rest_vars) == 0
then trav_vars rest_vars (bitpow+bitpow) sum_now
else
let coeff_j = get_coeff factor_n bitpow in
let var_j = get_coeff eval_lp bitpow in
trav_vars (rest_vars-bitpow) (bitpow+bitpow) (sum_now+/ (Int coeff_j) / (Int var_j))
in let val_factor_n = trav_vars factor_n.vars 1 (Int 0) in
trav (so_far/(power_num val_factor_n (Int factor_n.power))) (1+ n) in
trav nc 0;;

let check_substitution nc_v_lp =
let debug_lp_res = {power = -1; vars = 80; coeffs = [|3; -2|]} in
let v_eval = [|0;1000;1100;1110;-1000-1100-1110;101;3*(-1000-1100-1110)/2|] in
let var_bitpow = 1 lsl ((integer_length debug_lp_res.vars)-1) in
let substituted = v_lpower_subs_uniq nc_v_lp var_bitpow debug_lp_res in
let val_substituted = eval_nc_v_lp substituted v_eval in
let val_direct = eval_nc_v_lp nc_v_lp v_eval in
let difference = val_substituted-/val_direct in
(* List.map string_of_num [val_substituted;val_direct;difference] *)
string_of_num difference;;

let rec trav n =
if n == 10 then () else
(Printf.printf "%d ==> %s\n" n (check_substitution debug_v_lp);
trav (1+ n))
in trav 0;;

----- End forwarded message -----

vicuna · 2003-12-27T22:55:23Z

Comment author: administrator

On Friday, December 19, 2003, at 11:17 AM, Thomas Fischbacher wrote:

Please note that this bug concerns the "old" Num implementation that
was
removed in a recent update.

This bug is in fact in the ML part of the Num library, which is common
to both the old and new Num implementations. I think it dates back to
Caml Light.

The attachment will reproduce a condition where doing the same
calculation
(which is purely functional to the outside and does not use global
state
or side effects, but involves Num) many times in sequence will give
different results.

This "nondeterminism" comes from using a Nat.nat with an uninitialised
digit, using whatever was already in memory at that address.

arithmetics I am using? Yet to me it seems plausible that some serious
issue may have been overlooked here, since I suppose very few people
use
ratio arith (in contrast to bignum-only arith, which has applications
for,
e.g., RSA).

Indeed, the bug is in the big-int part of the library (arbitrary
precision
signed integers).

Thank you for your well-delimited bug report with a reproducible test
case. It was rather easy to find the bug with it.

I enclose a patch, which should work on 3.06 and on 3.07pl2.

-- Damien

*** otherlibs/num/big_int.ml.old Sat Dec 27 23:45:12 2003
--- otherlibs/num/big_int.ml Sat Dec 27 23:38:44 2003

*** 178,183 ****
--- 178,184 ----
if i = monster_int
then let res = create_nat size_res in
blit_nat res 0 (bi.abs_value) 0 size_bi;

        set_digit_nat res size_bi 0;
        mult_digit_nat res 0 size_res (bi.abs_value) 0 size_bi 
                       (nat_of_int biggest_int) 0;
        { sign = - (sign_big_int bi);

vicuna · 2003-12-29T19:26:30Z

Comment author: administrator

fixed DD 2003-12-29

vicuna closed this as completed Dec 29, 2003

vicuna added the bug label Mar 19, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[Thomas.Fischbacher@Physik.Uni-Muenchen.DE: Bug#224417: Problem with the old Num implementation of ocaml] #8423

[Thomas.Fischbacher@Physik.Uni-Muenchen.DE: Bug#224417: Problem with the old Num implementation of ocaml] #8423

vicuna commented Dec 18, 2003

vicuna commented Dec 27, 2003

vicuna commented Dec 29, 2003

[Thomas.Fischbacher@Physik.Uni-Muenchen.DE: Bug#224417: Problem with the old Num implementation of ocaml] #8423

[Thomas.Fischbacher@Physik.Uni-Muenchen.DE: Bug#224417: Problem with the old Num implementation of ocaml] #8423

Comments

vicuna commented Dec 18, 2003

Bug description

vicuna commented Dec 27, 2003

vicuna commented Dec 29, 2003