Re: [Rd] choose(n,k) when n is almost integer

From: Petr Savicky <savicky_at_cs.cas.cz>
Date: Wed, 23 Dec 2009 23:42:23 +0100

In a previous email, i suggested two patches A and B to choose(n, k), which solve some of its problems, but keep some of the inaccuracies of the original implementation. I would like to suggest another patch, which i will call C (patch-C.txt in an attachment), which eliminates the warnings obtained sometimes from the original implementation and which is more accurate in all ranges of the output.

For testing patch C a simpler script is sufficient, since we need not to take care of the warnings. Namely
  http://www.cs.cas.cz/~savicky/R-devel/test_choose_1.R which produces the output (the error is the relative error)

> source("test_choose_1.R")

  k <= 0 max err = 0 
  k <= 10 max err = 1.332268e-15 
  k <= 20 max err = 2.442491e-15 
  k <= 30 max err = 3.774758e-15 
  k <= 40 max err = 2.553513e-14 
  k <= 50 max err = 2.88658e-14 
  k <= 60 max err = 3.197442e-14 
  k <= 70 max err = 4.396483e-14 
  k <= 80 max err = 4.685141e-14 
  k <= 90 max err = 4.907186e-14 
  k <= 100 max err = 5.084821e-14 
  k <= 110 max err = 5.373479e-14 
  k <= 120 max err = 5.551115e-14 
  k <= 130 max err = 7.41629e-14 
  k <= 140 max err = 9.592327e-14 
  k <= 150 max err = 9.636736e-14 
  k <= 160 max err = 9.725554e-14 
  k <= 170 max err = 9.947598e-14 
  k <= 180 max err = 1.04583e-13 
  k <= 190 max err = 1.088019e-13 
  k <= 200 max err = 1.090239e-13 

  minimum log2() of a wrong result for integer n : 53.32468   maximum error for real n : 1.090239e-13

Increasing accuracy of choose(n, k) for n almost an integer needed to use additional transformations of it to those already used in the code. I will work out a description of these transformations and send a link to it. Similarly as patches A and B, patch C also does not modify lchoose().

It should be pointed out that choose(n, k) for non-integer n is mainly needed if n is a rational number like 1/2, 1/3, 2/3, .... However, making choose(n, k) accurate for all inputs seems to be not too complex as the patch C and its test results show.

I appreciate comments on patch C.

Petr Savicky.



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Wed 23 Dec 2009 - 22:50:15 GMT

This archive was generated by hypermail 2.2.0 : Thu 24 Dec 2009 - 00:41:13 GMT