← Back to branch summary

~ubuntu-branches/ubuntu/vivid/golang/vivid

« back to all changes in this revision

Viewing changes to src/pkg/math/big/rat.go

  • Committer: Package Import Robot
  • Author(s): James Page
  • Date: 2013-08-20 14:06:23 UTC
  • mfrom: (14.1.23 saucy-proposed)
  • Revision ID: package-import@ubuntu.com-20130820140623-b414jfxi3m0qkmrq
Tags: 2:1.1.2-2ubuntu1
https://launchpad.net/bugs/1211749
https://launchpad.net/bugs/1202027
* Merge from Debian unstable (LP: #1211749, #1202027). Remaining changes:
  - 016-armhf-elf-header.patch: Use correct ELF header for armhf binaries.
  - d/control,control.cross: Update Breaks/Replaces for Ubuntu
    versions to ensure smooth upgrades, regenerate control file.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pkg/math/big/rat.go
10
10
        "encoding/binary"
11
11
        "errors"
12
12
        "fmt"
 
13
        "math"
13
14
        "strings"
14
15
)
15
16
 
16
17
// A Rat represents a quotient a/b of arbitrary precision.
17
18
// The zero value for a Rat represents the value 0.
18
19
type Rat struct {
19
 
        a Int
20
 
        b nat // len(b) == 0 acts like b == 1
 
20
        // To make zero values for Rat work w/o initialization,
 
21
        // a zero value of b (len(b) == 0) acts like b == 1.
 
22
        // a.neg determines the sign of the Rat, b.neg is ignored.
 
23
        a, b Int
21
24
}
22
25
 
23
26
// NewRat creates a new Rat with numerator a and denominator b.
25
28
        return new(Rat).SetFrac64(a, b)
26
29
}
27
30
 
 
31
// SetFloat64 sets z to exactly f and returns z.
 
32
// If f is not finite, SetFloat returns nil.
 
33
func (z *Rat) SetFloat64(f float64) *Rat {
 
34
        const expMask = 1<<11 - 1
 
35
        bits := math.Float64bits(f)
 
36
        mantissa := bits & (1<<52 - 1)
 
37
        exp := int((bits >> 52) & expMask)
 
38
        switch exp {
 
39
        case expMask: // non-finite
 
40
                return nil
 
41
        case 0: // denormal
 
42
                exp -= 1022
 
43
        default: // normal
 
44
                mantissa |= 1 << 52
 
45
                exp -= 1023
 
46
        }
 
47
 
 
48
        shift := 52 - exp
 
49
 
 
50
        // Optimisation (?): partially pre-normalise.
 
51
        for mantissa&1 == 0 && shift > 0 {
 
52
                mantissa >>= 1
 
53
                shift--
 
54
        }
 
55
 
 
56
        z.a.SetUint64(mantissa)
 
57
        z.a.neg = f < 0
 
58
        z.b.Set(intOne)
 
59
        if shift > 0 {
 
60
                z.b.Lsh(&z.b, uint(shift))
 
61
        } else {
 
62
                z.a.Lsh(&z.a, uint(-shift))
 
63
        }
 
64
        return z.norm()
 
65
}
 
66
 
 
67
// isFinite reports whether f represents a finite rational value.
 
68
// It is equivalent to !math.IsNan(f) && !math.IsInf(f, 0).
 
69
func isFinite(f float64) bool {
 
70
        return math.Abs(f) <= math.MaxFloat64
 
71
}
 
72
 
 
73
// low64 returns the least significant 64 bits of natural number z.
 
74
func low64(z nat) uint64 {
 
75
        if len(z) == 0 {
 
76
                return 0
 
77
        }
 
78
        if _W == 32 && len(z) > 1 {
 
79
                return uint64(z[1])<<32 | uint64(z[0])
 
80
        }
 
81
        return uint64(z[0])
 
82
}
 
83
 
 
84
// quotToFloat returns the non-negative IEEE 754 double-precision
 
