Make RootOf.n() return the correct root

The way RootOf.n() works is the following: RootOf contains the interval (a, b)
where the (only) root is (as symbolically represented by RootOf). We then take
the center of the interval as an initial estimate for findroot() from mpmath.
However, sometimes findroot() finds a root that is not from the interval (a,
b), because the initial estimate was too far, and findroot() found a different
root (from some other interval) instead.

The solution is to check, that the findroot()'s result is from the interval
(a, b) and if it is not, just keep refining the interval. Eventually the
initial estimate will be close enough, so that findroot() converges to the
correct root.

Test (that previously failed and now works) added.
RootOf.eval_rational() implemented

Test added.
Add Gauss-Legendre and Gauss-Laguerre quadratures

Added gauss_legendre() and gauss_laguerre() functions and tests for calculating
the quadrature points and weights. The functions work for any order and any
precision.
Add Gaussian integrals into docs
Fix docstring formatting, added "Returns" section

Now it should follow SymPy's conventions. Also it should make clear that the
function returns numerical values (not symbolic).
Simplify condition to "a < root < b"
Use `` instead of '
Fix bisect's docstring
Fix RootOf.n() for complex roots

The problem was the same as for the real roots that was fixed in
5c58406b3be8fc6d75d0e258a67b6847d5274a4f. This commit provides the same fix for
complex roots by refining the complex interval.
Add a test for RootOf() for complex roots

These tests are correct, but currently fail because they have revealed a few
other bugs, which will be fixed in subsequent commits.  After the other bugs
are fixed, these tests will start passing on their own.
Do not store the refined interval in RootOf

This fixes a caching bug. A test has been added, that failed before.
What happened is that the evalf() method refines the interval and stores it
back into the RootOf() method, which uses global cache... Subsequent calls to
RootOf with the same polynomial and the *same* root index then assigns a wrong
interval for the root.

Ideally, the RootOf() implementation should be improved so that it is robust
against refining the intervals. For the time being, this patch fixes the issue.
Refine the complex interval properly

Now all tests pass. Previously a wrong root was returned due to numerical
inaccuracies.

In particular, the following case was happening: the imaginary part of the root
was supposed to be in [0, 3] and it was a complex root, so it couldn't be 0.
However, due to numerics, the returned imaginary part was 1e-69. Which belongs
into this interval technically, but in reality the returned root was simply a
real root, already returned with some other root index. (The real part happened
to fit into the same interval by accident.)

This patch makes sure that the imaginar part is always distinct from other
roots. That way things start working.
Make RootOf() work with gmpy

Now it works with gmpy as well as Python ground types.
Simplify a while loop
Merge pull request #1497 from certik/roots_bug

Calculate Gauss-Legendre and Gauss-Laguerre points and weights