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 ⊥
____________________
Racket Users list:
http://lists.racket-lang.org/users
