184 views

Skip to first unread message

Sep 24, 2014, 6:35:50 AM9/24/14

to sage-...@googlegroups.com

there is quite a bit of weirdness mentioned here:

e.g.

sage: hypergeometric([-2,-1],[2],-1).n(100)

---------------------------------------------------------------------------

ValueError Traceback (most recent call last)

<ipython-input-4-0bef797c6f48> in <module>()

----> 1 hypergeometric([-Integer(2),-Integer(1)],[Integer(2)],-Integer(1)).n(Integer(100))

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (build/cythonized/sage/symbolic/expression.cpp:26972)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:7667)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/functions/hypergeometric.pyc in _evalf_(self, a, b, z, parent, algorithm)

294 aa = [rational_param_as_tuple(c) for c in a]

295 bb = [rational_param_as_tuple(c) for c in b]

--> 296 return mpmath_utils.call(hyper, aa, bb, z, parent=parent)

297

298 def _tderivative_(self, a, b, z, *args, **kwargs):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (build/cythonized/sage/libs/mpmath/utils.c:6309)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in hyper(ctx, a_s, b_s, z, **kwargs)

223 elif q == 0: return ctx._hyp1f0(a_s[0][0], z)

224 elif p == 2:

--> 225 if q == 1: return ctx._hyp2f1(a_s, b_s, z, **kwargs)

226 elif q == 2: return ctx._hyp2f2(a_s, b_s, z, **kwargs)

227 elif q == 3: return ctx._hyp2f3(a_s, b_s, z, **kwargs)

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in _hyp2f1(ctx, a_s, b_s, z, **kwargs)

442 if absz <= 0.8 or (ctx.isint(a) and a <= 0 and a >= -1000) or \

443 (ctx.isint(b) and b <= 0 and b >= -1000):

--> 444 return ctx.hypsum(2, 1, (atype, btype, ctype), [a, b, c], z, **kwargs)

445

446 orig = ctx.prec

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/ctx_mp.pyc in hypsum(ctx, p, q, flags, coeffs, z, accurate_small, **kwargs)

686 while 1:

687 if extraprec > maxprec:

--> 688 raise ValueError(ctx._hypsum_msg % (prec, prec+extraprec))

689 wp = prec + extraprec

690 if magnitude_check:

ValueError: hypsum() failed to converge to the requested 100 bits of accuracy

using a working precision of 7135 bits. Try with a higher maxprec,

maxterms, or set zeroprec.

sage:

(and if I try .n(10000) the error is the same, only the last bit reads as

ValueError: hypsum() failed to converge to the requested 10000 bits of accuracy

using a working precision of 66315 bits. Try with a higher maxprec,

maxterms, or set zeroprec.)

and

sage: maxima_calculus.hgfpoly([1,2],[2],-1)

---------------------------------------------------------------------------

TypeError Traceback (most recent call last)

<ipython-input-6-a307d86456bf> in <module>()

----> 1 maxima_calculus.hgfpoly([Integer(1),Integer(2)],[Integer(2)],-Integer(1))

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, *args, **kwds)

561

562 def __call__(self, *args, **kwds):

--> 563 return self._parent.function_call(self._name, list(args), kwds)

564

565 def _sage_doc_(self):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in function_call(self, function, args, kwds)

487 [s.name() for s in args],

488 ['%s=%s'%(key,value.name()) for key, value in kwds.items()])

--> 489 return self.new(s)

490

491 def _function_call_string(self, function, args, kwds):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in new(self, code)

262

263 def new(self, code):

--> 264 return self(code)

265

266 ###################################################################

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)

197

198 if isinstance(x, basestring):

--> 199 return cls(self, x, name=name)

200 try:

201 return self._coerce_from_special_method(x)

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __init__(self, parent, value, is_name, name)

624 self._name = parent._create(value, name=name)

625 except (TypeError, RuntimeError, ValueError) as x:

--> 626 raise TypeError(x)

627

628 def _latex_(self):

