Source

mutated_ocaml / otherlibs / num / big_int.ml

Full commit
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
(***********************************************************************)
(*                                                                     *)
(*                                OCaml                                *)
(*                                                                     *)
(*    Valerie Menissier-Morain, projet Cristal, INRIA Rocquencourt     *)
(*                                                                     *)
(*  Copyright 1996 Institut National de Recherche en Informatique et   *)
(*  en Automatique.  All rights reserved.  This file is distributed    *)
(*  under the terms of the GNU Library General Public License, with    *)
(*  the special exception on linking described in file ../../LICENSE.  *)
(*                                                                     *)
(***********************************************************************)

(* $Id: big_int.ml 11156 2011-07-27 14:17:02Z doligez $ *)

open Int_misc
open Nat

type big_int =
   { sign : int;
     abs_value : nat }

let create_big_int sign nat =
 if sign = 1 || sign = -1 ||
    (sign = 0 &&
     is_zero_nat nat 0 (num_digits_nat nat 0 (length_nat nat)))
 then { sign = sign;
         abs_value = nat }
 else invalid_arg "create_big_int"

(* Sign of a big_int *)
let sign_big_int bi = bi.sign

let zero_big_int =
 { sign = 0;
   abs_value = make_nat 1 }

let unit_big_int =
  { sign = 1;
    abs_value = nat_of_int 1 }

(* Number of digits in a big_int *)
let num_digits_big_int bi =
 num_digits_nat (bi.abs_value) 0 (length_nat bi.abs_value)

(* Opposite of a big_int *)
let minus_big_int bi =
 { sign = - bi.sign;
   abs_value = copy_nat (bi.abs_value) 0 (num_digits_big_int bi)}

(* Absolute value of a big_int *)
let abs_big_int bi =
    { sign = if bi.sign = 0 then 0 else 1;
      abs_value = copy_nat (bi.abs_value) 0 (num_digits_big_int bi)}

(* Comparison operators on big_int *)

(*
   compare_big_int (bi, bi2) = sign of (bi-bi2)
   i.e. 1 if bi > bi2
        0 if bi = bi2
        -1 if bi < bi2
*)
let compare_big_int bi1 bi2 =
  if bi1.sign = 0 && bi2.sign = 0 then 0
  else if bi1.sign < bi2.sign then -1
  else if bi1.sign > bi2.sign then 1
  else if bi1.sign = 1 then
            compare_nat (bi1.abs_value) 0 (num_digits_big_int bi1)
                        (bi2.abs_value) 0 (num_digits_big_int bi2)
  else
            compare_nat (bi2.abs_value) 0 (num_digits_big_int bi2)
                        (bi1.abs_value) 0 (num_digits_big_int bi1)

let eq_big_int bi1 bi2 = compare_big_int bi1 bi2 = 0
and le_big_int bi1 bi2 = compare_big_int bi1 bi2 <= 0
and ge_big_int bi1 bi2 = compare_big_int bi1 bi2 >= 0
and lt_big_int bi1 bi2 = compare_big_int bi1 bi2 < 0
and gt_big_int bi1 bi2 = compare_big_int bi1 bi2 > 0

let max_big_int bi1 bi2 = if lt_big_int bi1 bi2 then bi2 else bi1
and min_big_int bi1 bi2 = if gt_big_int bi1 bi2 then bi2 else bi1

(* Operations on big_int *)

let pred_big_int bi =
 match bi.sign with
    0 -> { sign = -1; abs_value = nat_of_int 1}
  | 1 -> let size_bi = num_digits_big_int bi in
          let copy_bi = copy_nat (bi.abs_value) 0 size_bi in
            ignore (decr_nat copy_bi 0 size_bi 0);
            { sign = if is_zero_nat copy_bi 0 size_bi then 0 else 1;
              abs_value = copy_bi }
  | _ -> let size_bi = num_digits_big_int bi in
         let size_res = succ (size_bi) in
         let copy_bi = create_nat (size_res) in
          blit_nat copy_bi 0 (bi.abs_value) 0 size_bi;
          set_digit_nat copy_bi size_bi 0;
          ignore (incr_nat copy_bi 0 size_res 1);
          { sign = -1;
            abs_value = copy_bi }

