Following functions kept:
```pycon
>>> import inspect
... with_args = []
... with_kwargs = []
... for n in dir(mpmath):
... m = getattr(mpmath, n)
... try:
... s = inspect.signature(m)
... except:
... continue
... if any(_.kind == inspect._ParameterKind.VAR_POSITIONAL for _ in s.parameters.values()):
... for name in s.parameters:
... if s.parameters[name].kind == inspect._ParameterKind.VAR_POSITIONAL and name == 'args':
... with_args.append(n)
... break
... if any(_.kind == inspect._ParameterKind.VAR_KEYWORD for _ in s.parameters.values()):
... with_kwargs.append(n)
... print(with_args)
... print(with_kwargs)
...
['arange', 'ellipe', 'ellippi', 'linspace', 'matrix', 'ones', 'timing', 'zeros']
['multiplicity', 'timing']
```
We need support for multiple signatures in the first case. In the
second - it's impossible to implement these functions without kwargs.
Closes#1056
* Fix qr_solve() failure on well-conditioned matrices with zero pivot
In householder(), the sign convention
p[j] = -sign(Re(A[j,j])) * sqrt(s)
collapses to zero when A[j,j] is exactly zero, because ctx.sign(0) == 0.
A zero p[j] makes kappa = 1/(s - p[j]*A[j,j]) = 1/s instead of the
correct 2/||v||**2, corrupting the Householder reflection. The damage
propagates to subsequent columns and eventually trips the
"matrix is numerically singular" guard.
Default sign to ctx.one when A[j,j] is zero (matching LAPACK's dlarfg
convention) so the reflection is computed correctly. The existing
singularity check is left in place to catch genuinely zero column
slices.
Fixes#983.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
We use pytest for tests, it has everything that runtests.py has,
including support for profiling (pytest-profiling plugin).
So, let's drop this legacy stuff to avoid bugs like #579.
Make matrices type-stable, not allowing mixed types in one matrix. Always use the type of the context.
Therefore, specifying an element type by the "force_type" argument does not make sense.
It did not work properly before either, the data was not stored in the given type:
```
>>> M = mpmath.matrix(mpmath.matrix([[0,1,2]]), force_type=bool)
>>> M
matrix(
[['0.0', '1.0', '1.0']])
>>> type(M[0,0])
<class 'mpmath.ctx_mp_python.mpf'>
>>> M._matrix__data
{(0, 1): mpf('1.0'), (0, 2): mpf('1.0')}
```
Now the argument is removed and gives a warning instead:
>>> mpmath.matrix(mpmath.matrix([[0,1,2]]), force_type=bool)
/mpmath/mpmath/matrices/matrices.py:288: UserWarning: The force_type argument was removed, it did not work properly anyway. If you want to force floating-point or interval computations, use the respective methods from `fp` or `mp` instead, e.g., `fp.matrix()` or `iv.matrix()`. If you want to truncate values to integer, use .apply(int) instead.
warnings.warn("The force_type argument was removed, it did not work"
matrix([['0.0', '1.0', '2.0']])
Non-interval datatypes were not converted to intervals, resulting in unexpected behavior. For example, M[i, j] didn't return mpi interval objects, but mpf floats, so M[i, j] * M[k, l] and, under certain conditions, M * M, was not performed as interval multiplication.
This could have unexpected consequences, as outlined below.
"Exact" computation for reference:
>>> import mpmath
>>> x = mpmath.convert('1.00000000000001') # 1.0000...001 rounded to the next mpf floating point value
>>> mpmath.mp.dps=1000
>>> x*x # Good approximation of x*x
mpf('1.00000000000001998401444325291756783368705994138804689654360263219301518944348572404123842716217041015625')
Interval arithmetic should return an interval containing that value, with some uncertainty.
Restart python to reset mp.dps.
Before this commit, the following occurs:
>>> import mpmath
>>> x = mpmath.convert('1.00000000000001')
>>> x
mpf('1.00000000000001')
>>> A = mpmath.matrix([[x]])
>>> B = mpmath.iv.matrix(A)
>>> C = mpmath.iv.matrix([[x]])
>>> A*A
matrix(
[['1.00000000000002']])
>>> B*B
matrix(
[['[1.000000000000019984, 1.000000000000019984]']])
>>> (B*B)[0,0].delta
mpi('0.0', '0.0')
>>> C*C
matrix(
[['[1.000000000000019984, 1.0000000000000202061]']])
B*B is wrong, the interval width must be nonzero at the default precision of 15 digits.
B*B is different from C*C, although both were initialized from the same numerical value and computed the same way.
After this commit, the result is valid:
>>> B*B
matrix(
[['[1.000000000000019984, 1.0000000000000202061]']])
>>> (B*B)[0,0].delta
mpi('2.2204460492503131e-16', '2.2204460492503131e-16')
Some more insight:
Old incorrect behavior:
>>> mpmath.iv.matrix(mpmath.eye(1))[0,0]
mpf('1.0')
New correct behavior:
>>> mpmath.iv.matrix(mpmath.eye(1))[0,0]
mpf('1.0')
>>> mpmath.iv.eye(1)[0,0]
mpi('1.0', '1.0')
>>> import mpmath
>>> mpmath.iv.matrix(mpmath.eye(1))
matrix(
[['[1.0, 1.0]']])
>>> mpmath.iv.matrix(mpmath.eye(1))[0,0]
mpi('1.0', '1.0')
>>> mpmath.fp.matrix(mpmath.eye(1))[0,0]
1.0
>>> mpmath.matrix(mpmath.eye(1))[0,0]
mpf('1.0')
The type now exactly matches the type returned by the context's matrix functions such as eye():
>>> mpmath.matrix(mpmath.eye(1))[0,0]
mpf('1.0')
>>> mpmath.eye(1)[0,0]
mpf('1.0')
>>> mpmath.iv.matrix(mpmath.eye(1))[0,0]
mpi('1.0', '1.0')
>>> mpmath.iv.eye(1)[0,0]
mpi('1.0', '1.0')
>>> mpmath.fp.matrix(mpmath.eye(1))[0,0]
1.0
>>> mpmath.fp.eye(1)[0,0]
1.0
Sergey B Kirpichev