Mailman 3 September 2009 - Python-checkins

r75123 - in python/trunk: Lib/test/test_threading_local.py Misc/NEWS Modules/threadmodule.c
by philip.jenvey 29 Sep '09

29 Sep '09

Author: philip.jenvey Date: Tue Sep 29 06:32:44 2009 New Revision: 75123 Log: #6990: clear threading.local's key only after its thread state is removed: fixes local subclasses leaving old state around after a ref cycle GC which could be recycled by new locals Modified: python/trunk/Lib/test/test_threading_local.py python/trunk/Misc/NEWS python/trunk/Modules/threadmodule.c Modified: python/trunk/Lib/test/test_threading_local.py ============================================================================== --- python/trunk/Lib/test/test_threading_local.py (original) +++ python/trunk/Lib/test/test_threading_local.py Tue Sep 29 06:32:44 2009 @@ -67,6 +67,44 @@ for t in threads: t.join() + def test_derived_cycle_dealloc(self): + # http://bugs.python.org/issue6990 + class Local(threading.local): + pass + locals = None + passed = [False] + e1 = threading.Event() + e2 = threading.Event() + + def f(): + # 1) Involve Local in a cycle + cycle = [Local()] + cycle.append(cycle) + cycle[0].foo = 'bar' + + # 2) GC the cycle (triggers threadmodule.c::local_clear + # before local_dealloc) + del cycle + gc.collect() + e1.set() + e2.wait() + + # 4) New Locals should be empty + passed[0] = all(not hasattr(local, 'foo') for local in locals) + + t = threading.Thread(target=f) + t.start() + e1.wait() + + # 3) New Locals should recycle the original's address. Creating + # them in the thread overwrites the thread state and avoids the + # bug + locals = [Local() for i in range(10)] + e2.set() + t.join() + + self.assertTrue(passed[0]) + def test_main(): suite = unittest.TestSuite() Modified: python/trunk/Misc/NEWS ============================================================================== --- python/trunk/Misc/NEWS (original) +++ python/trunk/Misc/NEWS Tue Sep 29 06:32:44 2009 @@ -12,6 +12,9 @@ Core and Builtins ----------------- +- Issue #6990: Fix threading.local subclasses leaving old state around + after a reference cycle GC which could be recycled by new locals. + - Issue #6300: unicode.encode, unicode.decode, str.decode, and str.encode now take keyword arguments. Modified: python/trunk/Modules/threadmodule.c ============================================================================== --- python/trunk/Modules/threadmodule.c (original) +++ python/trunk/Modules/threadmodule.c Tue Sep 29 06:32:44 2009 @@ -244,7 +244,6 @@ static int local_clear(localobject *self) { - Py_CLEAR(self->key); Py_CLEAR(self->args); Py_CLEAR(self->kw); Py_CLEAR(self->dict); @@ -266,6 +265,7 @@ PyDict_DelItem(tstate->dict, self->key); } + Py_XDECREF(self->key); local_clear(self); Py_TYPE(self)->tp_free((PyObject*)self); }

1 0

Python Regression Test Failures refleak (1)
by Neal Norwitz 29 Sep '09

29 Sep '09

More important issues: ---------------------- test_ssl leaked [0, 0, -81] references, sum=-81 Less important issues: ---------------------- test_cmd_line leaked [0, 0, -25] references, sum=-25 test_file2k leaked [0, 77, -77] references, sum=0 test_popen2 leaked [-25, 29, 25] references, sum=29 test_smtplib leaked [0, 0, 88] references, sum=88 test_threadedtempfile leaked [-102, 0, 0] references, sum=-102 test_threading leaked [48, 48, 48] references, sum=144 test_threadsignals leaked [0, 0, -8] references, sum=-8 test_xmlrpc leaked [0, -6, 6] references, sum=0

1 0

r75122 - in python/branches/release26-maint: Lib/test/test_urllib2_localnet.py
by r.david.murray 29 Sep '09

29 Sep '09

Author: r.david.murray Date: Mon Sep 28 21:53:25 2009 New Revision: 75122 Log: Merged revisions 75111 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r75111 | r.david.murray | 2009-09-28 12:57:45 -0400 (Mon, 28 Sep 2009) | 5 lines Prevent test_bad_address failure when a domain in the dns search path implements a '*' default rule. Also update comment with a more complete explanation of the difficulties inherent in the test. ........ Modified: python/branches/release26-maint/ (props changed) python/branches/release26-maint/Lib/test/test_urllib2_localnet.py Modified: python/branches/release26-maint/Lib/test/test_urllib2_localnet.py ============================================================================== --- python/branches/release26-maint/Lib/test/test_urllib2_localnet.py (original) +++ python/branches/release26-maint/Lib/test/test_urllib2_localnet.py Mon Sep 28 21:53:25 2009 @@ -467,14 +467,24 @@ # Make sure proper exception is raised when connecting to a bogus # address. self.assertRaises(IOError, - # SF patch 809915: In Sep 2003, VeriSign started - # highjacking invalid .com and .net addresses to - # boost traffic to their own site. This test - # started failing then. One hopes the .invalid - # domain will be spared to serve its defined - # purpose. - # urllib2.urlopen, "http://www.sadflkjsasadf.com/") - urllib2.urlopen, "http://sadflkjsasf.i.nvali.d/") + # Given that both VeriSign and various ISPs have in + # the past or are presently hijacking various invalid + # domain name requests in an attempt to boost traffic + # to their own sites, finding a domain name to use + # for this test is difficult. RFC2606 leads one to + # believe that '.invalid' should work, but experience + # seemed to indicate otherwise. Single character + # TLDs are likely to remain invalid, so this seems to + # be the best choice. The trailing '.' prevents a + # related problem: The normal DNS resolver appends + # the domain names from the search path if there is + # no '.' the end and, and if one of those domains + # implements a '*' rule a result is returned. + # However, none of this will prevent the test from + # failing if the ISP hijacks all invalid domain + # requests. The real solution would be to be able to + # parameterize the framework with a mock resolver. + urllib2.urlopen, "http://sadflkjsasf.i.nvali.d./") def test_main():

1 0

Python Regression Test Failures refleak (1)
by Neal Norwitz 29 Sep '09

29 Sep '09

More important issues: ---------------------- test_ssl leaked [26, -26, -420] references, sum=-420 Less important issues: ---------------------- test_asynchat leaked [0, 126, -126] references, sum=0 test_docxmlrpc leaked [0, -3, 11] references, sum=8 test_popen2 leaked [0, 25, 4] references, sum=29 test_smtplib leaked [4, 0, 0] references, sum=4 test_sys leaked [42, 0, 0] references, sum=42 test_threading leaked [48, 48, 48] references, sum=144

1 0

r75121 - python/branches/release31-maint
by mark.dickinson 29 Sep '09

29 Sep '09

Author: mark.dickinson Date: Mon Sep 28 21:23:20 2009 New Revision: 75121 Log: Blocked revisions 75119 via svnmerge ................ r75119 | mark.dickinson | 2009-09-28 20:21:11 +0100 (Mon, 28 Sep 2009) | 10 lines Merged revisions 75117 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r75117 | mark.dickinson | 2009-09-28 19:54:55 +0100 (Mon, 28 Sep 2009) | 3 lines Issue #3366: Add gamma function to math module. (lgamma, erf and erfc to follow). ........ ................ Modified: python/branches/release31-maint/ (props changed)

