# [racket] q about code for partitions

 From: Neil Toronto (neil.toronto at gmail.com) Date: Wed Jun 4 17:11:09 EDT 2014 Previous message: [racket] q about code for partitions Next message: [racket] Adding a language to the Choose Language... menu Messages sorted by: [date] [thread] [subject] [author]

```On 06/04/2014 12:36 PM, Jos Koot wrote:
> Hi
>
> In share/pkgs/math-lib/math/private/number-theory I find the following two
> pieces of code:
>
> line 34: (define m (/ (+ 1.0 (flsqrt (+ 1.0 (* 24.0 n)))) 6.0))
> line 39: (exact-floor m)
>
> Obviously for finding the positive root of the equation n-k(3k-1)/2=0 for k
> in terms of n.
> (from wikipedia: http://en.wikipedia.org/wiki/Partition_(number_theory) )
>
> How can we be sure that (exact-floor m) never is too small?
> Because of the inexact operations, might line 34 not produce a value off by
> more than 1 for very large values of n?
> No criticism here, just wondering.

No problem! I don't think we've tested this before.

The following program tests the computation of `m` against bigfloat
versions of it with rounding up and down. (All the operations are
monotone increasing, so the bigfloat versions give upper and lower
bounds.) The plot suggests that `m` has less than 1 ulp error, at least
up to 1000, and usually has < 0.5 ulp error (i.e. just rounding error).

#lang racket

(require plot math)

(define (m n)
(/ (+ 1.0 (flsqrt (+ 1.0 (* 24.0 (fl n))))) 6.0))

(define (m-up n)
(parameterize ([bf-rounding-mode  'up])
(/ (+ 1 (bigfloat->real (bfsqrt (bf (+ 1 (* 24 n)))))) 6)))

(define (m-down n)
(parameterize ([bf-rounding-mode  'zero])
(/ (+ 1 (bigfloat->real (bfsqrt (bf (+ 1 (* 24 n)))))) 6)))

(plot (list
(function (λ (x)
(define n (floor x))
(flulp-error (m n) (m-up n))))
(function (λ (x)
(define n (floor x))
(flulp-error (m n) (m-down n)))
#:color 3))
#:x-min 0 #:x-max 1000)

This isn't too surprising: (+ 1.0 (* 24.0 (fl n))) is exact for
approximately n < 2^52, the remaining operations add up to 0.5 ulp error
each, adding 1.0 is exact except just below powers of 2, and some errors
cancel.

But that's not proof: there's about a 2^-52 probability that flooring
produces an off-by-one result. So here's an exhaustive proof for `n` in
[0..999999]:

(for ([n  (in-range 0 1000000)])
(when (= 0 (modulo (+ n 1) 1000))
(printf "n = ~v~n" n))
(define i (exact-floor (m n)))
(define i-up (exact-floor (m-up n)))
(define i-down (exact-floor (m-down n)))
(unless (= i i-up i-down)