set debug(jacobi) 0
proc jacobi {a b} {
  variable debug
  ## coded by rmelton 9/25/12
  # ref: http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx
  if {$debug(jacobi)} {
    puts stderr "${::STATUS}: jacobi $a $b"
  }
  if {$b<=0 || ($b&1)==0} {
    if {$debug(jacobi)} {
      puts stderr "${::STATUS}: jacobi=0 because b is negative, or even"
    }
    return 0;
  }
  set j 1
  if {$a<0} {
    set a [expr {0-$a}]
    if {($b & 3) == 3} {
      if {$debug(jacobi)} {
        puts stderr "${::STATUS}: use identity: Jacobi(-a,b) = -Jacobi(a,b) if b≡3 (mod 4)"
      }
      set j [expr {0-$j}]
    } else {
      if {$debug(jacobi)} {
        puts stderr "${::STATUS}: use identity: Jacobi(-a,b) = Jacobi(a,b) if b≡1 (mod 4)"
      }
    }
  }
  while {$a != 0} {
    while {($a&1) == 0} {
      ##/* Process factors of 2: Jacobi(2,b)=-1 if b=3,5 (mod 8) */
      set a [expr {$a>>1}]
      if {(($b & 7)==3) || (($b & 7)==5)} {
        set j [expr {0-$j}]
        if {$debug(jacobi)} {
          puts stderr "${::STATUS}: removed factors of 2 by identity: Jacobi(2a,b) = -Jacobi(a,b) if b≡3 (mod 8)"
        }
      } else {
        if {$debug(jacobi)} {
          puts stderr "${::STATUS}: removed factors of 2 by identity: Jacobi(2a,b) = Jacobi(a,b) if b≡1 (mod 8)"
        }
      }
    }
    ##/* Quadratic reciprocity: Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4) */
    lassign [list $a $b] b a
    if {(($a & 3)==3) && (($b & 3)==3)} {
      set j [expr {0-$j}]
      if {$debug(jacobi)} {
        puts stderr "${::STATUS}: use Identity Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4)"
      }
    } else {
      if {$debug(jacobi)} {
        puts stderr "${::STATUS}: use Identity Jacobi(a,b)=Jacobi(b,a) if a≡1 or b≡1 (mod 4)"
      }
    }
    if {$debug(jacobi)} {
      puts stderr "${::STATUS}: use Identity Jacobi(a Mod b,b) = Jacobi(a,b)"
    }
    set a [expr {$a % $b}]
    if {$debug(jacobi)} {
      puts stderr "${::STATUS}: Jacobi has been reduced to-> $j*Jacobi($a,$b)"
    }
  }
  if {$b==1} {
    if {$debug(jacobi)} {
      puts stderr "${::STATUS}: jacobi=$j (because we reduced to $j*Jacobi(0,1)=1)"
    }
    return $j
  } else {
    if {$debug(jacobi)} {
      puts stderr "${::STATUS}: jacobi=0 (because we reduced to $j*Jacobi(0,1)=0 n!=1)"
    }
    return 0
  }
}