let succ_big_int bi =
 match bi.sign with
    0 -> {sign = 1; abs_value = nat_of_int 1}
  | -1 -> let size_bi = num_digits_big_int bi in
           let copy_bi = copy_nat (bi.abs_value) 0 size_bi in
            ignore (decr_nat copy_bi 0 size_bi 0);
            { sign = if is_zero_nat copy_bi 0 size_bi then 0 else -1;
              abs_value = copy_bi }
  | _ -> let size_bi = num_digits_big_int bi in
         let size_res = succ (size_bi) in
         let copy_bi = create_nat (size_res) in
          blit_nat copy_bi 0 (bi.abs_value) 0 size_bi;
          set_digit_nat copy_bi size_bi 0;
          ignore (incr_nat copy_bi 0 size_res 1);
          { sign = 1;
            abs_value = copy_bi }

let add_big_int bi1 bi2 =
 let size_bi1 = num_digits_big_int bi1
 and size_bi2 = num_digits_big_int bi2 in
  if bi1.sign = bi2.sign
   then    (* Add absolute values if signs are the same *)
    { sign = bi1.sign;
      abs_value =
       match compare_nat (bi1.abs_value) 0 size_bi1
                         (bi2.abs_value) 0 size_bi2 with
        -1 -> let res = create_nat (succ size_bi2) in
                (blit_nat res 0 (bi2.abs_value) 0 size_bi2;
                 set_digit_nat res size_bi2 0;
                 ignore
                   (add_nat res 0 (succ size_bi2)
                      (bi1.abs_value) 0 size_bi1 0);
                 res)
       |_  -> let res = create_nat (succ size_bi1) in
               (blit_nat res 0 (bi1.abs_value) 0 size_bi1;
                set_digit_nat res size_bi1 0;
                ignore (add_nat res 0 (succ size_bi1)
                         (bi2.abs_value) 0 size_bi2 0);
                res)}

  else      (* Subtract absolute values if signs are different *)
    match compare_nat (bi1.abs_value) 0 size_bi1
                      (bi2.abs_value) 0 size_bi2 with
       0 -> zero_big_int
     | 1 -> { sign = bi1.sign;
               abs_value =
                let res = copy_nat (bi1.abs_value) 0 size_bi1 in
                 (ignore (sub_nat res 0 size_bi1
                            (bi2.abs_value) 0 size_bi2 1);
                  res) }
     | _ -> { sign = bi2.sign;
              abs_value =
               let res = copy_nat (bi2.abs_value) 0 size_bi2 in
                 (ignore (sub_nat res 0 size_bi2
                            (bi1.abs_value) 0 size_bi1 1);
                  res) }

(* Coercion with int type *)
let big_int_of_int i =
  { sign = sign_int i;
    abs_value =
      let res = (create_nat 1)
      in (if i = monster_int
             then (set_digit_nat res 0 biggest_int;
                   ignore (incr_nat res 0 1 1))
             else set_digit_nat res 0 (abs i));
      res }

let add_int_big_int i bi = add_big_int (big_int_of_int i) bi

let sub_big_int bi1 bi2 = add_big_int bi1 (minus_big_int bi2)

(* Returns i * bi *)
let mult_int_big_int i bi =
 let size_bi = num_digits_big_int bi in
 let size_res = succ size_bi in
  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;
            ignore (mult_digit_nat res 0 size_res (bi.abs_value) 0 size_bi
                      (nat_of_int biggest_int) 0);
            { sign = - (sign_big_int bi);
              abs_value = res }
     else let res = make_nat (size_res) in
          ignore (mult_digit_nat res 0 size_res (bi.abs_value) 0 size_bi
                    (nat_of_int (abs i)) 0);
          { sign = (sign_int i) * (sign_big_int bi);
            abs_value = res }