1 0

r75120 - python/branches/release26-maint
by mark.dickinson 29 Sep '09

29 Sep '09

Author: mark.dickinson Date: Mon Sep 28 21:22:34 2009 New Revision: 75120 Log: Blocked revisions 75117 via svnmerge ........ r75117 | mark.dickinson | 2009-09-28 19:54:55 +0100 (Mon, 28 Sep 2009) | 3 lines Issue #3366: Add gamma function to math module. (lgamma, erf and erfc to follow). ........ Modified: python/branches/release26-maint/ (props changed)

1 0

r75119 - in python/branches/py3k: Doc/library/math.rst Lib/test/math_testcases.txt Lib/test/test_math.py Misc/NEWS Modules/mathmodule.c
by mark.dickinson 29 Sep '09

29 Sep '09

Author: mark.dickinson Date: Mon Sep 28 21:21:11 2009 New Revision: 75119 Log: Merged revisions 75117 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r75117 | mark.dickinson | 2009-09-28 19:54:55 +0100 (Mon, 28 Sep 2009) | 3 lines Issue #3366: Add gamma function to math module. (lgamma, erf and erfc to follow). ........ Added: python/branches/py3k/Lib/test/math_testcases.txt Modified: python/branches/py3k/ (props changed) python/branches/py3k/Doc/library/math.rst python/branches/py3k/Lib/test/test_math.py python/branches/py3k/Misc/NEWS python/branches/py3k/Modules/mathmodule.c Modified: python/branches/py3k/Doc/library/math.rst ============================================================================== --- python/branches/py3k/Doc/library/math.rst (original) +++ python/branches/py3k/Doc/library/math.rst Mon Sep 28 21:21:11 2009 @@ -278,6 +278,16 @@ Return the hyperbolic tangent of *x*. +Special functions +----------------- + +.. function:: gamma(x) + + Return the Gamma function at *x*. + + .. versionadded:: 2.7 + + Constants --------- Added: python/branches/py3k/Lib/test/math_testcases.txt ============================================================================== --- (empty file) +++ python/branches/py3k/Lib/test/math_testcases.txt Mon Sep 28 21:21:11 2009 @@ -0,0 +1,146 @@ +-- Testcases for functions in math. +-- +-- Each line takes the form: +-- +-- <testid> <function> <input_value> -> <output_value> <flags> +-- +-- where: +-- +-- <testid> is a short name identifying the test, +-- +-- <function> is the function to be tested (exp, cos, asinh, ...), +-- +-- <input_value> is a string representing a floating-point value +-- +-- <output_value> is the expected (ideal) output value, again +-- represented as a string. +-- +-- <flags> is a list of the floating-point flags required by C99 +-- +-- The possible flags are: +-- +-- divide-by-zero : raised when a finite input gives a +-- mathematically infinite result. +-- +-- overflow : raised when a finite input gives a finite result that +-- is too large to fit in the usual range of an IEEE 754 double. +-- +-- invalid : raised for invalid inputs (e.g., sqrt(-1)) +-- +-- ignore-sign : indicates that the sign of the result is +-- unspecified; e.g., if the result is given as inf, +-- then both -inf and inf should be accepted as correct. +-- +-- Flags may appear in any order. +-- +-- Lines beginning with '--' (like this one) start a comment, and are +-- ignored. Blank lines, or lines containing only whitespace, are also +-- ignored. + +-- Many of the values below were computed with the help of +-- version 2.4 of the MPFR library for multiple-precision +-- floating-point computations with correct rounding. All output +-- values in this file are (modulo yet-to-be-discovered bugs) +-- correctly rounded, provided that each input and output decimal +-- floating-point value below is interpreted as a representation of +-- the corresponding nearest IEEE 754 double-precision value. See the +-- MPFR homepage at http://www.mpfr.org for more information about the +-- MPFR project. + +--------------------------- +-- gamma: Gamma function -- +--------------------------- + +-- special values +gam0000 gamma 0.0 -> inf divide-by-zero +gam0001 gamma -0.0 -> -inf divide-by-zero +gam0002 gamma inf -> inf +gam0003 gamma -inf -> nan invalid +gam0004 gamma nan -> nan + +-- negative integers inputs are invalid +gam0010 gamma -1 -> nan invalid +gam0011 gamma -2 -> nan invalid +gam0012 gamma -1e16 -> nan invalid +gam0013 gamma -1e300 -> nan invalid + +-- small positive integers give factorials +gam0020 gamma 1 -> 1 +gam0021 gamma 2 -> 1 +gam0022 gamma 3 -> 2 +gam0023 gamma 4 -> 6 +gam0024 gamma 5 -> 24 +gam0025 gamma 6 -> 120 + +-- half integers +gam0030 gamma 0.5 -> 1.7724538509055161 +gam0031 gamma 1.5 -> 0.88622692545275805 +gam0032 gamma 2.5 -> 1.3293403881791370 +gam0033 gamma 3.5 -> 3.3233509704478426 +gam0034 gamma -0.5 -> -3.5449077018110322 +gam0035 gamma -1.5 -> 2.3632718012073548 +gam0036 gamma -2.5 -> -0.94530872048294190 +gam0037 gamma -3.5 -> 0.27008820585226911 + +-- values near 0 +gam0040 gamma 0.1 -> 9.5135076986687306 +gam0041 gamma 0.01 -> 99.432585119150602 +gam0042 gamma 1e-8 -> 99999999.422784343 +gam0043 gamma 1e-16 -> 10000000000000000 +gam0044 gamma 1e-30 -> 9.9999999999999988e+29 +gam0045 gamma 1e-160 -> 1.0000000000000000e+160 +gam0046 gamma 1e-308 -> 1.0000000000000000e+308 +gam0047 gamma 5.6e-309 -> 1.7857142857142848e+308 +gam0048 gamma 5.5e-309 -> inf overflow +gam0049 gamma 1e-309 -> inf overflow +gam0050 gamma 1e-323 -> inf overflow +gam0051 gamma 5e-324 -> inf overflow +gam0060 gamma -0.1 -> -10.686287021193193 +gam0061 gamma -0.01 -> -100.58719796441078 +gam0062 gamma -1e-8 -> -100000000.57721567 +gam0063 gamma -1e-16 -> -10000000000000000 +gam0064 gamma -1e-30 -> -9.9999999999999988e+29 +gam0065 gamma -1e-160 -> -1.0000000000000000e+160 +gam0066 gamma -1e-308 -> -1.0000000000000000e+308 +gam0067 gamma -5.6e-309 -> -1.7857142857142848e+308 +gam0068 gamma -5.5e-309 -> -inf overflow +gam0069 gamma -1e-309 -> -inf overflow +gam0070 gamma -1e-323 -> -inf overflow +gam0071 gamma -5e-324 -> -inf overflow + +-- values near negative integers +gam0080 gamma -0.99999999999999989 -> -9007199254740992.0 +gam0081 gamma -1.0000000000000002 -> 4503599627370495.5 +gam0082 gamma -1.9999999999999998 -> 2251799813685248.5 +gam0083 gamma -2.0000000000000004 -> -1125899906842623.5 +gam0084 gamma -100.00000000000001 -> -7.5400833348831090e-145 +gam0085 gamma -99.999999999999986 -> 7.5400833348840962e-145 + +-- large inputs +gam0100 gamma 170 -> 4.2690680090047051e+304 +gam0101 gamma 171 -> 7.2574156153079990e+306 +gam0102 gamma 171.624 -> 1.7942117599248104e+308 +gam0103 gamma 171.625 -> inf overflow +gam0104 gamma 172 -> inf overflow +gam0105 gamma 2000 -> inf overflow +gam0106 gamma 1.7e308 -> inf overflow + +-- inputs for which gamma(x) is tiny +gam0120 gamma -100.5 -> -3.3536908198076787e-159 +gam0121 gamma -160.5 -> -5.2555464470078293e-286 +gam0122 gamma -170.5 -> -3.3127395215386074e-308 +gam0123 gamma -171.5 -> 1.9316265431711902e-310 +gam0124 gamma -176.5 -> -1.1956388629358166e-321 +gam0125 gamma -177.5 -> 4.9406564584124654e-324 +gam0126 gamma -178.5 -> -0.0 +gam0127 gamma -179.5 -> 0.0 +gam0128 gamma -201.0001 -> 0.0 +gam0129 gamma -202.9999 -> -0.0 +gam0130 gamma -1000.5 -> -0.0 +gam0131 gamma -1000000000.3 -> -0.0 +gam0132 gamma -4503599627370495.5 -> 0.0 + +-- inputs that cause problems for the standard reflection formula, +-- thanks to loss of accuracy in 1-x +gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88 +gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214 Modified: python/branches/py3k/Lib/test/test_math.py ============================================================================== --- python/branches/py3k/Lib/test/test_math.py (original) +++ python/branches/py3k/Lib/test/test_math.py Mon Sep 28 21:21:11 2009 @@ -7,6 +7,7 @@ import os import sys import random +import struct eps = 1E-05 NAN = float('nan') @@ -29,8 +30,50 @@ else: file = __file__ test_dir = os.path.dirname(file) or os.curdir +math_testcases = os.path.join(test_dir, 'math_testcases.txt') test_file = os.path.join(test_dir, 'cmath_testcases.txt') +def to_ulps(x): + """Convert a non-NaN float x to an integer, in such a way that + adjacent floats are converted to adjacent integers. Then + abs(ulps(x) - ulps(y)) gives the difference in ulps between two + floats. + + The results from this function will only make sense on platforms + where C doubles are represented in IEEE 754 binary64 format. + + """ + n = struct.unpack('q', struct.pack('<d', x))[0] + if n < 0: + n = ~(n+2**63) + return n + + +def parse_mtestfile(fname): + """Parse a file with test values + + -- starts a comment + blank lines, or lines containing only a comment, are ignored + other lines are expected to have the form + id fn arg -> expected [flag]* + + """ + with open(fname) as fp: + for line in fp: + # strip comments, and skip blank lines + if '--' in line: + line = line[:line.index('--')] + if not line.strip(): + continue + + lhs, rhs = line.split('->') + id, fn, arg = lhs.split() + rhs_pieces = rhs.split() + exp = rhs_pieces[0] + flags = rhs_pieces[1:] + + yield (id, fn, float(arg), float(exp), flags) + def parse_testfile(fname): """Parse a file with test values @@ -884,6 +927,51 @@ self.fail(message) self.ftest("%s:%s(%r)" % (id, fn, ar), result, er) + @unittest.skipUnless(float.__getformat__("double").startswith("IEEE"), + "test requires IEEE 754 doubles") + def test_mtestfile(self): + ALLOWED_ERROR = 20 # permitted error, in ulps + fail_fmt = "{}:{}({!r}): expected {!r}, got {!r}" + + failures = [] + for id, fn, arg, expected, flags in parse_mtestfile(math_testcases): + func = getattr(math, fn) + + if 'invalid' in flags or 'divide-by-zero' in flags: + expected = 'ValueError' + elif 'overflow' in flags: + expected = 'OverflowError' + + try: + got = func(arg) + except ValueError: + got = 'ValueError' + except OverflowError: + got = 'OverflowError' + + diff_ulps = None + if isinstance(got, float) and isinstance(expected, float): + if math.isnan(expected) and math.isnan(got): + continue + if not math.isnan(expected) and not math.isnan(got): + diff_ulps = to_ulps(expected) - to_ulps(got) + if diff_ulps <= ALLOWED_ERROR: + continue + + if isinstance(got, str) and isinstance(expected, str): + if got == expected: + continue + + fail_msg = fail_fmt.format(id, fn, arg, expected, got) + if diff_ulps is not None: + fail_msg += ' ({} ulps)'.format(diff_ulps) + failures.append(fail_msg) + + if failures: + self.fail('Failures in test_mtestfile:\n ' + + '\n '.join(failures)) + + def test_main(): from doctest import DocFileSuite suite = unittest.TestSuite() Modified: python/branches/py3k/Misc/NEWS ============================================================================== --- python/branches/py3k/Misc/NEWS (original) +++ python/branches/py3k/Misc/NEWS Mon Sep 28 21:21:11 2009 @@ -201,6 +201,9 @@ Extension Modules ----------------- + +- Issue #3366: Add gamma function to math module. + - Issue #6877: It is now possible to link the readline extension to the libedit readline emulation on OSX 10.5 or later. Modified: python/branches/py3k/Modules/mathmodule.c ============================================================================== --- python/branches/py3k/Modules/mathmodule.c (original) +++ python/branches/py3k/Modules/mathmodule.c Mon Sep 28 21:21:11 2009 @@ -60,44 +60,265 @@ extern double copysign(double, double); #endif -/* Call is_error when errno != 0, and where x is the result libm - * returned. is_error will usually set up an exception and return - * true (1), but may return false (0) without setting up an exception. - */ -static int -is_error(double x) +/* + sin(pi*x), giving accurate results for all finite x (especially x + integral or close to an integer). This is here for use in the + reflection formula for the gamma function. It conforms to IEEE + 754-2008 for finite arguments, but not for infinities or nans. +*/ + +static const double pi = 3.141592653589793238462643383279502884197; + +static double +sinpi(double x) { - int result = 1; /* presumption of guilt */ - assert(errno); /* non-zero errno is a precondition for calling */ - if (errno == EDOM) - PyErr_SetString(PyExc_ValueError, "math domain error"); + double y, r; + int n; + /* this function should only ever be called for finite arguments */ + assert(Py_IS_FINITE(x)); + y = fmod(fabs(x), 2.0); + n = (int)round(2.0*y); + assert(0 <= n && n <= 4); + switch (n) { + case 0: + r = sin(pi*y); + break; + case 1: + r = cos(pi*(y-0.5)); + break; + case 2: + /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give + -0.0 instead of 0.0 when y == 1.0. */ + r = sin(pi*(1.0-y)); + break; + case 3: + r = -cos(pi*(y-1.5)); + break; + case 4: + r = sin(pi*(y-2.0)); + break; + default: + assert(0); /* should never get here */ + r = -1.23e200; /* silence gcc warning */ + } + return copysign(1.0, x)*r; +} + +/* Implementation of the real gamma function. In extensive but non-exhaustive + random tests, this function proved accurate to within <= 10 ulps across the + entire float domain. Note that accuracy may depend on the quality of the + system math functions, the pow function in particular. Special cases + follow C99 annex F. The parameters and method are tailored to platforms + whose double format is the IEEE 754 binary64 format. + + Method: for x > 0.0 we use the Lanczos approximation with parameters N=13 + and g=6.024680040776729583740234375; these parameters are amongst those + used by the Boost library. Following Boost (again), we re-express the + Lanczos sum as a rational function, and compute it that way. The + coefficients below were computed independently using MPFR, and have been + double-checked against the coefficients in the Boost source code. + + For x < 0.0 we use the reflection formula. + + There's one minor tweak that deserves explanation: Lanczos' formula for + Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x + values, x+g-0.5 can be represented exactly. However, in cases where it + can't be represented exactly the small error in x+g-0.5 can be magnified + significantly by the pow and exp calls, especially for large x. A cheap + correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error + involved in the computation of x+g-0.5 (that is, e = computed value of + x+g-0.5 - exact value of x+g-0.5). Here's the proof: + + Correction factor + ----------------- + Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754 + double, and e is tiny. Then: + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e) + = pow(y, x-0.5)/exp(y) * C, + + where the correction_factor C is given by + + C = pow(1-e/y, x-0.5) * exp(e) + + Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so: + + C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y + + But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y), + + Note that for accuracy, when computing r*C it's better to do - else if (errno == ERANGE) { - /* ANSI C generally requires libm functions to set ERANGE - * on overflow, but also generally *allows* them to set - * ERANGE on underflow too. There's no consistency about - * the latter across platforms. - * Alas, C99 never requires that errno be set. - * Here we suppress the underflow errors (libm functions - * should return a zero on underflow, and +- HUGE_VAL on - * overflow, so testing the result for zero suffices to - * distinguish the cases). - * - * On some platforms (Ubuntu/ia64) it seems that errno can be - * set to ERANGE for subnormal results that do *not* underflow - * to zero. So to be safe, we'll ignore ERANGE whenever the - * function result is less than one in absolute value. - */ - if (fabs(x) < 1.0) - result = 0; - else - PyErr_SetString(PyExc_OverflowError, - "math range error"); + r + e*g/y*r; + + than + + r * (1 + e*g/y); + + since the addition in the latter throws away most of the bits of + information in e*g/y. +*/ + +#define LANCZOS_N 13 +static const double lanczos_g = 6.024680040776729583740234375; +static const double lanczos_g_minus_half = 5.524680040776729583740234375; +static const double lanczos_num_coeffs[LANCZOS_N] = { + 23531376880.410759688572007674451636754734846804940, + 42919803642.649098768957899047001988850926355848959, + 35711959237.355668049440185451547166705960488635843, + 17921034426.037209699919755754458931112671403265390, + 6039542586.3520280050642916443072979210699388420708, + 1439720407.3117216736632230727949123939715485786772, + 248874557.86205415651146038641322942321632125127801, + 31426415.585400194380614231628318205362874684987640, + 2876370.6289353724412254090516208496135991145378768, + 186056.26539522349504029498971604569928220784236328, + 8071.6720023658162106380029022722506138218516325024, + 210.82427775157934587250973392071336271166969580291, + 2.5066282746310002701649081771338373386264310793408 +}; + +/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ +static const double lanczos_den_coeffs[LANCZOS_N] = { + 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, + 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; + +/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ +#define NGAMMA_INTEGRAL 23 +static const double gamma_integral[NGAMMA_INTEGRAL] = { + 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, + 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, + 1307674368000.0, 20922789888000.0, 355687428096000.0, + 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, + 51090942171709440000.0, 1124000727777607680000.0, +}; + +/* Lanczos' sum L_g(x), for positive x */ + +static double +lanczos_sum(double x) +{ + double num = 0.0, den = 0.0; + int i; + assert(x > 0.0); + /* evaluate the rational function lanczos_sum(x). For large + x, the obvious algorithm risks overflow, so we instead + rescale the denominator and numerator of the rational + function by x**(1-LANCZOS_N) and treat this as a + rational function in 1/x. This also reduces the error for + larger x values. The choice of cutoff point (5.0 below) is + somewhat arbitrary; in tests, smaller cutoff values than + this resulted in lower accuracy. */ + if (x < 5.0) { + for (i = LANCZOS_N; --i >= 0; ) { + num = num * x + lanczos_num_coeffs[i]; + den = den * x + lanczos_den_coeffs[i]; + } } - else - /* Unexpected math error */ - PyErr_SetFromErrno(PyExc_ValueError); - return result; + else { + for (i = 0; i < LANCZOS_N; i++) { + num = num / x + lanczos_num_coeffs[i]; + den = den / x + lanczos_den_coeffs[i]; + } + } + return num/den; +} + +static double +m_tgamma(double x) +{ + double absx, r, y, z, sqrtpow; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x) || x > 0.0) + return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ + else { + errno = EDOM; + return Py_NAN; /* tgamma(-inf) = nan, invalid */ + } + } + if (x == 0.0) { + errno = EDOM; + return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */ + } + + /* integer arguments */ + if (x == floor(x)) { + if (x < 0.0) { + errno = EDOM; /* tgamma(n) = nan, invalid for */ + return Py_NAN; /* negative integers n */ + } + if (x <= NGAMMA_INTEGRAL) + return gamma_integral[(int)x - 1]; + } + absx = fabs(x); + + /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ + if (absx < 1e-20) { + r = 1.0/x; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; + } + + /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for + x > 200, and underflows to +-0.0 for x < -200, not a negative + integer. */ + if (absx > 200.0) { + if (x < 0.0) { + return 0.0/sinpi(x); + } + else { + errno = ERANGE; + return Py_HUGE_VAL; + } + } + + y = absx + lanczos_g_minus_half; + /* compute error in sum */ + if (absx > lanczos_g_minus_half) { + /* note: the correction can be foiled by an optimizing + compiler that (incorrectly) thinks that an expression like + a + b - a - b can be optimized to 0.0. This shouldn't + happen in a standards-conforming compiler. */ + double q = y - absx; + z = q - lanczos_g_minus_half; + } + else { + double q = y - lanczos_g_minus_half; + z = q - absx; + } + z = z * lanczos_g / y; + if (x < 0.0) { + r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx); + r -= z * r; + if (absx < 140.0) { + r /= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r /= sqrtpow; + r /= sqrtpow; + } + } + else { + r = lanczos_sum(absx) / exp(y); + r += z * r; + if (absx < 140.0) { + r *= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r *= sqrtpow; + r *= sqrtpow; + } + } + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; } /* @@ -188,6 +409,46 @@ } +/* Call is_error when errno != 0, and where x is the result libm + * returned. is_error will usually set up an exception and return + * true (1), but may return false (0) without setting up an exception. + */ +static int +is_error(double x) +{ + int result = 1; /* presumption of guilt */ + assert(errno); /* non-zero errno is a precondition for calling */ + if (errno == EDOM) + PyErr_SetString(PyExc_ValueError, "math domain error"); + + else if (errno == ERANGE) { + /* ANSI C generally requires libm functions to set ERANGE + * on overflow, but also generally *allows* them to set + * ERANGE on underflow too. There's no consistency about + * the latter across platforms. + * Alas, C99 never requires that errno be set. + * Here we suppress the underflow errors (libm functions + * should return a zero on underflow, and +- HUGE_VAL on + * overflow, so testing the result for zero suffices to + * distinguish the cases). + * + * On some platforms (Ubuntu/ia64) it seems that errno can be + * set to ERANGE for subnormal results that do *not* underflow + * to zero. So to be safe, we'll ignore ERANGE whenever the + * function result is less than one in absolute value. + */ + if (fabs(x) < 1.0) + result = 0; + else + PyErr_SetString(PyExc_OverflowError, + "math range error"); + } + else + /* Unexpected math error */ + PyErr_SetFromErrno(PyExc_ValueError); + return result; +} + /* math_1 is used to wrap a libm function f that takes a double arguments and returns a double. @@ -252,6 +513,26 @@ return (*from_double_func)(r); } +/* variant of math_1, to be used when the function being wrapped is known to + set errno properly (that is, errno = EDOM for invalid or divide-by-zero, + errno = ERANGE for overflow). */ + +static PyObject * +math_1a(PyObject *arg, double (*func) (double)) +{ + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + PyFPE_START_PROTECT("in math_1a", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + /* math_2 is used to wrap a libm function f that takes two double arguments and returns a double. @@ -330,6 +611,12 @@ }\ PyDoc_STRVAR(math_##funcname##_doc, docstring); +#define FUNC1A(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1a(args, func); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + #define FUNC2(funcname, func, docstring) \ static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ return math_2(args, func, #funcname); \ @@ -405,6 +692,8 @@ "floor(x)\n\nReturn the floor of x as an int.\n" "This is the largest integral value <= x."); +FUNC1A(gamma, m_tgamma, + "gamma(x)\n\nGamma function at x.") FUNC1(log1p, log1p, 1, "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\ The result is computed in a way which is accurate for x near zero.") @@ -1150,6 +1439,7 @@ {"fmod", math_fmod, METH_VARARGS, math_fmod_doc}, {"frexp", math_frexp, METH_O, math_frexp_doc}, {"fsum", math_fsum, METH_O, math_fsum_doc}, + {"gamma", math_gamma, METH_O, math_gamma_doc}, {"hypot", math_hypot, METH_VARARGS, math_hypot_doc}, {"isinf", math_isinf, METH_O, math_isinf_doc}, {"isnan", math_isnan, METH_O, math_isnan_doc},

