calculateLegendre.py 1.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. def isPrime(a):
  4. return all(a % i for i in xrange(2, a))
  5. # http://stackoverflow.com/a/14793082/562769
  6. def factorize(n):
  7. factors = []
  8. p = 2
  9. while True:
  10. while(n % p == 0 and n > 0): #while we can divide by smaller number, do so
  11. factors.append(p)
  12. n = n / p
  13. p += 1 #p is not necessary prime, but n%p == 0 only for prime numbers
  14. if p > n / p:
  15. break
  16. if n > 1:
  17. factors.append(n)
  18. return factors
  19. def calculateLegendre(a, p):
  20. """
  21. Calculate the legendre symbol (a, p) with p is prime.
  22. The result is either -1, 0 or 1
  23. >>> calculateLegendre(3, 29)
  24. -1
  25. >>> calculateLegendre(111, 41) # Beispiel aus dem Skript, S. 114
  26. -1
  27. >>> calculateLegendre(113, 41) # Beispiel aus dem Skript, S. 114
  28. 1
  29. """
  30. if a >= p or a < 0:
  31. return calculateLegendre(a % p, p)
  32. elif a == 0 or a == 1:
  33. return a
  34. elif a == 2:
  35. if a%8 == 1 or a%8 == 7:
  36. return 1
  37. else:
  38. return -1
  39. elif a == p-1:
  40. if p%4 == 1:
  41. return 1
  42. else:
  43. return -1
  44. elif not isPrime(a):
  45. factors = factorize(a)
  46. product = 1
  47. for pi in factors:
  48. product *= calculateLegendre(pi, p)
  49. return product
  50. else:
  51. if ((p-1)/2)%2==0:
  52. return calculateLegendre(p, a)
  53. else:
  54. return (-1)*calculateLegendre(p, a)
  55. if __name__ == "__main__":
  56. import doctest
  57. doctest.testmod()