let mult_big_int bi1 bi2 =
 let size_bi1 = num_digits_big_int bi1
 and size_bi2 = num_digits_big_int bi2 in
 let size_res = size_bi1 + size_bi2 in
 let res = make_nat (size_res) in
  { sign = bi1.sign * bi2.sign;
    abs_value =
         if size_bi2 > size_bi1
           then (ignore (mult_nat res 0 size_res (bi2.abs_value) 0 size_bi2
                           (bi1.abs_value) 0 size_bi1);res)
           else (ignore (mult_nat res 0 size_res (bi1.abs_value) 0 size_bi1
                           (bi2.abs_value) 0 size_bi2);res) }

(* (quotient, rest) of the euclidian division of 2 big_int *)
let quomod_big_int bi1 bi2 =
 if bi2.sign = 0 then raise Division_by_zero
 else
  let size_bi1 = num_digits_big_int bi1
  and size_bi2 = num_digits_big_int bi2 in
   match compare_nat (bi1.abs_value) 0 size_bi1
                     (bi2.abs_value) 0 size_bi2 with
      -1 -> (* 1/2  -> 0, reste 1, -1/2  -> -1, reste 1 *)
            (* 1/-2 -> 0, reste 1, -1/-2 -> 1, reste 1 *)
             if bi1.sign >= 0 then
               (big_int_of_int 0, bi1)
             else if bi2.sign >= 0 then
               (big_int_of_int(-1), add_big_int bi2 bi1)
             else
               (big_int_of_int 1, sub_big_int bi1 bi2)
    | 0 -> (big_int_of_int (bi1.sign * bi2.sign), zero_big_int)
    | _ -> let bi1_negatif = bi1.sign = -1 in
           let size_q =
            if bi1_negatif
             then succ (max (succ (size_bi1 - size_bi2)) 1)
             else max (succ (size_bi1 - size_bi2)) 1
           and size_r = succ (max size_bi1 size_bi2)
            (* r is long enough to contain both quotient and remainder *)
            (* of the euclidian division *)
           in
            (* set up quotient, remainder *)
            let q = create_nat size_q
            and r = create_nat size_r in
            blit_nat r 0 (bi1.abs_value) 0 size_bi1;
            set_to_zero_nat r size_bi1 (size_r - size_bi1);

            (* do the division of |bi1| by |bi2|
               - at the beginning, r contains |bi1|
               - at the end, r contains
                 * in the size_bi2 least significant digits, the remainder
                 * in the size_r-size_bi2 most significant digits, the quotient
              note the conditions for application of div_nat are verified here
             *)
            div_nat r 0 size_r (bi2.abs_value) 0 size_bi2;

            (* separate quotient and remainder *)
            blit_nat q 0 r size_bi2 (size_r - size_bi2);
            let not_null_mod = not (is_zero_nat r 0 size_bi2) in

            (* correct the signs, adjusting the quotient and remainder *)
            if bi1_negatif && not_null_mod
             then
              (* bi1<0, r>0, noting r for (r, size_bi2) the remainder,      *)
              (* we have |bi1|=q * |bi2| + r with 0 < r < |bi2|,            *)
              (* thus -bi1 = q * |bi2| + r                                  *)
              (* and bi1 = (-q) * |bi2| + (-r) with -|bi2| < (-r) < 0       *)
              (* thus bi1 = -(q+1) * |bi2| + (|bi2|-r)                      *)
              (* with 0 < (|bi2|-r) < |bi2|                                 *)
              (* so the quotient has for sign the opposite of the bi2'one   *)
              (*                 and for value q+1                          *)
              (* and the remainder is strictly positive                     *)
              (*                  has for value |bi2|-r                     *)
              (let new_r = copy_nat (bi2.abs_value) 0 size_bi2 in
                      (* new_r contains (r, size_bi2) the remainder *)
                { sign = - bi2.sign;
                  abs_value = (set_digit_nat q (pred size_q) 0;
                               ignore (incr_nat q 0 size_q 1); q) },
                { sign = 1;
                 abs_value =
                      (ignore (sub_nat new_r 0 size_bi2 r 0 size_bi2 1);
                      new_r) })
             else
              (if bi1_negatif then set_digit_nat q (pred size_q) 0;
                { sign = if is_zero_nat q 0 size_q
                          then 0
                          else bi1.sign * bi2.sign;
                  abs_value = q },
                { sign = if not_null_mod then 1 else 0;
                  abs_value = copy_nat r 0 size_bi2 })

let div_big_int bi1 bi2 = fst (quomod_big_int bi1 bi2)
and mod_big_int bi1 bi2 = snd (quomod_big_int bi1 bi2)

let gcd_big_int bi1 bi2 =
 let size_bi1 = num_digits_big_int bi1
 and size_bi2 = num_digits_big_int bi2 in
  if is_zero_nat (bi1.abs_value) 0 size_bi1 then abs_big_int bi2
  else if is_zero_nat (bi2.abs_value) 0 size_bi2 then
        { sign = 1;
          abs_value = bi1.abs_value }
  else
        { sign = 1;
          abs_value =
           match compare_nat (bi1.abs_value) 0 size_bi1
                             (bi2.abs_value) 0 size_bi2 with
           0 -> bi1.abs_value
         | 1 ->
            let res = copy_nat (bi1.abs_value) 0 size_bi1 in
            let len =
              gcd_nat res 0 size_bi1 (bi2.abs_value) 0 size_bi2 in
            copy_nat res 0 len
         | _ ->
            let res = copy_nat (bi2.abs_value) 0 size_bi2 in
            let len =
              gcd_nat res 0 size_bi2 (bi1.abs_value) 0 size_bi1 in
            copy_nat res 0 len
         }

(* Coercion operators *)

let monster_big_int = big_int_of_int monster_int;;

let monster_nat = monster_big_int.abs_value;;

let is_int_big_int bi =
  num_digits_big_int bi == 1 &&
  match compare_nat bi.abs_value 0 1 monster_nat 0 1 with
  | 0 -> bi.sign == -1
  | -1 -> true
  | _ -> false;;

let int_of_big_int bi =
  try let n = int_of_nat bi.abs_value in
      if bi.sign = -1 then - n else n
  with Failure _ ->
    if eq_big_int bi monster_big_int then monster_int
    else failwith "int_of_big_int";;

let big_int_of_nativeint i =
  if i = 0n then
    zero_big_int
  else if i > 0n then begin
    let res = create_nat 1 in
    set_digit_nat_native res 0 i;
    { sign = 1; abs_value = res }
  end else begin
    let res = create_nat 1 in
    set_digit_nat_native res 0 (Nativeint.neg i);
    { sign = -1; abs_value = res }
  end

let nativeint_of_big_int bi =
  if num_digits_big_int bi > 1 then failwith "nativeint_of_big_int";
  let i = nth_digit_nat_native bi.abs_value 0 in
  if bi.sign >= 0 then
    if i >= 0n then i else failwith "nativeint_of_big_int"
  else
    if i >= 0n || i = Nativeint.min_int
    then Nativeint.neg i
    else failwith "nativeint_of_big_int"

let big_int_of_int32 i = big_int_of_nativeint (Nativeint.of_int32 i)

let int32_of_big_int bi =
  let i = nativeint_of_big_int bi in
  if i <= 0x7FFF_FFFFn && i >= -0x8000_0000n
  then Nativeint.to_int32 i
  else failwith "int32_of_big_int"

let big_int_of_int64 i =
  if Sys.word_size = 64 then
    big_int_of_nativeint (Int64.to_nativeint i)
  else begin
    let (sg, absi) =
      if i = 0L then (0, 0L)
      else if i > 0L then (1, i)
      else (-1, Int64.neg i) in
    let res = create_nat 2 in
    set_digit_nat_native res 0 (Int64.to_nativeint absi);
    set_digit_nat_native res 1 (Int64.to_nativeint (Int64.shift_right absi 32));
    { sign = sg; abs_value = res }
  end

let int64_of_big_int bi =
  if Sys.word_size = 64 then
    Int64.of_nativeint (nativeint_of_big_int bi)
  else begin
    let i =
      match num_digits_big_int bi with
      | 1 -> Int64.logand
               (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 0))
               0xFFFFFFFFL
      | 2 -> Int64.logor
               (Int64.logand
                 (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 0))
                 0xFFFFFFFFL)
               (Int64.shift_left
                 (Int64.of_nativeint (nth_digit_nat_native bi.abs_value 1))
                 32)
      | _ -> failwith "int64_of_big_int" in
    if bi.sign >= 0 then
      if i >= 0L then i else failwith "int64_of_big_int"
    else
      if i >= 0L || i = Int64.min_int
      then Int64.neg i
      else failwith "int64_of_big_int"
  end