1 0

r75118 - python/branches/py3k
by r.david.murray 29 Sep '09

29 Sep '09

Author: r.david.murray Date: Mon Sep 28 20:55:36 2009 New Revision: 75118 Log: Blocked revisions 75115 via svnmerge ........ r75115 | r.david.murray | 2009-09-28 14:29:28 -0400 (Mon, 28 Sep 2009) | 2 lines Applying patches backported from 3.1, by Gregor Lingl. ........ Modified: python/branches/py3k/ (props changed)

1 0

r75117 - in python/trunk: Doc/library/math.rst Lib/test/math_testcases.txt Lib/test/test_math.py Misc/NEWS Modules/mathmodule.c
by mark.dickinson 29 Sep '09

29 Sep '09

Author: mark.dickinson Date: Mon Sep 28 20:54:55 2009 New Revision: 75117 Log: Issue #3366: Add gamma function to math module. (lgamma, erf and erfc to follow). Added: python/trunk/Lib/test/math_testcases.txt Modified: python/trunk/Doc/library/math.rst python/trunk/Lib/test/test_math.py python/trunk/Misc/NEWS python/trunk/Modules/mathmodule.c Modified: python/trunk/Doc/library/math.rst ============================================================================== --- python/trunk/Doc/library/math.rst (original) +++ python/trunk/Doc/library/math.rst Mon Sep 28 20:54:55 2009 @@ -308,6 +308,16 @@ Return the hyperbolic tangent of *x*. +Special functions +----------------- + +.. function:: gamma(x) + + Return the Gamma function at *x*. + + .. versionadded:: 2.7 + + Constants --------- Added: python/trunk/Lib/test/math_testcases.txt ============================================================================== --- (empty file) +++ python/trunk/Lib/test/math_testcases.txt Mon Sep 28 20:54:55 2009 @@ -0,0 +1,146 @@ +-- Testcases for functions in math. +-- +-- Each line takes the form: +-- +-- <testid> <function> <input_value> -> <output_value> <flags> +-- +-- where: +-- +-- <testid> is a short name identifying the test, +-- +-- <function> is the function to be tested (exp, cos, asinh, ...), +-- +-- <input_value> is a string representing a floating-point value +-- +-- <output_value> is the expected (ideal) output value, again +-- represented as a string. +-- +-- <flags> is a list of the floating-point flags required by C99 +-- +-- The possible flags are: +-- +-- divide-by-zero : raised when a finite input gives a +-- mathematically infinite result. +-- +-- overflow : raised when a finite input gives a finite result that +-- is too large to fit in the usual range of an IEEE 754 double. +-- +-- invalid : raised for invalid inputs (e.g., sqrt(-1)) +-- +-- ignore-sign : indicates that the sign of the result is +-- unspecified; e.g., if the result is given as inf, +-- then both -inf and inf should be accepted as correct. +-- +-- Flags may appear in any order. +-- +-- Lines beginning with '--' (like this one) start a comment, and are +-- ignored. Blank lines, or lines containing only whitespace, are also +-- ignored. + +-- Many of the values below were computed with the help of +-- version 2.4 of the MPFR library for multiple-precision +-- floating-point computations with correct rounding. All output +-- values in this file are (modulo yet-to-be-discovered bugs) +-- correctly rounded, provided that each input and output decimal +-- floating-point value below is interpreted as a representation of +-- the corresponding nearest IEEE 754 double-precision value. See the +-- MPFR homepage at http://www.mpfr.org for more information about the +-- MPFR project. + +--------------------------- +-- gamma: Gamma function -- +--------------------------- + +-- special values +gam0000 gamma 0.0 -> inf divide-by-zero +gam0001 gamma -0.0 -> -inf divide-by-zero +gam0002 gamma inf -> inf +gam0003 gamma -inf -> nan invalid +gam0004 gamma nan -> nan + +-- negative integers inputs are invalid +gam0010 gamma -1 -> nan invalid +gam0011 gamma -2 -> nan invalid +gam0012 gamma -1e16 -> nan invalid +gam0013 gamma -1e300 -> nan invalid + +-- small positive integers give factorials +gam0020 gamma 1 -> 1 +gam0021 gamma 2 -> 1 +gam0022 gamma 3 -> 2 +gam0023 gamma 4 -> 6 +gam0024 gamma 5 -> 24 +gam0025 gamma 6 -> 120 + +-- half integers +gam0030 gamma 0.5 -> 1.7724538509055161 +gam0031 gamma 1.5 -> 0.88622692545275805 +gam0032 gamma 2.5 -> 1.3293403881791370 +gam0033 gamma 3.5 -> 3.3233509704478426 +gam0034 gamma -0.5 -> -3.5449077018110322 +gam0035 gamma -1.5 -> 2.3632718012073548 +gam0036 gamma -2.5 -> -0.94530872048294190 +gam0037 gamma -3.5 -> 0.27008820585226911 + +-- values near 0 +gam0040 gamma 0.1 -> 9.5135076986687306 +gam0041 gamma 0.01 -> 99.432585119150602 +gam0042 gamma 1e-8 -> 99999999.422784343 +gam0043 gamma 1e-16 -> 10000000000000000 +gam0044 gamma 1e-30 -> 9.9999999999999988e+29 +gam0045 gamma 1e-160 -> 1.0000000000000000e+160 +gam0046 gamma 1e-308 -> 1.0000000000000000e+308 +gam0047 gamma 5.6e-309 -> 1.7857142857142848e+308 +gam0048 gamma 5.5e-309 -> inf overflow +gam0049 gamma 1e-309 -> inf overflow +gam0050 gamma 1e-323 -> inf overflow +gam0051 gamma 5e-324 -> inf overflow +gam0060 gamma -0.1 -> -10.686287021193193 +gam0061 gamma -0.01 -> -100.58719796441078 +gam0062 gamma -1e-8 -> -100000000.57721567 +gam0063 gamma -1e-16 -> -10000000000000000 +gam0064 gamma -1e-30 -> -9.9999999999999988e+29 +gam0065 gamma -1e-160 -> -1.0000000000000000e+160 +gam0066 gamma -1e-308 -> -1.0000000000000000e+308 +gam0067 gamma -5.6e-309 -> -1.7857142857142848e+308 +gam0068 gamma -5.5e-309 -> -inf overflow +gam0069 gamma -1e-309 -> -inf overflow +gam0070 gamma -1e-323 -> -inf overflow +gam0071 gamma -5e-324 -> -inf overflow + +-- values near negative integers +gam0080 gamma -0.99999999999999989 -> -9007199254740992.0 +gam0081 gamma -1.0000000000000002 -> 4503599627370495.5 +gam0082 gamma -1.9999999999999998 -> 2251799813685248.5 +gam0083 gamma -2.0000000000000004 -> -1125899906842623.5 +gam0084 gamma -100.00000000000001 -> -7.5400833348831090e-145 +gam0085 gamma -99.999999999999986 -> 7.5400833348840962e-145 + +-- large inputs +gam0100 gamma 170 -> 4.2690680090047051e+304 +gam0101 gamma 171 -> 7.2574156153079990e+306 +gam0102 gamma 171.624 -> 1.7942117599248104e+308 +gam0103 gamma 171.625 -> inf overflow +gam0104 gamma 172 -> inf overflow +gam0105 gamma 2000 -> inf overflow +gam0106 gamma 1.7e308 -> inf overflow + +-- inputs for which gamma(x) is tiny +gam0120 gamma -100.5 -> -3.3536908198076787e-159 +gam0121 gamma -160.5 -> -5.2555464470078293e-286 +gam0122 gamma -170.5 -> -3.3127395215386074e-308 +gam0123 gamma -171.5 -> 1.9316265431711902e-310 +gam0124 gamma -176.5 -> -1.1956388629358166e-321 +gam0125 gamma -177.5 -> 4.9406564584124654e-324 +gam0126 gamma -178.5 -> -0.0 +gam0127 gamma -179.5 -> 0.0 +gam0128 gamma -201.0001 -> 0.0 +gam0129 gamma -202.9999 -> -0.0 +gam0130 gamma -1000.5 -> -0.0 +gam0131 gamma -1000000000.3 -> -0.0 +gam0132 gamma -4503599627370495.5 -> 0.0 + +-- inputs that cause problems for the standard reflection formula, +-- thanks to loss of accuracy in 1-x +gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88 +gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214 Modified: python/trunk/Lib/test/test_math.py ============================================================================== --- python/trunk/Lib/test/test_math.py (original) +++ python/trunk/Lib/test/test_math.py Mon Sep 28 20:54:55 2009 @@ -7,6 +7,7 @@ import os import sys import random +import struct eps = 1E-05 NAN = float('nan') @@ -29,8 +30,50 @@ else: file = __file__ test_dir = os.path.dirname(file) or os.curdir +math_testcases = os.path.join(test_dir, 'math_testcases.txt') test_file = os.path.join(test_dir, 'cmath_testcases.txt') +def to_ulps(x): + """Convert a non-NaN float x to an integer, in such a way that + adjacent floats are converted to adjacent integers. Then + abs(ulps(x) - ulps(y)) gives the difference in ulps between two + floats. + + The results from this function will only make sense on platforms + where C doubles are represented in IEEE 754 binary64 format. + + """ + n = struct.unpack('q', struct.pack('<d', x))[0] + if n < 0: + n = ~(n+2**63) + return n + + +def parse_mtestfile(fname): + """Parse a file with test values + + -- starts a comment + blank lines, or lines containing only a comment, are ignored + other lines are expected to have the form + id fn arg -> expected [flag]* + + """ + with open(fname) as fp: + for line in fp: + # strip comments, and skip blank lines + if '--' in line: + line = line[:line.index('--')] + if not line.strip(): + continue + + lhs, rhs = line.split('->') + id, fn, arg = lhs.split() + rhs_pieces = rhs.split() + exp = rhs_pieces[0] + flags = rhs_pieces[1:] + + yield (id, fn, float(arg), float(exp), flags) + def parse_testfile(fname): """Parse a file with test values @@ -887,6 +930,51 @@ self.fail(message) self.ftest("%s:%s(%r)" % (id, fn, ar), result, er) + @unittest.skipUnless(float.__getformat__("double").startswith("IEEE"), + "test requires IEEE 754 doubles") + def test_mtestfile(self): + ALLOWED_ERROR = 20 # permitted error, in ulps + fail_fmt = "{}:{}({!r}): expected {!r}, got {!r}" + + failures = [] + for id, fn, arg, expected, flags in parse_mtestfile(math_testcases): + func = getattr(math, fn) + + if 'invalid' in flags or 'divide-by-zero' in flags: + expected = 'ValueError' + elif 'overflow' in flags: + expected = 'OverflowError' + + try: + got = func(arg) + except ValueError: + got = 'ValueError' + except OverflowError: + got = 'OverflowError' + + diff_ulps = None + if isinstance(got, float) and isinstance(expected, float): + if math.isnan(expected) and math.isnan(got): + continue + if not math.isnan(expected) and not math.isnan(got): + diff_ulps = to_ulps(expected) - to_ulps(got) + if diff_ulps <= ALLOWED_ERROR: + continue + + if isinstance(got, str) and isinstance(expected, str): + if got == expected: + continue + + fail_msg = fail_fmt.format(id, fn, arg, expected, got) + if diff_ulps is not None: + fail_msg += ' ({} ulps)'.format(diff_ulps) + failures.append(fail_msg) + + if failures: + self.fail('Failures in test_mtestfile:\n ' + + '\n '.join(failures)) + + def test_main(): from doctest import DocFileSuite suite = unittest.TestSuite() Modified: python/trunk/Misc/NEWS ============================================================================== --- python/trunk/Misc/NEWS (original) +++ python/trunk/Misc/NEWS Mon Sep 28 20:54:55 2009 @@ -1341,6 +1341,8 @@ Extension Modules ----------------- +- Issue #3366: Add gamma function to math module. + - Issue #6823: Allow time.strftime() to accept a tuple with a isdst field outside of the range of [-1, 1] by normalizing the value to within that range. Modified: python/trunk/Modules/mathmodule.c ============================================================================== --- python/trunk/Modules/mathmodule.c (original) +++ python/trunk/Modules/mathmodule.c Mon Sep 28 20:54:55 2009 @@ -60,44 +60,265 @@ extern double copysign(double, double); #endif -/* Call is_error when errno != 0, and where x is the result libm - * returned. is_error will usually set up an exception and return - * true (1), but may return false (0) without setting up an exception. - */ -static int -is_error(double x) +/* + sin(pi*x), giving accurate results for all finite x (especially x + integral or close to an integer). This is here for use in the + reflection formula for the gamma function. It conforms to IEEE + 754-2008 for finite arguments, but not for infinities or nans. +*/ + +static const double pi = 3.141592653589793238462643383279502884197; + +static double +sinpi(double x) { - int result = 1; /* presumption of guilt */ - assert(errno); /* non-zero errno is a precondition for calling */ - if (errno == EDOM) - PyErr_SetString(PyExc_ValueError, "math domain error"); + double y, r; + int n; + /* this function should only ever be called for finite arguments */ + assert(Py_IS_FINITE(x)); + y = fmod(fabs(x), 2.0); + n = (int)round(2.0*y); + assert(0 <= n && n <= 4); + switch (n) { + case 0: + r = sin(pi*y); + break; + case 1: + r = cos(pi*(y-0.5)); + break; + case 2: + /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give + -0.0 instead of 0.0 when y == 1.0. */ + r = sin(pi*(1.0-y)); + break; + case 3: + r = -cos(pi*(y-1.5)); + break; + case 4: + r = sin(pi*(y-2.0)); + break; + default: + assert(0); /* should never get here */ + r = -1.23e200; /* silence gcc warning */ + } + return copysign(1.0, x)*r; +} + +/* Implementation of the real gamma function. In extensive but non-exhaustive + random tests, this function proved accurate to within <= 10 ulps across the + entire float domain. Note that accuracy may depend on the quality of the + system math functions, the pow function in particular. Special cases + follow C99 annex F. The parameters and method are tailored to platforms + whose double format is the IEEE 754 binary64 format. + + Method: for x > 0.0 we use the Lanczos approximation with parameters N=13 + and g=6.024680040776729583740234375; these parameters are amongst those + used by the Boost library. Following Boost (again), we re-express the + Lanczos sum as a rational function, and compute it that way. The + coefficients below were computed independently using MPFR, and have been + double-checked against the coefficients in the Boost source code. + + For x < 0.0 we use the reflection formula. + + There's one minor tweak that deserves explanation: Lanczos' formula for + Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x + values, x+g-0.5 can be represented exactly. However, in cases where it + can't be represented exactly the small error in x+g-0.5 can be magnified + significantly by the pow and exp calls, especially for large x. A cheap + correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error + involved in the computation of x+g-0.5 (that is, e = computed value of + x+g-0.5 - exact value of x+g-0.5). Here's the proof: + + Correction factor + ----------------- + Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754 + double, and e is tiny. Then: + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e) + = pow(y, x-0.5)/exp(y) * C, + + where the correction_factor C is given by + + C = pow(1-e/y, x-0.5) * exp(e) + + Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so: + + C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y + + But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and + + pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y), + + Note that for accuracy, when computing r*C it's better to do - else if (errno == ERANGE) { - /* ANSI C generally requires libm functions to set ERANGE - * on overflow, but also generally *allows* them to set - * ERANGE on underflow too. There's no consistency about - * the latter across platforms. - * Alas, C99 never requires that errno be set. - * Here we suppress the underflow errors (libm functions - * should return a zero on underflow, and +- HUGE_VAL on - * overflow, so testing the result for zero suffices to - * distinguish the cases). - * - * On some platforms (Ubuntu/ia64) it seems that errno can be - * set to ERANGE for subnormal results that do *not* underflow - * to zero. So to be safe, we'll ignore ERANGE whenever the - * function result is less than one in absolute value. - */ - if (fabs(x) < 1.0) - result = 0; - else - PyErr_SetString(PyExc_OverflowError, - "math range error"); + r + e*g/y*r; + + than + + r * (1 + e*g/y); + + since the addition in the latter throws away most of the bits of + information in e*g/y. +*/ + +#define LANCZOS_N 13 +static const double lanczos_g = 6.024680040776729583740234375; +static const double lanczos_g_minus_half = 5.524680040776729583740234375; +static const double lanczos_num_coeffs[LANCZOS_N] = { + 23531376880.410759688572007674451636754734846804940, + 42919803642.649098768957899047001988850926355848959, + 35711959237.355668049440185451547166705960488635843, + 17921034426.037209699919755754458931112671403265390, + 6039542586.3520280050642916443072979210699388420708, + 1439720407.3117216736632230727949123939715485786772, + 248874557.86205415651146038641322942321632125127801, + 31426415.585400194380614231628318205362874684987640, + 2876370.6289353724412254090516208496135991145378768, + 186056.26539522349504029498971604569928220784236328, + 8071.6720023658162106380029022722506138218516325024, + 210.82427775157934587250973392071336271166969580291, + 2.5066282746310002701649081771338373386264310793408 +}; + +/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ +static const double lanczos_den_coeffs[LANCZOS_N] = { + 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, + 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; + +/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ +#define NGAMMA_INTEGRAL 23 +static const double gamma_integral[NGAMMA_INTEGRAL] = { + 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, + 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, + 1307674368000.0, 20922789888000.0, 355687428096000.0, + 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, + 51090942171709440000.0, 1124000727777607680000.0, +}; + +/* Lanczos' sum L_g(x), for positive x */ + +static double +lanczos_sum(double x) +{ + double num = 0.0, den = 0.0; + int i; + assert(x > 0.0); + /* evaluate the rational function lanczos_sum(x). For large + x, the obvious algorithm risks overflow, so we instead + rescale the denominator and numerator of the rational + function by x**(1-LANCZOS_N) and treat this as a + rational function in 1/x. This also reduces the error for + larger x values. The choice of cutoff point (5.0 below) is + somewhat arbitrary; in tests, smaller cutoff values than + this resulted in lower accuracy. */ + if (x < 5.0) { + for (i = LANCZOS_N; --i >= 0; ) { + num = num * x + lanczos_num_coeffs[i]; + den = den * x + lanczos_den_coeffs[i]; + } } - else - /* Unexpected math error */ - PyErr_SetFromErrno(PyExc_ValueError); - return result; + else { + for (i = 0; i < LANCZOS_N; i++) { + num = num / x + lanczos_num_coeffs[i]; + den = den / x + lanczos_den_coeffs[i]; + } + } + return num/den; +} + +static double +m_tgamma(double x) +{ + double absx, r, y, z, sqrtpow; + + /* special cases */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_NAN(x) || x > 0.0) + return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ + else { + errno = EDOM; + return Py_NAN; /* tgamma(-inf) = nan, invalid */ + } + } + if (x == 0.0) { + errno = EDOM; + return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */ + } + + /* integer arguments */ + if (x == floor(x)) { + if (x < 0.0) { + errno = EDOM; /* tgamma(n) = nan, invalid for */ + return Py_NAN; /* negative integers n */ + } + if (x <= NGAMMA_INTEGRAL) + return gamma_integral[(int)x - 1]; + } + absx = fabs(x); + + /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ + if (absx < 1e-20) { + r = 1.0/x; + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; + } + + /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for + x > 200, and underflows to +-0.0 for x < -200, not a negative + integer. */ + if (absx > 200.0) { + if (x < 0.0) { + return 0.0/sinpi(x); + } + else { + errno = ERANGE; + return Py_HUGE_VAL; + } + } + + y = absx + lanczos_g_minus_half; + /* compute error in sum */ + if (absx > lanczos_g_minus_half) { + /* note: the correction can be foiled by an optimizing + compiler that (incorrectly) thinks that an expression like + a + b - a - b can be optimized to 0.0. This shouldn't + happen in a standards-conforming compiler. */ + double q = y - absx; + z = q - lanczos_g_minus_half; + } + else { + double q = y - lanczos_g_minus_half; + z = q - absx; + } + z = z * lanczos_g / y; + if (x < 0.0) { + r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx); + r -= z * r; + if (absx < 140.0) { + r /= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r /= sqrtpow; + r /= sqrtpow; + } + } + else { + r = lanczos_sum(absx) / exp(y); + r += z * r; + if (absx < 140.0) { + r *= pow(y, absx - 0.5); + } + else { + sqrtpow = pow(y, absx / 2.0 - 0.25); + r *= sqrtpow; + r *= sqrtpow; + } + } + if (Py_IS_INFINITY(r)) + errno = ERANGE; + return r; } /* @@ -188,6 +409,46 @@ } +/* Call is_error when errno != 0, and where x is the result libm + * returned. is_error will usually set up an exception and return + * true (1), but may return false (0) without setting up an exception. + */ +static int +is_error(double x) +{ + int result = 1; /* presumption of guilt */ + assert(errno); /* non-zero errno is a precondition for calling */ + if (errno == EDOM) + PyErr_SetString(PyExc_ValueError, "math domain error"); + + else if (errno == ERANGE) { + /* ANSI C generally requires libm functions to set ERANGE + * on overflow, but also generally *allows* them to set + * ERANGE on underflow too. There's no consistency about + * the latter across platforms. + * Alas, C99 never requires that errno be set. + * Here we suppress the underflow errors (libm functions + * should return a zero on underflow, and +- HUGE_VAL on + * overflow, so testing the result for zero suffices to + * distinguish the cases). + * + * On some platforms (Ubuntu/ia64) it seems that errno can be + * set to ERANGE for subnormal results that do *not* underflow + * to zero. So to be safe, we'll ignore ERANGE whenever the + * function result is less than one in absolute value. + */ + if (fabs(x) < 1.0) + result = 0; + else + PyErr_SetString(PyExc_OverflowError, + "math range error"); + } + else + /* Unexpected math error */ + PyErr_SetFromErrno(PyExc_ValueError); + return result; +} + /* math_1 is used to wrap a libm function f that takes a double arguments and returns a double. @@ -247,6 +508,26 @@ return PyFloat_FromDouble(r); } +/* variant of math_1, to be used when the function being wrapped is known to + set errno properly (that is, errno = EDOM for invalid or divide-by-zero, + errno = ERANGE for overflow). */ + +static PyObject * +math_1a(PyObject *arg, double (*func) (double)) +{ + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + PyFPE_START_PROTECT("in math_1a", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + /* math_2 is used to wrap a libm function f that takes two double arguments and returns a double. @@ -313,6 +594,12 @@ }\ PyDoc_STRVAR(math_##funcname##_doc, docstring); +#define FUNC1A(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1a(args, func); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + #define FUNC2(funcname, func, docstring) \ static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ return math_2(args, func, #funcname); \ @@ -350,6 +637,8 @@ FUNC1(floor, floor, 0, "floor(x)\n\nReturn the floor of x as a float.\n" "This is the largest integral value <= x.") +FUNC1A(gamma, m_tgamma, + "gamma(x)\n\nGamma function at x.") FUNC1(log1p, log1p, 1, "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\ The result is computed in a way which is accurate for x near zero.") @@ -1077,6 +1366,7 @@ {"fmod", math_fmod, METH_VARARGS, math_fmod_doc}, {"frexp", math_frexp, METH_O, math_frexp_doc}, {"fsum", math_fsum, METH_O, math_fsum_doc}, + {"gamma", math_gamma, METH_O, math_gamma_doc}, {"hypot", math_hypot, METH_VARARGS, math_hypot_doc}, {"isinf", math_isinf, METH_O, math_isinf_doc}, {"isnan", math_isnan, METH_O, math_isnan_doc},