TypeError: ECL says: In function -, the value of the only argument is

NIL

which is not of the expected type NUMBER

sage:

-----------------------------------------------------------------------------

is any of these known? In the 1st case, looks like numeric evaluation of qFp by mpmath is

seriously broken...

Sep 24, 2014, 7:11:57 PM9/24/14

to sage-...@googlegroups.com

On Wednesday, September 24, 2014 12:35:50 PM UTC+2, Dima Pasechnik wrote:

there is quite a bit of weirdness mentioned here:e.g.sage: hypergeometric([-2,-1],[2],-1).n(100)---------------------------------------------------------------------------

-----------------------------------------------------------------------------is any of these known? In the 1st case, looks like numeric evaluation of qFp by mpmath isseriously broken...

hyp2f1(-2,-1,2,-1) doesn't converge because the value is zero, and the hypergeometric function evaluation code in mpmath always tries to achieve full *relative* accuracy. One solution is to do:

>>> hyp2f1(-2,-1,2,-1,zeroprec=1000)

mpf('0.0')

This is saying that if the result cannot be distinguished from zero when using 1000 bits of working precision, it is probably actually zero and not just merely something very small (so don't be fussy about it).

In fact, hyp2f1(-2,-1,2,z) is just the polynomial 1+z. mpmath has no way to distinguish 1 + (-1) from 1 + (-1 + eps) in which the eps disappears due rounding, so it has to assume that something may have disappeared. Indeed:

>>> mp.dps = 1000

>>> z = mpf(-1) + mpf("1e-900")

>>> mp.dps = 15

>>> hyper([-2,-1],[2],z) # conservatively throws an error

Traceback (most recent call last):

...

ValueError: hypsum() failed to converge to the requested 53 bits of accuracy

using a working precision of 3568 bits. Try with a higher maxprec,

maxterms, or set zeroprec.

>>> hyper([-2,-1],[2],z,zeroprec=1000) # WRONG!!! (relatively speaking)

mpf('0.0')

>>> hyper([-2,-1],[2],z,maxprec=10000) # correct

mpf('9.9999999999999998e-901')

Exact zero detection could be done at least in trivial cases by identifying easy factors such as z-1 or z+1. In general, I think one basically has to fall back to using exact rational arithmetic for the evaluation. This would be worth implementing in mpmath, because the current behavior is definitely annoying.

For Sage, fixing the problem is actually trivial: when the hypergeometric function is a polynomial (and at least when the inputs are exact), don't call mpmath; just evaluate the polynomial directly and then call .n() on the result.

Fredrik

Sep 25, 2014, 3:02:27 AM9/25/14

to sage-...@googlegroups.com

For Sage, fixing the problem is actually trivial: when the hypergeometric function is a polynomial (and at least when the inputs are exact), don't call mpmath; just evaluate the polynomial directly and then call .n() on the result.

Except then Sage would have to know when it is a polynomial, and probably we would need to ask Maxima for that (assuming it knows). So maybe not completely trivial to make sure it works.

- kcrisman

Sep 25, 2014, 3:05:44 AM9/25/14

to sage-...@googlegroups.com

Fredrik

Sep 25, 2014, 3:51:36 PM9/25/14

to sage-...@googlegroups.com

>> For Sage, fixing the problem is actually trivial: when the hypergeometric

>> function is a polynomial (and at least when the inputs are exact), don't

>> call mpmath; just evaluate the polynomial directly and then call .n() on the

>> result.

>>

>

> Except then Sage would have to know when it is a polynomial, and probably we

> would need to ask Maxima for that (assuming it knows). So maybe not

> completely trivial to make sure it works.

It's a polynomial when any of the first parameters is a nonpositive integer.

Is that "if and only if"? That would certainly be convenient.

Sep 26, 2014, 4:10:26 AM9/26/14

to SAGE devel

expansion of (1-x)^{-n}, there are rising factorials and these hit

zero if and only if they start at a nonpositive integer. But you

should check the definition of hypergeometric function to be 100%

sure.

John

>

> --

> You received this message because you are subscribed to the Google Groups

> "sage-devel" group.

> To unsubscribe from this group and stop receiving emails from it, send an

> email to sage-devel+...@googlegroups.com.

> To post to this group, send email to sage-...@googlegroups.com.

> Visit this group at http://groups.google.com/group/sage-devel.

> For more options, visit https://groups.google.com/d/optout.

Sep 26, 2014, 4:17:17 AM9/26/14

to sage-...@googlegroups.com

Dima

Sep 26, 2014, 5:34:52 AM9/26/14

to sage-...@googlegroups.com

I should clarify that the hypergeometric function can be zero for

rational input even if it is not a polynomial in z. For example

1F1(3/2, 1/2, z) = (1 + 2*z) * exp(z) (you can verify that mpmath does

not converge with z = -0.5). You are less likely to run into such

zeros, though.

Fredrik

Sep 26, 2014, 8:30:52 AM9/26/14

to sage-...@googlegroups.com

>> >> For Sage, fixing the problem is actually trivial: when the

>> >> hypergeometric

>> >> function is a polynomial (and at least when the inputs are exact),

>> >> don't

>> >> call mpmath; just evaluate the polynomial directly and then call .n()

>> >> on the

>> >> result.

>> >>

>> >

>> > Except then Sage would have to know when it is a polynomial, and

>> > probably we

>> > would need to ask Maxima for that (assuming it knows). So maybe not

>> > completely trivial to make sure it works.

>>

>> It's a polynomial when any of the first parameters is a nonpositive

>> integer.

>>

>

> Is that "if and only if"? That would certainly be convenient.

Yes, that follows immediately from the definition.

Well, then I guess this shouldn't be too hard to implement... I will admit that I will not be the one doing it, but sounds like a no-brainer.

Sep 26, 2014, 9:55:39 AM9/26/14

to sage-...@googlegroups.com

(AFAIK, hypergeometric_U is vulnerable to the very same problem, so this

would mean the workaround would either be needed in several places...)

Dima

Sep 26, 2014, 10:25:58 AM9/26/14

to sage-...@googlegroups.com

On Friday, September 26, 2014 6:55:39 AM UTC-7, Dima Pasechnik wrote:

IMHO it would be great to have this in mpmath rather than a Sage workaround.

(AFAIK, hypergeometric_U is vulnerable to the very same problem, so this

would mean the workaround would either be needed in several places...)

It would be good if mpmath used a different approach when it's known a hypergeometric function is a polynomial, but it's a little orthogonal to what happens in sage. By the time we hit mpmath, we're doing multiprecision float computation. Polynomial hypergeometric functions can easily be evaluated at arbitrary inputs (rational, algebraic, symbolic). We should not be bothering with multiprecision floats if we don't have to, or at least only do so if there's a benefit and we can control precision issues.

Sep 30, 2014, 12:51:12 PM9/30/14

to sage-...@googlegroups.com

On Friday, September 26, 2014 4:25:58 PM UTC+2, Nils Bruin wrote:

It would be good if mpmath used a different approach when it's known a hypergeometric function is a polynomial, but it's a little orthogonal to what happens in sage. By the time we hit mpmath, we're doing multiprecision float computation. Polynomial hypergeometric functions can easily be evaluated at arbitrary inputs (rational, algebraic, symbolic). We should not be bothering with multiprecision floats if we don't have to, or at least only do so if there's a benefit and we can control precision issues.

In most symbolic special functions we check for special values and, if found, evaluate immediately. This is not so in hypergeometric.py but would be trivial to add. I can only speculate that the reason was the then forced call of maxima. Since this doesn't seem a problem with commenters please review

and this can be simply emulated until inclusion of the ticket by saying

`sage: hypergeometric([-2,-1],[2],-1).simplify_hypergeometric()`

0

No need to explicitly call Maxima here.

Regards,

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu