Source file
src/math/big/rat.go
1
2
3
4
5
6
7 package big
8
9 import (
10 "fmt"
11 "math"
12 )
13
14
15
16 type Rat struct {
17
18
19
20 a, b Int
21 }
22
23
24 func NewRat(a, b int64) *Rat {
25 return new(Rat).SetFrac64(a, b)
26 }
27
28
29
30 func (z *Rat) SetFloat64(f float64) *Rat {
31 const expMask = 1<<11 - 1
32 bits := math.Float64bits(f)
33 mantissa := bits & (1<<52 - 1)
34 exp := int((bits >> 52) & expMask)
35 switch exp {
36 case expMask:
37 return nil
38 case 0:
39 exp -= 1022
40 default:
41 mantissa |= 1 << 52
42 exp -= 1023
43 }
44
45 shift := 52 - exp
46
47
48 for mantissa&1 == 0 && shift > 0 {
49 mantissa >>= 1
50 shift--
51 }
52
53 z.a.SetUint64(mantissa)
54 z.a.neg = f < 0
55 z.b.Set(intOne)
56 if shift > 0 {
57 z.b.Lsh(&z.b, uint(shift))
58 } else {
59 z.a.Lsh(&z.a, uint(-shift))
60 }
61 return z.norm()
62 }
63
64
65
66
67
68 func quotToFloat32(a, b nat) (f float32, exact bool) {
69 const (
70
71 Fsize = 32
72
73
74 Msize = 23
75 Msize1 = Msize + 1
76 Msize2 = Msize1 + 1
77
78
79 Esize = Fsize - Msize1
80 Ebias = 1<<(Esize-1) - 1
81 Emin = 1 - Ebias
82 Emax = Ebias
83 )
84
85
86 alen := a.bitLen()
87 if alen == 0 {
88 return 0, true
89 }
90 blen := b.bitLen()
91 if blen == 0 {
92 panic("division by zero")
93 }
94
95
96
97
98
99
100
101 exp := alen - blen
102 var a2, b2 nat
103 a2 = a2.set(a)
104 b2 = b2.set(b)
105 if shift := Msize2 - exp; shift > 0 {
106 a2 = a2.shl(a2, uint(shift))
107 } else if shift < 0 {
108 b2 = b2.shl(b2, uint(-shift))
109 }
110
111
112
113
114 var q nat
115 q, r := q.div(a2, a2, b2)
116 mantissa := low32(q)
117 haveRem := len(r) > 0
118
119
120
121 if mantissa>>Msize2 == 1 {
122 if mantissa&1 == 1 {
123 haveRem = true
124 }
125 mantissa >>= 1
126 exp++
127 }
128 if mantissa>>Msize1 != 1 {
129 panic(fmt.Sprintf("expected exactly %d bits of result", Msize2))
130 }
131
132
133 if Emin-Msize <= exp && exp <= Emin {
134
135 shift := uint(Emin - (exp - 1))
136 lostbits := mantissa & (1<<shift - 1)
137 haveRem = haveRem || lostbits != 0
138 mantissa >>= shift
139 exp = 2 - Ebias
140 }
141
142 exact = !haveRem
143 if mantissa&1 != 0 {
144 exact = false
145 if haveRem || mantissa&2 != 0 {
146 if mantissa++; mantissa >= 1<<Msize2 {
147
148 mantissa >>= 1
149 exp++
150 }
151 }
152 }
153 mantissa >>= 1
154
155 f = float32(math.Ldexp(float64(mantissa), exp-Msize1))
156 if math.IsInf(float64(f), 0) {
157 exact = false
158 }
159 return
160 }
161
162
163
164
165
166 func quotToFloat64(a, b nat) (f float64, exact bool) {
167 const (
168
169 Fsize = 64
170
171
172 Msize = 52
173 Msize1 = Msize + 1
174 Msize2 = Msize1 + 1
175
176
177 Esize = Fsize - Msize1
178 Ebias = 1<<(Esize-1) - 1
179 Emin = 1 - Ebias
180 Emax = Ebias
181 )
182
183
184 alen := a.bitLen()
185 if alen == 0 {
186 return 0, true
187 }
188 blen := b.bitLen()
189 if blen == 0 {
190 panic("division by zero")
191 }
192
193
194
195
196
197
198
199 exp := alen - blen
200 var a2, b2 nat
201 a2 = a2.set(a)
202 b2 = b2.set(b)
203 if shift := Msize2 - exp; shift > 0 {
204 a2 = a2.shl(a2, uint(shift))
205 } else if shift < 0 {
206 b2 = b2.shl(b2, uint(-shift))
207 }
208
209
210
211
212 var q nat
213 q, r := q.div(a2, a2, b2)
214 mantissa := low64(q)
215 haveRem := len(r) > 0
216
217
218
219 if mantissa>>Msize2 == 1 {
220 if mantissa&1 == 1 {
221 haveRem = true
222 }
223 mantissa >>= 1
224 exp++
225 }
226 if mantissa>>Msize1 != 1 {
227 panic(fmt.Sprintf("expected exactly %d bits of result", Msize2))
228 }
229
230
231 if Emin-Msize <= exp && exp <= Emin {
232
233 shift := uint(Emin - (exp - 1))
234 lostbits := mantissa & (1<<shift - 1)
235 haveRem = haveRem || lostbits != 0
236 mantissa >>= shift
237 exp = 2 - Ebias
238 }
239
240 exact = !haveRem
241 if mantissa&1 != 0 {
242 exact = false
243 if haveRem || mantissa&2 != 0 {
244 if mantissa++; mantissa >= 1<<Msize2 {
245
246 mantissa >>= 1
247 exp++
248 }
249 }
250 }
251 mantissa >>= 1
252
253 f = math.Ldexp(float64(mantissa), exp-Msize1)
254 if math.IsInf(f, 0) {
255 exact = false
256 }
257 return
258 }
259
260
261
262
263
264 func (x *Rat) Float32() (f float32, exact bool) {
265 b := x.b.abs
266 if len(b) == 0 {
267 b = b.set(natOne)
268 }
269 f, exact = quotToFloat32(x.a.abs, b)
270 if x.a.neg {
271 f = -f
272 }
273 return
274 }
275
276
277
278
279
280 func (x *Rat) Float64() (f float64, exact bool) {
281 b := x.b.abs
282 if len(b) == 0 {
283 b = b.set(natOne)
284 }
285 f, exact = quotToFloat64(x.a.abs, b)
286 if x.a.neg {
287 f = -f
288 }
289 return
290 }
291
292
293 func (z *Rat) SetFrac(a, b *Int) *Rat {
294 z.a.neg = a.neg != b.neg
295 babs := b.abs
296 if len(babs) == 0 {
297 panic("division by zero")
298 }
299 if &z.a == b || alias(z.a.abs, babs) {
300 babs = nat(nil).set(babs)
301 }
302 z.a.abs = z.a.abs.set(a.abs)
303 z.b.abs = z.b.abs.set(babs)
304 return z.norm()
305 }
306
307
308 func (z *Rat) SetFrac64(a, b int64) *Rat {
309 z.a.SetInt64(a)
310 if b == 0 {
311 panic("division by zero")
312 }
313 if b < 0 {
314 b = -b
315 z.a.neg = !z.a.neg
316 }
317 z.b.abs = z.b.abs.setUint64(uint64(b))
318 return z.norm()
319 }
320
321
322 func (z *Rat) SetInt(x *Int) *Rat {
323 z.a.Set(x)
324 z.b.abs = z.b.abs[:0]
325 return z
326 }
327
328
329 func (z *Rat) SetInt64(x int64) *Rat {
330 z.a.SetInt64(x)
331 z.b.abs = z.b.abs[:0]
332 return z
333 }
334
335
336 func (z *Rat) Set(x *Rat) *Rat {
337 if z != x {
338 z.a.Set(&x.a)
339 z.b.Set(&x.b)
340 }
341 return z
342 }
343
344
345 func (z *Rat) Abs(x *Rat) *Rat {
346 z.Set(x)
347 z.a.neg = false
348 return z
349 }
350
351
352 func (z *Rat) Neg(x *Rat) *Rat {
353 z.Set(x)
354 z.a.neg = len(z.a.abs) > 0 && !z.a.neg
355 return z
356 }
357
358
359 func (z *Rat) Inv(x *Rat) *Rat {
360 if len(x.a.abs) == 0 {
361 panic("division by zero")
362 }
363 z.Set(x)
364 a := z.b.abs
365 if len(a) == 0 {
366 a = a.set(natOne)
367 }
368 b := z.a.abs
369 if b.cmp(natOne) == 0 {
370 b = b[:0]
371 }
372 z.a.abs, z.b.abs = a, b
373 return z
374 }
375
376
377
378
379
380
381
382 func (x *Rat) Sign() int {
383 return x.a.Sign()
384 }
385
386
387 func (x *Rat) IsInt() bool {
388 return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
389 }
390
391
392
393
394
395 func (x *Rat) Num() *Int {
396 return &x.a
397 }
398
399
400
401
402 func (x *Rat) Denom() *Int {
403 x.b.neg = false
404 if len(x.b.abs) == 0 {
405 x.b.abs = x.b.abs.set(natOne)
406 }
407 return &x.b
408 }
409
410 func (z *Rat) norm() *Rat {
411 switch {
412 case len(z.a.abs) == 0:
413
414 z.a.neg = false
415 z.b.abs = z.b.abs[:0]
416 case len(z.b.abs) == 0:
417
418 case z.b.abs.cmp(natOne) == 0:
419
420 z.b.abs = z.b.abs[:0]
421 default:
422 neg := z.a.neg
423 z.a.neg = false
424 z.b.neg = false
425 if f := NewInt(0).lehmerGCD(nil, nil, &z.a, &z.b); f.Cmp(intOne) != 0 {
426 z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
427 z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
428 if z.b.abs.cmp(natOne) == 0 {
429
430 z.b.abs = z.b.abs[:0]
431 }
432 }
433 z.a.neg = neg
434 }
435 return z
436 }
437
438
439
440
441 func mulDenom(z, x, y nat) nat {
442 switch {
443 case len(x) == 0:
444 return z.set(y)
445 case len(y) == 0:
446 return z.set(x)
447 }
448 return z.mul(x, y)
449 }
450
451
452
453 func scaleDenom(x *Int, f nat) *Int {
454 var z Int
455 if len(f) == 0 {
456 return z.Set(x)
457 }
458 z.abs = z.abs.mul(x.abs, f)
459 z.neg = x.neg
460 return &z
461 }
462
463
464
465
466
467
468
469 func (x *Rat) Cmp(y *Rat) int {
470 return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
471 }
472
473
474 func (z *Rat) Add(x, y *Rat) *Rat {
475 a1 := scaleDenom(&x.a, y.b.abs)
476 a2 := scaleDenom(&y.a, x.b.abs)
477 z.a.Add(a1, a2)
478 z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
479 return z.norm()
480 }
481
482
483 func (z *Rat) Sub(x, y *Rat) *Rat {
484 a1 := scaleDenom(&x.a, y.b.abs)
485 a2 := scaleDenom(&y.a, x.b.abs)
486 z.a.Sub(a1, a2)
487 z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
488 return z.norm()
489 }
490
491
492 func (z *Rat) Mul(x, y *Rat) *Rat {
493 if x == y {
494
495 z.a.neg = false
496 z.a.abs = z.a.abs.sqr(x.a.abs)
497 z.b.abs = z.b.abs.sqr(x.b.abs)
498 return z
499 }
500 z.a.Mul(&x.a, &y.a)
501 z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
502 return z.norm()
503 }
504
505
506
507 func (z *Rat) Quo(x, y *Rat) *Rat {
508 if len(y.a.abs) == 0 {
509 panic("division by zero")
510 }
511 a := scaleDenom(&x.a, y.b.abs)
512 b := scaleDenom(&y.a, x.b.abs)
513 z.a.abs = a.abs
514 z.b.abs = b.abs
515 z.a.neg = a.neg != b.neg
516 return z.norm()
517 }
518
View as plain text