1 0

r75116 - in python/branches/release26-maint: Demo/turtle/turtleDemo.py Lib/lib-tk/turtle.py
by r.david.murray 29 Sep '09

29 Sep '09

Author: r.david.murray Date: Mon Sep 28 20:45:13 2009 New Revision: 75116 Log: Merged revisions 75115 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r75115 | r.david.murray | 2009-09-28 14:29:28 -0400 (Mon, 28 Sep 2009) | 2 lines Applying patches backported from 3.1, by Gregor Lingl. ........ Modified: python/branches/release26-maint/ (props changed) python/branches/release26-maint/Demo/turtle/turtleDemo.py python/branches/release26-maint/Lib/lib-tk/turtle.py Modified: python/branches/release26-maint/Demo/turtle/turtleDemo.py ============================================================================== --- python/branches/release26-maint/Demo/turtle/turtleDemo.py (original) +++ python/branches/release26-maint/Demo/turtle/turtleDemo.py Mon Sep 28 20:45:13 2009 @@ -98,6 +98,7 @@ turtle._Screen._canvas = turtle.ScrolledCanvas(g_frame, 800, 600, 1000, 800) #xturtle.Screen._canvas.pack(expand=1, fill="both") self.screen = _s_ = turtle.Screen() + turtle.TurtleScreen.__init__(_s_, _s_._canvas) self.scanvas = _s_._canvas #xturtle.RawTurtle.canvases = [self.scanvas] turtle.RawTurtle.screens = [_s_] Modified: python/branches/release26-maint/Lib/lib-tk/turtle.py ============================================================================== --- python/branches/release26-maint/Lib/lib-tk/turtle.py (original) +++ python/branches/release26-maint/Lib/lib-tk/turtle.py Mon Sep 28 20:45:13 2009 @@ -1,8 +1,8 @@ # # turtle.py: a Tkinter based turtle graphics module for Python -# Version 1.0b1 - 31. 5. 2008 +# Version 1.0.1 - 24. 9. 2009 # -# Copyright (C) 2006 - 2008 Gregor Lingl +# Copyright (C) 2006 - 2009 Gregor Lingl # email: glingl(a)aon.at # # This software is provided 'as-is', without any express or implied @@ -1244,9 +1244,12 @@ def update(self): """Perform a TurtleScreen update. """ + tracing = self._tracing + self._tracing = True for t in self.turtles(): t._update_data() t._drawturtle() + self._tracing = tracing self._update() def window_width(self): @@ -1337,7 +1340,8 @@ ### consequently drawing a hexagon """ if fun == None: - self._keys.remove(key) + if key in self._keys: + self._keys.remove(key) elif key not in self._keys: self._keys.append(key) self._onkey(fun, key) @@ -3574,8 +3578,8 @@ topbottom = _CFG["topbottom"] self._root.setupcanvas(width, height, canvwidth, canvheight) _Screen._canvas = self._root._getcanvas() + TurtleScreen.__init__(self, _Screen._canvas) self.setup(width, height, leftright, topbottom) - TurtleScreen.__init__(self, _Screen._canvas) def setup(self, width=_CFG["width"], height=_CFG["height"], startx=_CFG["leftright"], starty=_CFG["topbottom"]): @@ -3615,6 +3619,7 @@ if starty is None: starty = (sh - height) / 2 self._root.set_geometry(width, height, startx, starty) + self.update() def title(self, titlestring): """Set title of turtle-window

1 0