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