When arithmetic on a computer bite back

92 views
Skip to first unread message

kawas

unread,
Jun 1, 2012, 12:19:53 PM6/1/12
to Clojure
Hi,

Can someone explain to me the behavior of this function when applied
to different kind of numbers.

(defn prime-factors [n]
(loop [f 2 n n res []]
(cond
(<= n 1) res
(zero? (rem n f)) (recur f (quot n f) (conj res f))
:else (recur (inc f) n res))))

Problem 1 (solved): If you use (= n 1) in the first cond clause, the
function may not terminate because (= 1 1.0) may be false. But we all
know that :)

Problem 2: Now compare theses function calls
(prime-factors (Math/pow 2 58))
(prime-factors (bigdec (Math/pow 2 58)))
(prime-factors (bigint (bigdec (Math/pow 2 58))))

The last two function calls are really really slow.
- Is this only a problem of rem/quot slow on big numbers ?

-
Laurent

Tassilo Horn

unread,
Jun 1, 2012, 1:21:24 PM6/1/12
to clo...@googlegroups.com
kawas <laurent...@gmail.com> writes:

Hi Laurent,

> (defn prime-factors [n]
> (loop [f 2 n n res []]
> (cond
> (<= n 1) res
> (zero? (rem n f)) (recur f (quot n f) (conj res f))
> :else (recur (inc f) n res))))
>
> Problem 1 (solved): If you use (= n 1) in the first cond clause, the
> function may not terminate because (= 1 1.0) may be false. But we all
> know that :)

Use (== n 1) instead.

> Problem 2: Now compare theses function calls
> (prime-factors (Math/pow 2 58))
> (prime-factors (bigdec (Math/pow 2 58)))
> (prime-factors (bigint (bigdec (Math/pow 2 58))))
>
> The last two function calls are really really slow.
> - Is this only a problem of rem/quot slow on big numbers ?

No, debugging it a bit it seems to be a bug in
clojure.lang.Numbers.quotient() applied to BigDecs:

user> (quot 1.4411518807585587E17 2) ;; correct with doubles
7.2057594037927936E16
user> (quot 1.4411518807585587E+17M 2) ;; wrong with BigDecs
72057594037927935M

Please file a bug report about that.

Bye,
Tassilo

kawas

unread,
Jun 1, 2012, 7:53:10 PM6/1/12
to Clojure
I've checked with a python repl, the correct value seems to be the one
using BigDecs

>>> from decimal import *
>>> Decimal('1.4411518807585587E17') / Decimal(2)
Decimal('72057594037927935')

I've submitted a bug and commented on it about BigDecimal constructors
but I may have made a mistake about which value is the correct one...

Bye,
Laurent

On 1 juin, 19:21, Tassilo Horn <tass...@member.fsf.org> wrote:

Tassilo Horn

unread,
Jun 2, 2012, 5:31:24 AM6/2/12
to clo...@googlegroups.com
kawas <laurent...@gmail.com> writes:

> I've checked with a python repl, the correct value seems to be the one
> using BigDecs
>
>>>> from decimal import *
>>>> Decimal('1.4411518807585587E17') / Decimal(2)
> Decimal('72057594037927935')
>
> I've submitted a bug and commented on it about BigDecimal constructors
> but I may have made a mistake about which value is the correct one...

Yes, you did. How can a power of two divided by two be *odd* (well,
except for 2/2 = 1, of course)? ;-)

But strangely, clojure is wrong, python is wrong, and ruby is wrong,
too.

--8<---------------cut here---------------start------------->8---
irb(main):012:0> x = BigDecimal.new("1.4411518807585587E17")
=> #<BigDecimal:7fabd400d958,'0.1441151880 7585587E18',18(27)>
irb(main):013:0> x
=> #<BigDecimal:7fabd400d958,'0.1441151880 7585587E18',18(27)>
irb(main):014:0> x / 2
=> #<BigDecimal:7fabd4008368,'0.7205759403 7927935E17',18(81)>
--8<---------------cut here---------------end--------------->8---

So maybe that's a general bigdecimal representation issue?

Bye,
Tassilo

Mark Engelberg

unread,
Jun 2, 2012, 7:01:16 AM6/2/12
to clo...@googlegroups.com
1.4411518807585587E17 ends in 0, and therefore when you divide by 2, it should end in 5.
It's not a power of 2, it is a merely an inexact approximation of one.

Tassilo Horn

unread,
Jun 3, 2012, 4:29:22 AM6/3/12
to clo...@googlegroups.com
Mark Engelberg <mark.en...@gmail.com> writes:

> 1.4411518807585587E17 ends in 0,

Oh, my counting was bad yesterday.

> and therefore when you divide by 2, it should end in 5. It's not a
> power of 2, it is a merely an inexact approximation of one.

Yes.

Bye,
Tassilo
Reply all
Reply to author
Forward
0 new messages