* Bug 15058 fixed: gcd and lcm result could be puzzingly <0
[scilab.git] / scilab / modules / polynomials / tests / unit_tests / gcd.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2017 - Samuel GOUGEON
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 //
8 // <-- CLI SHELL MODE -->
9 // <-- NO CHECK REF -->
10 //
11 // unit tests for gcd()
12 // ====================
13
14 // Checking arguments
15 // ------------------
16 [g, f] = gcd([]);
17 assert_checkequal(g, []);
18 assert_checkequal(f, 1);
19 d = 9*49*13*23*29*37*47;    // 6.650D+09   2^63 ~ 9.223D+18
20 v = [d*2*25*31  d*29*37  d*31];
21 [g, f] = gcd(v);
22 assert_checkequal(size(g), [1 1]);
23 assert_checkequal(size(f), [3 3]);
24
25 // Checking errors
26 // ---------------
27 assert_checkfalse(execstr("gcd()", "errcatch")==0);
28 assert_checkfalse(execstr("gcd(%z,%s)", "errcatch")==0);
29 assert_checkfalse(execstr("[a,b,c]=gcd(v)", "errcatch")==0);
30 assert_checkfalse(execstr("gcd(sparse([25 21 110]))", "errcatch")==0);
31 assert_checkfalse(execstr("gcd([25 21 110]+0.1)", "errcatch")==0);
32
33 // Checking results
34 // ================
35 // With encoded integers
36 // ---------------------
37 //--> pr = primes(50)
38 // pr  =
39 //   2.  3.  5.  7.  11.  13.  17.  19.  23.  29.  31.  37.  41.  43.  47.
40 // int8
41 Gcd = 10;     // 2^7 = 128
42 v = [10*3 ; 10*11];
43 // int16
44 d =  3*7*13;  // 273.  2^15 = 32768
45 Gcd = [Gcd  d];
46 v = [v  [d*3*17 ; d*5*23]];
47 // int32
48 d =  3*7*13*29*37;  // 292929.  2^31 ~ 2.147D+09
49 Gcd = [Gcd  d];
50 v = [v  [d*2*19*27 ; d*5*11*23]];
51 // int64
52 d = 9*49*13*23*29*37*47;  // 6.650D+09   2^63 ~ 9.223D+18
53 Gcd = [Gcd  d];
54 v = [v  [d*2*25*31 ; d*7*29]];  // 1.031D+13 ; 1.350D+12
55 for j = 1:length(Gcd)
56     for i = [0 10]
57         ty = 2^(j-1)+i;  // inttype: [1 11 2 12 4 14 8 18]
58         [g, m] = gcd(iconvert(v(:,j), ty));
59         assert_checkequal(inttype(g), ty);
60         assert_checkequal(g, iconvert(Gcd(j), ty));
61         assert_checkequal(inttype(m), modulo(ty,10));
62         if ty~=18
63             assert_checkequal(v(:,j)'*m, iconvert([0 Gcd(j)],inttype(m)));
64             // Otherwise there are expected round-off errors
65         end
66     end
67 end
68
69 // With decimal integers
70 // ---------------------
71 // Decimals  2^31 ~ 2.147D+09 < X < 2^52 ~ 4.504D+15
72 v = [d*2*25*31 ; d*29*37];  //    1.031D+13 ; 7.135D+12
73 [g, f] = gcd(v);
74 assert_checkequal(g, d);
75 assert_checkequal(v'*f, [0 g]);
76
77
78 // With negative numbers
79 // ---------------------
80 // http://bugzilla.scilab.org/15058
81 v = v(:)';
82 v(1) = -v(1);
83 [g, f] = gcd(v);
84 assert_checktrue(g>0);
85 assert_checkequal(g, d);
86 assert_checkequal(v*f, [0 g]);
87 v = -v;
88 [g, f] = gcd(v);
89 assert_checktrue(g>0);
90 assert_checkequal(g, d);
91 assert_checkequal(v*f, [0 g]);
92 v = -abs(v);
93 [g, f] = gcd(v);
94 assert_checktrue(g>0);
95 assert_checkequal(g, d);
96 assert_checkequal(v*f, [0 g]);
97
98
99 // With polynomials
100 // ----------------
101 d = (1-%z)*(2+%z)*(5*%z -4);
102 v = d*[%z (%z+1)];
103 [g, f] = gcd(v);
104 assert_checkequal(type(g), 2);
105 assert_checktrue(clean(g+d/5)==0);
106 assert_checkequal(clean(f - [-5-5*%z, 0.2 ; 5*%z, -0.2]),zeros(2,2)*%z);