85
// value nearest to the quotient a/b, using round-to-even in halfway
 
86
// cases.  It does not mutate its arguments.
 
87
// Preconditions: b is non-zero; a and b have no common factors.
 
88
func quotToFloat(a, b nat) (f float64, exact bool) {
 
89
        // TODO(adonovan): specialize common degenerate cases: 1.0, integers.
 
90
        alen := a.bitLen()
 
91
        if alen == 0 {
 
92
                return 0, true
 
93
        }
 
94
        blen := b.bitLen()
 
95
        if blen == 0 {
 
96
                panic("division by zero")
 
97
        }
 
98
 
 
99
        // 1. Left-shift A or B such that quotient A/B is in [1<<53, 1<<55).
 
100
        // (54 bits if A<B when they are left-aligned, 55 bits if A>=B.)
 
101
        // This is 2 or 3 more than the float64 mantissa field width of 52:
 
102
        // - the optional extra bit is shifted away in step 3 below.
 
103
        // - the high-order 1 is omitted in float64 "normal" representation;
 
104
        // - the low-order 1 will be used during rounding then discarded.
 
105
        exp := alen - blen
 
106
        var a2, b2 nat
 
107
        a2 = a2.set(a)
 
108
        b2 = b2.set(b)
 
109
        if shift := 54 - exp; shift > 0 {
 
110
                a2 = a2.shl(a2, uint(shift))
 
111
        } else if shift < 0 {
 
112
                b2 = b2.shl(b2, uint(-shift))
 
113
        }
 
114
 
 
115
        // 2. Compute quotient and remainder (q, r).  NB: due to the
 
116
        // extra shift, the low-order bit of q is logically the
 
117
        // high-order bit of r.
 
118
        var q nat
 
119
        q, r := q.div(a2, a2, b2) // (recycle a2)
 
120
        mantissa := low64(q)
 
121
        haveRem := len(r) > 0 // mantissa&1 && !haveRem => remainder is exactly half
 
122
 
 
123
        // 3. If quotient didn't fit in 54 bits, re-do division by b2<<1
 
124
        // (in effect---we accomplish this incrementally).
 
125
        if mantissa>>54 == 1 {
 
126
                if mantissa&1 == 1 {
 
127
                        haveRem = true
 
128
                }
 
129
                mantissa >>= 1
 
130
                exp++
 
131
        }
 
132
        if mantissa>>53 != 1 {
 
133
                panic("expected exactly 54 bits of result")
 
134
        }
 
135
 
 
136
        // 4. Rounding.
 
137
        if -1022-52 <= exp && exp <= -1022 {
 
138
                // Denormal case; lose 'shift' bits of precision.
 
139
                shift := uint64(-1022 - (exp - 1)) // [1..53)
 
140
                lostbits := mantissa & (1<<shift - 1)
 
141
                haveRem = haveRem || lostbits != 0
 
142
                mantissa >>= shift
 
143
                exp = -1023 + 2
 
144
        }
 
145
        // Round q using round-half-to-even.
 
146
        exact = !haveRem
 
147
        if mantissa&1 != 0 {
 
148
                exact = false
 
149
                if haveRem || mantissa&2 != 0 {
 
150
                        if mantissa++; mantissa >= 1<<54 {
 
151
                                // Complete rollover 11...1 => 100...0, so shift is safe
 
152
                                mantissa >>= 1
 
153
                                exp++
 
154
                        }
 
155
                }
 
156
        }
 
157
        mantissa >>= 1 // discard rounding bit.  Mantissa now scaled by 2^53.
 
158
 
 
159
        f = math.Ldexp(float64(mantissa), exp-53)
 
160
        if math.IsInf(f, 0) {
 
161
                exact = false
 
162
        }
 
163
        return
 
164
}
 
165
 
 
166
// Float64 returns the nearest float64 value for x and a bool indicating
 
167
// whether f represents x exactly. The sign of f always matches the sign
 