(* Coercion with nat type *)
let nat_of_big_int bi =
 if bi.sign = -1
 then failwith "nat_of_big_int"
 else copy_nat (bi.abs_value) 0 (num_digits_big_int bi)

let sys_big_int_of_nat nat off len =
 let length = num_digits_nat nat off len in
    { sign = if is_zero_nat nat off  length then 0 else 1;
      abs_value = copy_nat nat off length }

let big_int_of_nat nat =
 sys_big_int_of_nat nat 0 (length_nat nat)

(* Coercion with string type *)

let string_of_big_int bi =
  if bi.sign = -1
  then "-" ^ string_of_nat bi.abs_value
  else string_of_nat bi.abs_value


let sys_big_int_of_string_aux s ofs len sgn =
  if len < 1 then failwith "sys_big_int_of_string";
  let n = sys_nat_of_string 10 s ofs len in
  if is_zero_nat n 0 (length_nat n) then zero_big_int
  else {sign = sgn; abs_value = n}
;;

let sys_big_int_of_string s ofs len =
  if len < 1 then failwith "sys_big_int_of_string";
  match s.[ofs] with
  | '-' -> sys_big_int_of_string_aux s (ofs+1) (len-1) (-1)
  | '+' -> sys_big_int_of_string_aux s (ofs+1) (len-1) 1
  | _ -> sys_big_int_of_string_aux s ofs len 1
