1 /*++
2 Copyright (c) 2006 Microsoft Corporation
3
4 Module Name:
5
6 inf_rational.cpp
7
8 Abstract:
9
10 Rational numbers with infenitesimals
11
12 Author:
13
14 Nikolaj Bjorner (nbjorner) 2006-12-05.
15
16 Revision History:
17
18 --*/
19 #include "util/inf_rational.h"
20
21 inf_rational inf_rational::m_zero;
22 inf_rational inf_rational::m_one;
23 inf_rational inf_rational::m_minus_one;
24
inf_mult(inf_rational const & r1,inf_rational const & r2)25 inf_rational inf_mult(inf_rational const& r1, inf_rational const& r2)
26 {
27 inf_rational result;
28 result.m_first = r1.m_first * r2.m_first;
29 result.m_second = (r1.m_first * r2.m_second) + (r1.m_second * r2.m_first);
30
31 if (r1.m_second.is_pos() && r2.m_second.is_neg()) {
32 --result.m_second;
33 }
34 else if (r1.m_second.is_neg() && r2.m_second.is_pos()) {
35 --result.m_second;
36 }
37 return result;
38 }
39
sup_mult(inf_rational const & r1,inf_rational const & r2)40 inf_rational sup_mult(inf_rational const& r1, inf_rational const& r2)
41 {
42 inf_rational result;
43 result.m_first = r1.m_first * r2.m_first;
44 result.m_second = (r1.m_first * r2.m_second) + (r1.m_second * r2.m_first);
45
46 if (r1.m_second.is_pos() && r2.m_second.is_pos()) {
47 ++result.m_second;
48 }
49 else if (r1.m_second.is_neg() && r2.m_second.is_neg()) {
50 ++result.m_second;
51 }
52 return result;
53 }
54
55 //
56 // Find rationals c, x, such that c + epsilon*x <= r1/r2
57 //
58 // let r1 = a + d_1
59 // let r2 = b + d_2
60 //
61 // suppose b != 0:
62 //
63 // r1/b <= r1/r2
64 // <=> { if b > 0, then r2 > 0, and cross multiplication does not change the sign }
65 // { if b < 0, then r2 < 0, and cross multiplication changes sign twice }
66 // r1 * r2 <= b * r1
67 // <=>
68 // r1 * (b + d_2) <= r1 * b
69 // <=>
70 // r1 * d_2 <= 0
71 //
72 // if r1 * d_2 > 0, then r1/(b + sign_of(r1)*1/2*|b|) <= r1/r2
73 //
74 // Not handled here:
75 // if b = 0, then d_2 != 0
76 // if r1 * d_2 = 0 then it's 0.
77 // if r1 * d_2 > 0, then result is +oo
78 // if r1 * d_2 < 0, then result is -oo
79 //
inf_div(inf_rational const & r1,inf_rational const & r2)80 inf_rational inf_div(inf_rational const& r1, inf_rational const& r2)
81 {
82 SASSERT(!r2.m_first.is_zero());
83 inf_rational result;
84
85 if (r2.m_second.is_neg() && r1.is_neg()) {
86 result = r1 / (r2.m_first - (abs(r2.m_first)/rational(2)));
87 }
88 else if (r2.m_second.is_pos() && r1.is_pos()) {
89 result = r1 / (r2.m_first + (abs(r2.m_first)/rational(2)));
90 }
91 else {
92 result = r1 / r2.m_first;
93 }
94 return result;
95 }
96
sup_div(inf_rational const & r1,inf_rational const & r2)97 inf_rational sup_div(inf_rational const& r1, inf_rational const& r2)
98 {
99 SASSERT(!r2.m_first.is_zero());
100 inf_rational result;
101
102 if (r2.m_second.is_pos() && r1.is_neg()) {
103 result = r1 / (r2.m_first + (abs(r2.m_first)/rational(2)));
104 }
105 else if (r2.m_second.is_neg() && r1.is_pos()) {
106 result = r1 / (r2.m_first - (abs(r2.m_first)/rational(2)));
107 }
108 else {
109 result = r1 / r2.m_first;
110 }
111 return result;
112
113 }
114
inf_power(inf_rational const & r,unsigned n)115 inf_rational inf_power(inf_rational const& r, unsigned n)
116 {
117 bool is_even = (0 == (n & 0x1));
118 inf_rational result;
119 if (n == 1) {
120 result = r;
121 }
122 else if ((r.m_second.is_zero()) ||
123 (r.m_first.is_pos() && r.m_second.is_pos()) ||
124 (r.m_first.is_neg() && r.m_second.is_neg() && is_even)) {
125 result.m_first = r.m_first.expt(n);
126 }
127 else if (is_even) {
128 // 0 will work.
129 }
130 else if (r.m_first.is_zero()) {
131 result.m_first = rational::minus_one();
132 }
133 else if (r.m_first.is_pos()) {
134 result.m_first = rational(r.m_first - r.m_first/rational(2)).expt(n);
135 }
136 else {
137 result.m_first = rational(r.m_first + r.m_first/rational(2)).expt(n);
138 }
139 return result;
140 }
141
sup_power(inf_rational const & r,unsigned n)142 inf_rational sup_power(inf_rational const& r, unsigned n)
143 {
144 bool is_even = (0 == (n & 0x1));
145 inf_rational result;
146 if (n == 1) {
147 result = r;
148 }
149 else if (r.m_second.is_zero() ||
150 (r.m_first.is_pos() && r.m_second.is_neg()) ||
151 (r.m_first.is_neg() && r.m_second.is_pos() && is_even)) {
152 result.m_first = r.m_first.expt(n);
153 }
154 else if (r.m_first.is_zero() || (n == 0)) {
155 result.m_first = rational::one();
156 }
157 else if (r.m_first.is_pos() || is_even) {
158 result.m_first = rational(r.m_first + r.m_first/rational(2)).expt(n);
159 }
160 else {
161 // r (r.m_first) is negative, n is odd.
162 result.m_first = rational(r.m_first - r.m_first/rational(2)).expt(n);
163 }
164 return result;
165 }
166
inf_root(inf_rational const & r,unsigned n)167 inf_rational inf_root(inf_rational const& r, unsigned n)
168 {
169 SASSERT(!r.is_neg());
170 // use 0
171 return inf_rational();
172 }
173
sup_root(inf_rational const & r,unsigned n)174 inf_rational sup_root(inf_rational const& r, unsigned n)
175 {
176 SASSERT(!r.is_neg());
177 // use r.
178 return r;
179 }
180
initialize_inf_rational()181 void initialize_inf_rational() {
182 inf_rational::init();
183 }
184
init()185 void inf_rational::init() {
186 m_zero.m_first = rational::zero();
187 m_one.m_first = rational::one();
188 m_minus_one.m_first = rational::minus_one();
189 }
190
finalize_inf_rational()191 void finalize_inf_rational() {
192 inf_rational::finalize();
193 }
194
finalize()195 void inf_rational::finalize() {
196 m_zero.~inf_rational();
197 m_one.~inf_rational();
198 m_minus_one.~inf_rational();
199 }
200