inapy/chapter04/interval.py at master · tkouya/inapy · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
# interval.py: 区間演算クラス
# R.T.Kneusel, "Numbers and Computers", Springer, 2015.
# Interval -> [left, right]
import rmode
import math


class Interval:

    # 開始前のデフォルト丸めモード
    default_rmode = rmode.get_rmode()

    # コンストラクタ
    def __init__(self, left, right=None):
        if right is None:
            self.left = left
            self.right = left
        else:
            self.left = left
            self.right = right

    # + : 加算
    def __add__(self, y):
        rmode.set_rmode(rmode.FE_DOWNWARD)
        left = self.left + y.left
        rmode.set_rmode(rmode.FE_UPWARD)
        right = self.right + y.right
        rmode.set_rmode(Interval.default_rmode)

        return Interval(left, right)

    # - : 減算
    def __sub__(self, y):
        rmode.set_rmode(rmode.FE_DOWNWARD)
        left = self.left - y.right
        rmode.set_rmode(rmode.FE_UPWARD)
        right = self.right - y.left
        rmode.set_rmode(Interval.default_rmode)

        return Interval(left, right)

    # - :マイナス
    def __neg__(self):
        left = -self.right
        right = -self.left

        return Interval(left, right)

    # * : 乗算
    def __mul__(self, y):
        if not isinstance(y, Interval):
            y = Interval(float(y))
        rmode.set_rmode(rmode.FE_DOWNWARD)
        left = min([
            self.left * y.left,
            self.right * y.left,
            self.left * y.right,
            self.right * y.right
        ])
        rmode.set_rmode(rmode.FE_UPWARD)
        right = max([
            self.left * y.left,
            self.right * y.left,
            self.left * y.right,
            self.right * y.right
        ])
        rmode.set_rmode(Interval.default_rmode)

        return Interval(left, right)

    __rmul__ = __mul__

    # recip(x) = 1/x: 逆数
    def recip(self):
        if (self.left < 0.0) and (self.right > 0.0):
            return Interval(-math.inf, math.inf)

        rmode.set_rmode(rmode.FE_DOWNWARD)
        left = 1.0 / self.right
        rmode.set_rmode(rmode.FE_UPWARD)
        right = 1.0 / self.left
        rmode.set_rmode(Interval.default_rmode)

        return Interval(left, right)

    # / : 除算
    def __truediv__(self, y):
        return self.__mul__(y.recip())

    # 絶対値
    def abs(self):
        return max([math.abs(self.left), math.abs(self.right)])

    # 平方根
    def sqrt(self):
        if (self.left < 0) or (self.right < 0):
            return Interval(-math.nan, -math.nan)

        return Interval(math.sqrt(self.left), math.sqrt(self.right))

    # ** : べき乗
    def __pow__(self, n):
        if (n % 2) == 1:
            rmode.set_rmode(rmode.FE_DOWNWARD)
            left = self.left ** n
            rmode.set_rmode(rmode.FE_UPWARD)
            right = self.right ** n
        elif self.left >= 0:
            rmode.set_rmode(rmode.FE_DOWNWARD)
            left = self.left ** n
            rmode.set_rmode(rmode.FE_UPWARD)
            right = self.right ** n
        elif self.right < 0:
            rmode.set_rmode(rmode.FE_DOWNWARD)
            left = self.right ** n
            rmode.set_rmode(rmode.FE_UPWARD)
            right = self.left ** n
        else:
            left = 0.0
            rmode.set_rmode(rmode.FE_UPWARD)
            right = self.left ** n
            temp = self.right ** n
            if temp > right:
                right = temp

        rmode.set_rmode(Interval.default_rmode)
        return Interval(left, right)

    # 出力用
    def __str__(self):
        return "[ %g, %g ]" % (self.left, self.right)

    __repr__ = __str__


# -------------------------------------
# Copyright (c) 2021 Tomonori Kouya
# All rights reserved.
# -------------------------------------