[racket] q about code for partitions

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)
     (printf "bad: ~v~n" n)))


I haven't had the patience to wait for (partitions 100000) to finish, 
let alone (partitions 1000000). If your implementation is faster on 
large inputs, we're interested in a patch.

Neil ⊥