calculateLegendre.py 1.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576
  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. >>> calculateLegendre(2, 31)
  30. 1
  31. >>> calculateLegendre(5, 31)
  32. 1
  33. >>> calculateLegendre(150, 1009) # http://math.stackexchange.com/q/221223/6876
  34. 1
  35. >>> calculateLegendre(25, 1009) # http://math.stackexchange.com/q/221223/6876
  36. 1
  37. >>> calculateLegendre(2, 1009) # http://math.stackexchange.com/q/221223/6876
  38. 1
  39. >>> calculateLegendre(3, 1009) # http://math.stackexchange.com/q/221223/6876
  40. 1
  41. """
  42. if a >= p or a < 0:
  43. return calculateLegendre(a % p, p)
  44. elif a == 0 or a == 1:
  45. return a
  46. elif a == 2:
  47. if p%8 == 1 or p%8 == 7:
  48. return 1
  49. else:
  50. return -1
  51. elif a == p-1:
  52. if p%4 == 1:
  53. return 1
  54. else:
  55. return -1
  56. elif not isPrime(a):
  57. factors = factorize(a)
  58. product = 1
  59. for pi in factors:
  60. product *= calculateLegendre(pi, p)
  61. return product
  62. else:
  63. if ((p-1)/2)%2==0 or ((a-1)/2)%2==0:
  64. return calculateLegendre(p, a)
  65. else:
  66. return (-1)*calculateLegendre(p, a)
  67. if __name__ == "__main__":
  68. import doctest
  69. doctest.testmod()