← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/pkg/math/big/nat.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/nat.go
236
236
// Operands that are shorter than karatsubaThreshold are multiplied using
237
237
// "grade school" multiplication; for longer operands the Karatsuba algorithm
238
238
// is used.
239
 
var karatsubaThreshold int = 32 // computed by calibrate.go
 
239
var karatsubaThreshold int = 40 // computed by calibrate.go
240
240
 
241
241
// karatsuba multiplies x and y and leaves the result in z.
242
242
// Both x and y must have the same length n and n must be a
342
342
        return cap(x) > 0 && cap(y) > 0 && &x[0:cap(x)][cap(x)-1] == &y[0:cap(y)][cap(y)-1]
343
343
}
344
344
 
345
 
// addAt implements z += x*(1<<(_W*i)); z must be long enough.
 
345
// addAt implements z += x<<(_W*i); z must be long enough.
346
346
// (we don't use nat.add because we need z to stay the same
347
347
// slice, and we don't need to normalize z after each addition)
348
348
func addAt(z, x nat, i int) {
396
396
        }
397
397
 
398
398
        // use basic multiplication if the numbers are small
399
 
        if n < karatsubaThreshold || n < 2 {
 
399
        if n < karatsubaThreshold {
400
400
                z = z.make(m + n)
401
401
                basicMul(z, x, y)
402
402
                return z.norm()
405
405
 
406
406
        // determine Karatsuba length k such that
407
407
        //
408
 
        //   x = x1*b + x0
409
 
        //   y = y1*b + y0  (and k <= len(y), which implies k <= len(x))
 
408
        //   x = xh*b + x0  (0 <= x0 < b)
 
409
        //   y = yh*b + y0  (0 <= y0 < b)
410
410
        //   b = 1<<(_W*k)  ("base" of digits xi, yi)
411
411
        //
412
412
        k := karatsubaLen(n)
417
417
        y0 := y[0:k]              // y0 is not normalized
418
418
        z = z.make(max(6*k, m+n)) // enough space for karatsuba of x0*y0 and full result of x*y
419
419
        karatsuba(z, x0, y0)
420
 
        z = z[0 : m+n] // z has final length but may be incomplete, upper portion is garbage
 
420
        z = z[0 : m+n]  // z has final length but may be incomplete
 
421
        z[2*k:].clear() // upper portion of z is garbage (and 2*k <= m+n since k <= n <= m)
421
422
 
422
 
        // If x1 and/or y1 are not 0, add missing terms to z explicitly:
423
 
        //
424
 
        //     m+n       2*k       0
425
 
        //   z = [   ...   | x0*y0 ]
426
 
        //     +   [ x1*y1 ]
427
 
        //     +   [ x1*y0 ]
428
 
        //     +   [ x0*y1 ]
 
423
        // If xh != 0 or yh != 0, add the missing terms to z. For
 
424
        //
 
425
        //   xh = xi*b^i + ... + x2*b^2 + x1*b (0 <= xi < b)
 
426
        //   yh =                         y1*b (0 <= y1 < b)
 
427
        //
 
428
        // the missing terms are
 
429
        //
 
430
        //   x0*y1*b and xi*y0*b^i, xi*y1*b^(i+1) for i > 0
 
431
        //
 
432
        // since all the yi for i > 1 are 0 by choice of k: If any of them
 
433
        // were > 0, then yh >= b^2 and thus y >= b^2. Then k' = k*2 would
 
434
        // be a larger valid threshold contradicting the assumption about k.
429
435
        //
430
436
        if k < n || m != n {
431
 
                x1 := x[k:] // x1 is normalized because x is
432
 
                y1 := y[k:] // y1 is normalized because y is
433
437
                var t nat
434
 
                t = t.mul(x1, y1)
435
 
                copy(z[2*k:], t)
436
 
                z[2*k+len(t):].clear() // upper portion of z is garbage
437
 
                t = t.mul(x1, y0.norm())
438
 
                addAt(z, t, k)
439
 
                t = t.mul(x0.norm(), y1)
440
 
                addAt(z, t, k)
 
438
 
 
439
                // add x0*y1*b
 
440
                x0 := x0.norm()
 
441
                y1 := y[k:]       // y1 is normalized because y is
 
442
                t = t.mul(x0, y1) // update t so we don't lose t's underlying array
 
443
                addAt(z, t, k)
 
444
 
 
445
                // add xi*y0<<i, xi*y1*b<<(i+k)
 
446
                y0 := y0.norm()
 
447
                for i := k; i < len(x); i += k {
 
448
                        xi := x[i:]
 
449
                        if len(xi) > k {
 
450
                                xi = xi[:k]
 
451
                        }
 
452
                        xi = xi.norm()
 
453
                        t = t.mul(xi, y0)
 
454
                        addAt(z, t, i)
 
455
                        t = t.mul(xi, y1)
 
456
                        addAt(z, t, i+k)
 
457
                }
441
458
        }
442
459
 
443
460
        return z.norm()
493
510
        }
494
511
 
495
512
        if len(v) == 1 {
496
 
                var rprime Word
497
 
                q, rprime = z.divW(u, v[0])
498
 
                if rprime > 0 {
499
 
                        r = z2.make(1)
500
 
                        r[0] = rprime
501
 
                } else {
502
 
                        r = z2.make(0)
503
 
                }
 
513
                var r2 Word
 
514
                q, r2 = z.divW(u, v[0])
 
515
                r = z2.setWord(r2)
504
516
                return
505
517
        }
506
518
 
740
752
        // convert power of two and non power of two bases separately
741
753
        if b == b&-b {
742
754
                // shift is base-b digit size in bits
743
 
                shift := uint(trailingZeroBits(b)) // shift > 0 because b >= 2
 
755
                shift := trailingZeroBits(b) // shift > 0 because b >= 2
744
756
                mask := Word(1)<<shift - 1
745
757
                w := x[0]
746
758
                nbits := uint(_W) // number of unprocessed bits in w
814
826
 
815
827
// Convert words of q to base b digits in s. If q is large, it is recursively "split in half"
816
828
// by nat/nat division using tabulated divisors. Otherwise, it is converted iteratively using
817
 
// repeated nat/Word divison.
 
829
// repeated nat/Word division.
818
830
//
819
 
// The iterative method processes n Words by n divW() calls, each of which visits every Word in the 
820
 
// incrementally shortened q for a total of n + (n-1) + (n-2) ... + 2 + 1, or n(n+1)/2 divW()'s. 
821
 
// Recursive conversion divides q by its approximate square root, yielding two parts, each half 
 
831
// The iterative method processes n Words by n divW() calls, each of which visits every Word in the
 
832
// incrementally shortened q for a total of n + (n-1) + (n-2) ... + 2 + 1, or n(n+1)/2 divW()'s.
 
833
// Recursive conversion divides q by its approximate square root, yielding two parts, each half
822
834
// the size of q. Using the iterative method on both halves means 2 * (n/2)(n/2 + 1)/2 divW()'s
823
835
// plus the expensive long div(). Asymptotically, the ratio is favorable at 1/2 the divW()'s, and
824
 
// is made better by splitting the subblocks recursively. Best is to split blocks until one more 
825
 
// split would take longer (because of the nat/nat div()) than the twice as many divW()'s of the 
826
 
// iterative approach. This threshold is represented by leafSize. Benchmarking of leafSize in the 
827
 
// range 2..64 shows that values of 8 and 16 work well, with a 4x speedup at medium lengths and 
828
 
// ~30x for 20000 digits. Use nat_test.go's BenchmarkLeafSize tests to optimize leafSize for 
 
836
// is made better by splitting the subblocks recursively. Best is to split blocks until one more
 
837
// split would take longer (because of the nat/nat div()) than the twice as many divW()'s of the
 
838
// iterative approach. This threshold is represented by leafSize. Benchmarking of leafSize in the
 
839
// range 2..64 shows that values of 8 and 16 work well, with a 4x speedup at medium lengths and
 
840
// ~30x for 20000 digits. Use nat_test.go's BenchmarkLeafSize tests to optimize leafSize for
829
841
// specific hardware.
830
842
//
831
843
func (q nat) convertWords(s []byte, charset string, b Word, ndigits int, bb Word, table []divisor) {
908
920
        ndigits int // digit length of divisor in terms of output base digits
909
921
}
910
922
 
911
 
var cacheBase10 [64]divisor // cached divisors for base 10
912
 
var cacheLock sync.Mutex    // protects cacheBase10
 
923
var cacheBase10 struct {
 
924
        sync.Mutex
 
925
        table [64]divisor // cached divisors for base 10
 
926
}
913
927
 
914
928
// expWW computes x**y
915
929
func (z nat) expWW(x, y Word) nat {
925
939
 
926
940
        // determine k where (bb**leafSize)**(2**k) >= sqrt(x)
927
941
        k := 1
928
 
        for words := leafSize; words < m>>1 && k < len(cacheBase10); words <<= 1 {
 
942
        for words := leafSize; words < m>>1 && k < len(cacheBase10.table); words <<= 1 {
929
943
                k++
930
944
        }
931
945
 
932
 
        // create new table of divisors or extend and reuse existing table as appropriate
933
 
        var table []divisor
934
 
        var cached bool
935
 
        switch b {
936
 
        case 10:
937
 
                table = cacheBase10[0:k] // reuse old table for this conversion
938
 
                cached = true
939
 
        default:
940
 
                table = make([]divisor, k) // new table for this conversion
 
946
        // reuse and extend existing table of divisors or create new table as appropriate
 
947
        var table []divisor // for b == 10, table overlaps with cacheBase10.table
 
948
        if b == 10 {
 
949
                cacheBase10.Lock()
 
950
                table = cacheBase10.table[0:k] // reuse old table for this conversion
 
951
        } else {
 
952
                table = make([]divisor, k) // create new table for this conversion
941
953
        }
942
954
 
943
955
        // extend table
944
956
        if table[k-1].ndigits == 0 {
945
 
                if cached {
946
 
                        cacheLock.Lock() // begin critical section
947
 
                }
948
 
 
949
957
                // add new entries as needed
950
958
                var larger nat
951
959
                for i := 0; i < k; i++ {
952
960
                        if table[i].ndigits == 0 {
953
961
                                if i == 0 {
954
 
                                        table[i].bbb = nat(nil).expWW(bb, Word(leafSize))
955
 
                                        table[i].ndigits = ndigits * leafSize
 
962
                                        table[0].bbb = nat(nil).expWW(bb, Word(leafSize))
 
963
                                        table[0].ndigits = ndigits * leafSize
956
964
                                } else {
957
965
                                        table[i].bbb = nat(nil).mul(table[i-1].bbb, table[i-1].bbb)
958
966
                                        table[i].ndigits = 2 * table[i-1].ndigits
968
976
                                table[i].nbits = table[i].bbb.bitLen()
969
977
                        }
970
978
                }
 
979
        }
971
980
 
972
 
                if cached {
973
 
                        cacheLock.Unlock() // end critical section
974
 
                }
 
981
        if b == 10 {
 
982
                cacheBase10.Unlock()
975
983
        }
976
984
 
977
985
        return table
993
1001
        54, 26, 40, 15, 34, 20, 31, 10, 25, 14, 19, 9, 13, 8, 7, 6,
994
1002
}
995
1003
 
996
 
// trailingZeroBits returns the number of consecutive zero bits on the right
997
 
// side of the given Word.
998
 
// See Knuth, volume 4, section 7.3.1
999
 
func trailingZeroBits(x Word) int {
 
1004
// trailingZeroBits returns the number of consecutive least significant zero
 
1005
// bits of x.
 
1006
func trailingZeroBits(x Word) uint {
1000
1007
        // x & -x leaves only the right-most bit set in the word. Let k be the
1001
1008
        // index of that bit. Since only a single bit is set, the value is two
1002
1009
        // to the power of k. Multiplying by a power of two is equivalent to
1005
1012
        // Therefore, if we have a left shifted version of this constant we can
1006
1013
        // find by how many bits it was shifted by looking at which six bit
1007
1014
        // substring ended up at the top of the word.
 
1015
        // (Knuth, volume 4, section 7.3.1)
1008
1016
        switch _W {
1009
1017
        case 32:
1010
 
                return int(deBruijn32Lookup[((x&-x)*deBruijn32)>>27])
 
1018
                return uint(deBruijn32Lookup[((x&-x)*deBruijn32)>>27])
1011
1019
        case 64:
1012
 
                return int(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58])
 
1020
                return uint(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58])
1013
1021
        default:
1014
 
                panic("Unknown word size")
 
1022
                panic("unknown word size")
1015
1023
        }
 
1024
}
1016
1025
 
1017
 
        return 0
 
1026
// trailingZeroBits returns the number of consecutive least significant zero
 
1027
// bits of x.
 
1028
func (x nat) trailingZeroBits() uint {
 
1029
        if len(x) == 0 {
 
1030
                return 0
 
1031
        }
 
1032
        var i uint
 
1033
        for x[i] == 0 {
 
1034
                i++
 
1035
        }
 
1036
        // x[i] != 0
 
1037
        return i*_W + trailingZeroBits(x[i])
1018
1038
}
1019
1039
 
1020
1040
// z = x << s
1169
1189
        return divWVW(q, 0, x, d)
1170
1190
}
1171
1191
 
1172
 
// powersOfTwoDecompose finds q and k with x = q * 1<<k and q is odd, or q and k are 0.
1173
 
func (x nat) powersOfTwoDecompose() (q nat, k int) {
1174
 
        if len(x) == 0 {
1175
 
                return x, 0
1176
 
        }
1177
 
 
1178
 
        // One of the words must be non-zero by definition,
1179
 
        // so this loop will terminate with i < len(x), and
1180
 
        // i is the number of 0 words.
1181
 
        i := 0
1182
 
        for x[i] == 0 {
1183
 
                i++
1184
 
        }
1185
 
        n := trailingZeroBits(x[i]) // x[i] != 0
1186
 
 
1187
 
        q = make(nat, len(x)-i)
1188
 
        shrVU(q, x[i:], uint(n))
1189
 
 
1190
 
        q = q.norm()
1191
 
        k = i*_W + n
1192
 
        return
1193
 
}
1194
 
 
1195
1192
// random creates a random integer in [0..limit), using the space in z if
1196
1193
// possible. n is the bit length of limit.
1197
1194
func (z nat) random(rand *rand.Rand, limit nat, n int) nat {
1207
1204
        mask := Word((1 << bitLengthOfMSW) - 1)
1208
1205
 
1209
1206
        for {
1210
 
                for i := range z {
1211
 
                        switch _W {
1212
 
                        case 32:
 
1207
                switch _W {
 
1208
                case 32:
 
1209
                        for i := range z {
1213
1210
                                z[i] = Word(rand.Uint32())
1214
 
                        case 64:
 
1211
                        }
 
1212
                case 64:
 
1213
                        for i := range z {
1215
1214
                                z[i] = Word(rand.Uint32()) | Word(rand.Uint32())<<32
1216
1215
                        }
 
1216
                default:
 
1217
                        panic("unknown word size")
1217
1218
                }
1218
 
 
1219
1219
                z[len(limit)-1] &= mask
1220
 
 
1221
1220
                if z.cmp(limit) < 0 {
1222
1221
                        break
1223
1222
                }
1226
1225
        return z.norm()
1227
1226
}
1228
1227
 
1229
 
// If m != nil, expNN calculates x**y mod m. Otherwise it calculates x**y. It
1230
 
// reuses the storage of z if possible.
 
1228
// If m != 0 (i.e., len(m) != 0), expNN sets z to x**y mod m;
 
1229
// otherwise it sets z to x**y. The result is the value of z.
1231
1230
func (z nat) expNN(x, y, m nat) nat {
1232
1231
        if alias(z, x) || alias(z, y) {
1233
 
                // We cannot allow in place modification of x or y.
 
1232
                // We cannot allow in-place modification of x or y.
1234
1233
                z = nil
1235
1234
        }
1236
1235
 
1239
1238
                z[0] = 1
1240
1239
                return z
1241
1240
        }
 
1241
        // y > 0
1242
1242
 
1243
 
        if m != nil {
 
1243
        if len(m) != 0 {
1244
1244
                // We likely end up being as long as the modulus.
1245
1245
                z = z.make(len(m))
1246
1246
        }
1247
1247
        z = z.set(x)
1248
 
        v := y[len(y)-1]
1249
 
        // It's invalid for the most significant word to be zero, therefore we
1250
 
        // will find a one bit.
 
1248
 
 
1249
        // If the base is non-trivial and the exponent is large, we use
 
1250
        // 4-bit, windowed exponentiation. This involves precomputing 14 values
 
1251
        // (x^2...x^15) but then reduces the number of multiply-reduces by a
 
1252
        // third. Even for a 32-bit exponent, this reduces the number of
 
1253
        // operations.
 
1254
        if len(x) > 1 && len(y) > 1 && len(m) > 0 {
 
1255
                return z.expNNWindowed(x, y, m)
 
1256
        }
 
1257
 
 
1258
        v := y[len(y)-1] // v > 0 because y is normalized and y > 0
1251
1259
        shift := leadingZeros(v) + 1
1252
1260
        v <<= shift
1253
1261
        var q nat
1259
1267
        // we also multiply by x, thus adding one to the power.
1260
1268
 
1261
1269
        w := _W - int(shift)
 
1270
        // zz and r are used to avoid allocating in mul and div as
 
1271
        // otherwise the arguments would alias.
 
1272
        var zz, r nat
1262
1273
        for j := 0; j < w; j++ {
1263
 
                z = z.mul(z, z)
 
1274
                zz = zz.mul(z, z)
 
1275
                zz, z = z, zz
1264
1276
 
1265
1277
                if v&mask != 0 {
1266
 
                        z = z.mul(z, x)
 
1278
                        zz = zz.mul(z, x)
 
1279
                        zz, z = z, zz
1267
1280
                }
1268
1281
 
1269
 
                if m != nil {
1270
 
                        q, z = q.div(z, z, m)
 
1282
                if len(m) != 0 {
 
1283
                        zz, r = zz.div(r, z, m)
 
1284
                        zz, r, q, z = q, z, zz, r
1271
1285
                }
1272
1286
 
1273
1287
                v <<= 1
1277
1291
                v = y[i]
1278
1292
 
1279
1293
                for j := 0; j < _W; j++ {
1280
 
                        z = z.mul(z, z)
 
1294
                        zz = zz.mul(z, z)
 
1295
                        zz, z = z, zz
1281
1296
 
1282
1297
                        if v&mask != 0 {
1283
 
                                z = z.mul(z, x)
 
1298
                                zz = zz.mul(z, x)
 
1299
                                zz, z = z, zz
1284
1300
                        }
1285
1301
 
1286
 
                        if m != nil {
1287
 
                                q, z = q.div(z, z, m)
 
1302
                        if len(m) != 0 {
 
1303
                                zz, r = zz.div(r, z, m)
 
1304
                                zz, r, q, z = q, z, zz, r
1288
1305
                        }
1289
1306
 
1290
1307
                        v <<= 1
1294
1311
        return z.norm()
1295
1312
}
1296
1313
 
 
1314
// expNNWindowed calculates x**y mod m using a fixed, 4-bit window.
 
1315
func (z nat) expNNWindowed(x, y, m nat) nat {
 
1316
        // zz and r are used to avoid allocating in mul and div as otherwise
 
1317
        // the arguments would alias.
 
1318
        var zz, r nat
 
1319
 
 
1320
        const n = 4
 
1321
        // powers[i] contains x^i.
 
1322
        var powers [1 << n]nat
 
1323
        powers[0] = natOne
 
1324
        powers[1] = x
 
1325
        for i := 2; i < 1<<n; i += 2 {
 
1326
                p2, p, p1 := &powers[i/2], &powers[i], &powers[i+1]
 
1327
                *p = p.mul(*p2, *p2)
 
1328
                zz, r = zz.div(r, *p, m)
 
1329
                *p, r = r, *p
 
1330
                *p1 = p1.mul(*p, x)
 
1331
                zz, r = zz.div(r, *p1, m)
 
1332
                *p1, r = r, *p1
 
1333
        }
 
1334
 
 
1335
        z = z.setWord(1)
 
1336
 
 
1337
        for i := len(y) - 1; i >= 0; i-- {
 
1338
                yi := y[i]
 
1339
                for j := 0; j < _W; j += n {
 
1340
                        if i != len(y)-1 || j != 0 {
 
1341
                                // Unrolled loop for significant performance
 
1342
                                // gain.  Use go test -bench=".*" in crypto/rsa
 
1343
                                // to check performance before making changes.
 
1344
                                zz = zz.mul(z, z)
 
1345
                                zz, z = z, zz
 
1346
                                zz, r = zz.div(r, z, m)
 
1347
                                z, r = r, z
 
1348
 
 
1349
                                zz = zz.mul(z, z)
 
1350
                                zz, z = z, zz
 
1351
                                zz, r = zz.div(r, z, m)
 
1352
                                z, r = r, z
 
1353
 
 
1354
                                zz = zz.mul(z, z)
 
1355
                                zz, z = z, zz
 
1356
                                zz, r = zz.div(r, z, m)
 
1357
                                z, r = r, z
 
1358
 
 
1359
                                zz = zz.mul(z, z)
 
1360
                                zz, z = z, zz
 
1361
                                zz, r = zz.div(r, z, m)
 
1362
                                z, r = r, z
 
1363
                        }
 
1364
 
 
1365
                        zz = zz.mul(z, powers[yi>>(_W-n)])
 
1366
                        zz, z = z, zz
 
1367
                        zz, r = zz.div(r, z, m)
 
1368
                        z, r = r, z
 
1369
 
 
1370
                        yi <<= n
 
1371
                }
 
1372
        }
 
1373
 
 
1374
        return z.norm()
 
1375
}
 
1376
 
1297
1377
// probablyPrime performs reps Miller-Rabin tests to check whether n is prime.
1298
1378
// If it returns true, n is prime with probability 1 - 1/4^reps.
1299
1379
// If it returns false, n is not prime.
1343
1423
        }
1344
1424
 
1345
1425
        nm1 := nat(nil).sub(n, natOne)
1346
 
        // 1<<k * q = nm1;
1347
 
        q, k := nm1.powersOfTwoDecompose()
 
1426
        // determine q, k such that nm1 = q << k
 
1427
        k := nm1.trailingZeroBits()
 
1428
        q := nat(nil).shr(nm1, k)
1348
1429
 
1349
1430
        nm3 := nat(nil).sub(nm1, natTwo)
1350
1431
        rand := rand.New(rand.NewSource(int64(n[0])))
1360
1441
                if y.cmp(natOne) == 0 || y.cmp(nm1) == 0 {
1361
1442
                        continue
1362
1443
                }
1363
 
                for j := 1; j < k; j++ {
 
1444
                for j := uint(1); j < k; j++ {
1364
1445
                        y = y.mul(y, y)
1365
1446
                        quotient, y = quotient.div(y, y, n)
1366
1447
                        if y.cmp(nm1) == 0 {