/*
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 2.1 of the License, or
(at your option) any later version. See .
*/
#include
#include "flint.h"
#include "ulong_extras.h"
/* return (x|y)*(-1)^tstbit(r,1) */
int _n_jacobi_unsigned(mp_limb_t x, mp_limb_t y, unsigned int r)
{
mp_limb_t t, st;
int e;
FLINT_ASSERT(y & 1);
r ^= 2;
while (y > 1)
{
if (x == 0)
return 0;
/* x = odd part of x */
count_trailing_zeros(e, 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(mp_limb_t x, mp_limb_t y)
{
return _n_jacobi_unsigned(x, y, 0);
}
int n_jacobi(mp_limb_signed_t x, mp_limb_t y)
{
return _n_jacobi_unsigned(FLINT_ABS(x), y, FLINT_SIGN_EXT(x) & y);
}