168
// of x, even if f == 0.
 
169
func (x *Rat) Float64() (f float64, exact bool) {
 
170
        b := x.b.abs
 
171
        if len(b) == 0 {
 
172
                b = b.set(natOne) // materialize denominator
 
173
        }
 
174
        f, exact = quotToFloat(x.a.abs, b)
 
175
        if x.a.neg {
 
176
                f = -f
 
177
        }
 
178
        return
 
179
}
 
180
 
28
181
// SetFrac sets z to a/b and returns z.
29
182
func (z *Rat) SetFrac(a, b *Int) *Rat {
30
183
        z.a.neg = a.neg != b.neg
36
189
                babs = nat(nil).set(babs) // make a copy
37
190
        }
38
191
        z.a.abs = z.a.abs.set(a.abs)
39
 
        z.b = z.b.set(babs)
 
192
        z.b.abs = z.b.abs.set(babs)
40
193
        return z.norm()
41
194
}
42
195
 
50
203
                b = -b
51
204
                z.a.neg = !z.a.neg
52
205
        }
53
 
        z.b = z.b.setUint64(uint64(b))
 
206
        z.b.abs = z.b.abs.setUint64(uint64(b))
54
207
        return z.norm()
55
208
}
56
209
 
57
210
// SetInt sets z to x (by making a copy of x) and returns z.
58
211
func (z *Rat) SetInt(x *Int) *Rat {
59
212
        z.a.Set(x)
60
 
        z.b = z.b.make(0)
 
213
        z.b.abs = z.b.abs.make(0)
61
214
        return z
62
215
}
63
216
 
64
217
// SetInt64 sets z to x and returns z.
65
218
func (z *Rat) SetInt64(x int64) *Rat {
66
219
        z.a.SetInt64(x)
67
 
        z.b = z.b.make(0)
 
220
        z.b.abs = z.b.abs.make(0)
68
221
        return z
69
222
}
70
223
 
72
225
func (z *Rat) Set(x *Rat) *Rat {
73
226
        if z != x {
74
227
                z.a.Set(&x.a)
75
 
                z.b = z.b.set(x.b)
 
228
                z.b.Set(&x.b)
76
229
        }
77
230
        return z
78
231
}
97
250
                panic("division by zero")
98
251
        }
99
252
        z.Set(x)
100
 
        a := z.b
 
253
        a := z.b.abs
101
254
        if len(a) == 0 {
102
 
                a = a.setWord(1) // materialize numerator
 
255
                a = a.set(natOne) // materialize numerator
103
256
        }
104
257
        b := z.a.abs
105
258
        if b.cmp(natOne) == 0 {
106
259
                b = b.make(0) // normalize denominator
107
260
        }
108
 
        z.a.abs, z.b = a, b // sign doesn't change
 
261
        z.a.abs, z.b.abs = a, b // sign doesn't change
109
262
        return z
110
263
}
111
264
 
121
274
 
122
275
// IsInt returns true if the denominator of x is 1.
123
276
func (x *Rat) IsInt() bool {
124
 
        return len(x.b) == 0 || x.b.cmp(natOne) == 0
 
277
        return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
125
278
}
126
279
 
127
280
// Num returns the numerator of x; it may be <= 0.
128
281
// The result is a reference to x's numerator; it
129
 
// may change if a new value is assigned to x.
 
282
// may change if a new value is assigned to x, and vice versa.
 
283
// The sign of the numerator corresponds to the sign of x.
130
284
func (x *Rat) Num() *Int {
131
285
        return &x.a
132
286
}
133
287
 
134
288
// Denom returns the denominator of x; it is always > 0.
135
289
// The result is a reference to x's denominator; it
136
 
// may change if a new value is assigned to x.
 
290
// may change if a new value is assigned to x, and vice versa.
137
291
func (x *Rat) Denom() *Int {
138
 
        if len(x.b) == 0 {
139
 
                return &Int{abs: nat{1}}
140
 
        }
141
 
        return &Int{abs: x.b}
142
 
}
143
 
 
144
 
