File: cr.ml

package info (click to toggle)
ocamlcreal 0.7-6
  • links: PTS, VCS
  • area: main
  • in suites: buster, jessie, jessie-kfreebsd, sid, squeeze, stretch, wheezy
  • size: 384 kB
  • ctags: 751
  • sloc: ml: 2,699; ansic: 1,067; makefile: 15
file content (860 lines) | stat: -rw-r--r-- 26,632 bytes parent folder | download | duplicates (3)
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
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
(* This module is a translation in ocaml of Hans's Boehm Java library CR *)

(*i
// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED 
// 
// Permission is granted free of charge to copy, modify, use and distribute
// this software  provided you include the entirety of this notice in all
// copies made.
// 
// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
// FOR A PARTICULAR PURPOSE OR NON-INFRINGING.   SGI ASSUMES NO RISK AS TO THE
// QUALITY AND PERFORMANCE OF THE SOFTWARE.   SHOULD THE SOFTWARE PROVE
// DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY
// SERVICING, REPAIR OR CORRECTION.  THIS DISCLAIMER OF WARRANTY CONSTITUTES
// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
// 
// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
// OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF
// THE POSSIBILITY OF SUCH DAMAGES.  THIS LIMITATION OF LIABILITY SHALL NOT
// APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE
// LAW PROHIBITS SUCH LIMITATION.  SOME JURISDICTIONS DO NOT ALLOW THE
// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
// 
// These license terms shall be governed by and construed in accordance with
// the laws of the United States and the State of California as applied to
// agreements entered into and to be performed entirely within California
// between California residents.  Any litigation relating to these license
// terms shall be subject to the exclusive jurisdiction of the Federal Courts
// of the Northern District of California (or, absent subject matter
// jurisdiction in such courts, the courts of the State of California), with
// venue lying exclusively in Santa Clara County, California. 

// Copyright (c) 2001-2004, Hewlett-Packard Development Company, L.P. 
// 
// Permission is granted free of charge to copy, modify, use and distribute
// this software  provided you include the entirety of this notice in all
// copies made.
// 
// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
// FOR A PARTICULAR PURPOSE OR NON-INFRINGING.   HEWLETT-PACKARD ASSUMES
// NO RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE.
// SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, 
// HEWLETT-PACKARD ASSUMES NO COST OR LIABILITY FOR ANY
// SERVICING, REPAIR OR CORRECTION.  THIS DISCLAIMER OF WARRANTY CONSTITUTES
// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
// 
// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
// OTHERWISE, SHALL HEWLETT-PACKARD BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF HEWLETT-PACKARD SHALL
// HAVE BEEN INFORMED OF THE POSSIBILITY OF SUCH DAMAGES.
// THIS LIMITATION OF LIABILITY SHALL NOT APPLY TO LIABILITY RESULTING
// FROM HEWLETT-PACKARD's NEGLIGENCE TO THE EXTENT APPLICABLE
// LAW PROHIBITS SUCH LIMITATION.  SOME JURISDICTIONS DO NOT ALLOW THE
// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
// 
i*)

open Gmp

type t = {
  mutable cache : (int * Z.t) option;
  approximate : int -> Z.t 
}

let create f = { cache = None; approximate = f }

let bound_log2 n =
  let np1 = float (abs n + 1) in
  truncate (ceil (log np1 /. log 2.0))

let z_zero = Z.from_int 0
let z_one = Z.from_int 1
let z_mone = Z.from_int (-1)
let z_two = Z.from_int 2
let z_three = Z.from_int 3
let z_four = Z.from_int 4
let z_six = Z.from_int 6

let z_shift_left = Z.mul2exp
let z_shift_right = Z.fdiv_q_2exp

let shift k n =
  if n == 0 then k
  else if n < 0 then z_shift_right k (-n)
  else z_shift_left k n

let scale k n =
  if n >= 0 then
    z_shift_left k n
  else
    let adj_k = Z.add_ui (shift k (n+1)) 1 in
    z_shift_right adj_k 1

exception PrecisionOverflow

let check_prec n = 
  let high = n lsr 28 in
  let high_shifted = n lsr 29 in
  (* TODO *)
  ()

let approx x p =
  check_prec p;
  match x.cache with
    | Some (min_p, ma) when p >= min_p -> 
	scale ma (min_p - p)
    | _ -> 
	let r = x.approximate p in
	x.cache <- Some (p, r);
	r

(* ceil(log2(if z < 0 then -z else z+1)) *)
let z_bit_length z =
  let rec loop k two_k =
    if Z.cmp z two_k <= 0 then k else loop (succ k) (Z.mul2exp two_k 1)
  in
  let z = if Z.sgn z < 0 then Z.neg z else Z.add_ui z 1 in
  if Z.cmp z z_one = 0 then 0 else loop 1 z_two

let known_msd x = match x.cache with
  | Some (mp, ma) ->
      let length = 
	if Z.sgn ma >= 0 then z_bit_length ma else z_bit_length (Z.neg ma)
      in
      mp + length - 1
  | None -> 
      assert false

let msd_n x n =
  let close_to_0 = match x.cache with
    | None -> true
    | Some (_, ma) -> Z.cmp ma z_one <= 0 && Z.cmp ma z_mone >= 0
  in
  if close_to_0 then begin
    ignore (approx x (n-1));
    match x.cache with
      | Some (_, ma) when Z.cmp (Z.abs ma) z_one <= 0 -> min_int
      | Some _ -> known_msd x
      | None -> assert false
  end else
    known_msd x

let iter_msd x n =
  let rec iter prec =
    if prec > n + 30 then
      let msd = msd_n x prec in
      if msd != min_int then msd
      else begin
	check_prec prec;
	iter ((prec * 3)/2 - 16)
      end
    else
      msd_n x n
  in
  iter 0

let msd x = iter_msd x min_int

let shifted x n = create (fun p -> approx x (p - n))
let shift_left x n = shifted x n
let shift_right x n = shifted x (-n)

let of_z z = create (fun p -> scale z (-p))

let of_int n = of_z (Z.from_int n)

let of_int64 n = let z = Z.from_string (Int64.to_string n) in of_z z

let zero = of_int 0
let one = of_int 1
let two = of_int 2
let four = of_int 4

(* Operations *)

let add x y =
  create (fun p -> scale (Z.add (approx x (p-2)) (approx y (p-2))) (-2))

let (+!) = add

let neg x = create (fun p -> Z.neg (approx x p))

let sub x y = add x (neg y)

let (-!) = sub

exception Zero

let mul x y =
  create 
    (fun p ->
       let half_prec = (p asr 1) - 1 in
       try
	 let x,y,msd_x =
	   let msd_x = msd_n x half_prec in
	   if msd_x == min_int then
	     let msd_y = msd_n y half_prec in
	     if msd_y == min_int then raise Zero else y,x,msd_y
	   else
	     x,y,msd_x
	 in
	 let prec2 = p - msd_x - 3 in
	 let appry = approx y prec2 in
	 if Z.sgn appry == 0 then raise Zero;
	 let msd_y = known_msd y in
	 let prec1 = p - msd_y - 3 in
	 let apprx = approx x prec1 in
	 let scale_digits = prec1 + prec2 - p in
	 scale (Z.mul apprx appry) scale_digits
       with Zero ->
	 z_zero)

let ( *! ) = mul

let inv x =
  create
    (fun p ->
       let msd = msd x in
       let inv_msd = 1 - msd in
       let digits_needed = inv_msd - p + 3 in
       let prec_needed = msd - digits_needed in
       let log_scale_factor = -p - prec_needed in
       if log_scale_factor < 0 then
	 z_zero
       else
	 let dividend = z_shift_left z_one log_scale_factor in
	 let scaled_divisor = approx x prec_needed in
	 let abs_scaled_divisor = Z.abs scaled_divisor in
	 let adj_dividend = 
	   Z.add dividend (z_shift_right abs_scaled_divisor 1) 
	 in
	 let result = Z.fdiv_q adj_dividend abs_scaled_divisor in
	 if Z.sgn scaled_divisor < 0 then Z.neg result else result)

let div x y = mul x (inv y)

let (/!) = div

let sgn_n x a = 
  let quick_try = match x.cache with
    | Some (_, ma) -> Z.sgn ma
    | None -> 0
  in
  if quick_try != 0 then 
    quick_try
  else
    let needed_prec = a - 1 in
    let appr = approx x needed_prec in
    Z.sgn appr

let sgn x =
  let rec loop a =
    check_prec a;
    let result = sgn_n x a in
    if result != 0 then result else loop (2 * a)
  in
  loop (-20)

let select s x y =
  let selector_sign = ref (Z.sgn (approx s (-20))) in
  create
    (fun p ->
       if !selector_sign < 0 then
	 approx x p
       else if !selector_sign > 0 then
	 approx y p
       else
	 let x_appr = approx x (p-1) in
	 let y_appr = approx y (p-1) in
	 let diff = Z.abs (Z.sub x_appr y_appr) in
	 if Z.cmp diff z_one <= 0 then
	   scale x_appr (-1)
	 else if sgn s < 0 then begin
	   selector_sign := -1;
	   scale x_appr (-1)
	 end else begin
	   selector_sign := 1;
	   scale y_appr (-1)
	 end)

let abs x = select x (neg x) x

let max x y = select (sub x y) y x
let min x y = select (sub x y) x y

let compare_a x y a =
  let needed_prec = a - 1 in
  let x_appr = approx x needed_prec in
  let y_appr = approx y needed_prec in
  let comp1 = Z.cmp x_appr (Z.add_ui y_appr 1) in
  if comp1 > 0 then 
    1
  else 
    let comp2 = Z.cmp x_appr (Z.add_ui y_appr (-1)) in
    if comp2 < 0 then -1 else 0

let compare x y =
  let rec loop a =
    check_prec a;
    let r = compare_a x y a in
    if r <> 0 then r else loop (2 * a)
  in
  loop (-20)

let rec pow_int x n =
  if n == 0 then
    one
  else if n < 0 then
    inv (pow_int x (-n))
  else 
    let y = pow_int (mul x x) (n / 2) in
    if n mod 2 == 0 then y else mul y x

let prescaled_exp x =
  create
    (fun p ->
       if p >= 1 then
	 z_zero
       else
	 let iterations_needed = -p/2 + 2 in
	 let calc_precision = p - bound_log2(2*iterations_needed) - 4 in
	 let op_prec = p - 3 in
	 let op_appr = approx x op_prec in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let current_term = ref scaled_1 in
	 let current_sum = ref scaled_1 in
	 let n = ref 0 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   incr n;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   current_term := Z.fdiv_q_ui !current_term !n;
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let rec exp x =
  let low_prec = -10 in
  let rough_appr = approx x low_prec in
  if Z.sgn rough_appr < 0 then
    inv (exp (neg x))
  else if Z.cmp rough_appr z_two > 0 then
    let square_root = exp (shift_right x 1) in
    mul square_root square_root
  else
    prescaled_exp x

let e = prescaled_exp one

let sqrt x =
  let fp_prec = 50 in
  let fp_op_prec = 60 in
  let rec sqrt_rec p =
    let max_prec_needed = 2*p - 1 in
    let msd = msd_n x max_prec_needed in
    if msd <= max_prec_needed then
      z_zero
    else
      let result_msd = msd / 2 in
      let result_digits = result_msd - p in
      if result_digits > fp_prec then
	let appr_digits = result_digits/2 + 6 in
	let appr_prec = result_msd - appr_digits in
	let last_appr = sqrt_rec appr_prec in
	let prod_prec = 2 * appr_prec in
	let op_appr = approx x prod_prec in
	let prod_prec_scaled_numerator = 
	  Z.add (Z.mul last_appr last_appr) op_appr
	in
	let scaled_numerator = 
	  scale prod_prec_scaled_numerator (appr_prec - p)
	in
	assert (Z.cmp last_appr z_zero != 0);
	let shifted_result = Z.fdiv_q scaled_numerator last_appr in
	z_shift_right (Z.add shifted_result z_one) 1
      else begin
	let op_prec = (msd - fp_op_prec) land (lnot 1) in
	let working_prec = op_prec - fp_op_prec in
	let scaled_bi_appr = z_shift_left (approx x op_prec) fp_op_prec in
	let scaled_appr = Z.float_from scaled_bi_appr in
	if scaled_appr < 0.0 then invalid_arg "Cr.sqrt";
	let scaled_fp_sqrt = sqrt scaled_appr in
	let scaled_sqrt = 
	  Z.from_string (Int64.to_string (Int64.of_float scaled_fp_sqrt)) 
	in
	let shift_count = working_prec/2 - p in
	shift scaled_sqrt shift_count
      end
  in
  create sqrt_rec


let prescaled_ln x = 
  create
    (fun p ->
       if p >= 0 then
	 z_zero
       else
	 let iterations_needed = -p in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 4 in
	 let op_prec = p - 3 in
	 let op_appr = approx x op_prec in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let x_nth = ref (scale op_appr (op_prec - calc_precision)) in
	 let current_term = ref !x_nth in
	 let current_sum = ref !current_term in
	 let n = ref 1 in
	 let current_sign = ref 1 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   incr n;
	   current_sign := - !current_sign;
	   x_nth := scale (Z.mul !x_nth op_appr) op_prec;
	   current_term := Z.fdiv_q !x_nth (Z.of_int (!n * !current_sign));
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let simple_ln x = prescaled_ln (sub x one)

(* ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) *)
let ten_ninths = of_int 10 /! of_int 9
let ln2_1 = of_int 7 *! simple_ln ten_ninths
let twentyfive_twentyfourths = of_int 25 /! of_int 24
let ln2_2 = of_int 2 *! simple_ln twentyfive_twentyfourths
let eightyone_eightyeths = of_int 81 /! of_int 80
let ln2_3 = of_int 3 *! simple_ln eightyone_eightyeths
let ln2 = ln2_1 -! ln2_2 +! ln2_3

let low_ln_limit = Z.of_int 8
let high_ln_limit = Z.of_int (16 + 8)
let scaled_4 = Z.of_int (4 * 16)

let rec ln x =
  let low_prec = -4 in
  let rough_appr = approx x low_prec in
  if Z.cmp rough_appr z_zero < 0 then invalid_arg "Cr.ln";
  if Z.cmp rough_appr low_ln_limit <= 0 then
    neg (ln (inv x))
  else if Z.cmp rough_appr high_ln_limit >= 0 then
    if Z.cmp rough_appr scaled_4 <= 0 then
      let quarter = ln (sqrt (sqrt x)) in
      shift_left quarter 2
    else
      let extra_bits = z_bit_length rough_appr - 3 in
      let scaled_result = ln (shift_right x extra_bits) in
      add scaled_result (mul (of_int extra_bits) ln2)
  else
    simple_ln x

let log ~base:x y = ln y /! ln x

let pow x y = exp (y *! ln x)

let root n x = pow x (inv (of_int n))

let arctan_reciproqual op =
  create
    (fun p ->
       if p >= 1 then 
	 z_zero
       else
	 let iterations_needed = -p/2 + 2 in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 2 in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let big_op = Z.of_int op in
	 let big_op_squared = Z.of_int (op*op) in
	 let op_inverse = Z.fdiv_q scaled_1 big_op in
	 let current_power = ref op_inverse in
	 let current_term = ref op_inverse in
	 let current_sum = ref op_inverse in
	 let current_sign = ref 1 in
	 let n = ref 1 in
	 let max_trunc_error = z_shift_left z_one (p - 2 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   n := !n + 2;
	   current_power := Z.fdiv_q !current_power big_op_squared;
	   current_sign := - !current_sign;
	   current_term := 
	     Z.fdiv_q !current_power (Z.of_int (!current_sign * !n));
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let pi = 
     (of_int 48 *! arctan_reciproqual 18)
  +! (of_int 32 *! arctan_reciproqual 57)
  -! (of_int 20 *! arctan_reciproqual 239)

let half_pi = shift_right pi 1

let prescaled_cos x =
  create
    (fun p ->
       if p >= 1 then
	 z_zero
       else begin
	 let iterations_needed = -p/2 + 4 in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 4 in
	 let op_prec = p - 2 in
	 let op_appr = approx x op_prec in
	 let current_term = ref (z_shift_left z_one (-calc_precision)) in
	 let n = ref 0 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 let current_sum = ref !current_term in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   n := !n + 2;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   let divisor = Z.mul (Z.of_int (- !n)) (Z.of_int (!n-1)) in
	   current_term := Z.fdiv_q !current_term divisor;
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p)
       end)

let rec cos x =
  let rough_appr = approx x (-1) in
  let abs_rough_appr = Z.abs rough_appr in
  if Z.cmp abs_rough_appr z_six >= 0 then
    let multiplier = Z.fdiv_q_ui rough_appr 6 in
    let adjustment = mul pi (of_z multiplier) in
    if Z.sgn (Z.band multiplier z_one) != 0 then
      neg (cos (x -! adjustment))
    else
      cos (x -! adjustment)
  else if Z.cmp abs_rough_appr z_two >= 0 then
    let cos_half = cos (shift_right x 1) in
    (shift_left (cos_half *! cos_half) 1) -! one
  else
    prescaled_cos x

let sin x = cos (half_pi -! x)

let tan x = sin x /! cos x

(*s Hyperbolic functions. *)

let sinh x = let expx = exp x in (expx -! inv expx) /! two

let cosh x = let expx = exp x in (expx +! inv expx) /! two

let tanh x = 
  let expx = exp x in
  let exp_minus_x = inv expx in
  (expx -! exp_minus_x) /! (expx +! exp_minus_x)

let arcsinh x = ln (x +! sqrt (x *! x +! one))
let arccosh x = ln (x +! sqrt (x *! x -! one))
let arctanh x = ln ((one +! x) /! (one -! x)) /! two

let of_float n =
  begin match classify_float n with
    | FP_nan | FP_infinite -> invalid_arg "Cr.of_float"
    | _ -> ()
  end;
  let negative = n < 0.0 in
  let bits = Int64.bits_of_float (abs_float n) in
  let mantissa = Int64.logand bits 0xfffffffffffffL in
  let biased_exp = Int64.to_int (Int64.shift_right_logical bits 52) in
  let exp = biased_exp - 1075 in
  let mantissa = 
    if biased_exp != 0 then
      Int64.add mantissa (Int64.shift_left Int64.one 52)
    else
      Int64.shift_left mantissa 1
  in
  let result = shift_left (of_int64 mantissa) exp in
  if negative then neg result else result

let to_string ?(radix=10) x n =
  if n < 0 then invalid_arg "Cr.to_string";
  let scaled_x =
    if radix == 16 then 
      shift_left x (4 * n)
    else
      let scale_factor = Z.ui_pow_ui radix n in
      mul x (of_z scale_factor)
  in
  let scaled_int = approx scaled_x 0 in
  let scaled_string = Z.to_string_base ~base:radix (Z.abs scaled_int) in
  let result = 
    if n == 0 then
      scaled_string
    else
      let len = String.length scaled_string in
      let len, scaled_string =
	if len <= n then
	  n + 1, String.make (n + 1 - len) '0' ^ scaled_string
	else
	  len, scaled_string
      in
      let whole = String.sub scaled_string 0 (len - n) in
      let fraction = String.sub scaled_string (len - n) n in
      whole ^ "." ^ fraction
  in
  if Z.sgn scaled_int < 0 then "-" ^ result else result

let of_string ?(radix=10) s =
  try
    begin
      try
	let n = String.length s in
	let p = String.index s '.' in
	let dec = n - p - 1 in
	let s' = (String.sub s 0 p) ^ (String.sub s (p + 1) dec) in
	div (of_z (Z.from_string_base radix s')) (of_z (Z.pow_ui_ui radix dec))
      with Not_found -> 
	of_z (Z.from_string_base radix s)
    end
  with Invalid_argument _ -> invalid_arg "Cr.of_string"

let to_q x n =
  let xn = approx x (-n) in
  Q.div (Q.from_z xn) (Q.from_z (z_shift_left z_one n))

let to_float x n = Q.float_from (to_q x n)

(* Inverse of a monotone function (see file UnaryCRFunction.java) *)

let sloppy_compare x y =
  let difference = Z.sub x y in
  if Z.cmp difference z_one > 0 then 1 
  else if Z.cmp difference z_mone < 0 then -1 else 0

let trace = false
let printf = Format.printf

let inverse_monotone f ~low ~high =
  let f_low = f low in
  let f_high = f high in
  let f,negated,f_low,f_high =
    if compare f_low f_high > 0 then
      (fun x -> neg (f x)), true, neg f_low, neg f_high
    else
      f, false, f_low, f_high
  in
  let max_msd = msd (max (abs low) (abs high)) in
  let max_arg_prec = msd (high -! low) - 4 in
  let deriv_msd = msd ((f_high -! f_low) /! (high -! low)) in
  fun x ->
    let arg = if negated then neg x else x in
    let rec r = { 
      cache = None;
      approximate = fun p ->
	let digits_needed = max_msd - p in
	if digits_needed < 0 then 
	  z_zero
	else
	  let extra_arg_prec = 4 in
	  let working_arg_prec = 
	    Pervasives.min (p - extra_arg_prec) max_arg_prec 
	  in
	  let working_eval_prec = ref (working_arg_prec + deriv_msd - 20) in
	  let low_appr = Z.add_ui (approx low working_arg_prec) 1 in
	  let high_appr = Z.sub_ui (approx high working_arg_prec) 1 in
	  let arg_appr = ref (approx arg !working_eval_prec) in
	  let have_good_appr = match r.cache with
	    | Some (min_prec, _) -> min_prec < max_msd
	    | None -> false
	  in
	  let small_steps = ref 0 in
	  let l,f_l,h,f_h,at_left,at_right =
	    if digits_needed < 30 && not have_good_appr then begin
	      if trace then printf "Setting interval to entire domain@.";
	      small_steps := 2;
	      low_appr, approx f_low !working_eval_prec,
	      high_appr, approx f_high !working_eval_prec,
	      true, true
	    end else
	      let rough_prec = p + digits_needed/2 in
	      let rough_prec = match r.cache with
		| Some (min_prec, _) when
                  digits_needed < 30 || min_prec<p+3*digits_needed/4 ->
		    min_prec
		| _ -> rough_prec
	      in
	      let rough_appr = approx r rough_prec in
	      if trace then begin
		printf "Setting interval based on prev. appr@.";
		printf "prev. prec = %d appr = %a@." 
		  rough_prec Z.print rough_appr
	      end;
	      let h = z_shift_left (Z.add_ui rough_appr 1) 
		(rough_prec - working_arg_prec)
	      in
	      let l = z_shift_left (Z.sub_ui rough_appr 1)
		(rough_prec - working_arg_prec)
	      in
	      let h,f_h,at_right = 
		if Z.cmp h high_appr > 0 then
		  high_appr, approx f_high !working_eval_prec, true
		else
		  let h_cr = shift_left (of_z h) working_arg_prec in
		  let f_h = approx (f h_cr) !working_eval_prec in
		  h, f_h, false
	      in
	      let l,f_l,at_left =
		if Z.cmp l low_appr < 0 then
		  low_appr, approx f_low !working_eval_prec, true
		else
		  let l_cr = shift_left (of_z l) working_arg_prec in
		  let f_l = approx (f l_cr) !working_eval_prec in
		  l, f_l, false
	      in
	      l,f_l,h,f_h,at_left,at_right
	  in
	  let l,f_l,h,f_h,at_left,at_right =
	    ref l, ref f_l, ref h, ref f_h, ref at_left, ref at_right
	  in
	  let difference = ref (Z.sub !h !l) in
	  let rec loop i =
	    if trace then begin
	      printf "***Iteration: %d@." i;
	      printf "Arg prec = %d eval prec = %d arg appr. = %a@."
		working_arg_prec !working_eval_prec Z.print !arg_appr;
	      printf "l = %a h = %a@." Z.print !l Z.print !h;
	      printf "f(l) = %a; f(h) = %a@." Z.print !f_l Z.print !f_h
	    end;
	    if Z.cmp !difference z_six < 0 then
	      scale !h (-extra_arg_prec)
	    else
	      let f_difference = Z.sub !f_h !f_l in
	      let guess = 
		if !small_steps >= 2 || Z.sgn f_difference = 0 then
		  z_shift_right (Z.add !l !h) 1 
		else
		  let arg_difference = Z.sub !arg_appr !f_l in
		  let t = Z.mul arg_difference !difference in
		  let adj = Z.fdiv_q t f_difference in
		  let adj =
		    if Z.cmp adj (z_shift_right !difference 2) < 0 then
		      z_shift_left adj 1
		    else if Z.cmp adj 
		      (z_shift_right (Z.mul_ui !difference 3) 2) > 0 then
			Z.sub !difference 
			  (z_shift_left (Z.sub !difference adj) 1)
		    else
		      adj
		  in
		  let adj = if Z.sgn adj <= 0 then z_two else adj in
		  let adj = 
		    if Z.cmp adj !difference >= 0 then
		      Z.sub_ui !difference 2 
		    else 
		      adj
		  in
		  if Z.sgn adj <= 0 then Z.add_ui !l 2 else Z.add !l adj 
	      in
	      let guess = ref guess in
	      let tweak = ref z_two in
	      let rec loop2 adj_prec =
		let guess_cr = shift_left (of_z !guess) working_arg_prec in
		if trace then begin
		  printf "Evaluating at %s with precision %d@."
		    (to_string guess_cr 10) !working_eval_prec;
		end;
		let f_guess_cr = f guess_cr in
		if trace then begin
		  printf "fn value = %s@." (to_string f_guess_cr 10)
		end;
		let f_guess = approx f_guess_cr !working_eval_prec in
		let outcome = sloppy_compare f_guess !arg_appr in
		if outcome <> 0 then begin (* break *) 
		  if outcome > 0 then begin
		    h := !guess; f_h := f_guess; at_right := false
		  end else begin
		    l := !guess; f_l := f_guess; at_left := false
		  end;
		  let new_difference = Z.sub !h !l in
		  if Z.cmp new_difference (z_shift_right !difference 1) >= 0
		  then incr small_steps
		  else small_steps := 0;
 		  difference := new_difference;
		  loop (i+1)
		end else begin
		  if adj_prec then begin
		    let adjustment = 
		      if deriv_msd > 0 then -20 else deriv_msd - 20
		    in
		    let l_cr = shift_left (of_z !l) working_arg_prec in
		    let h_cr = shift_left (of_z !h) working_arg_prec in
		    working_eval_prec := !working_eval_prec + adjustment;
		    if trace then begin
		      printf "New eval prec = %d %s%s@." !working_eval_prec
			(if !at_left then "(at left)" else "")
			(if !at_right then "(at right)" else "")
		    end;
		    if !at_left then
		      f_l := approx f_low !working_eval_prec
		    else
		      f_l := approx (f l_cr) !working_eval_prec;
		    if !at_right then
		      f_h := approx f_high !working_eval_prec
		    else
		      f_h := approx (f h_cr) !working_eval_prec;
		    arg_appr := approx arg !working_eval_prec;
		  end else begin
		    if trace then printf "tweaking guess@.";
		    let new_guess = Z.add !guess !tweak in
		    if Z.cmp new_guess !h >= 0 then
		      guess := Z.sub !guess !tweak
		    else
		      guess := new_guess;
		    tweak := Z.neg !tweak
		  end;
		  loop2 (not adj_prec)
		end
	      in
	      loop2 false
	  in
	  loop 0
    }
    in
    r

(* Application to the inverse trigonometric functions *)

let arcsin = inverse_monotone sin ~low:(neg half_pi) ~high:half_pi

let arccos x = half_pi -! arcsin x

(* uses the identity [(sin x)^2 = (tan x)^2/(1 + (tan x)^2)] *)
let arctan x =
  let x2 = x *! x in
  let abs_sin_atan = sqrt (x2 /! (one +! x2)) in
  let sin_atan = select x (neg abs_sin_atan) abs_sin_atan in
  arcsin sin_atan

(*s Format pretty-printer. *)

let print_precision = ref 10
let set_print_precision = (:=) print_precision
let print fmt x = Format.fprintf fmt "%s" (to_string x !print_precision)

(* Infix notations *)

module Infixes = struct
  let (+!) = add
  let (-!) = sub
  let ( *! ) = mul
  let (/!) = div
end