acos获取误差数学域误差

Question

想弄清楚为什么我会犯错。我的数字介于-1和1之间，但仍然是错误的。

ValueError:数学域错误

有什么想法吗？

谢谢

from math import sqrt, acos, pi
from decimal import Decimal, getcontext

getcontext().prec = 30


class Vector(object):
    CANNOT_NORMALIZE_ZERO_VECTOR_MSG = 'Cannot normalize the zero vector'

    def __init__(self, coordinates):
        try:
            if not coordinates:
                raise ValueError
            self.coordinates = tuple([Decimal(x) for x in coordinates])
            self.dimension = len(self.coordinates)

        except ValueError:
            raise ValueError('The coordinates must be nonempty')

        except TypeError:
            raise TypeError('The coordinates must be an iterable')

    def __str__(self):
        return 'Vector: {}'.format(self.coordinates)

    def __eq__(self, v):
        return self.coordinates == v.coordinates

    def magnitude(self):
        coordinates_squared = [x ** 2 for x in self.coordinates]
        return sqrt(sum(coordinates_squared))

    def normalized(self):
        try:
            magnitude = self.magnitude()
            return self.times_scalar(Decimal(1.0 / magnitude))

        except ZeroDivisionError:
            raise Exception('Cannot normalize the zero vector')

    def plus(self, v):
        new_coordinates = [x + y for x, y in zip(self.coordinates, v.coordinates)]
        return Vector(new_coordinates)

    def minus(self, v):
        new_coordinates = [x - y for x, y in zip(self.coordinates, v.coordinates)]
        return Vector(new_coordinates)

    def times_scalar(self, c):
        new_coordinates = [Decimal(c) * x for x in self.coordinates]
        return Vector(new_coordinates)

    def dot(self, v):
        return sum([x * y for x, y in zip(self.coordinates, v.coordinates)])

    def angle_with(self, v, in_degrees=False):
        try:
            u1 = self.normalized()
            u2 = v.normalized()
            angle_in_radians = acos(u1.dot(u2))

            if in_degrees:
                degrees_per_radian = 180. / pi
                return angle_in_radians * degrees_per_radian
            else:
                return angle_in_radians

        except Exception as e:
            if str(e) == self.CANNOT_NORMALIZE_ZERO_VECTOR_MSG:
                raise Exception('Cannot comput an angle with a zero vector')
            else:
                raise e

    def is_orthogonal_to(self, v, tolerance=1e-10):
        return abs(self.dot(v)) < tolerance

    def is_parallel_to(self, v):
        return self.is_zero() or v.is_zero() or self.angle_with(v) == 0 or self.angle_with(v) == pi

    def is_zero(self, tolerance=1e-10):
        return self.magnitude() < tolerance


print('first pair...')
v = Vector(['-7.579', '-7.88'])
w = Vector(['22.737', '23.64'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('second pair...')
v = Vector(['-2.029', '9.97', '4.172'])
w = Vector(['-9.231', '-6.639', '-7.245'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('third pair...')
v = Vector(['-2.328', '-7.284', '-1.214'])
w = Vector(['-1.821', '1.072', '-2.94'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

print('fourth pair...')
v = Vector(['2.118', '4.827'])
w = Vector(['0', '0'])
print('is parallel:', v.is_parallel_to(w))
print('is orthogonal:', v.is_orthogonal_to(w))

Answer 1

难道u1.dot(u2)等于-1.00000000000000018058942747512吗

print(u2)
print(u1.dot(u2))
angle_in_radians = acos(u1.dot(u2))

这是在60号线附近。

更新，并进行进一步测试：

getcontext().prec = 16

......

def dot(self, v):
    print(self.coordinates, v.coordinates)
    print("asf")
    result = 0
    for x, y in zip(self.coordinates, v.coordinates):
        print("=================")
        print("x: ", x)
        print("y: ", y)
        print("x*y: ", x*y)
        result += (x*y)
        print("=================")
    print("Result: ", result)
    print(sum([x * y for x, y in zip(self.coordinates, v.coordinates)]))
    return sum([x * y for x, y in zip(self.coordinates, v.coordinates)])

在以下方面的成果：

=================
x:  -0.6932074151971374
y:  0.6932074151971375
x*y:  -0.4805365204842965
=================
=================
x:  -0.7207381490636552
y:  0.7207381490636553
x*y:  -0.5194634795157037
=================
Result:  -1.000000000000000
-1.000000000000000

但有：

getcontext().prec = 30

小数开始漂移。

=================
x:  -0.693207415197137377521618972764
y:  0.693207415197137482701372768190
x*y:  -0.480536520484296481693529594664
=================
=================
x:  -0.720738149063655170190045851086
y:  0.720738149063655279547013776664
x*y:  -0.519463479515703698895897880460
=================
Result:  -1.00000000000000018058942747512

这使得结果小于-1，破坏了acos()函数。

在发现浮点数退出后，我查看了您的代码，注意到了一些返回浮点数的函数。罪魁祸首是sqrt()函数，它没有足够高的精度。

def magnitude(self):
    coordinates_squared = [x ** 2 for x in self.coordinates]
    return Decimal(sum(coordinates_squared)).sqrt()

def normalized(self):
    try:
        magnitude = self.magnitude()
        return self.times_scalar(Decimal(1.0) / magnitude)

使用Decimal(x).sqrt()函数将解决您的问题。然后，您还需要更新normalized()函数。