1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
/*
Copyright (C) 2008 Peter Shrimpton
Copyright (C) 2009 William Hart
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "ulong_extras.h"
/* return (x|y)*(-1)^tstbit(r,1) */
int _n_jacobi_unsigned(ulong x, ulong y, unsigned int r)
{
ulong t, st;
int e;
FLINT_ASSERT(y & 1);
r ^= 2;
while (y > 1)
{
if (x == 0)
return 0;
/* x = odd part of x */
e = flint_ctz(x);
x >>= e;
r ^= ((y ^ (y>>1)) & (2*e)); /* (2|y) = (-1)^((y^2-1)/8) */
/* (x, y) = (|x - y|, min(x, y)) */
sub_ddmmss(st, t, UWORD(0), x, UWORD(0), y);
r ^= (x & y & st); /* if y > x, (x|y) = (y|x)*(-1)^((x-1)(y-1)/4) */
y += (st & t);
x = (t ^ st) - st;
}
return (int)(r & 2) - 1;
}
int n_jacobi_unsigned(ulong x, ulong y)
{
return _n_jacobi_unsigned(x, y, 0);
}
int n_jacobi(slong x, ulong y)
{
return _n_jacobi_unsigned(FLINT_UABS(x), y, FLINT_SIGN_EXT(x) & y);
}