authorgravatar for joachim.schmidt557@outlook.comJoachim Schmidt <joachim.schmidt557@outlook.com> 2020-04-06 23:22:03+02:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2020-04-07 13:43:15-04:00
log1ee59c5c31efca394d7e6d9da3f64c289a996b99
tree87c81299ed1adfc22f3e34901bbfcd3062dcf752
parenta20f3e3f02b24219f43c2da4568b794051e7d348

move big.rational.gcd to big.int.gcd

2 files changed, 188 insertions(+), 188 deletions(-)

lib/std/math/big/int.zig+187
...@@ -10,6 +10,7 @@ const minInt = std.math.minInt;...@@ -10,6 +10,7 @@ const minInt = std.math.minInt;
1010
11pub const Limb = usize;11pub const Limb = usize;
12pub const DoubleLimb = std.meta.IntType(false, 2 * Limb.bit_count);12pub const DoubleLimb = std.meta.IntType(false, 2 * Limb.bit_count);
13pub const SignedDoubleLimb = std.meta.IntType(true, DoubleLimb.bit_count);
13pub const Log2Limb = math.Log2Int(Limb);14pub const Log2Limb = math.Log2Int(Limb);
1415
15comptime {16comptime {
...@@ -1359,8 +1360,129 @@ pub const Int = struct {...@@ -1359,8 +1360,129 @@ pub const Int = struct {
1359 r[i] = a[i];1360 r[i] = a[i];
1360 }1361 }
1361 }1362 }
1363
1364 pub fn gcd(rma: *Int, x: Int, y: Int) !void {
1365 rma.assertWritable();
1366 var r = rma;
1367 var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;
1368
1369 var sr: Int = undefined;
1370 if (aliased) {
1371 sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len()));
1372 r = &sr;
1373 aliased = true;
1374 }
1375 defer if (aliased) {
1376 rma.swap(r);
1377 r.deinit();
1378 };
1379
1380 try gcdLehmer(r, x, y);
1381 }
1382
1383 fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
1384 var x = try xa.clone();
1385 x.abs();
1386 defer x.deinit();
1387
1388 var y = try ya.clone();
1389 y.abs();
1390 defer y.deinit();
1391
1392 if (x.cmp(y) == .lt) {
1393 x.swap(&y);
1394 }
1395
1396 var T = try Int.init(r.allocator.?);
1397 defer T.deinit();
1398
1399 while (y.len() > 1) {
1400 debug.assert(x.isPositive() and y.isPositive());
1401 debug.assert(x.len() >= y.len());
1402
1403 var xh: SignedDoubleLimb = x.limbs[x.len() - 1];
1404 var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1];
1405
1406 var A: SignedDoubleLimb = 1;
1407 var B: SignedDoubleLimb = 0;
1408 var C: SignedDoubleLimb = 0;
1409 var D: SignedDoubleLimb = 1;
1410
1411 while (yh + C != 0 and yh + D != 0) {
1412 const q = @divFloor(xh + A, yh + C);
1413 const qp = @divFloor(xh + B, yh + D);
1414 if (q != qp) {
1415 break;
1416 }
1417
1418 var t = A - q * C;
1419 A = C;
1420 C = t;
1421 t = B - q * D;
1422 B = D;
1423 D = t;
1424
1425 t = xh - q * yh;
1426 xh = yh;
1427 yh = t;
1428 }
1429
1430 if (B == 0) {
1431 // T = x % y, r is unused
1432 try Int.divTrunc(r, &T, x, y);
1433 debug.assert(T.isPositive());
1434
1435 x.swap(&y);
1436 y.swap(&T);
1437 } else {
1438 var storage: [8]Limb = undefined;
1439 const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]);
1440 const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]);
1441 const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]);
1442 const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]);
1443
1444 // T = Ax + By
1445 try r.mul(x, Ap);
1446 try T.mul(y, Bp);
1447 try T.add(r.*, T);
1448
1449 // u = Cx + Dy, r as u
1450 try x.mul(x, Cp);
1451 try r.mul(y, Dp);
1452 try r.add(x, r.*);
1453
1454 x.swap(&T);
1455 y.swap(r);
1456 }
1457 }
1458
1459 // euclidean algorithm
1460 debug.assert(x.cmp(y) != .lt);
1461
1462 while (!y.eqZero()) {
1463 try Int.divTrunc(&T, r, x, y);
1464 x.swap(&y);
1465 y.swap(r);
1466 }
1467
1468 r.swap(&x);
1469 }
1470
1362};1471};
13631472
1473// Storage must live for the lifetime of the returned value
1474fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
1475 std.debug.assert(storage.len >= 2);
1476
1477 var A_is_positive = A >= 0;
1478 const Au = @intCast(DoubleLimb, if (A < 0) -A else A);
1479 storage[0] = @truncate(Limb, Au);
1480 storage[1] = @truncate(Limb, Au >> Limb.bit_count);
1481 var Ap = Int.initFixed(storage[0..2]);
1482 Ap.setSign(A_is_positive);
1483 return Ap;
1484}
1485
1364// NOTE: All the following tests assume the max machine-word will be 64-bit.1486// NOTE: All the following tests assume the max machine-word will be 64-bit.
1365//1487//
1366// They will still run on larger than this and should pass, but the multi-limb code-paths1488// They will still run on larger than this and should pass, but the multi-limb code-paths
...@@ -2738,3 +2860,68 @@ test "big.int var args" {...@@ -2738,3 +2860,68 @@ test "big.int var args" {
2738 defer d.deinit();2860 defer d.deinit();
2739 testing.expect(a.cmp(d) != .gt);2861 testing.expect(a.cmp(d) != .gt);
2740}2862}
2863
2864test "big.int gcd non-one small" {
2865 var a = try Int.initSet(testing.allocator, 17);
2866 defer a.deinit();
2867 var b = try Int.initSet(testing.allocator, 97);
2868 defer b.deinit();
2869 var r = try Int.init(testing.allocator);
2870 defer r.deinit();
2871
2872 try r.gcd(a, b);
2873
2874 testing.expect((try r.to(u32)) == 1);
2875}
2876
2877test "big.int gcd non-one small" {
2878 var a = try Int.initSet(testing.allocator, 4864);
2879 defer a.deinit();
2880 var b = try Int.initSet(testing.allocator, 3458);
2881 defer b.deinit();
2882 var r = try Int.init(testing.allocator);
2883 defer r.deinit();
2884
2885 try r.gcd(a, b);
2886
2887 testing.expect((try r.to(u32)) == 38);
2888}
2889
2890test "big.int gcd non-one large" {
2891 var a = try Int.initSet(testing.allocator, 0xffffffffffffffff);
2892 defer a.deinit();
2893 var b = try Int.initSet(testing.allocator, 0xffffffffffffffff7777);
2894 defer b.deinit();
2895 var r = try Int.init(testing.allocator);
2896 defer r.deinit();
2897
2898 try r.gcd(a, b);
2899
2900 testing.expect((try r.to(u32)) == 4369);
2901}
2902
2903test "big.int gcd large multi-limb result" {
2904 var a = try Int.initSet(testing.allocator, 0x12345678123456781234567812345678123456781234567812345678);
2905 defer a.deinit();
2906 var b = try Int.initSet(testing.allocator, 0x12345671234567123456712345671234567123456712345671234567);
2907 defer b.deinit();
2908 var r = try Int.init(testing.allocator);
2909 defer r.deinit();
2910
2911 try r.gcd(a, b);
2912
2913 testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
2914}
2915
2916test "big.int gcd one large" {
2917 var a = try Int.initSet(testing.allocator, 1897056385327307);
2918 defer a.deinit();
2919 var b = try Int.initSet(testing.allocator, 2251799813685248);
2920 defer b.deinit();
2921 var r = try Int.init(testing.allocator);
2922 defer r.deinit();
2923
2924 try r.gcd(a, b);
2925
2926 testing.expect((try r.to(u64)) == 1);
2927}
lib/std/math/big/rational.zig+1-188
...@@ -447,7 +447,7 @@ pub const Rational = struct {...@@ -447,7 +447,7 @@ pub const Rational = struct {
447447
448 const sign = r.p.isPositive();448 const sign = r.p.isPositive();
449 r.p.abs();449 r.p.abs();
450 try gcd(&a, r.p, r.q);450 try a.gcd(r.p, r.q);
451 r.p.setSign(sign);451 r.p.setSign(sign);
452452
453 const one = Int.initFixed(([_]Limb{1})[0..]);453 const one = Int.initFixed(([_]Limb{1})[0..]);
...@@ -463,193 +463,6 @@ pub const Rational = struct {...@@ -463,193 +463,6 @@ pub const Rational = struct {
463 }463 }
464};464};
465465
466const SignedDoubleLimb = std.meta.IntType(true, DoubleLimb.bit_count);
467
468fn gcd(rma: *Int, x: Int, y: Int) !void {
469 rma.assertWritable();
470 var r = rma;
471 var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;
472
473 var sr: Int = undefined;
474 if (aliased) {
475 sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len()));
476 r = &sr;
477 aliased = true;
478 }
479 defer if (aliased) {
480 rma.swap(r);
481 r.deinit();
482 };
483
484 try gcdLehmer(r, x, y);
485}
486
487// Storage must live for the lifetime of the returned value
488fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
489 std.debug.assert(storage.len >= 2);
490
491 var A_is_positive = A >= 0;
492 const Au = @intCast(DoubleLimb, if (A < 0) -A else A);
493 storage[0] = @truncate(Limb, Au);
494 storage[1] = @truncate(Limb, Au >> Limb.bit_count);
495 var Ap = Int.initFixed(storage[0..2]);
496 Ap.setSign(A_is_positive);
497 return Ap;
498}
499
500fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
501 var x = try xa.clone();
502 x.abs();
503 defer x.deinit();
504
505 var y = try ya.clone();
506 y.abs();
507 defer y.deinit();
508
509 if (x.cmp(y) == .lt) {
510 x.swap(&y);
511 }
512
513 var T = try Int.init(r.allocator.?);
514 defer T.deinit();
515
516 while (y.len() > 1) {
517 debug.assert(x.isPositive() and y.isPositive());
518 debug.assert(x.len() >= y.len());
519
520 var xh: SignedDoubleLimb = x.limbs[x.len() - 1];
521 var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1];
522
523 var A: SignedDoubleLimb = 1;
524 var B: SignedDoubleLimb = 0;
525 var C: SignedDoubleLimb = 0;
526 var D: SignedDoubleLimb = 1;
527
528 while (yh + C != 0 and yh + D != 0) {
529 const q = @divFloor(xh + A, yh + C);
530 const qp = @divFloor(xh + B, yh + D);
531 if (q != qp) {
532 break;
533 }
534
535 var t = A - q * C;
536 A = C;
537 C = t;
538 t = B - q * D;
539 B = D;
540 D = t;
541
542 t = xh - q * yh;
543 xh = yh;
544 yh = t;
545 }
546
547 if (B == 0) {
548 // T = x % y, r is unused
549 try Int.divTrunc(r, &T, x, y);
550 debug.assert(T.isPositive());
551
552 x.swap(&y);
553 y.swap(&T);
554 } else {
555 var storage: [8]Limb = undefined;
556 const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]);
557 const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]);
558 const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]);
559 const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]);
560
561 // T = Ax + By
562 try r.mul(x, Ap);
563 try T.mul(y, Bp);
564 try T.add(r.*, T);
565
566 // u = Cx + Dy, r as u
567 try x.mul(x, Cp);
568 try r.mul(y, Dp);
569 try r.add(x, r.*);
570
571 x.swap(&T);
572 y.swap(r);
573 }
574 }
575
576 // euclidean algorithm
577 debug.assert(x.cmp(y) != .lt);
578
579 while (!y.eqZero()) {
580 try Int.divTrunc(&T, r, x, y);
581 x.swap(&y);
582 y.swap(r);
583 }
584
585 r.swap(&x);
586}
587
588test "big.rational gcd non-one small" {
589 var a = try Int.initSet(testing.allocator, 17);
590 defer a.deinit();
591 var b = try Int.initSet(testing.allocator, 97);
592 defer b.deinit();
593 var r = try Int.init(testing.allocator);
594 defer r.deinit();
595
596 try gcd(&r, a, b);
597
598 testing.expect((try r.to(u32)) == 1);
599}
600
601test "big.rational gcd non-one small" {
602 var a = try Int.initSet(testing.allocator, 4864);
603 defer a.deinit();
604 var b = try Int.initSet(testing.allocator, 3458);
605 defer b.deinit();
606 var r = try Int.init(testing.allocator);
607 defer r.deinit();
608
609 try gcd(&r, a, b);
610
611 testing.expect((try r.to(u32)) == 38);
612}
613
614test "big.rational gcd non-one large" {
615 var a = try Int.initSet(testing.allocator, 0xffffffffffffffff);
616 defer a.deinit();
617 var b = try Int.initSet(testing.allocator, 0xffffffffffffffff7777);
618 defer b.deinit();
619 var r = try Int.init(testing.allocator);
620 defer r.deinit();
621
622 try gcd(&r, a, b);
623
624 testing.expect((try r.to(u32)) == 4369);
625}
626
627test "big.rational gcd large multi-limb result" {
628 var a = try Int.initSet(testing.allocator, 0x12345678123456781234567812345678123456781234567812345678);
629 defer a.deinit();
630 var b = try Int.initSet(testing.allocator, 0x12345671234567123456712345671234567123456712345671234567);
631 defer b.deinit();
632 var r = try Int.init(testing.allocator);
633 defer r.deinit();
634
635 try gcd(&r, a, b);
636
637 testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
638}
639
640test "big.rational gcd one large" {
641 var a = try Int.initSet(testing.allocator, 1897056385327307);
642 defer a.deinit();
643 var b = try Int.initSet(testing.allocator, 2251799813685248);
644 defer b.deinit();
645 var r = try Int.init(testing.allocator);
646 defer r.deinit();
647
648 try gcd(&r, a, b);
649
650 testing.expect((try r.to(u64)) == 1);
651}
652
653fn extractLowBits(a: Int, comptime T: type) T {466fn extractLowBits(a: Int, comptime T: type) T {
654 testing.expect(@typeInfo(T) == .Int);467 testing.expect(@typeInfo(T) == .Int);
655468