1 // This file is part of the FXT library.
2 // Copyright (C) 2010, 2012 Joerg Arndt
3 // License: GNU General Public License version 3 or later,
4 // see the file COPYING.txt in the main directory.
5
6 #include "mod/mtypes.h"
7 //#include "mod/factor.h"
8
9
10 int
kronecker(umod_t a,umod_t b)11 kronecker(umod_t a, umod_t b)
12 // Return Kronecker symbol (a/b).
13 // Equal to Legendre symbol (a/b) if b is an odd prime.
14 //.
15 // Cf. Cohen p.29
16 {
17 // if ( a>=b ) a %= b;
18 static const int tab2[] = {0, 1, 0, -1, 0, -1, 0, 1};
19 // tab2[ a & 7 ] := (-1)^((a^2-1)/8)
20
21 if ( 0==b ) return (1==a);
22
23 if ( 0==((a|b)&1) ) return 0; // a and b both even ?
24
25 int v = 0;
26 while ( 0==(b&1) ) { ++v; b>>=1; }
27
28 int k;
29 if ( 0==(v&1) ) k = 1;
30 else k = tab2[ a & 7 ];
31
32 // signed: if ( b<0 ) { b=-b; if (a<0) {k=-k}; }
33
34 // step3:
35 while ( 1 )
36 {
37 // here b is odd
38 if ( 0==a ) return ( b>1 ? 0 : k );
39
40 v = 0;
41 while ( 0==(a&1) ) { ++v; a>>=1; }
42
43 if ( 1==(v&1) ) k *= tab2[ b & 7 ]; // k *= (-1)**((b*b-1)/8)
44
45 // step4:
46 if ( a & b & 2 ) k = -k; // k = k*(-1)**((a-1)*(b-1)/4)
47
48 umod_t r = a; // signed: r = abs(a)
49 a = b % r;
50 b = r;
51 // goto step3;
52 }
53 }
54 // -------------------------
55
56
57
58 //int
59 //reciprocity(umod_t a, umod_t b)
60 ////
61 //// (a/b) = (b/a) * (-1)^((a1-1)*(b1-1)/4 + (as-1)*(bs-1)/4)
62 //// where a==a1*2^ae, b=b1*2^be
63 //// and as=sign(a), bs=sign(b)
64 //// cf. Cohen p.43
65 ////
66 //// UNTESTED
67 //{
68 // if ( (0==a) || (0==b) ) return 0;
69 // while ( 0==(a&1) ) a >>= 1;
70 // while ( 0==(b&1) ) b >>= 1;
71 // int k = 1;
72 // if ( a & b & 2 ) k = -k;
73 // return k;
74 //}
75 //// -------------------------
76