Calculating Euler's constant with OCaml

Question

I started learning OCaml today. I already knew Scheme so I thought it would be good to try to translate some example Scheme programs to ML. I have some code below for calculating Euler's number which works, but it does not work in OCaml. I get this error: Exception: Division_by_zero. I think that this might be a problem with mixing floats and ints somewhere. I have not been able to figure out how to set breakpoints for functions with ocamldebug. Can anyone identify where my mistake is happening? Thank you.

(define (factorial n)
    (if (zero? n) 1 (* n (factorial (sub1 n)))))

(define (find-e accuracy sum)
    (if (zero? accuracy) (add1 sum)
        (find-e (sub1 accuracy) (+ sum (/ 1 (factorial accuracy))))))

(display (format "~f" (find-e 100 0)))

let rec factorial n = if n == 0 then 1 else n * factorial (n - 1) ;;

let rec find_e accuracy sum =
    if accuracy == 0
        then (sum + 1)
        else find_e (accuracy - 1) (sum + (1 / factorial accuracy)) ;;

let result = find_e 100 0 ;;

Jeffrey Scofield · Accepted Answer · 2020-11-24T18:37:05.307

As I recall, scheme has a "numerical tower" that tries to keep you from losing accuracy in numeric computations.

OCaml doesn't have any fancy automatic handling of numbers. Your code is using type int, which is a fixed size integer (31 or 63 bits in the usual implementations). Thus, your expression 1 / factorial accuracy will be 0 in almost all cases and values for factorial accuracy will be unrepresentable for all but the smallest values. The value of factorial 100 will be 0 because it is a multiple of 2^63:

# let rec fact n = if n < 2 then 1 else n * fact (n - 1);;
val fact : int -> int = <fun>
# fact 100;;
- : int = 0

There are no floats in your code, and hence there couldn't possibly be any mixing. But OCaml doesn't allow mixing in the first place. It's a strongly typed language where int and float are two different types.

Here is your code converted to use floats:

let rec factorial n =
    if n = 0.0 then 1.0 else n *. factorial (n -. 1.0) 

let rec find_e accuracy sum =
    if accuracy = 0.0 then
        sum +. 1.0
    else
        find_e (accuracy -. 1.0) (sum +. 1.0 /. factorial accuracy)

let result = find_e 100.0 0.0

If I copy/paste this into the OCaml REPL (also known as "toplevel") I see this:

val factorial : float -> float = <fun>
val find_e : float -> float -> float = <fun>
val result : float = 2.71828182845904509

As a side comment, the equality comparison operator in OCaml is =. Don't use == for comparisons. It's a completely different operator:

# 1.0 == 1.0;;
- : bool = false

do you mean that I should change each instance of `int` to `float`? — Caspian Ahlberg, Nov 22 '20 at 19:12
This will work better, but 100 factorial isn't representable as a float either. You would also need to change your operators to be ones that work with float. I don't think this is homework :-) so I'll write up a version that uses float. — Jeffrey Scofield, Nov 22 '20 at 19:14
Does that mean that there are alternative operators just for floats? — Caspian Ahlberg, Nov 22 '20 at 19:16
Yes, OCaml is strongly typed. A given function (like `+`) works only with values of one type (`int` in this case). The operator for adding floats is `+.` (with trailing period or decimal point). — Jeffrey Scofield, Nov 22 '20 at 19:21
(I can't edit my comment above, but 100 factorial can actually be represented pretty well as a float. OCaml floats (which are actually double precision) can represent a wider range than I was thinking.) — Jeffrey Scofield, Nov 22 '20 at 19:47

coredump · Answer 2 · 2020-11-24T14:59:11.270

This answer is meant as a complement to the currently accepted one.

Since factorials are integer values, it is somewhat unfortunate to use floats to represent the computation, as it introduces a loss of precision for large enough values. Scheme and Lisp have a transparent support for bignums and ratios, but OCaml has a library for that.

opam install zarith

Then in your ocaml toplevel:

# #use "topfind";;
# #require "zarith";;
/home/user/.opam/4.11.0/lib/zarith: added to search path
/home/user/.opam/4.11.0/lib/zarith/zarith.cma: loaded

There are two main modules, Z and Q (integers and rationals), as well as a compatibility module for an older big integers module. You probably want to have a look at the documentation: https://antoinemine.github.io/Zarith/doc/latest/Z.html

The library define its own type for integers:

# Z.of_int(3);;
- : Z.t = <abstr>

If you install the pretty-printer:

# #install_printer Z.pp_print;;
# Z.of_int(3);;
- : Z.t = 3

Likewise, the Q module can help when computing fractions. You probably can implement an equivalent version of your code in OCaml that way.

Calculating Euler's constant with OCaml

2 Answers2