func gcd(x, y nat) nat {
145
 
        // Euclidean algorithm.
146
 
        var a, b nat
147
 
        a = a.set(x)
148
 
        b = b.set(y)
149
 
        for len(b) != 0 {
150
 
                var q, r nat
151
 
                _, r = q.div(r, a, b)
152
 
                a = b
153
 
                b = r
154
 
        }
155
 
        return a
 
292
        x.b.neg = false // the result is always >= 0
 
293
        if len(x.b.abs) == 0 {
 
294
                x.b.abs = x.b.abs.set(natOne) // materialize denominator
 
295
        }
 
296
        return &x.b
156
297
}
157
298
 
158
299
func (z *Rat) norm() *Rat {
160
301
        case len(z.a.abs) == 0:
161
302
                // z == 0 - normalize sign and denominator
162
303
                z.a.neg = false
163
 
                z.b = z.b.make(0)
164
 
        case len(z.b) == 0:
 
304
                z.b.abs = z.b.abs.make(0)
 
305
        case len(z.b.abs) == 0:
165
306
                // z is normalized int - nothing to do
166
 
        case z.b.cmp(natOne) == 0:
 
307
        case z.b.abs.cmp(natOne) == 0:
167
308
                // z is int - normalize denominator
168
 
                z.b = z.b.make(0)
 
309
                z.b.abs = z.b.abs.make(0)
169
310
        default:
170
 
                if f := gcd(z.a.abs, z.b); f.cmp(natOne) != 0 {
171
 
                        z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f)
172
 
                        z.b, _ = z.b.div(nil, z.b, f)
 
311
                neg := z.a.neg
 
312
                z.a.neg = false
 
313
                z.b.neg = false
 
314
                if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
 
315
                        z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
 
316
                        z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
 
317
                        if z.b.abs.cmp(natOne) == 0 {
 
318
                                // z is int - normalize denominator
 
319
                                z.b.abs = z.b.abs.make(0)
 
320
                        }
173
321
                }
 
322
                z.a.neg = neg
174
323
        }
175
324
        return z
176
325
}
207
356
//   +1 if x >  y
208
357
//
209
358
func (x *Rat) Cmp(y *Rat) int {
210
 
        return scaleDenom(&x.a, y.b).Cmp(scaleDenom(&y.a, x.b))
 
359
        return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
211
360
}
212
361
 
213
362
// Add sets z to the sum x+y and returns z.
214
363
func (z *Rat) Add(x, y *Rat) *Rat {
215
 
        a1 := scaleDenom(&x.a, y.b)
216
 
        a2 := scaleDenom(&y.a, x.b)
 
364
        a1 := scaleDenom(&x.a, y.b.abs)
 
365
        a2 := scaleDenom(&y.a, x.b.abs)
217
366
        z.a.Add(a1, a2)
218
 
        z.b = mulDenom(z.b, x.b, y.b)
 
367
        z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
219
368
        return z.norm()
220
369
}
221
370
 
222
371
// Sub sets z to the difference x-y and returns z.
223
372
func (z *Rat) Sub(x, y *Rat) *Rat {
224
 
        a1 := scaleDenom(&x.a, y.b)
225
 
        a2 := scaleDenom(&y.a, x.b)
 
373
        a1 := scaleDenom(&x.a, y.b.abs)
 
374
        a2 := scaleDenom(&y.a, x.b.abs)
226
375
        z.a.Sub(a1, a2)
227
 
        z.b = mulDenom(z.b, x.b, y.b)
 
376
        z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
228
377
        return z.norm()
229
378
}
230
379
 
231
380
// Mul sets z to the product x*y and returns z.
232
381
func (z *Rat) Mul(x, y *Rat) *Rat {
233
382
        z.a.Mul(&x.a, &y.a)
234
 
        z.b = mulDenom(z.b, x.b, y.b)
 
383
        z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
235
384
        return z.norm()
236
385
}
237
386
 
