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 ;;

Answer 1

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

Answer 2

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.