NZMATH  1.2.0 About: NZMATH is a Python based number theory oriented calculation system.   Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)
arith1.py
Go to the documentation of this file.
1 """
2 Miscellaneous arithmetic functions
3 """
4
5 import itertools
6 import random
7 import bisect
8 import nzmath.gcd as gcd
9 from nzmath.plugins import MATHMODULE as math
10
11
12 def floorsqrt(a):
13  """
14  Return the floor of square root of the given integer.
15  """
16  if a < (1 << 59):
17  return int(math.sqrt(a))
18  else:
19  # Newton method
20  x = pow(10, (log(a, 10) >> 1) + 1) # compute initial value
21  while True:
22  x_new = (x + a//x) >> 1
23  if x <= x_new:
24  return x
25  x = x_new
26
27 def floorpowerroot(n, k, return_power = False):
28  """
29  Return the floor of k-th power root of the given integer n.
30  Use Newton method.
31  """
32  if k == 1:
33  if return_power:
34  return n, n
35  else:
36  return n
37  elif k == 2:
38  rslt = floorsqrt(n)
39  if return_power:
40  return rslt, rslt ** 2
41  else:
42  return rslt
43
44  if n <= 0:
45  if 0 == n:
46  if return_power:
47  return 0, 0
48  else:
49  return 0
50  elif 0 == k & 1:
51  raise ValueError("%d has no real %d-th root." % (n, k))
52  else:
53  sign = -1
54  n = -n
55  else:
56  sign = 1
57
58  # compute initial values
59  exp, rem = divmod(log(n) + 1, k)
60  if 0 != rem: # not ceiling
61  exp += 1
62  q = n >> (exp * (k - 1))
63  x = 1 << exp
64
65  # iteration by tangent approximation
66  while True:
67  x += (q - x) // k
68  z = x ** (k - 1)
69  q = n // z
70  if x <= q:
71  break
72
73  if sign < 0:
74  x = -x
75
76  if return_power:
77  return x, x * z
78  else:
79  return x
80
81 def powerDetection(n, largest_exp=False):
82  """
83  param positive integer n
84  param boolean largest_exp
85  return integer x, k s.t. n == x ** k
86  (2 <= k if exist else x, k == n, 1)
87  if largest_exp is true then return largest k
88  """
89  from nzmath.prime import generator_eratosthenes as generator_eratosthenes
90
91  ge = generator_eratosthenes(log(n, 2))
92  for exp in ge:
93  power_root, power = floorpowerroot(n, exp, True)
94  if power == n:
95  if largest_exp:
96  x, k = powerDetection(power_root, True)
97  return x, k * exp
98  else:
99  return power_root, exp
100
101  return n, 1
102
103 def legendre(a, m):
104  """
105  This function returns the Legendre symbol (a/m).
106  If m is an odd composite then this is the Jacobi symbol.
107  """
108  a = a % m
109  symbol = 1
110  while a != 0:
111  while a & 1 == 0:
112  a >>= 1
113  if m & 7 == 3 or m & 7 == 5:
114  symbol = -symbol
115  a, m = m, a
116  if a & 3 == 3 and m & 3 == 3:
117  symbol = -symbol
118  a = a % m
119  if m == 1:
120  return symbol
121  return 0
122
123 def modsqrt(n, p, e=1):
124  """
125  This function returns one of the square roots of n for mod p**e.
126  p must be an odd prime.
127  e must be a positive integer.
128  If 1 < e then n must be relatively prime with p.
129  """
130  if 1 < e:
131  x = modsqrt(n, p)
132  if 0 == x:
133  raise ValueError("if 1 < e then n must be relatively prime with p")
134  ppower = p
135  z = inverse(x << 1, p)
136  for i in range(e - 1):
137  x += (n - x ** 2) // ppower * z % p * ppower
138  ppower *= p
139  return x
140
141  symbol = legendre(n, p)
142  if symbol == 1:
143  pmod8 = p & 7
144  if pmod8 != 1:
145  n %= p
146  if pmod8 == 3 or pmod8 == 7:
147  x = pow(n, (p >> 2) + 1, p)
148  else: # p & 7==5
149  x = pow(n, (p >> 3) + 1, p)
150  c = pow(x, 2, p)
151  if c != n:
152  x = (x * pow(2, p >> 2, p)) % p
153  else: #p & 7==1
154  d = 2
155  while legendre(d, p) != -1:
156  d = random.randrange(3, p)
157  s, t = vp(p-1, 2)
158  A = pow(n, t, p)
159  D = pow(d, t, p)
160  m = 0
161  for i in range(1, s):
162  if pow(A*(D**m), 1 << (s-1-i), p) == (p-1):
163  m += 1 << i
164  x = (pow(n, (t+1) >> 1, p) * pow(D, m >> 1, p)) % p
165  return x
166  elif symbol == 0:
167  return 0
168  else:
169  raise ValueError("There is no solution")
170
171 def expand(n, m):
172  """
173  This function returns m-adic expansion of n.
174  n and m should satisfy 0 <= n, 2 <= m.
175  """
176  k = []
177  while n >= m:
178  k.append(n % m)
179  n //= m
180  k.append(n)
181  return k
182
183 def inverse(x, n):
184  """
185  This function returns inverse of x for modulo n.
186  """
187  x = x % n
188  y = gcd.extgcd(n, x)
189  if y[2] == 1:
190  if y[1] < 0:
191  r = n + y[1]
192  return r
193  else:
194  return y[1]
195  raise ZeroDivisionError("There is no inverse for %d modulo %d." % (x, n))
196
197 def CRT(nlist):
198  """
199  This function is Chinese Remainder Theorem using Algorithm 2.1.7
200  of C.Pomerance and R.Crandall's book.
201
202  For example:
203  >>> CRT([(1,2),(2,3),(3,5)])
204  23
205  """
206  r = len(nlist)
207  if r == 1 :
208  return nlist [ 0 ] [ 0 ]
209
210  product = []
211  prodinv = []
212  m = 1
213  for i in range(1, r):
214  m = m*nlist[i-1][1]
215  c = inverse(m, nlist[i][1])
216  product.append(m)
217  prodinv.append(c)
218
219  M = product[r-2]*nlist[r-1][1]
220  n = nlist[0][0]
221  for i in range(1, r):
222  u = ((nlist[i][0]-n)*prodinv[i-1]) % nlist[i][1]
223  n += u*product[i-1]
224  return n % M
225
226 def AGM(a, b):
227  """
228  Arithmetic-Geometric Mean.
229  """
230  x = (a+b) * 0.5
231  y = math.sqrt(a*b)
232  while abs(x-y) > y*1e-15:
233  x, y = (x+y) * 0.5, math.sqrt(x*y)
234  return x
235
236 #def _BhaskaraBrouncker(n):
237 # """
238 #
239 # _BhaskaraBrouncker returns the minimum tuple (p,q) such that:
240 # p ** 2 - n * q ** 2 = 1 or -1,
241 # for positive integer n, which is not a square.
242 #
243 # A good approximation for square root of n is given by the ratio
244 # p/q; the error is at most 1/2*q**2.
245 #
246 # """
247 # floorOfSqrt = floorsqrt(n)
248 # a = floorOfSqrt
249 # b0, b1 = 0, floorOfSqrt
250 # c0, c1 = 1, n - floorOfSqrt * floorOfSqrt
251 # p0, p1 = 1, floorOfSqrt
252 # q0, q1 = 0, 1
253 # while c1 != 1:
254 # a = (floorOfSqrt + b1)//c1
255 # b0, b1 = b1, a * c1 - b1
256 # c0, c1 = c1, c0 + a * (b0 - b1)
257 # p0, p1 = p1, p0 + a * p1
258 # q0, q1 = q1, q0 + a * q1
259 # return (p1, q1)
260
261 def vp(n, p, k=0):
262  """
263  Return p-adic valuation and indivisible part of given integer.
264
265  For example:
266  >>> vp(100, 2)
267  (2, 25)
268
269  That means, 100 is 2 times divisible by 2, and the factor 25 of
270  100 is indivisible by 2.
271
272  The optional argument k will be added to the valuation.
273  """
274  q = p
275  while not (n % q):
276  k += 1
277  q *= p
278
279  return (k, n // (q // p))
280
281 class _Issquare:
282  """
283  A class for testing whether a number is square or not.
284  The function issquare is an instance of the class, indeed.
285  """
286  q11 = [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0]
287  q63 = [1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
288  0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0,
289  1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]
290  q64 = [1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
291  0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,
292  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
293  q65 = [1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
294  0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0,
295  0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1]
296
297  def __call__(self, c):
298  """
299  Test whether a given number is a square number or not. If
300  the number is a square number, the function returns its square
301  root. Otherwise zero is returned.
302  """
303  t = c & 63
304  if not self.q64[t]:
305  return 0
306  r = c % 45045 # 45045 = 63 * 65 * 11
307  if not self.q63[r % 63]:
308  return 0
309  if not self.q65[r % 65]:
310  return 0
311  if not self.q11[r % 11]:
312  return 0
313  t = floorsqrt(c)
314  if t * t == c:
315  return t
316  else:
317  return 0
318
319 # test whether a given number is a square number or not.
320 issquare = _Issquare()
321 """test whether a given number is a square number or not."""
322
323
324 # choosing 4096 to cover at least 1000 digits
325 POWERS_LIST_LENGTH = 4096
326 POWERS_OF_2 = [2**i for i in range(POWERS_LIST_LENGTH)]
327
328 def log(n, base=2):
329  """
330  Return the integer part of logarithm of the given natural number
331  'n' to the 'base'. The default value for 'base' is 2.
332  """
333  if n < base:
334  return 0
335  if base == 10:
336  return _log10(n)
337  if base == 2 and n < POWERS_OF_2[-1]:
338  # a shortcut
339  return bisect.bisect_right(POWERS_OF_2, n) - 1
340  fit = base
341  result = 1
342  stock = [(result, fit)]
343  while fit < n:
344  next = fit ** 2
345  if next <= n:
346  fit = next
347  result += result
348  stock.append((result, fit))
349  else:
350  break
351  else: # just fit
352  return result
353  stock.reverse()
354  for index, power in stock:
355  prefit = fit * power
356  if prefit == n:
357  result += index
358  break
359  elif prefit < n:
360  fit = prefit
361  result += index
362  return result
363
364 def _log10(n):
365  return len(str(n))-1
366
367 def product(iterable, init=None):
368  """
369  product(iterable) is a product of all elements in iterable. If
370  init is given, the multiplication starts with init instead of the
371  first element in iterable. If the iterable is empty, then init or
372  1 will be returned.
373
374  If iterable is an iterator, it will be exhausted.
375  """
376  iterator = iter(iterable)
377  if init is None:
378  try:
379  result = iterator.next()
380  except StopIteration:
381  return 1 # empty product
382  else:
383  result = init
384  for item in iterator:
385  result *= item
386  return result
nzmath.arith1.expand
def expand(n, m)
Definition: arith1.py:171
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.arith1._log10
def _log10(n)
Definition: arith1.py:364
nzmath.arith1._Issquare.q11
list q11
Definition: arith1.py:286
nzmath.arith1.product
def product(iterable, init=None)
Definition: arith1.py:367
nzmath.gcd
Definition: gcd.py:1
nzmath.prime.generator_eratosthenes
def generator_eratosthenes(n)
Definition: prime.py:272
nzmath.arith1._Issquare.q65
list q65
Definition: arith1.py:293
nzmath.arith1.floorpowerroot
def floorpowerroot(n, k, return_power=False)
Definition: arith1.py:27
nzmath.arith1.powerDetection
def powerDetection(n, largest_exp=False)
Definition: arith1.py:81
nzmath.arith1._Issquare.q64
list q64
Definition: arith1.py:290
nzmath.arith1._Issquare
Definition: arith1.py:281
nzmath.arith1.floorsqrt
def floorsqrt(a)
Definition: arith1.py:12
nzmath.arith1.inverse
def inverse(x, n)
Definition: arith1.py:183
nzmath.plugins
Definition: plugins.py:1
nzmath.prime
Definition: prime.py:1
nzmath.arith1._Issquare.q63
list q63
Definition: arith1.py:287
nzmath.arith1.log
def log(n, base=2)
Definition: arith1.py:328
nzmath.arith1._Issquare.__call__
def __call__(self, c)
Definition: arith1.py:297
nzmath.arith1.vp
def vp(n, p, k=0)
Definition: arith1.py:261
nzmath.arith1.AGM
def AGM(a, b)
Definition: arith1.py:226
nzmath.arith1.CRT
def CRT(nlist)
Definition: arith1.py:197
nzmath.arith1.legendre
def legendre(a, m)
Definition: arith1.py:103
nzmath.arith1.modsqrt
def modsqrt(n, p, e=1)
Definition: arith1.py:123