;;

let big_int_of_string s =
  sys_big_int_of_string s 0 (String.length s)

let power_base_nat base nat off len =
  if base = 0 then nat_of_int 0 else
  if is_zero_nat nat off len || base = 1 then nat_of_int 1 else
  let power_base = make_nat (succ length_of_digit) in
  let (pmax, pint) = make_power_base base power_base in
  let (n, rem) =
      let (x, y) = quomod_big_int (sys_big_int_of_nat nat off len)
                                  (big_int_of_int (succ pmax)) in
        (int_of_big_int x, int_of_big_int y) in
  if n = 0 then copy_nat power_base (pred rem) 1 else
   begin
    let res = make_nat n
    and res2 = make_nat (succ n)
    and l = num_bits_int n - 2 in
    let p = ref (1 lsl l) in
      blit_nat res 0 power_base pmax 1;
      for i = l downto 0 do
        let len = num_digits_nat res 0 n in
        let len2 = min n (2 * len) in
        let succ_len2 = succ len2 in
          ignore (square_nat res2 0 len2 res 0 len);
          begin
           if n land !p > 0
              then (set_to_zero_nat res 0 len;
                    ignore (mult_digit_nat res 0 succ_len2
                              res2 0 len2 power_base pmax))
              else blit_nat res 0 res2 0 len2
          end;
          set_to_zero_nat res2 0 len2;
          p := !p lsr 1
      done;
    if rem > 0
     then (ignore (mult_digit_nat res2 0 (succ n)
                     res 0 n power_base (pred rem));
           res2)
     else res
  end

let power_int_positive_int i n =
  match sign_int n with
    0 -> unit_big_int
  | -1 -> invalid_arg "power_int_positive_int"
  | _ -> let nat = power_base_int (abs i) n in
           { sign = if i >= 0
                       then sign_int i
                       else if n land 1 = 0
                               then 1
                               else -1;
             abs_value = nat}

