authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-03-27 22:28:38+13:00
committergravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-04-11 19:36:35+12:00
logd3e1f323626ba5955aa98628c9c817a19baa35fb
tree37b96965ffb27bdb13df6e8af50efc92540e2a91
parentb3ecdfd7bfb9c6b3bd085b489fe20dd0ca1e12d8

Small fixes for big.Rational and corrections for gcdLehmer

The int div code still causes some edge cases to fail right now.

2 files changed, 91 insertions(+), 76 deletions(-)

std/math/big/int.zig+5-2
...@@ -93,6 +93,7 @@ pub const Int = struct {...@@ -93,6 +93,7 @@ pub const Int = struct {
93 }93 }
9494
95 pub fn clone(other: Int) !Int {95 pub fn clone(other: Int) !Int {
96 other.assertWritable();
96 return Int{97 return Int{
97 .allocator = other.allocator,98 .allocator = other.allocator,
98 .positive = other.positive,99 .positive = other.positive,
...@@ -804,11 +805,13 @@ pub const Int = struct {...@@ -804,11 +805,13 @@ pub const Int = struct {
804 rem.positive = true;805 rem.positive = true;
805 } else {806 } else {
806 // x and y are modified during division807 // x and y are modified during division
807 var x = try a.clone();808 var x = try Int.initCapacity(quo.allocator.?, a.len);
808 defer x.deinit();809 defer x.deinit();
810 try x.copy(a);
809811
810 var y = try b.clone();812 var y = try Int.initCapacity(quo.allocator.?, b.len);
811 defer y.deinit();813 defer y.deinit();
814 try y.copy(b);
812815
813 // x may grow one limb during normalization816 // x may grow one limb during normalization
814 try quo.ensureCapacity(a.len + y.len);817 try quo.ensureCapacity(a.len + y.len);
std/math/big/rational.zig+86-74
...@@ -3,6 +3,7 @@ const builtin = @import("builtin");...@@ -3,6 +3,7 @@ const builtin = @import("builtin");
3const debug = std.debug;3const debug = std.debug;
4const math = std.math;4const math = std.math;
5const mem = std.mem;5const mem = std.mem;
6const testing = std.testing;
6const Allocator = mem.Allocator;7const Allocator = mem.Allocator;
7const ArrayList = std.ArrayList;8const ArrayList = std.ArrayList;
89
...@@ -425,6 +426,7 @@ var al = debug.global_allocator;...@@ -425,6 +426,7 @@ var al = debug.global_allocator;
425const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);426const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);
426427
427fn gcd(rma: *Int, x: Int, y: Int) !void {428fn gcd(rma: *Int, x: Int, y: Int) !void {
429 rma.assertWritable();
428 var r = rma;430 var r = rma;
429 var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;431 var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;
430432
...@@ -463,19 +465,23 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {...@@ -463,19 +465,23 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
463//465//
464// r = gcd(x, y) where x, y > 0466// r = gcd(x, y) where x, y > 0
465fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {467fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
466 debug.assert(xa.positive and ya.positive);
467 debug.assert(xa.cmp(ya) >= 0);
468
469 var x = try xa.clone();468 var x = try xa.clone();
469 x.abs();
470 defer x.deinit();470 defer x.deinit();
471471
472 var y = try ya.clone();472 var y = try ya.clone();
473 y.abs();
473 defer y.deinit();474 defer y.deinit();
474475
476 if (x.cmp(y) < 0) {
477 x.swap(&y);
478 }
479
475 var T = try Int.init(r.allocator.?);480 var T = try Int.init(r.allocator.?);
476 defer T.deinit();481 defer T.deinit();
477482
478 while (y.len > 1) {483 while (y.len > 1) {
484 debug.assert(x.positive and y.positive);
479 debug.assert(x.len >= y.len);485 debug.assert(x.len >= y.len);
480486
481 // chop the leading zeros of the limbs and normalize487 // chop the leading zeros of the limbs and normalize
...@@ -540,7 +546,13 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {...@@ -540,7 +546,13 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
540 try r.add(x, r.*);546 try r.add(x, r.*);
541547
542 x.swap(&T);548 x.swap(&T);
549 x.abs();
543 y.swap(r);550 y.swap(r);
551 y.abs();
552
553 if (x.cmp(y) < 0) {
554 x.swap(&y);
555 }
544 }556 }
545 }557 }
546558
...@@ -563,7 +575,7 @@ test "big.rational gcd non-one small" {...@@ -563,7 +575,7 @@ test "big.rational gcd non-one small" {
563575
564 try gcd(&r, a, b);576 try gcd(&r, a, b);
565577
566 debug.assert((try r.to(u32)) == 1);578 testing.expect((try r.to(u32)) == 1);
567}579}
568580
569test "big.rational gcd non-one small" {581test "big.rational gcd non-one small" {
...@@ -573,7 +585,7 @@ test "big.rational gcd non-one small" {...@@ -573,7 +585,7 @@ test "big.rational gcd non-one small" {
573585
574 try gcd(&r, a, b);586 try gcd(&r, a, b);
575587
576 debug.assert((try r.to(u32)) == 38);588 testing.expect((try r.to(u32)) == 38);
577}589}
578590
579test "big.rational gcd non-one large" {591test "big.rational gcd non-one large" {
...@@ -583,7 +595,7 @@ test "big.rational gcd non-one large" {...@@ -583,7 +595,7 @@ test "big.rational gcd non-one large" {
583595
584 try gcd(&r, a, b);596 try gcd(&r, a, b);
585597
586 debug.assert((try r.to(u32)) == 4369);598 testing.expect((try r.to(u32)) == 4369);
587}599}
588600
589test "big.rational gcd large multi-limb result" {601test "big.rational gcd large multi-limb result" {
...@@ -593,11 +605,11 @@ test "big.rational gcd large multi-limb result" {...@@ -593,11 +605,11 @@ test "big.rational gcd large multi-limb result" {
593605
594 try gcd(&r, a, b);606 try gcd(&r, a, b);
595607
596 debug.assert((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);608 testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
597}609}
598610
599fn extractLowBits(a: Int, comptime T: type) T {611fn extractLowBits(a: Int, comptime T: type) T {
600 debug.assert(@typeId(T) == builtin.TypeId.Int);612 testing.expect(@typeId(T) == builtin.TypeId.Int);
601613
602 if (T.bit_count <= Limb.bit_count) {614 if (T.bit_count <= Limb.bit_count) {
603 return @truncate(T, a.limbs[0]);615 return @truncate(T, a.limbs[0]);
...@@ -619,71 +631,71 @@ test "big.rational extractLowBits" {...@@ -619,71 +631,71 @@ test "big.rational extractLowBits" {
619 var a = try Int.initSet(al, 0x11112222333344441234567887654321);631 var a = try Int.initSet(al, 0x11112222333344441234567887654321);
620632
621 const a1 = extractLowBits(a, u8);633 const a1 = extractLowBits(a, u8);
622 debug.assert(a1 == 0x21);634 testing.expect(a1 == 0x21);
623635
624 const a2 = extractLowBits(a, u16);636 const a2 = extractLowBits(a, u16);
625 debug.assert(a2 == 0x4321);637 testing.expect(a2 == 0x4321);
626638
627 const a3 = extractLowBits(a, u32);639 const a3 = extractLowBits(a, u32);
628 debug.assert(a3 == 0x87654321);640 testing.expect(a3 == 0x87654321);
629641
630 const a4 = extractLowBits(a, u64);642 const a4 = extractLowBits(a, u64);
631 debug.assert(a4 == 0x1234567887654321);643 testing.expect(a4 == 0x1234567887654321);
632644
633 const a5 = extractLowBits(a, u128);645 const a5 = extractLowBits(a, u128);
634 debug.assert(a5 == 0x11112222333344441234567887654321);646 testing.expect(a5 == 0x11112222333344441234567887654321);
635}647}
636648
637test "big.rational set" {649test "big.rational set" {
638 var a = try Rational.init(al);650 var a = try Rational.init(al);
639651
640 try a.setInt(5);652 try a.setInt(5);
641 debug.assert((try a.p.to(u32)) == 5);653 testing.expect((try a.p.to(u32)) == 5);
642 debug.assert((try a.q.to(u32)) == 1);654 testing.expect((try a.q.to(u32)) == 1);
643655
644 try a.setRatio(7, 3);656 try a.setRatio(7, 3);
645 debug.assert((try a.p.to(u32)) == 7);657 testing.expect((try a.p.to(u32)) == 7);
646 debug.assert((try a.q.to(u32)) == 3);658 testing.expect((try a.q.to(u32)) == 3);
647659
648 try a.setRatio(9, 3);660 try a.setRatio(9, 3);
649 debug.assert((try a.p.to(i32)) == 3);661 testing.expect((try a.p.to(i32)) == 3);
650 debug.assert((try a.q.to(i32)) == 1);662 testing.expect((try a.q.to(i32)) == 1);
651663
652 try a.setRatio(-9, 3);664 try a.setRatio(-9, 3);
653 debug.assert((try a.p.to(i32)) == -3);665 testing.expect((try a.p.to(i32)) == -3);
654 debug.assert((try a.q.to(i32)) == 1);666 testing.expect((try a.q.to(i32)) == 1);
655667
656 try a.setRatio(9, -3);668 try a.setRatio(9, -3);
657 debug.assert((try a.p.to(i32)) == -3);669 testing.expect((try a.p.to(i32)) == -3);
658 debug.assert((try a.q.to(i32)) == 1);670 testing.expect((try a.q.to(i32)) == 1);
659671
660 try a.setRatio(-9, -3);672 try a.setRatio(-9, -3);
661 debug.assert((try a.p.to(i32)) == 3);673 testing.expect((try a.p.to(i32)) == 3);
662 debug.assert((try a.q.to(i32)) == 1);674 testing.expect((try a.q.to(i32)) == 1);
663}675}
664676
665test "big.rational setFloat" {677test "big.rational setFloat" {
666 var a = try Rational.init(al);678 var a = try Rational.init(al);
667679
668 try a.setFloat(f64, 2.5);680 try a.setFloat(f64, 2.5);
669 debug.assert((try a.p.to(i32)) == 5);681 testing.expect((try a.p.to(i32)) == 5);
670 debug.assert((try a.q.to(i32)) == 2);682 testing.expect((try a.q.to(i32)) == 2);
671683
672 try a.setFloat(f32, -2.5);684 try a.setFloat(f32, -2.5);
673 debug.assert((try a.p.to(i32)) == -5);685 testing.expect((try a.p.to(i32)) == -5);
674 debug.assert((try a.q.to(i32)) == 2);686 testing.expect((try a.q.to(i32)) == 2);
675687
676 try a.setFloat(f32, 3.141593);688 try a.setFloat(f32, 3.141593);
677689
678 // = 3.14159297943115234375690 // = 3.14159297943115234375
679 debug.assert((try a.p.to(u32)) == 3294199);691 testing.expect((try a.p.to(u32)) == 3294199);
680 debug.assert((try a.q.to(u32)) == 1048576);692 testing.expect((try a.q.to(u32)) == 1048576);
681693
682 try a.setFloat(f64, 72.141593120712409172417410926841290461290467124);694 try a.setFloat(f64, 72.141593120712409172417410926841290461290467124);
683695
684 // = 72.1415931207124145885245525278151035308837890625696 // = 72.1415931207124145885245525278151035308837890625
685 debug.assert((try a.p.to(u128)) == 5076513310880537);697 testing.expect((try a.p.to(u128)) == 5076513310880537);
686 debug.assert((try a.q.to(u128)) == 70368744177664);698 testing.expect((try a.q.to(u128)) == 70368744177664);
687}699}
688700
689test "big.rational setFloatString" {701test "big.rational setFloatString" {
...@@ -692,8 +704,8 @@ test "big.rational setFloatString" {...@@ -692,8 +704,8 @@ test "big.rational setFloatString" {
692 try a.setFloatString("72.14159312071241458852455252781510353");704 try a.setFloatString("72.14159312071241458852455252781510353");
693705
694 // = 72.1415931207124145885245525278151035308837890625706 // = 72.1415931207124145885245525278151035308837890625
695 debug.assert((try a.p.to(u128)) == 7214159312071241458852455252781510353);707 testing.expect((try a.p.to(u128)) == 7214159312071241458852455252781510353);
696 debug.assert((try a.q.to(u128)) == 100000000000000000000000000000000000);708 testing.expect((try a.q.to(u128)) == 100000000000000000000000000000000000);
697}709}
698710
699test "big.rational toFloat" {711test "big.rational toFloat" {
...@@ -701,11 +713,11 @@ test "big.rational toFloat" {...@@ -701,11 +713,11 @@ test "big.rational toFloat" {
701713
702 // = 3.14159297943115234375714 // = 3.14159297943115234375
703 try a.setRatio(3294199, 1048576);715 try a.setRatio(3294199, 1048576);
704 debug.assert((try a.toFloat(f64)) == 3.14159297943115234375);716 testing.expect((try a.toFloat(f64)) == 3.14159297943115234375);
705717
706 // = 72.1415931207124145885245525278151035308837890625718 // = 72.1415931207124145885245525278151035308837890625
707 try a.setRatio(5076513310880537, 70368744177664);719 try a.setRatio(5076513310880537, 70368744177664);
708 debug.assert((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124);720 testing.expect((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124);
709}721}
710722
711test "big.rational set/to Float round-trip" {723test "big.rational set/to Float round-trip" {
...@@ -720,7 +732,7 @@ test "big.rational set/to Float round-trip" {...@@ -720,7 +732,7 @@ test "big.rational set/to Float round-trip" {
720 while (i < 512) : (i += 1) {732 while (i < 512) : (i += 1) {
721 const r = prng.random.float(f64);733 const r = prng.random.float(f64);
722 try a.setFloat(f64, r);734 try a.setFloat(f64, r);
723 debug.assert((try a.toFloat(f64)) == r);735 testing.expect((try a.toFloat(f64)) == r);
724 }736 }
725}737}
726738
...@@ -730,54 +742,54 @@ test "big.rational copy" {...@@ -730,54 +742,54 @@ test "big.rational copy" {
730 const b = try Int.initSet(al, 5);742 const b = try Int.initSet(al, 5);
731743
732 try a.copyInt(b);744 try a.copyInt(b);
733 debug.assert((try a.p.to(u32)) == 5);745 testing.expect((try a.p.to(u32)) == 5);
734 debug.assert((try a.q.to(u32)) == 1);746 testing.expect((try a.q.to(u32)) == 1);
735747
736 const c = try Int.initSet(al, 7);748 const c = try Int.initSet(al, 7);
737 const d = try Int.initSet(al, 3);749 const d = try Int.initSet(al, 3);
738750
739 try a.copyRatio(c, d);751 try a.copyRatio(c, d);
740 debug.assert((try a.p.to(u32)) == 7);752 testing.expect((try a.p.to(u32)) == 7);
741 debug.assert((try a.q.to(u32)) == 3);753 testing.expect((try a.q.to(u32)) == 3);
742754
743 const e = try Int.initSet(al, 9);755 const e = try Int.initSet(al, 9);
744 const f = try Int.initSet(al, 3);756 const f = try Int.initSet(al, 3);
745757
746 try a.copyRatio(e, f);758 try a.copyRatio(e, f);
747 debug.assert((try a.p.to(u32)) == 3);759 testing.expect((try a.p.to(u32)) == 3);
748 debug.assert((try a.q.to(u32)) == 1);760 testing.expect((try a.q.to(u32)) == 1);
749}761}
750762
751test "big.rational negate" {763test "big.rational negate" {
752 var a = try Rational.init(al);764 var a = try Rational.init(al);
753765
754 try a.setInt(-50);766 try a.setInt(-50);
755 debug.assert((try a.p.to(i32)) == -50);767 testing.expect((try a.p.to(i32)) == -50);
756 debug.assert((try a.q.to(i32)) == 1);768 testing.expect((try a.q.to(i32)) == 1);
757769
758 a.negate();770 a.negate();
759 debug.assert((try a.p.to(i32)) == 50);771 testing.expect((try a.p.to(i32)) == 50);
760 debug.assert((try a.q.to(i32)) == 1);772 testing.expect((try a.q.to(i32)) == 1);
761773
762 a.negate();774 a.negate();
763 debug.assert((try a.p.to(i32)) == -50);775 testing.expect((try a.p.to(i32)) == -50);
764 debug.assert((try a.q.to(i32)) == 1);776 testing.expect((try a.q.to(i32)) == 1);
765}777}
766778
767test "big.rational abs" {779test "big.rational abs" {
768 var a = try Rational.init(al);780 var a = try Rational.init(al);
769781
770 try a.setInt(-50);782 try a.setInt(-50);
771 debug.assert((try a.p.to(i32)) == -50);783 testing.expect((try a.p.to(i32)) == -50);
772 debug.assert((try a.q.to(i32)) == 1);784 testing.expect((try a.q.to(i32)) == 1);
773785
774 a.abs();786 a.abs();
775 debug.assert((try a.p.to(i32)) == 50);787 testing.expect((try a.p.to(i32)) == 50);
776 debug.assert((try a.q.to(i32)) == 1);788 testing.expect((try a.q.to(i32)) == 1);
777789
778 a.abs();790 a.abs();
779 debug.assert((try a.p.to(i32)) == 50);791 testing.expect((try a.p.to(i32)) == 50);
780 debug.assert((try a.q.to(i32)) == 1);792 testing.expect((try a.q.to(i32)) == 1);
781}793}
782794
783test "big.rational swap" {795test "big.rational swap" {
...@@ -787,19 +799,19 @@ test "big.rational swap" {...@@ -787,19 +799,19 @@ test "big.rational swap" {
787 try a.setRatio(50, 23);799 try a.setRatio(50, 23);
788 try b.setRatio(17, 3);800 try b.setRatio(17, 3);
789801
790 debug.assert((try a.p.to(u32)) == 50);802 testing.expect((try a.p.to(u32)) == 50);
791 debug.assert((try a.q.to(u32)) == 23);803 testing.expect((try a.q.to(u32)) == 23);
792804
793 debug.assert((try b.p.to(u32)) == 17);805 testing.expect((try b.p.to(u32)) == 17);
794 debug.assert((try b.q.to(u32)) == 3);806 testing.expect((try b.q.to(u32)) == 3);
795807
796 a.swap(&b);808 a.swap(&b);
797809
798 debug.assert((try a.p.to(u32)) == 17);810 testing.expect((try a.p.to(u32)) == 17);
799 debug.assert((try a.q.to(u32)) == 3);811 testing.expect((try a.q.to(u32)) == 3);
800812
801 debug.assert((try b.p.to(u32)) == 50);813 testing.expect((try b.p.to(u32)) == 50);
802 debug.assert((try b.q.to(u32)) == 23);814 testing.expect((try b.q.to(u32)) == 23);
803}815}
804816
805test "big.rational cmp" {817test "big.rational cmp" {
...@@ -808,11 +820,11 @@ test "big.rational cmp" {...@@ -808,11 +820,11 @@ test "big.rational cmp" {
808820
809 try a.setRatio(500, 231);821 try a.setRatio(500, 231);
810 try b.setRatio(18903, 8584);822 try b.setRatio(18903, 8584);
811 debug.assert((try a.cmp(b)) < 0);823 testing.expect((try a.cmp(b)) < 0);
812824
813 try a.setRatio(890, 10);825 try a.setRatio(890, 10);
814 try b.setRatio(89, 1);826 try b.setRatio(89, 1);
815 debug.assert((try a.cmp(b)) == 0);827 testing.expect((try a.cmp(b)) == 0);
816}828}
817829
818test "big.rational add single-limb" {830test "big.rational add single-limb" {
...@@ -821,11 +833,11 @@ test "big.rational add single-limb" {...@@ -821,11 +833,11 @@ test "big.rational add single-limb" {
821833
822 try a.setRatio(500, 231);834 try a.setRatio(500, 231);
823 try b.setRatio(18903, 8584);835 try b.setRatio(18903, 8584);
824 debug.assert((try a.cmp(b)) < 0);836 testing.expect((try a.cmp(b)) < 0);
825837
826 try a.setRatio(890, 10);838 try a.setRatio(890, 10);
827 try b.setRatio(89, 1);839 try b.setRatio(89, 1);
828 debug.assert((try a.cmp(b)) == 0);840 testing.expect((try a.cmp(b)) == 0);
829}841}
830842
831test "big.rational add" {843test "big.rational add" {
...@@ -838,7 +850,7 @@ test "big.rational add" {...@@ -838,7 +850,7 @@ test "big.rational add" {
838 try a.add(a, b);850 try a.add(a, b);
839851
840 try r.setRatio(984786924199, 290395044174);852 try r.setRatio(984786924199, 290395044174);
841 debug.assert((try a.cmp(r)) == 0);853 testing.expect((try a.cmp(r)) == 0);
842}854}
843855
844test "big.rational sub" {856test "big.rational sub" {
...@@ -851,7 +863,7 @@ test "big.rational sub" {...@@ -851,7 +863,7 @@ test "big.rational sub" {
851 try a.sub(a, b);863 try a.sub(a, b);
852864
853 try r.setRatio(979040510045, 290395044174);865 try r.setRatio(979040510045, 290395044174);
854 debug.assert((try a.cmp(r)) == 0);866 testing.expect((try a.cmp(r)) == 0);
855}867}
856868
857test "big.rational mul" {869test "big.rational mul" {
...@@ -864,7 +876,7 @@ test "big.rational mul" {...@@ -864,7 +876,7 @@ test "big.rational mul" {
864 try a.mul(a, b);876 try a.mul(a, b);
865877
866 try r.setRatio(571481443, 17082061422);878 try r.setRatio(571481443, 17082061422);
867 debug.assert((try a.cmp(r)) == 0);879 testing.expect((try a.cmp(r)) == 0);
868}880}
869881
870test "big.rational div" {882test "big.rational div" {
...@@ -877,7 +889,7 @@ test "big.rational div" {...@@ -877,7 +889,7 @@ test "big.rational div" {
877 try a.div(a, b);889 try a.div(a, b);
878890
879 try r.setRatio(75531824394, 221015929);891 try r.setRatio(75531824394, 221015929);
880 debug.assert((try a.cmp(r)) == 0);892 testing.expect((try a.cmp(r)) == 0);
881}893}
882894
883test "big.rational div" {895test "big.rational div" {
...@@ -888,11 +900,11 @@ test "big.rational div" {...@@ -888,11 +900,11 @@ test "big.rational div" {
888 a.invert();900 a.invert();
889901
890 try r.setRatio(23341, 78923);902 try r.setRatio(23341, 78923);
891 debug.assert((try a.cmp(r)) == 0);903 testing.expect((try a.cmp(r)) == 0);
892904
893 try a.setRatio(-78923, 23341);905 try a.setRatio(-78923, 23341);
894 a.invert();906 a.invert();
895907
896 try r.setRatio(-23341, 78923);908 try r.setRatio(-23341, 78923);
897 debug.assert((try a.cmp(r)) == 0);909 testing.expect((try a.cmp(r)) == 0);
898}910}