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:
54:
55:
56:
57:
58:
59:
60:
61:

{---------------------------------------------------------------------------}
function jacobi32(a,b: longint): longint;
  {-Compute the Jacobi/Legendre symbol (a|b), b: odd and > 2}
var
  j,m8: integer;
  t: longint;
begin
  {This is the classic method using only mod see e.g. [10] Alg. 2.3.5.}
  {Binary Pascal algorithm as used in kron_intern is NOT faster.}

  jacobi32 := 0;
  if (b<3) or (b and 1 = 0) then begin
    {$ifdef MPC_HaltOnError}
      {$ifdef MPC_UseExceptions}
        raise MPXRange.Create('jacobi32: b<3 or even');
      {$else}
        RunError(MP_RTE_RANGE);
      {$endif}
    {$else}
      set_mp_error(MP_RANGE);
      exit;
    {$endif}
  end;

  {initialize accumulated result}
  j := 1;

  {if a<0 use property of Jacobi symbol to make it positive. Don't}
  {rely on Pascal mod function to return positive result.}
  if a<0 then begin
    {(-a|b) = (a|b)*(-1|b); change sign if b=3 (mod 4)}
    if b and 3 = 3 then j := -j;
    a := -a;
  end;

  {do initial reduction to force a < b}
  if a>=b then a := a mod b;

  {loop invariant: a < b}
  while a<>0 do begin
    {divide out powers of 4}
    while a and 3 = 0 do a := a shr 2;
    {if a is even, divide by 2 and adjust result}
    if a and 1 = 0 then begin
      a  := a shr 1;
      m8 := b and 7;
      { (2|b) = -1  if b = +-3 (mod 8)}
      if (m8=3) or (m8=5) then j:=-j;
    end;
    if a=1 then begin
      jacobi32 := j;
      exit;
    end;
    {swap variables, reduce and adjust result using quadratic reciprocity}
    if (a and b and 3) = 3 then j := -j;
    t := b;
    b := a;
    a := t mod a;
  end;
  if b=1 then jacobi32 := j;
end;