let power_big_int_positive_int bi n =
  match sign_int n with
    0 -> unit_big_int
  | -1 -> invalid_arg "power_big_int_positive_int"
  | _ -> let bi_len = num_digits_big_int bi in
         let res_len = bi_len * n in
         let res = make_nat res_len
         and res2 = make_nat res_len
         and l = num_bits_int n - 2 in
         let p = ref (1 lsl l) in
         blit_nat res 0 bi.abs_value 0 bi_len;
         for i = l downto 0 do
           let len = num_digits_nat res 0 res_len in
           let len2 = min res_len (2 * len) in
           set_to_zero_nat res2 0 len2;
           ignore (square_nat res2 0 len2 res 0 len);
           if n land !p > 0 then begin
             let lenp = min res_len (len2 + bi_len) in
             set_to_zero_nat res 0 lenp;
             ignore(mult_nat res 0 lenp res2 0 len2 (bi.abs_value) 0 bi_len)
           end else begin
             blit_nat res 0 res2 0 len2
           end;
           p := !p lsr 1
         done;
         {sign = if bi.sign >=  0 then bi.sign
                 else if n land 1 = 0 then 1 else -1;
            abs_value = res}

let power_int_positive_big_int i bi =
  match sign_big_int bi with
    0 -> unit_big_int
  | -1 -> invalid_arg "power_int_positive_big_int"
  | _ -> let nat = power_base_nat
                     (abs i) (bi.abs_value) 0 (num_digits_big_int bi) in
           { sign = if i >= 0
                       then sign_int i
                       else if is_digit_odd (bi.abs_value) 0
                               then -1
                               else 1;
             abs_value = nat }

let power_big_int_positive_big_int bi1 bi2 =
  match sign_big_int bi2 with
    0 -> unit_big_int
  | -1 -> invalid_arg "power_big_int_positive_big_int"
  | _ -> try
           power_big_int_positive_int bi1 (int_of_big_int bi2)
         with Failure _ ->
         try
           power_int_positive_big_int (int_of_big_int bi1) bi2
         with Failure _ ->
           raise Out_of_memory
           (* If neither bi1 nor bi2 is a small integer, bi1^bi2 is not
              representable.  Indeed, on a 32-bit platform,
              |bi1| >= 2 and |bi2| >= 2^30, hence bi1^bi2 has at least
              2^30 bits = 2^27 bytes, greater than the max size of
              allocated blocks.  On a 64-bit platform,
              |bi1| >= 2 and |bi2| >= 2^62, hence bi1^bi2 has at least
              2^62 bits = 2^59 bytes, greater than the max size of
              allocated blocks. *)

(* base_power_big_int compute bi*base^n *)
let base_power_big_int base n bi =
  match sign_int n with
    0 -> bi
  | -1 -> let nat = power_base_int base (-n) in
           let len_nat = num_digits_nat nat 0 (length_nat nat)
           and len_bi = num_digits_big_int bi in
             if len_bi < len_nat then
               invalid_arg "base_power_big_int"
             else if len_bi = len_nat &&
                     compare_digits_nat (bi.abs_value) len_bi nat len_nat = -1
               then invalid_arg "base_power_big_int"
             else
               let copy = create_nat (succ len_bi) in
                      blit_nat copy 0 (bi.abs_value) 0 len_bi;
                      set_digit_nat copy len_bi 0;
                      div_nat copy 0 (succ len_bi)
                              nat 0 len_nat;
                      if not (is_zero_nat copy 0 len_nat)
                         then invalid_arg "base_power_big_int"
                         else { sign = bi.sign;
                                abs_value = copy_nat copy len_nat 1 }
  | _ -> let nat = power_base_int base n in
         let len_nat = num_digits_nat nat 0 (length_nat nat)
         and len_bi = num_digits_big_int bi in
         let new_len = len_bi + len_nat in
         let res = make_nat new_len in
         ignore
           (if len_bi > len_nat
               then mult_nat res 0 new_len
                              (bi.abs_value) 0 len_bi
                              nat 0 len_nat
               else mult_nat res 0 new_len
                              nat 0 len_nat
                              (bi.abs_value) 0 len_bi)
          ; if is_zero_nat res 0 new_len
               then zero_big_int
               else create_big_int (bi.sign) res

(* Coercion with float type *)

let float_of_big_int bi =
  float_of_string (string_of_big_int bi)

(* XL: suppression de big_int_of_float et nat_of_float. *)

(* Other functions needed *)

(* Integer part of the square root of a big_int *)
let sqrt_big_int bi =
 match bi.sign with
 | 0 -> zero_big_int
 | -1 -> invalid_arg "sqrt_big_int"
 | _ -> {sign = 1;
         abs_value = sqrt_nat (bi.abs_value) 0 (num_digits_big_int bi)}