241
390
        if len(y.a.abs) == 0 {
242
391
                panic("division by zero")
243
392
        }
244
 
        a := scaleDenom(&x.a, y.b)
245
 
        b := scaleDenom(&y.a, x.b)
 
393
        a := scaleDenom(&x.a, y.b.abs)
 
394
        b := scaleDenom(&y.a, x.b.abs)
246
395
        z.a.abs = a.abs
247
 
        z.b = b.abs
 
396
        z.b.abs = b.abs
248
397
        z.a.neg = a.neg != b.neg
249
398
        return z.norm()
250
399
}
286
435
                }
287
436
                s = s[sep+1:]
288
437
                var err error
289
 
                if z.b, _, err = z.b.scan(strings.NewReader(s), 10); err != nil {
 
438
                if z.b.abs, _, err = z.b.abs.scan(strings.NewReader(s), 10); err != nil {
290
439
                        return nil, false
291
440
                }
292
441
                return z.norm(), true
317
466
        }
318
467
        powTen := nat(nil).expNN(natTen, exp.abs, nil)
319
468
        if exp.neg {
320
 
                z.b = powTen
 
469
                z.b.abs = powTen
321
470
                z.norm()
322
471
        } else {
323
472
                z.a.abs = z.a.abs.mul(z.a.abs, powTen)
324
 
                z.b = z.b.make(0)
 
473
                z.b.abs = z.b.abs.make(0)
325
474
        }
326
475
 
327
476
        return z, true
330
479
// String returns a string representation of z in the form "a/b" (even if b == 1).
331
480
func (x *Rat) String() string {
332
481
        s := "/1"
333
 
        if len(x.b) != 0 {
334
 
                s = "/" + x.b.decimalString()
 
482
        if len(x.b.abs) != 0 {
 
483
                s = "/" + x.b.abs.decimalString()
335
484
        }
336
485
        return x.a.String() + s
337
486
}
355
504
                }
356
505
                return s
357
506
        }
358
 
        // x.b != 0
 
507
        // x.b.abs != 0
359
508
 
360
 
        q, r := nat(nil).div(nat(nil), x.a.abs, x.b)
 
509
        q, r := nat(nil).div(nat(nil), x.a.abs, x.b.abs)
361
510
 
362
511
        p := natOne
363
512
        if prec > 0 {
365
514
        }
366
515
 
367
516
        r = r.mul(r, p)
368
 
        r, r2 := r.div(nat(nil), r, x.b)
 
517
        r, r2 := r.div(nat(nil), r, x.b.abs)
369
518
 
370
519
        // see if we need to round up
371
520
        r2 = r2.add(r2, r2)
372
 
        if x.b.cmp(r2) <= 0 {
 
521
        if x.b.abs.cmp(r2) <= 0 {
373
522
                r = r.add(r, natOne)
374
523
                if r.cmp(p) >= 0 {
375
524
                        q = nat(nil).add(q, natOne)
396
545
 
397
546
// GobEncode implements the gob.GobEncoder interface.
398
547
func (x *Rat) GobEncode() ([]byte, error) {
399
 
        buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
400
 
        i := x.b.bytes(buf)
 
548
        buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b.abs))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
 
549
        i := x.b.abs.bytes(buf)
401
550
        j := x.a.abs.bytes(buf[0:i])
402
551
        n := i - j
403
552
        if int(uint32(n)) != n {
427
576
        i := j + binary.BigEndian.Uint32(buf[j-4:j])
428
577
        z.a.neg = b&1 != 0
429
578
        z.a.abs = z.a.abs.setBytes(buf[j:i])
430
 
        z.b = z.b.setBytes(buf[i:])
 
579
        z.b.abs = z.b.abs.setBytes(buf[i:])
431
580
        return nil
432
581
}