let square_big_int bi =
  if bi.sign == 0 then zero_big_int else
  let len_bi = num_digits_big_int bi in
  let len_res = 2 * len_bi in
  let res = make_nat len_res in
  ignore (square_nat res 0 len_res (bi.abs_value) 0 len_bi);
  {sign = 1; abs_value = res}

(* round off of the futur last digit (of the integer represented by the string
   argument of the function) that is now the previous one.
   if s contains an integer of the form (10^n)-1
    then s <- only 0 digits and the result_int is true
   else s <- the round number and the result_int is false *)
let round_futur_last_digit s off_set length =
 let l = pred (length + off_set) in
  if Char.code(String.get s l) >= Char.code '5'
    then
     let rec round_rec l =
       if l < off_set then true else begin
         let current_char = String.get s l in
         if current_char = '9' then
           (String.set s l '0'; round_rec (pred l))
         else
           (String.set s l (Char.chr (succ (Char.code current_char)));
            false)
       end
     in round_rec (pred l)
   else false


(* Approximation with floating decimal point a` la approx_ratio_exp *)
let approx_big_int prec bi =
  let len_bi = num_digits_big_int bi in
  let n =
    max 0
        (int_of_big_int (
          add_int_big_int
            (-prec)
            (div_big_int (mult_big_int (big_int_of_int (pred len_bi))
                                      (big_int_of_string "963295986"))
                        (big_int_of_string "100000000")))) in
  let s =
    string_of_big_int (div_big_int bi (power_int_positive_int 10 n)) in
  let (sign, off, len) =
    if String.get s 0 = '-'
       then ("-", 1, succ prec)
       else ("", 0, prec) in
  if (round_futur_last_digit s off (succ prec))
       then (sign^"1."^(String.make prec '0')^"e"^
             (string_of_int (n + 1 - off + String.length s)))
       else (sign^(String.sub s off 1)^"."^
             (String.sub s (succ off) (pred prec))
             ^"e"^(string_of_int (n - succ off + String.length s)))

(* Logical operations *)

(* Shift left by N bits *)

let shift_left_big_int bi n =
  if n < 0 then invalid_arg "shift_left_big_int"
  else if n = 0 then bi
  else if bi.sign = 0 then bi
  else begin
    let size_bi = num_digits_big_int bi in
    let size_res = size_bi + ((n + length_of_digit - 1) / length_of_digit) in
    let res = create_nat size_res in
    let ndigits = n / length_of_digit in
    set_to_zero_nat res 0 ndigits;
    blit_nat res ndigits bi.abs_value 0 size_bi;
    let nbits = n mod length_of_digit in
    if nbits > 0 then
      shift_left_nat res ndigits size_bi res (ndigits + size_bi) nbits;
    { sign = bi.sign; abs_value = res }
  end

(* Shift right by N bits (rounds toward zero) *)

let shift_right_towards_zero_big_int bi n =
  if n < 0 then invalid_arg "shift_right_towards_zero_big_int"
  else if n = 0 then bi
  else if bi.sign = 0 then bi
  else begin
    let size_bi = num_digits_big_int bi in
    let ndigits = n / length_of_digit in
    let nbits = n mod length_of_digit in
    if ndigits >= size_bi then zero_big_int else begin
      let size_res = size_bi - ndigits in
      let res = create_nat size_res in
      blit_nat res 0 bi.abs_value ndigits size_res;
      if nbits > 0 then begin
        let tmp = create_nat 1 in
        shift_right_nat res 0 size_res tmp 0 nbits
      end;
      if is_zero_nat res 0 size_res
      then zero_big_int
      else { sign = bi.sign; abs_value = res }
    end
  end

(* Compute 2^n - 1 *)

let two_power_m1_big_int n =
  if n < 0 then invalid_arg "two_power_m1_big_int"
  else if n = 0 then zero_big_int
  else begin
    let size_res = (n + length_of_digit - 1) / length_of_digit in
    let res = make_nat size_res in
    set_digit_nat_native res (n / length_of_digit)
                             (Nativeint.shift_left 1n (n mod length_of_digit));
    ignore (decr_nat res 0 size_res 0);
    { sign = 1; abs_value = res }
  end

(* Shift right by N bits (rounds toward minus infinity) *)

let shift_right_big_int bi n =
  if n < 0 then invalid_arg "shift_right_big_int"
  else if bi.sign >= 0 then shift_right_towards_zero_big_int bi n
  else shift_right_towards_zero_big_int (sub_big_int bi (two_power_m1_big_int n)) n

(* Extract N bits starting at ofs.
   Treats bi in two's complement.
   Result is always positive. *)

let extract_big_int bi ofs n =
  if ofs < 0 || n < 0 then invalid_arg "extract_big_int"
  else if bi.sign = 0 then bi
  else begin
    let size_bi = num_digits_big_int bi in
    let size_res = (n + length_of_digit - 1) / length_of_digit in
    let ndigits = ofs / length_of_digit in
    let nbits = ofs mod length_of_digit in
    let res = make_nat size_res in
    if ndigits < size_bi then
      blit_nat res 0 bi.abs_value ndigits (min size_res (size_bi - ndigits));
    if bi.sign < 0 then begin
      (* Two's complement *)
      complement_nat res 0 size_res;
      ignore (incr_nat res 0 size_res 1)
    end;
    if nbits > 0 then begin
      let tmp = create_nat 1 in
      shift_right_nat res 0 size_res tmp 0 nbits
    end;
    let n' = n mod length_of_digit in
    if n' > 0 then begin
      let tmp = create_nat 1 in
      set_digit_nat_native tmp 0
          (Nativeint.shift_right_logical (-1n) (length_of_digit - n'));
      land_digit_nat res (size_res - 1) tmp 0
    end;
    if is_zero_nat res 0 size_res
    then zero_big_int
    else { sign = 1; abs_value = res }
  end

(* Bitwise logical operations.  Arguments must be >= 0. *)

let and_big_int a b =
  if a.sign < 0 || b.sign < 0 then invalid_arg "and_big_int"
  else if a.sign = 0 || b.sign = 0 then zero_big_int
  else begin
    let size_a = num_digits_big_int a
    and size_b = num_digits_big_int b in
    let size_res = min size_a size_b in
    let res = create_nat size_res in
    blit_nat res 0 a.abs_value 0 size_res;
    for i = 0 to size_res - 1 do
      land_digit_nat res i b.abs_value i
    done;
    if is_zero_nat res 0 size_res
    then zero_big_int
    else { sign = 1; abs_value = res }
  end

let or_big_int a b =
  if a.sign < 0 || b.sign < 0 then invalid_arg "or_big_int"
  else if a.sign = 0 then b
  else if b.sign = 0 then a
  else begin
    let size_a = num_digits_big_int a
    and size_b = num_digits_big_int b in
    let size_res = max size_a size_b in
    let res = create_nat size_res in
    let or_aux a' b' size_b' =
      blit_nat res 0 a'.abs_value 0 size_res;
      for i = 0 to size_b' - 1 do
        lor_digit_nat res i b'.abs_value i
      done in
    if size_a >= size_b
    then or_aux a b size_b
    else or_aux b a size_a;
    if is_zero_nat res 0 size_res
    then zero_big_int
    else { sign = 1; abs_value = res }
  end

let xor_big_int a b =
  if a.sign < 0 || b.sign < 0 then invalid_arg "xor_big_int"
  else if a.sign = 0 then b
  else if b.sign = 0 then a
  else begin
    let size_a = num_digits_big_int a
    and size_b = num_digits_big_int b in
    let size_res = max size_a size_b in
    let res = create_nat size_res in
    let xor_aux a' b' size_b' =
      blit_nat res 0 a'.abs_value 0 size_res;
      for i = 0 to size_b' - 1 do
        lxor_digit_nat res i b'.abs_value i
      done in
    if size_a >= size_b
    then xor_aux a b size_b
    else xor_aux b a size_a;
    if is_zero_nat res 0 size_res
    then zero_big_int
    else { sign = 1; abs_value = res }
  end