DSPython  00.03.03 — 25 juin 2012
 Tout Classes Espaces de nommage Fichiers Fonctions Variables Pages
factors.py
Aller à la documentation de ce fichier.
1 #!/usr/bin/env python
2 # -*- coding: latin-1 -*-
3 ##\package DSPython.factors Facteurs premiers
4 
5 ##\file
6 # Facteurs premiers
7 
8 # (c) Olivier Pirson --- DragonSoft
9 # http://www.opimedia.be/DS/
10 # Débuté le 8 décembre 2006
11 ####################################
12 from __future__ import print_function
13 
14 ## Date du dernier changement pour ce module
15 VERSION = 'factors --- 2010 April 12'
16 
17 import bisect, fractions, types
18 
19 import DSPython
20 import DSPython.nat32 as nat32
21 import DSPython.natural as natural
22 
23 
24 
25 # ###########
26 # Fonctions #
27 #############
28 ##\brief Pour s = \htmlonly(&hellip; (P<sub>i</sub>, e<sub>i</sub>) &hellip;)\endhtmlonly
29 ## où P<sub>i</sub> est le i<sup>e</sup> nombre premier
30 ## renvoie \htmlonly&sum;<sub>i</sub> e<sub>i</sub>.2<sup>i - 1</sup>\endhtmlonly
31 def bef(s):
32  """Pour s = [... (P_i, e_i) ...] où P_i est le ième nombre premier
33  renvoie Sum_i e_i * 2**(i - 1)
34  [OEIS A048675] (Binary Encoding of Factorizations)
35 
36  Pre: s: séquence strictement ordonnée de primaries
37 
38  Result: natural
39 
40  O(s) = ..."""
41  assert primaries_is(s), s
42 
43  n = 0
44  for p in s:
45  i = bisect.bisect_left(nat32.PRIME_16_ARRAY, p[0])
46  if i >= len(nat32.PRIME_16_ARRAY):
47  raise NotImplementedError
48  n += p[1] * (2**i)
49  return n
50 
51 
52 ## s et t sont premiers entre eux ?
53 def coprime_is(s, t):
54  """Renvoie True si s et t sont premiers entre eux,
55  False sinon
56 
57  Pre: s: séquence strictement ordonnée de primaries
58  t: séquence strictement ordonnée de primaries
59 
60  Result: boolean
61 
62  O(s, t) = ..."""
63  assert primaries_is(s), s
64  assert primaries_is(t), t
65 
66  return gcd(s, t) == []
67 
68 
69 ##\brief Distance de Dominici à partir des primaries
70 #
71 # (cf. <i lang="en">An Arithmetic Metric</i>,
72 # Diego <span style="font-variant:small-caps">Dominici</span>,
73 # <a href="http://arxiv.org/abs/0906.0632/"
74 # target="_blank"><tt>http://arxiv.org/abs/0906.0632/</tt></a>)
75 def distance_dominici(s, t):
76  """Renvoie la distance de Dominici entre les naturels
77  ayant s et t pour liste de leurs primaries
78 
79  Pre: s: séquence strictement ordonnée de primaries
80  t: séquence strictement ordonnée de primaries
81 
82  Result: naturel
83  O(s) = ..."""
84  assert primaries_is(s), s
85  assert primaries_is(t), t
86 
87  d = 0
88  for p in lcm(s, t):
89  d += p[1]
90  for p in gcd(s, t):
91  d -= p[1]
92  return d
93 
94 
95 ## Nombre des diviseurs à partir des primaries (fonction \htmlonly &nu;\endhtmlonly)
96 def divisors_nb(s):
97  """Renvoie le nombre des diviseurs du naturel
98  ayant s pour liste des primaries
99  == divisors_sum_pow(s, 0)
100 
101  Pre: s: séquence strictement ordonnée de primaries
102 
103  Result: naturel >= 1
104 
105  O(s) = ..."""
106  assert primaries_is(s), s
107 
108  n = 1
109  for i in s:
110  n *= i[1] + 1
111  return n
112 
113 
114 ## Produit des diviseurs à partir des primaries (fonction \htmlonly &piv;\endhtmlonly)
115 def divisors_prod(s):
116  """Renvoie le produit des diviseurs du naturel
117  ayant s pour liste des primaries
118 
119  Pre: s: séquence strictement ordonnée de primaries
120 
121  Result: naturel >= 1
122 
123  O(s) = ..."""
124  assert primaries_is(s), s
125 
126  nb = divisors_nb(s)
127  if nb&1: # nb impair (ssi le naturel correspondant à s est un carré)
128  p = 1
129  for i in s:
130  assert i[1]%2 == 0, (i, s)
131 
132  p *= i[0]**(i[1]//2)
133  return p**nb
134  else: # nb pair (ssi le naturel correspondant à s n'est pas un carré)
135  return primaries_to_n(s)**(nb//2)
136 
137 
138 ## Somme des diviseurs à partir des primaries (fonction \htmlonly &sigma;\endhtmlonly)
139 def divisors_sum(s):
140  """Renvoie la somme des diviseurs du naturel
141  ayant s pour liste des primaries
142  == divisors_sum_pow(s, 1)
143 
144  Pre: s: séquence strictement ordonnée de primaries
145 
146  Result: naturel >= 1
147 
148  O(s) = ..."""
149  assert primaries_is(s), s
150 
151  n = divisors_sum_odd(s)
152  return (n if (len(s) == 0) or (s[0][0] != 2) # n impair
153  else n * natural.mersenne(s[0][1] + 1)) # n pair : M_(k+1) * divisors_sum_odd
154 
155 
156 ## Somme des diviseurs pairs à partir des primaries
157 def divisors_sum_even(s):
158  """Renvoie la somme des diviseurs pairs du naturel
159  ayant s pour liste des primaries
160 
161  Pre: s: séquence strictement ordonnée de primaries
162 
163  Result: naturel >= 1
164 
165  O(s) = ..."""
166  assert primaries_is(s), s
167 
168  if (len(s) == 0) or (s[0][0] != 2): # nombre impair
169  return 0
170  else: # nombre pair
171  n = 1
172  for i in s[1:]:
173  n *= (i[0]**(i[1] + 1) - 1) // (i[0] - 1)
174  return (n * natural.mersenne(s[0][1])) << 1 # == 2 * M_k * divisors_sum_odd
175 
176 
177 ## Somme des diviseurs impairs à partir des primaries
178 def divisors_sum_odd(s):
179  """Renvoie la somme des diviseurs pairs du naturel
180  ayant s pour liste des primaries
181 
182  Pre: s: séquence strictement ordonnée de primaries
183 
184  Result: naturel >= 1
185 
186  O(s) = ..."""
187  assert primaries_is(s), s
188 
189  n = 1
190  if (len(s) == 0) or (s[0][0] != 2): # nombre impair
191  for i in s:
192  n *= (i[0]**(i[1] + 1) - 1) // (i[0] - 1)
193  else: # nombre pair
194  for i in s[1:]:
195  n *= (i[0]**(i[1] + 1) - 1) // (i[0] - 1)
196  return n
197 
198 
199 ## Somme des diviseurs à partir des primaries (fonction \htmlonly &sigma;\endhtmlonly<sub>k</sub>)
200 def divisors_sum_pow(s, k):
201  """Renvoie la somme des diviseurs ** k du naturel
202  ayant s pour liste des primaries
203 
204  Pre: s: séquence strictement ordonnée de primaries
205  k: naturel
206 
207  Result: naturel >= 1
208 
209  O(s) = ..."""
210  assert primaries_is(s), s
211  assert DSPython.natural_is(k), k
212 
213  if k > 1:
214  n = 1
215  for i in s:
216  n *= (i[0]**(k*(i[1] + 1)) - 1) // (i[0]**k - 1)
217  return n
218  elif k == 1:
219  return divisors_sum(s)
220  else:
221  return divisors_nb(s)
222 
223 
224 ##\brief Pour n =
225 ## \htmlonly&sum;<sub>i=0</sub><sup>&infin;</sup> b<sub>i</sub>.2<sup>i</sup>\endhtmlonly
226 ## renvoie les primaries de
227 ## \htmlonly&prod;<sub>i=1</sub><sup>&infin;</sup>
228 ## P<sub>i</sub><sup>b<sub>i-1</sub></sup>\endhtmlonly
229 ## où P<sub>i</sub> est le i<sup>e</sup> nombre premier
230 def feb_primaries(n):
231  """Pour n = Sum_(i=0)^infini b_i * 2**i
232  renvoie les primaries de Prod_(i=1)^infini P_i**(b_(i-1))
233  où P_i est le ième nombre premier
234  [OEIS A019565]
235  bef(feb_primaries(n)) == n
236 
237  Pre: n: natural
238 
239  Result: séquence strictement ordonnée de primaries (naturel premier, 1)
240 
241  O(n) = ..."""
242  assert DSPython.natural_is(n), n
243 
244  i = 0
245  s = []
246  while (n > 0) and (i < len(nat32.PRIME_16_ARRAY)):
247  if n&1:
248  s.append((nat32.PRIME_16_ARRAY[i], 1))
249  i += 1
250  n >>= 1
251  if n > 0:
252  raise NotImplementedError
253  return s
254 
255 
256 ## PGCD de s et t (Plus Grand Commun Diviseur/ Greatest Common Divisor)
257 def gcd(s, t):
258  """Renvoie le PGCD de s et t
259  (Plus Grand Commun Diviseur/ Greatest Common Divisor)
260 
261  Pre: s: séquence strictement ordonnée de primaries
262  t: séquence strictement ordonnée de primaries
263 
264  Result: séquence strictement ordonnée de primaries
265 
266  O(s, t) = ..."""
267  assert primaries_is(s), s
268  assert primaries_is(t), t
269 
270  if (len(s) > 0) and (len(t) > 0):
271  s_i = 0
272  t_i = 0
273  l = []
274  while (len(s) > s_i) and (len(t) > t_i):
275  while s[s_i][0] < t[t_i][0]:
276  s_i += 1
277  if len(s) <= s_i:
278  return l
279  while s[s_i][0] > t[t_i][0]:
280  t_i += 1
281  if len(t) <= t_i:
282  return l
283  if s[s_i][0] != t[t_i][0]:
284  continue
285  l.append((s[s_i][0], min(s[s_i][1], t[t_i][1])))
286  s_i += 1
287  t_i += 1
288  return l
289  else:
290  return []
291 
292 
293 ##\brief Renvoie le "nombre de Gödel" de la séquence s :
294 ## \htmlonly
295 ## pour s == (a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>, &hellip;.., a<sub>k</sub>),
296 ## renvoie &prod;<sub>i=1</sub><sup>k</sup> P<sub>i</sub><sup>a<sub>i</sub></sup>
297 ## \endhtmlonly
298 def godelnumber(s):
299  """Renvoie le "nombre de Gödel" de la séquence s
300  Pour s == (a_1, a_2, a_3, ..., a_k), renvoie prod_i=1^k (P_i)^(a_i)
301  où P_i est le ième nombre premier
302 
303  Pre: s: séquence de naturels
304 
305  Result: naturel >= 1
306 
307  O() = ..."""
308  assert isinstance(s, list) or isinstance(s, tuple), type(s)
309 
310  n = 1
311  if len(s) > len(nat32.PRIME_16_ARRAY):
312  raise NotImplementedError
313  for i in range(len(s)):
314  assert DSPython.natural_is(s[i]), (i, type(s[i]), s[i])
315 
316  n *= nat32.PRIME_16_ARRAY[i]**s[i]
317  return n
318 
319 
320 ##\brief
321 ## \htmlonly
322 ## Pour n == &prod;<sub>i=1</sub><sup>k</sup> P<sub>i</sub><sup>a<sub>i</sub></sup>,
323 ## renvoie [a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>, &hellip;.., a<sub>k</sub>]
324 ## \endhtmlonly
325 def godelnumber_to_list(n):
326  """Pour n == prod_i=1^k (P_i)^(a_i) renvoie [a_1, a_2, a_3, ..., a_k]
327  Si n == 1 alors renvoie []
328 
329  Pre: n: naturel >= 1
330 
331  Result: [] ou list de naturels dont le dernier élément est != 0
332 
333  O() = ..."""
334  assert DSPython.natural_is(n), n
335  assert n > 0, n
336 
337  if n == 1:
338  return []
339  l = [natural.scan1(n)]
340  n >>= l[0]
341  for p in nat32.PRIME_16_ARRAY[1:]: # tant qu'on sait diviser n
342  if n == 1:
343  return l
344  a = 0
345  q, r = divmod(n, p)
346  while r == 0: # tant que divisible
347  a += 1
348  n = q
349  q, r = divmod(n, p)
350  l.append(a)
351  if n == 1:
352  return l
353  raise NotImplementedError # plus de nombre premier dans la table
354 
355 
356 ## PPCM de s et t (Plus Petit Commun Multiple/ Least Common Multiple)
357 def lcm(s, t):
358  """Renvoie PPCM de s et t
359  (Plus Petit Commun Multiple/ Least Common Multiple)
360 
361  Pre: s: séquence strictement ordonnée de primaries
362  t: séquence strictement ordonnée de primaries
363 
364  Result: séquence strictement ordonnée de primaries
365 
366  O(s, t) = ..."""
367  assert primaries_is(s), s
368  assert primaries_is(t), t
369 
370  if (len(s) > 0) and (len(t) > 0):
371  s_i = 0
372  t_i = 0
373  l = []
374  while (len(s) > s_i) and (len(t) > t_i):
375  if s[s_i][0] < t[t_i][0]:
376  l.append((s[s_i][0], s[s_i][1]))
377  s_i += 1
378  elif s[s_i][0] > t[t_i][0]:
379  l.append((t[t_i][0], t[t_i][1]))
380  t_i += 1
381  else:
382  l.append((s[s_i][0], max(s[s_i][1], t[t_i][1])))
383  s_i += 1
384  t_i += 1
385  if len(s) > s_i:
386  return l + s[s_i:]
387  return (l + t[t_i:] if len(t) > t_i
388  else l)
389  elif len(s) > 0: # t == []
390  return s
391  else: # s == []
392  return t
393 
394 
395 ## Fonction de Möbius à partir des primaries (fonction \htmlonly &mu;\endhtmlonly)
396 def mobius(s):
397  """Renvoie la fonction de Möbius du naturel
398  ayant s pour liste des primaries
399 
400  Pre: s: séquence strictement ordonnée de primaries
401 
402  Result: -1, 0 ou 1
403 
404  O(s) = ..."""
405  assert primaries_is(s), s
406 
407  k = 0
408  for i in s:
409  if i[1] > 1:
410  return 0
411  k += 1
412  return (1 if k&1 == 0
413  else -1)
414 
415 
416 ## Produit de s et t
417 def mul(s, t):
418  """Renvoie le produit s et t
419 
420  Pre: s: séquence strictement ordonnée de primaries
421  t: séquence strictement ordonnée de primaries
422 
423  Result: séquence strictement ordonnée de primaries
424 
425  O(s, t) = ..."""
426  assert primaries_is(s), s
427  assert primaries_is(t), t
428 
429  if (len(s) > 0) and (len(t) > 0):
430  s_i = 0
431  t_i = 0
432  l = []
433  while (len(s) > s_i) and (len(t) > t_i):
434  while s[s_i][0] < t[t_i][0]:
435  l.append((s[s_i][0], s[s_i][1]))
436  s_i += 1
437  if len(s) <= s_i:
438  return l + t[t_i:]
439  while s[s_i][0] > t[t_i][0]:
440  l.append((t[t_i][0], t[t_i][1]))
441  t_i += 1
442  if len(t) <= t_i:
443  return l + s[s_i:]
444  if s[s_i][0] == t[t_i][0]:
445  l.append((s[s_i][0], s[s_i][1] + t[t_i][1]))
446  s_i += 1
447  t_i += 1
448  if len(s) > s_i:
449  return l + s[s_i:]
450  return (l + t[t_i:] if len(t) > t_i
451  else l)
452  elif len(s) > 0: # t == []
453  return s
454  else: # s == []
455  return t
456 
457 
458 ##\brief Nombre de décompositions (en tenant compte de l'ordre) en somme de 4 carrés d'entiers
459 ## (fonction GR)
460 def nb_in_integers_4sqr(s):
461  """Renvoie le nombre de décompositions (en tenant compte de l'ordre)
462  en somme de 4 carrés d'entiers du naturel
463  ayant s pour liste des primaries
464 
465  Pre: s: séquence strictement ordonnée de primaries
466 
467  Result: naturel multiple de 8
468 
469  O(s) = ..."""
470  assert primaries_is(s), s
471 
472  return ((divisors_sum_odd(s) * 3) << 3 if (len(s) > 0) and (s[0][0] == 2) # nombre pair
473  else divisors_sum_odd(s) << 3) # nombre impair
474 
475 
476 ## Produit des facteurs non communs de s et t (Not common Factors Product)
477 def nfp(s, t):
478  """Renvoie le produit des facteurs non communs de s et t
479  (Not common Factors Product)
480 
481  Pre: s: séquence strictement ordonnée de primaries
482  t: séquence strictement ordonnée de primaries
483 
484  Result: séquence strictement ordonnée de primaries
485 
486  O(s, t) = ..."""
487  assert primaries_is(s), s
488  assert primaries_is(t), t
489 
490  if (len(s) > 0) and (len(t) > 0):
491  s_i = 0
492  t_i = 0
493  l = []
494  while (len(s) > s_i) and (len(t) > t_i):
495  while s[s_i][0] < t[t_i][0]:
496  l.append((s[s_i][0], s[s_i][1]))
497  s_i += 1
498  if len(s) <= s_i:
499  return l + t[t_i:]
500  while s[s_i][0] > t[t_i][0]:
501  l.append((t[t_i][0], t[t_i][1]))
502  t_i += 1
503  if len(t) <= t_i:
504  return l + s[s_i:]
505  if s[s_i][0] == t[t_i][0]:
506  if s[s_i][1] > t[t_i][1]:
507  l.append((s[s_i][0], s[s_i][1] - t[t_i][1]))
508  elif s[s_i][1] < t[t_i][1]:
509  l.append((s[s_i][0], t[t_i][1] - s[s_i][1]))
510  s_i += 1
511  t_i += 1
512  if len(s) > s_i:
513  return l + s[s_i:]
514  return (l + t[t_i:] if len(t) > t_i
515  else l)
516  elif len(s) > 0: # t == []
517  return s
518  else: # s == []
519  return t
520 
521 
522 ## Nontotient à partir des primaries
523 def nontotient(s):
524  """Renvoie le nombre de nontotatives du naturel
525  ayant s pour liste des primaries
526  (Si s == [] alors renvoie 0)
527 
528  Pre: s: séquence strictement ordonnée de primaries
529 
530  Result: naturel
531 
532  O(n) = ..."""
533  assert primaries_is(s), s
534 
535  return primaries_to_n(s) - totient(s)
536 
537 
538 ##\brief
539 ## \htmlonly
540 ## Pour s = ((p<sub>1</sub>, &alpha;<sub>1</sub>), (p<sub>2</sub>, &alpha;<sub>2</sub>),
541 ## (p<sub>3</sub>, &alpha;<sub>3</sub>), &hellip;, (p<sub>k</sub>, &alpha;<sub>k</sub>))
542 ## et t = ((q<sub>1</sub>, &beta;<sub>1</sub>), (q<sub>2</sub>, &beta;<sub>2</sub>),
543 ## (q<sub>3</sub>, &beta;<sub>3</sub>), &hellip;, (q<sub>l</sub>, &beta;<sub>l</sub>)),<br>
544 ## renvoie [(r<sub>1</sub>, &gamma;<sub>1</sub>), (r<sub>2</sub>, &gamma;<sub>2</sub>),
545 ## (r<sub>3</sub>, &gamma;<sub>3</sub>), &hellip;, (r<sub>m</sub>, &gamma;<sub>m</sub>)]<br>
546 ## tel que pour les p<sub>i</sub> = q<sub>j</sub>
547 ## on ai r<sub>n</sub> = p<sub>i</sub>
548 ## et &gamma;<sub>n</sub> = f(&alpha;<sub>i</sub>, &beta;<sub>j</sub>).<br>
549 ## Pour les p<sub>i</sub> &ne; q<sub>j</sub>
550 ## on prend &alpha;<sub>i</sub> = 0 ou &beta;<sub>j</sub> = 0.
551 ## On s'arrête dès il n'y a plus de primary.
552 ## \endhtmlonly
553 def ope_exp(s, t, f):
554  """Renvoie []
555 
556  Pre: s: séquence strictement ordonnée de primaries
557  t: séquence strictement ordonnée de primaries
558  f: fonction : (naturel >= 1, naturel >= 1) -> naturel
559 
560  Result: séquence strictement ordonnée de primaries
561 
562  O(s, t) = ..."""
563  assert primaries_is(s), s
564  assert primaries_is(t), t
565  assert isinstance(f, types.FunctionType), type(f)
566 
567  s_i = 0 # position dans s
568  t_i = 0 # position dans t
569  l = [] # résultat
570  while (len(s) > s_i) and (len(t) > t_i): # il reste des p_i et des q_j
571  if s[s_i][0] == t[t_i][0]: # p_i = q_j
572  gamma = f(s[s_i][1], t[t_i][1])
573  if gamma > 0:
574  l.append((s[s_i][0], gamma))
575  s_i += 1
576  t_i += 1
577  elif s[s_i][0] > t[t_i][0]: # p_i > q_j
578  gamma = f(0, t[t_i][1])
579  if gamma > 0:
580  l.append((t[t_i][0], gamma))
581  t_i += 1
582  else: # p_i < q_j
583  gamma = f(s[s_i][1], 0)
584  if gamma > 0:
585  l.append((s[s_i][0], gamma))
586  s_i += 1
587  while len(s) > s_i: # il ne reste que des p_i
588  gamma = f(s[s_i][1], 0)
589  if gamma > 0:
590  l.append((s[s_i][0], gamma))
591  s_i += 1
592  while len(t) > t_i: # il ne reste que des q_j
593  gamma = f(0, t[t_i][1])
594  if gamma > 0:
595  l.append((t[t_i][0], gamma))
596  t_i += 1
597 
598  assert primaries_is(l), l
599 
600  return l
601 
602 
603 ##\brief Liste des primaries de n,\n
604 ## \htmlonly
605 ## c.-à-d. [(p<sub>1</sub>, &alpha;<sub>1</sub>), (p<sub>2</sub>, &alpha;<sub>2</sub>),
606 ## (p<sub>3</sub>, &alpha;<sub>3</sub>), &hellip;, (p<sub>k</sub>, &alpha;<sub>k</sub>)]
607 ## tel que n = &prod;<sub>i=1</sub><sup>k</sup> p<sub>i</sub><sup>&alpha;<sub>i</sub></sup>
608 ## avec p<sub>i</sub> premier,
609 ## p<sub>i</sub> &lt; p<sub>i + 1</sub>
610 ## et &alpha;<sub>i</sub> &ge; 1
611 ## \endhtmlonly
612 def primaries(n):
613  """Renvoie la liste des primaries (dans l'ordre strictement croissant) de n
614  c.-à-d. les couples (facteur premier, exposant) décomposant n
615  (Si n == 1 alors renvoie []).
616  Le nombre de primaries (la longueur de la liste renvoyée) correspond
617  au nombre de facteurs premiers distincts de n
618  (fonction omega)
619 
620  Pre: n: naturel >= 1
621 
622  Result: séquence strictement ordonnée (sur le premier élément) de couples
623  (naturel premier, naturel >= 1)
624 
625  O(n) = ..."""
626  assert DSPython.natural_is(n), n
627  assert n > 0, n
628 
629  # Facteurs 2
630  nb = natural.scan1(n)
631  if nb > 0:
632  n >>= nb
633  l = [(2, nb)]
634  else:
635  l = []
636 
637  # Facteurs premiers tenant sur 16 bits
638  for i in range(1, len(nat32.PRIME_16_ARRAY)):
639  nb = 0
640  d = nat32.PRIME_16_ARRAY[i]
641  q, r = divmod(n, d)
642  while r == 0:
643  nb += 1
644  n = q
645  q, r = divmod(n, d)
646  if nb > 0:
647  l.append((d, nb))
648  if n == 1:
649  return l
650 
651  assert d%nat32.PRIMORIAL32_TUPLE[4] == 1, \
652  (d, nat32.PRIMORIAL32_TUPLE[4], d%nat32.PRIMORIAL32_TUPLE[4])
653 
654  # Facteurs premiers ne tenant pas sur 16 bits
655  while True:
656  for m in nat32.MODS_PRIMORIAL_4_DIFF_TUPLE:
657  # Parcourt les restes modulo 210 premiers avec 210
658  nb = 0
659  q, r = divmod(n, d)
660  while r == 0:
661  nb += 1
662  n = q
663  q, r = divmod(n, d)
664  if nb > 0:
665  l.append((d, nb))
666  if n == 1:
667  return l
668  d += m
669 
670 
671 ## Renvoie True si p est une séquence strictement ordonnée de primaries, False sinon
672 def primaries_is(s):
673  """Renvoie True si p est une séquence strictement ordonnée de primaries,
674  False sinon
675 
676  Pre: s: quelconque
677 
678  Result: boolean"""
679  if isinstance(s, list) or isinstance(s, tuple):
680  prev_n = -1
681  prev_e = -1
682  for i in s:
683  if ((not isinstance(i, tuple)) or (len(i) != 2)
684  or (not DSPython.natural_is(i[0])) or (not DSPython.natural_is(i[1]))
685  or (i[1] == 0) or (prev_n >= i[0]) or (not natural.prime_is(i[0]))):
686  return False
687  prev_n = i[0]
688  prev_e = i[1]
689  return True
690  else:
691  return False
692 
693 
694 ## Renvoie True si s représente un nombre premier, False sinon
695 def primaries_prime_is(s):
696  """Renvoie True si s représente un nombre premier, False sinon
697 
698  Pre: s: séquence strictement ordonnée de primaries
699 
700  Result: boolean"""
701  assert primaries_is(s), s
702 
703  return (len(s) == 1) and (s[0][1] == 1)
704 
705 
706 ## String représentant les primaries
707 def primaries_str(s, times='.', exp='^'):
708  """Renvoie un string représentant les primaries
709  (Si s est vide alors renvoie '1')
710 
711  Pre: s: séquence strictement ordonnée de primaries
712  times: string représentant l'opérateur de multiplication
713  exp: string représentant l'opérateur d'exponentiation
714 
715  Result: string
716 
717  O(s) = ..."""
718  assert primaries_is(s), s
719 
720  if len(s) > 0:
721  l = []
722  for i in s:
723  l.append('{0}{1}{2}'.format(i[0], exp, i[1]) if i[1] > 1
724  else str(i[0]))
725  return times.join(l)
726  else:
727  return '1'
728 
729 
730 ## Produit des primaries
731 def primaries_to_n(s):
732  """Renvoie le naturel ayant s pour liste des primaries
733  (Si s est vide alors renvoie 1)
734 
735  Pre: s: séquence strictement ordonnée de primaries
736 
737  Result: naturel >= 1
738 
739  O(s) = ..."""
740  assert primaries_is(s), s
741 
742  p = 1
743  for i in s:
744  p *= i[0]**i[1]
745  return p
746 
747 
748 ## Liste des primes à partir des primaries
749 def primaries_to_primes(s):
750  """Renvoie la liste des primes à partir de la liste des primaries
751  (Si s est vide alors renvoie [])
752 
753  Pre: s: séquence strictement ordonnée de primaries
754 
755  Result: séquence ordonnée de primes
756 
757  O(s) = ..."""
758  assert primaries_is(s), s
759 
760  r = []
761  for c in s:
762  r += [c[0]] * c[1]
763  return r
764 
765 
766 ## Liste des primes de n
767 def primes(n):
768  """Renvoie la liste des primes (dans l'ordre croissant) de n
769  c.-à-d. les facteurs premiers de n
770  (Si n == 1 alors renvoie []).
771  Le nombre de primes (la longueur de la liste renvoyée) correspond
772  au nombre de facteurs premiers (non nécessairement distincts) de n
773  (fonction Omega)
774 
775  Pre: n: naturel >= 1
776 
777  Result: séquence ordonnée de naturels premiers
778 
779  O(n) = ..."""
780  assert DSPython.natural_is(n), n
781  assert n > 0, n
782 
783  # Facteurs 2
784  nb = natural.scan1(n)
785  n >>= nb
786  l = [2] * nb
787 
788  # Facteurs premiers tenant sur 16 bits
789  for i in range(1, len(nat32.PRIME_16_ARRAY)):
790  nb = 0
791  d = nat32.PRIME_16_ARRAY[i]
792  q, r = divmod(n, d)
793  while r == 0:
794  nb += 1
795  n = q
796  q, r = divmod(n, d)
797  if nb > 0:
798  l += [d] * nb
799  if n == 1:
800  return l
801 
802  assert d%nat32.PRIMORIAL32_TUPLE[4] == 1, \
803  (d, nat32.PRIMORIAL32_TUPLE[4], d%nat32.PRIMORIAL32_TUPLE[4])
804 
805  # Facteurs premiers ne tenant pas sur 16 bits
806  while True:
807  for m in nat32.MODS_PRIMORIAL_4_DIFF_TUPLE:
808  # Parcourt les restes modulo 210 premiers avec 210
809  nb = 0
810  q, r = divmod(n, d)
811  while r == 0:
812  nb += 1
813  n = q
814  q, r = divmod(n, d)
815  if nb > 0:
816  l += [d] * nb
817  if n == 1:
818  return l
819  d += m
820 
821 
822 ## Renvoie True si p est une séquence ordonnée de primes, False sinon
823 def primes_is(s):
824  """Renvoie True si p est une séquence ordonnée de primes,
825  False sinon
826 
827  Pre: s: quelconque
828 
829  Result: boolean"""
830  if isinstance(s, list) or isinstance(s, tuple):
831  prev_n = 1
832  for n in s:
833  if (not DSPython.natural_is(n)) or (prev_n > n) or (not natural.prime_is(n)):
834  return False
835  prev_n = n
836  return True
837  else:
838  return False
839 
840 
841 ## String représentant les primes
842 def primes_str(s, times='.'):
843  """Renvoie un string représentant les primes
844  (Si s est vide alors renvoie '1')
845 
846  Pre: s: séquence ordonnée de primes
847  times: string représentant l'opérateur de multiplication
848 
849  Result: string
850 
851  O(s) = ..."""
852  assert primes_is(s), s
853 
854  if len(s) > 0:
855  l = []
856  for i in s:
857  l.append(str(i))
858  return times.join(l)
859  else:
860  return '1'
861 
862 
863 ## Produit des primes
864 def primes_to_n(s):
865  """Renvoie le naturel ayant s pour liste des primes
866  (Si s est vide alors renvoie 1)
867 
868  Pre: s: séquence ordonnée de primes
869 
870  Result: naturel >= 1
871 
872  O(s) = ..."""
873  assert primes_is(s), s
874 
875  p = 1
876  for i in s:
877  p *= i
878  return p
879 
880 
881 ## Liste des primaries à partir des primes
882 def primes_to_primaries(s):
883  """Renvoie la liste des primaries à partir de la liste des primes
884  (Si s est vide alors renvoie [])
885 
886  Pre: s: séquence ordonnée de primes
887 
888  Result: séquence strictement ordonnée de primaries
889 
890  O(s) = ..."""
891  assert primes_is(s), s
892 
893  if len(s) > 0:
894  r = []
895  prev = None
896  nb = 1
897  for n in s:
898  if n == prev:
899  nb += 1
900  else:
901  if prev != None:
902  r.append((prev, nb))
903  nb = 1
904  prev = n
905  r.append((prev, nb))
906  return r
907  else:
908  return []
909 
910 
911 ## Somme des diviseurs à partir des primaries (fonction s)
912 def properdivisors_sum(s):
913  """Renvoie la somme des diviseurs propres du naturel
914  ayant s pour liste des primaries
915  == divisors_sum(s) - primaries_to_n(s)
916 
917  Pre: s: séquence strictement ordonnée de primaries
918 
919  Result: naturel >= 1
920 
921  O(s) = ..."""
922  assert primaries_is(s), s
923 
924  p = 1
925  n = 1
926  for i in s:
927  t = i[0]**i[1]
928  p *= t
929  n *= (t*i[0] - 1) // (i[0] - 1)
930  return n - p
931 
932 
933 ## Produit des facteurs premiers distincts à partir des primaries (fonction radical)
934 def rad(s):
935  """Pour s = [... (P_i, e_i) ...] où P_i est le ième nombre premier
936  renvoie [... (P_i, 1) ...]
937  Renvoie le produit des facteurs premiers distincts
938  du naturel ayant s pour liste des primaries
939  [OEIS A007947] Largest square-free number dividing n (the square-free kernel of n)
940 
941  Pre: s: séquence strictement ordonnée de primaries
942 
943  Result: séquence strictement ordonnée de primaries (naturel premier, 1)
944 
945  O(s) = ..."""
946  assert primaries_is(s), s
947 
948  l = []
949  for p in s:
950  l.append((p[0], 1))
951  return l
952 
953 
954 ## Renvoie True si p représente un nombre sans carré, False sinon
955 def squarefree_is(s):
956  """Renvoie True si p représente un nombre sans carré [OEIS A005117],
957  False sinon
958 
959  Pre: s: séquence strictement ordonnée de primaries
960 
961  Result: boolean"""
962  assert primaries_is(s), s
963 
964  for p in s:
965  if p[1] > 1:
966  return False
967  return True
968 
969 
970 ## Totient à partir des primaries
971 def totient(s):
972  """Renvoie le nombre de totatives du naturel
973  ayant s pour liste des primaries
974  (Si s == [] alors renvoie 1)
975 
976  Pre: s: séquence strictement ordonnée de primaries
977 
978  Result: naturel
979 
980  O(n) = ..."""
981  assert primaries_is(s), s
982 
983  nb = 1
984  for i in s:
985  nb *= (i[0] - 1) * i[0]**(i[1] - 1)
986  return nb
987 
988 
989 
990 # ######\cond MAINTEST
991 # Main #
992 ########
993 if __name__ == '__main__':
994  def main_test():
995  """Test du module"""
996  import sys
997 
998  try:
999  if not 'profile' in dir():
1000  import psyco
1001  psyco.full()
1002  except ImportError:
1003  pass
1004 
1005  import DSPython.debug as debug
1006 
1007  import DSPython.natseq as natseq
1008 
1009  debug.test_begin(VERSION, __debug__)
1010 
1011  print('bef()...', end=''); sys.stdout.flush()
1012  assert bef([]) == 0, bef([])
1013  assert bef([(2, 1)]) == 1, bef([(2, 1)])
1014  assert bef([(2, 3)]) == 3, bef([(2, 3)])
1015  assert bef([(5, 1)]) == 4, bef([(5, 1)])
1016  assert bef([(5, 3)]) == 12, bef([(5, 3)])
1017  assert bef([(2, 10), (5, 3)]) == 22, bef([(2, 10), (5, 3)])
1018  assert bef([(2, 10), (5, 3), (11, 2)]) == 54, bef([(2, 10), (5, 3), (11, 2)])
1019  for i in range(100):
1020  assert bef([(nat32.PRIME_16_ARRAY[i], 1)]) == 2**i, \
1021  (i, nat32.PRIME_16_ARRAY[i], bef([(nat32.PRIME_16_ARRAY[i], 1)]), 2**i)
1022  print('ok'); sys.stdout.flush()
1023 
1024 
1025  print('coprime_is()...', end=''); sys.stdout.flush()
1026  assert coprime_is([], []), coprime_is([], [])
1027  for m in range(1, 50 if debug.assertspeed >= debug.ASSERT_NORMAL else 30):
1028  f = primaries(m)
1029  assert coprime_is(f, []), (m, f, coprime_is(primaries(m), []))
1030  for p in nat32.PRIME_16_ARRAY[:1000 if debug.assertspeed >= debug.ASSERT_NORMAL else
1031  100]:
1032  if m%p == 0:
1033  assert not coprime_is(f, [(p, 1)]), (m, f, p, coprime_is(f, [(p, 1)]))
1034  else:
1035  assert coprime_is(f, [(p, 1)]), (m, f, p, coprime_is(f, [(p, 1)]))
1036  assert coprime_is(f, primaries(m + 1)), m
1037  for m in range(1, 50 if debug.assertspeed >= debug.ASSERT_NORMAL else 30):
1038  f_m = primaries(m)
1039  for n in range(1, 10):
1040  f_n = primaries(n)
1041  assert coprime_is(f_m, f_n) == coprime_is(f_n, f_m), \
1042  (m, n, coprime_is(f_m, f_n), gcd(f_n, f_m))
1043  if debug.assertspeed >= debug.ASSERT_NORMAL:
1044  for m in range(2**32 - 1, 2**32 + 1):
1045  assert coprime_is(primaries(m), primaries(m + 1)), m
1046  else:
1047  print(debug.assertspeed_str(), end='')
1048  print('ok'); sys.stdout.flush()
1049 
1050 
1051  print('distance_dominici()...', end=''); sys.stdout.flush()
1052  f = primaries(11)
1053  g = primaries(12)
1054  assert distance_dominici(f, g) == 4, distance_dominici(f, g)
1055  for a in range(1, 100):
1056  f = primaries(a)
1057  assert distance_dominici(f, f) == 0, distance_dominici(f, f)
1058  for b in range(1, 50):
1059  g = primaries(b)
1060  assert distance_dominici(f, g) == distance_dominici(g, f), \
1061  (distance_dominici(f, g), distance_dominici(g, f))
1062  print('???', end='')
1063  print('ok'); sys.stdout.flush()
1064 
1065 
1066  print('divisors_nb()...', end=''); sys.stdout.flush()
1067  assert divisors_nb([]) == 1, divisors_nb([])
1068  assert divisors_nb([(2, 1)]) == 2, divisors_nb([(2, 1)])
1069  assert divisors_nb([(2, 2)]) == 3, divisors_nb([(2, 2)])
1070  assert divisors_nb([(2, 1), (3, 1)]) == 4, divisors_nb([(2, 1), (3, 1)])
1071  assert divisors_nb([(2, 5), (3, 2)]) == 18, divisors_nb([(2, 5), (3, 2)])
1072 
1073  for n in range(1, 1000):
1074  f = primaries(n)
1075  d = natural.divisors(n)
1076  assert divisors_nb(f) == len(d), (n, divisors_nb(f), len(d), f, d)
1077  print('ok'); sys.stdout.flush()
1078 
1079 
1080  print('divisors_prod()...', end=''); sys.stdout.flush()
1081  assert divisors_prod([]) == 1, divisors_prod([])
1082  assert divisors_prod([(2, 1)]) == 2, divisors_prod([(2, 1)])
1083  assert divisors_prod([(2, 2)]) == 8, divisors_prod([(2, 2)])
1084  assert divisors_prod([(2, 1), (3, 1)]) == 36, divisors_prod([(2, 1), (3, 1)])
1085  assert divisors_prod([(2, 5), (3, 2)]) == 13631146639813244878848, \
1086  divisors_prod([(2, 5), (3, 2)])
1087 
1088  for n in range(1, 1000):
1089  f = primaries(n)
1090  d = natural.divisors(n)
1091  assert divisors_prod(f) == natseq.prod(d), (n, divisors_prod(f), natseq.prod(d), f, d)
1092  print('ok'); sys.stdout.flush()
1093 
1094 
1095  print('divisors_sum()...', end=''); sys.stdout.flush()
1096  assert divisors_sum([]) == 1, divisors_sum([])
1097  assert divisors_sum([(2, 1)]) == 3, divisors_sum([(2, 1)])
1098  assert divisors_sum([(2, 2)]) == 7, divisors_sum([(2, 2)])
1099  assert divisors_sum([(2, 1), (3, 1)]) == 12, divisors_sum([(2, 1), (3, 1)])
1100  assert divisors_sum([(2, 5), (3, 2)]) == 63*13, divisors_sum([(2, 5), (3, 2)])
1101 
1102  for n in range(1, 1000):
1103  f = primaries(n)
1104  d = natural.divisors(n)
1105  assert divisors_sum(f) == natseq.sum(d), (n, divisors_sum(f), natseq.sum(d), f, d)
1106  print('ok'); sys.stdout.flush()
1107 
1108 
1109  print('divisors_sum_odd()...', end=''); sys.stdout.flush()
1110  assert divisors_sum_odd([]) == 1, divisors_sum_odd([])
1111  assert divisors_sum_odd([(2, 1)]) == 1, divisors_sum_odd([(2, 1)])
1112  assert divisors_sum_odd([(2, 2)]) == 1, divisors_sum_odd([(2, 2)])
1113  assert divisors_sum_odd([(2, 1), (3, 1)]) == 4, divisors_sum_odd([(2, 1), (3, 1)])
1114  assert divisors_sum_odd([(2, 5), (3, 2)]) == 13, divisors_sum_odd([(2, 5), (3, 2)])
1115 
1116  for n in range(1, 1000):
1117  f = primaries(n)
1118  d = natural.divisors_cond(n, lambda d: d&1 != 0)
1119  assert divisors_sum_odd(f) == natseq.sum(d), \
1120  (n, divisors_sum_odd(f), natseq.sum(d), f, d)
1121  print('ok'); sys.stdout.flush()
1122 
1123 
1124  print('divisors_sum_even()...', end=''); sys.stdout.flush()
1125  assert divisors_sum_even([]) == 0, divisors_sum_even([])
1126  assert divisors_sum_even([(2, 1)]) == 2, divisors_sum_even([(2, 1)])
1127  assert divisors_sum_even([(2, 2)]) == 6, divisors_sum_even([(2, 2)])
1128  assert divisors_sum_even([(2, 1), (3, 1)]) == 8, divisors_sum_even([(2, 1), (3, 1)])
1129  assert divisors_sum_even([(2, 5), (3, 2)]) == 62*13, divisors_sum_even([(2, 5), (3, 2)])
1130 
1131  for n in range(1, 1000):
1132  f = primaries(n)
1133  d = natural.divisors_cond(n, lambda d: d&1 == 0)
1134  assert divisors_sum_even(f) == natseq.sum(d), \
1135  (n, divisors_sum_even(f), natseq.sum(d), f, d)
1136  assert divisors_sum_even(f) + divisors_sum_odd(f) == divisors_sum(f), \
1137  (n, d, divisors_sum_even(f), divisors_sum_odd(f), divisors_sum(f), f)
1138  print('ok'); sys.stdout.flush()
1139 
1140 
1141  print('divisors_sum_pow()...', end=''); sys.stdout.flush()
1142  assert divisors_sum_pow([], 0) == 1, divisors_sum_pow([], 0)
1143  assert divisors_sum_pow([(2, 1)], 0) == 2, divisors_sum_pow([(2, 1)], 0)
1144  assert divisors_sum_pow([(2, 2)], 0) == 3, divisors_sum_pow([(2, 2)], 0)
1145  assert divisors_sum_pow([(2, 1), (3, 1)], 0) == 4, divisors_sum_pow([(2, 1), (3, 1)], 0)
1146  assert divisors_sum_pow([(2, 5), (3, 2)], 0) == 18, divisors_sum_pow([(2, 5), (3, 2)], 0)
1147 
1148  assert divisors_sum_pow([], 1) == 1, divisors_sum_pow([], 1)
1149  assert divisors_sum_pow([(2, 1)], 1) == 3, divisors_sum_pow([(2, 1)], 1)
1150  assert divisors_sum_pow([(2, 2)], 1) == 7, divisors_sum_pow([(2, 2)], 1)
1151  assert divisors_sum_pow([(2, 1), (3, 1)], 1) == 12, divisors_sum_pow([(2, 1), (3, 1)], 1)
1152  assert divisors_sum_pow([(2, 5), (3, 2)], 1) == 63*13, divisors_sum_pow([(2, 5), (3, 2)], 1)
1153 
1154  assert divisors_sum_pow([], 2) == 1, divisors_sum_pow([], 2)
1155  assert divisors_sum_pow([(2, 1)], 2) == 5, divisors_sum_pow([(2, 1)], 2)
1156  assert divisors_sum_pow([(2, 2)], 2) == 21, divisors_sum_pow([(2, 2)], 2)
1157  assert divisors_sum_pow([(2, 1), (3, 1)], 2) == 50, divisors_sum_pow([(2, 1), (3, 1)], 2)
1158  assert divisors_sum_pow([(2, 5), (3, 2)], 2) == 124215, \
1159  divisors_sum_pow([(2, 5), (3, 2)], 2)
1160 
1161  for n in range(1, 1000):
1162  f = primaries(n)
1163  d = natural.divisors(n)
1164  assert divisors_sum_pow(f, 0) == len(d), \
1165  (n, divisors_sum_pow(f, 0), len(d), f, d)
1166  assert divisors_sum_pow(f, 1) == natseq.sum(d), \
1167  (n, divisors_sum_pow(f, 1), natseq.sum(d), f, d)
1168  for k in range(10):
1169  assert divisors_sum_pow(f, k) == natseq.sum_pow(d, k), \
1170  (n, k, divisors_sum_pow(f, k), natseq.sum_pow(d, k), f, d)
1171  print('ok'); sys.stdout.flush()
1172 
1173 
1174  print('feb_primaries()...', end=''); sys.stdout.flush()
1175  assert feb_primaries(0) == [], feb_primaries(0)
1176  assert feb_primaries(1) == [(2, 1)], feb_primaries(1)
1177  assert feb_primaries(2) == [(3, 1)], feb_primaries(2)
1178  assert feb_primaries(3) == [(2, 1), (3, 1)], feb_primaries(3)
1179  assert feb_primaries(4) == [(5, 1)], feb_primaries(4)
1180  assert feb_primaries(5) == [(2, 1), (5, 1)], feb_primaries(5)
1181  assert feb_primaries(6) == [(3, 1), (5, 1)], feb_primaries(6)
1182  assert feb_primaries(7) == [(2, 1), (3, 1), (5, 1)], feb_primaries(7)
1183  for i in range(100):
1184  assert feb_primaries(2**i) == [(nat32.PRIME_16_ARRAY[i], 1)], (i, feb_primaries(2**i))
1185  for j in range(100):
1186  if i != j:
1187  assert feb_primaries(2**i + 2**j) == \
1188  [(nat32.PRIME_16_ARRAY[min(i, j)], 1),
1189  (nat32.PRIME_16_ARRAY[max(i, j)], 1)], \
1190  (i, j, feb_primaries(2**i + 2**j))
1191  for n in range(1000):
1192  for p in feb_primaries(n):
1193  assert nat32.prime_is(p[0]), (n, p)
1194  assert p[1] == 1, (n, p)
1195  assert bef(feb_primaries(n)) == n , (n, bef(feb_primaries(n)), feb_primaries(n))
1196  print('ok'); sys.stdout.flush()
1197 
1198 
1199  print('gcd()...', end=''); sys.stdout.flush()
1200  assert gcd([], []) == [], gcd([], [])
1201  for m in range(1, 1000 if debug.assertspeed >= debug.ASSERT_NORMAL else 100):
1202  f = primaries(m)
1203  assert gcd(f, []) == [], (m, gcd(f, []), f)
1204  assert gcd([], f) == [], (m, gcd([], f), f)
1205  assert gcd(f, f) == f, (m, gcd(f, f), f)
1206  for p in nat32.PRIME_16_ARRAY[:200]:
1207  if m%p == 0:
1208  assert gcd(f, [(p, 1)]) == [(p, 1)], (m, p, gcd(f, [(p, 1)]), f)
1209  else:
1210  assert gcd(f, [(p, 1)]) == [], (m, p, gcd(f, [(p, 1)]), f)
1211  for m in range(1, 100):
1212  m_f = primaries(m)
1213  for n in range(1, 100):
1214  n_f = primaries(n)
1215  d_f = primaries(fractions.gcd(m, n))
1216  assert gcd(m_f, n_f) == d_f, (m, n, gcd(m_f, n_f),
1217  m_f, n_f, fractions.gcd(m, n), d_f)
1218  assert gcd(n_f, m_f) == d_f, (m, n, gcd(n_f, m_f),
1219  m_f, n_f, fractions.gcd(m, n), d_f)
1220  if debug.assertspeed < debug.ASSERT_NORMAL:
1221  print(debug.assertspeed_str(), end='')
1222  print('ok'); sys.stdout.flush()
1223 
1224 
1225  print('godelnumber()...', end=''); sys.stdout.flush()
1226  assert godelnumber([]) == 1, godelnumber([])
1227  assert godelnumber(()) == 1, godelnumber(())
1228 
1229  for a in range(10):
1230  assert godelnumber([a]) == 2**a, (a, godelnumber([a]), 2**a)
1231  assert godelnumber((a, )) == 2**a, (a, godelnumber((a, )), 2**a)
1232  for a in range(10):
1233  for b in range(10):
1234  assert godelnumber((a, b)) == 2**a * 3**b, (a, b, godelnumber((a, b)), 2**a * 3**b)
1235  p = 1
1236  for n in range(10):
1237  assert godelnumber((1, )*n) == p, (n, godelnumber((1, )*n), p)
1238  p *= nat32.PRIME_16_ARRAY[n]
1239  for a in range(10):
1240  assert godelnumber((0, )*n + (a, )) == nat32.PRIME_16_ARRAY[n]**a, \
1241  (n, a, godelnumber((0, )*n + (a, )),
1242  nat32.PRIME_16_ARRAY[n]**a, (0, )*n + (a, ))
1243  print('ok'); sys.stdout.flush()
1244 
1245 
1246  print('godelnumber_to_list()...', end=''); sys.stdout.flush()
1247  assert godelnumber_to_list(1) == [], godelnumber_to_list(1)
1248  assert godelnumber_to_list(2) == [1], godelnumber_to_list(2)
1249  assert godelnumber_to_list(3) == [0, 1], godelnumber_to_list(3)
1250  assert godelnumber_to_list(4) == [2], godelnumber_to_list(4)
1251  assert godelnumber_to_list(5) == [0, 0, 1], godelnumber_to_list(5)
1252  for n in range(1, 1000):
1253  for k in range(10):
1254  assert godelnumber(godelnumber_to_list(n) + [0]*k) == n, \
1255  (n, k, godelnumber_to_list(n) + [0]*k,
1256  godelnumber(godelnumber_to_list(n) + [0]*k))
1257  print('ok'); sys.stdout.flush()
1258 
1259 
1260  print('lcm()...', end=''); sys.stdout.flush()
1261  assert lcm([], []) == [], lcm([], [])
1262  for m in range(1, 1000 if debug.assertspeed >= debug.ASSERT_NORMAL else 100):
1263  f = primaries(m)
1264  assert lcm(f, []) == f, (m, lcm(f, []), f)
1265  assert lcm([], f) == f, (m, lcm([], f), f)
1266  assert lcm(f, f) == f, (m, lcm(f, f), f)
1267  for p in nat32.PRIME_16_ARRAY[:30]:
1268  if m%p == 0:
1269  assert lcm(f, [(p, 1)]) == f, (m, p, lcm(f, [(p, 1)]), f)
1270  else:
1271  assert lcm(f, [(p, 1)]) == primaries(m * p), (m, p, lcm(f, [(p, 1)]), f)
1272  for m in range(1, 100):
1273  m_f = primaries(m)
1274  for n in range(1, 100):
1275  n_f = primaries(n)
1276  d_f = primaries(nat32.lcm(m, n))
1277  assert lcm(m_f, n_f) == d_f, (m, n, lcm(m_f, n_f), m_f, n_f, nat32.lcm(m, n), d_f)
1278  assert lcm(n_f, m_f) == d_f, (m, n, lcm(n_f, m_f), m_f, n_f, nat32.lcm(m, n), d_f)
1279  if debug.assertspeed < debug.ASSERT_NORMAL:
1280  print(debug.assertspeed_str(), end='')
1281  print('ok'); sys.stdout.flush()
1282 
1283 
1284  print('mobius()...', end=''); sys.stdout.flush()
1285  assert mobius([]) == 1, mobius([])
1286  assert mobius([(2, 1)]) == -1, mobius([(2, 1)])
1287  assert mobius([(3, 1)]) == -1, mobius([(3, 1)])
1288  assert mobius([(2, 2)]) == 0, mobius([(2, 2)])
1289  assert mobius([(5, 1)]) == -1, mobius([(5, 1)])
1290  assert mobius([(2, 1), (5, 1)]) == 1, mobius([(2, 1), (5, 1)])
1291  for p in nat32.PRIME_16_ARRAY[2:]:
1292  assert mobius([(p, 1)]) == -1, (p, mobius([(p, 1)]))
1293  assert mobius([(p, 2)]) == 0, (p, mobius([(p, 2)]))
1294  assert mobius([(2, 1), (p, 1)]) == 1, (p, mobius([(2, 1), (p, 1)]))
1295  assert mobius([(2, 1), (3, 1), (p, 1)]) == -1, (p, mobius([(2, 1), (3, 1), (p, 1)]))
1296  print('ok'); sys.stdout.flush()
1297 
1298 
1299  print('mul()...', end=''); sys.stdout.flush()
1300  assert mul([], []) == [], mul([], [])
1301  for n in range(1, 1000):
1302  f = primaries(m)
1303  assert mul(f, []) == f, (m, mul(f, []), f)
1304  assert mul([], f) == f, (m, mul([], f), f)
1305  assert mul(f, f) == primaries(m*m), (m, mul(f, f), primaries(m*m))
1306  for m in range(1, 100):
1307  m_f = primaries(m)
1308  for n in range(1, 100):
1309  n_f = primaries(n)
1310  d_f = primaries(m*n)
1311  assert mul(m_f, n_f) == d_f, (m, n, mul(m_f, n_f), m_f, n_f, m*n, d_f)
1312  assert mul(n_f, m_f) == d_f, (m, n, mul(n_f, m_f), m_f, n_f, m*n, d_f)
1313  print('ok'); sys.stdout.flush()
1314 
1315 
1316  print('nb_in_integers_4sqr()...', end=''); sys.stdout.flush()
1317  assert nb_in_integers_4sqr([]) == 8, nb_in_integers_4sqr([])
1318  assert nb_in_integers_4sqr([(2, 1)]) == 24, nb_in_integers_4sqr([(2, 1)])
1319  assert nb_in_integers_4sqr([(3, 1)]) == 32, nb_in_integers_4sqr([(3, 1)])
1320  assert nb_in_integers_4sqr([(2, 2)]) == 24, nb_in_integers_4sqr([(2, 2)])
1321  assert nb_in_integers_4sqr([(5, 1)]) == 48, nb_in_integers_4sqr([(5, 1)])
1322  for n in range(1, 1000):
1323  f = primaries(n)
1324  if n&1 == 0:
1325  assert nb_in_integers_4sqr(f) == 24 * divisors_sum(f[1:]), \
1326  (n, nb_in_integers_4sqr(f), divisors_sum(f[1:]), f)
1327  else:
1328  assert nb_in_integers_4sqr(f) == 8 * divisors_sum(f), \
1329  (n, nb_in_integers_4sqr(f), divisors_sum(f), f)
1330  print('ok'); sys.stdout.flush()
1331 
1332 
1333  print('nfp()...', end=''); sys.stdout.flush()
1334  assert nfp([], []) == [], nfp([], [])
1335  for n in range(1, 1000 if debug.assertspeed >= debug.ASSERT_NORMAL else 100):
1336  f = primaries(m)
1337  assert nfp(f, []) == f, (m, nfp(f, []), f)
1338  assert nfp([], f) == f, (m, nfp([], f), f)
1339  assert nfp(f, f) == [], (m, nfp(f, f), f)
1340  for p in nat32.PRIME_16_ARRAY[:50]:
1341  if m%p == 0:
1342  assert nfp(f, [(p, 1)]) == primaries(m//p), (m, p, nfp(f, [(p, 1)]), f)
1343  else:
1344  assert nfp(f, [(p, 1)]) == primaries(m*p), (m, p, nfp(f, [(p, 1)]), f)
1345  for m in range(1, 100):
1346  m_f = primaries(m)
1347  for n in range(1, 100):
1348  n_f = primaries(n)
1349  d_f = primaries(nat32.nfp(m, n))
1350  assert nfp(m_f, n_f) == d_f, (m, n, nfp(m_f, n_f), m_f, n_f, nat32.nfp(m, n), d_f)
1351  assert nfp(n_f, m_f) == d_f, (m, n, nfp(n_f, m_f), m_f, n_f, nat32.nfp(m, n), d_f)
1352  if debug.assertspeed < debug.ASSERT_NORMAL:
1353  print(debug.assertspeed_str(), end='')
1354  print('ok'); sys.stdout.flush()
1355 
1356 
1357  print('nontotient()...', end=''); sys.stdout.flush()
1358  assert nontotient([]) == 0, nontotient([])
1359  assert nontotient([(2, 1)]) == 1, nontotient([(2, 1)])
1360  assert nontotient([(3, 1)]) == 1, nontotient([(3, 1)])
1361  assert nontotient([(2, 2)]) == 2, nontotient([(2, 2)])
1362  assert nontotient([(5, 1)]) == 1, nontotient([(5, 1)])
1363  assert nontotient([(2, 1), (3, 1)]) == 4, nontotient([(2, 1), (3, 1)])
1364  assert nontotient([(2, 3), (3, 1)]) == 16, nontotient([(2, 3), (3, 1)])
1365  for n in range(1, 500 if debug.assertspeed >= debug.ASSERT_NORMAL else 100):
1366  assert nontotient(primaries(n)) == len(natural.nontotatives(n)), \
1367  (n, nontotient(primaries(n)), len(natural.nontotatives(n)),
1368  natural.nontotatives(n))
1369  for p in nat32.PRIME_16_ARRAY:
1370  assert nontotient([(p, 1)]) == 1, (p, nontotient([(p, 1)]))
1371  if debug.assertspeed < debug.ASSERT_NORMAL:
1372  print(debug.assertspeed_str(), end='')
1373  print('ok'); sys.stdout.flush()
1374 
1375 
1376  print('ope_exp()...', end=''); sys.stdout.flush()
1377  assert ope_exp([], [], lambda a, b: a + b) == [], ope_exp([], [], lambda a, b: a + b)
1378  for m in range(1, 200):
1379  m_f = primaries(m)
1380  for n in range(1, 200 if debug.assertspeed >= debug.ASSERT_NORMAL else 20):
1381  n_f = primaries(n)
1382  r = ope_exp(m_f, n_f, lambda a, b: a + b)
1383  assert r == mul(m_f, n_f), (m, n, r, mul(m_f, n_f))
1384  if m%n == 0:
1385  r = ope_exp(m_f, n_f, lambda a, b: a - b)
1386  assert r == primaries(m//n), (m, n, r, primaries(m//n))
1387  r = ope_exp(m_f, n_f, lambda a, b: min(a, b))
1388  assert r == gcd(m_f, n_f), (m, n, r, gcd(m_f, n_f))
1389  r = ope_exp(m_f, n_f, lambda a, b: max(a, b))
1390  assert r == lcm(m_f, n_f), (m, n, r, lcm(m_f, n_f))
1391  r = ope_exp(m_f, n_f, lambda a, b: abs(a - b))
1392  assert r == nfp(m_f, n_f), (m, n, r, nfp(m_f, n_f))
1393  if debug.assertspeed < debug.ASSERT_NORMAL:
1394  print(debug.assertspeed_str(), end='')
1395  print('ok'); sys.stdout.flush()
1396 
1397 
1398  print('primaries()...', end=''); sys.stdout.flush()
1399  assert primaries(1) == [], primaries(1)
1400  assert primaries(2) == [(2, 1)], primaries(2)
1401  assert primaries(3) == [(3, 1)], primaries(3)
1402  assert primaries(4) == [(2, 2)], primaries(4)
1403  assert primaries(5) == [(5, 1)], primaries(5)
1404  assert primaries(6) == [(2, 1), (3, 1)], primaries(6)
1405 
1406  q = nat32.PRIME_16_ARRAY[100]
1407  for p in nat32.PRIME_16_ARRAY[:100]:
1408  for k in range(1, 10):
1409  assert primaries(p**k) == [(p, k)], (p, k, primaries(p**k))
1410  assert primaries(p**k * q**(k + 1)) == [(p, k), (q, k + 1)], \
1411  (p, k, primaries(p**k * q**(k + 1)))
1412 
1413  for n in range(1, 1000):
1414  f = primaries(n)
1415  for p in f:
1416  assert nat32.prime_is(p[0]), (n, p)
1417  assert DSPython.natural_is(p[1]), (n, p)
1418  assert p[1] >= 1, (n, p[1])
1419  print('ok'); sys.stdout.flush()
1420 
1421 
1422  print('primaries_is()...', end=''); sys.stdout.flush()
1423  assert primaries_is([])
1424  assert primaries_is(())
1425  assert primaries_is([(2, 1)])
1426  assert not primaries_is([(2, 0)])
1427  assert not primaries_is([(1, 1)])
1428  assert primaries_is([(2, 1), (3, 1)])
1429  assert not primaries_is([(2, 1), (2, 1)])
1430  assert primaries_is([(2, 4), (3, 1)])
1431  assert primaries_is([(2, 24), (3, 19), (17, 3)])
1432  for n in range(1, 1000):
1433  assert primaries_is(primaries(n)), (n, primaries(n))
1434  print('ok'); sys.stdout.flush()
1435 
1436 
1437  print('primaries_prime_is()...', end=''); sys.stdout.flush()
1438  for p in nat32.PRIME_16_ARRAY[:1000]:
1439  f = primaries(p)
1440  assert primaries_prime_is(f), (p, f)
1441  for n in range(1, 100):
1442  f = primaries(n)
1443  assert primaries_prime_is(f) == nat32.prime_is(n), \
1444  (n, primaries_prime_is(f), nat32.prime_is(n))
1445  print('ok'); sys.stdout.flush()
1446 
1447 
1448  print('primaries_str()...', end=''); sys.stdout.flush()
1449  assert primaries_str([]) == '1', primaries_str([])
1450  assert primaries_str([(3, 1)]) == '3', primaries_str([(3, 1)])
1451  assert primaries_str([(3, 2)]) == '3^2', primaries_str([(3, 2)])
1452  assert primaries_str([(3, 1), (5, 7)]) == '3.5^7', primaries_str([(3, 1), (5, 7)])
1453  assert primaries_str([(3, 2), (5, 1)]) == '3^2.5', primaries_str([(3, 2), (5, 1)])
1454  assert primaries_str([(3, 2), (5, 7)]) == '3^2.5^7', primaries_str([(3, 2), (5, 7)])
1455 
1456  assert primaries_str([(3, 2), (5, 7)], times=' x ', exp='**') == '3**2 x 5**7', \
1457  primaries_str([(3, 2), (5, 7)])
1458  print('ok'); sys.stdout.flush()
1459 
1460 
1461  print('primaries_to_n()...', end=''); sys.stdout.flush()
1462  assert primaries_to_n([]) == 1, primaries([])
1463  assert primaries_to_n([(2, 1)]) == 2, primaries([(2, 1)])
1464  assert primaries_to_n([(2, 2)]) == 4, primaries([(2, 2)])
1465  assert primaries_to_n([(2, 3), (5, 1)]) == 40, primaries([(2, 3), (5, 1)])
1466 
1467  for n in range(1, 1000):
1468  assert primaries_to_n(primaries(n)) == n, \
1469  (n, primaries_to_n(primaries(n)), primaries(n))
1470  if debug.assertspeed >= debug.ASSERT_NORMAL:
1471  for n in range(2**16 - 50, 2**16 + 100):
1472  assert primaries_to_n(primaries(n)) == n, \
1473  (n, primaries_to_n(primaries(n)), primaries(n))
1474  for n in range(2**24 - 5, 2**24 + 10):
1475  assert primaries_to_n(primaries(n)) == n, \
1476  (n, primaries_to_n(primaries(n)), primaries(n))
1477  else:
1478  print(debug.assertspeed_str(), end='')
1479  print('ok'); sys.stdout.flush()
1480 
1481 
1482  print('primaries_to_primes()...', end=''); sys.stdout.flush()
1483  assert primaries_to_primes([]) == [], primaries_to_primes([])
1484  assert primaries_to_primes([(2, 1)]) == [2], primaries_to_primes([(2, 1)])
1485  assert primaries_to_primes([(2, 2)]) == [2, 2], primaries_to_primes([(2, 2)])
1486  assert primaries_to_primes([(2, 2), (3, 1)]) == [2, 2, 3], \
1487  primaries_to_primes([(2, 2), (3, 1)])
1488  for n in range(1, 1000):
1489  assert primaries_to_primes(primaries(n)) == primes(n), \
1490  (n, primaries_to_primes(primaries(n)), primes(n), primaries(n))
1491  print('ok'); sys.stdout.flush()
1492 
1493 
1494  print('primes()...', end=''); sys.stdout.flush()
1495  assert primes(1) == [], primes(1)
1496  assert primes(2) == [2], primes(2)
1497  assert primes(3) == [3], primes(3)
1498  assert primes(4) == [2, 2], primes(4)
1499  assert primes(5) == [5], primes(5)
1500  assert primes(6) == [2, 3], primes(6)
1501 
1502  q = nat32.PRIME_16_ARRAY[100]
1503  for p in nat32.PRIME_16_ARRAY[:100]:
1504  for k in range(1, 10):
1505  assert primes(p**k) == [p]*k, (p, k, primes(p**k))
1506  assert primes(p**k * q**(k + 1)) == [p]*k + [q]*(k + 1), \
1507  (p, k, primes(p**k * q**(k + 1)))
1508 
1509  for n in range(1, 1000):
1510  f = primes(n)
1511  for p in f:
1512  assert nat32.prime_is(p), (n, p)
1513 
1514  fp = primaries(n)
1515  nb = 0
1516  for p in fp:
1517  nb += p[1]
1518  assert len(f) == nb, (n, len(f), nb)
1519  print('ok'); sys.stdout.flush()
1520 
1521 
1522  print('primes_is()...', end=''); sys.stdout.flush()
1523  assert primes_is([])
1524  assert primes_is(())
1525  assert primes_is([2])
1526  assert not primes_is([0])
1527  assert not primes_is([1])
1528  assert primes_is([2, 3])
1529  assert primes_is([2, 2])
1530  assert primes_is([2, 3])
1531  assert primes_is([2, 3, 3])
1532  assert not primes_is([2, 3, 3, 6])
1533  for n in range(1, 1000):
1534  assert primes_is(primes(n)), (n, primes(n))
1535  print('ok'); sys.stdout.flush()
1536 
1537 
1538  print('primes_str()...', end=''); sys.stdout.flush()
1539  assert primes_str([]) == '1', primes_str([])
1540  assert primes_str([2]) == '2', primes_str([2])
1541  assert primes_str([2, 2]) == '2.2', primes_str([2, 2])
1542  assert primes_str([2, 3]) == '2.3', primes_str([2, 3])
1543 
1544  assert primes_str([], times=' x ') == '1', primes_str([], times=' x ')
1545  assert primes_str([2], times=' x ') == '2', primes_str([2], times=' x ')
1546  assert primes_str([2, 2], times=' x ') == '2 x 2', primes_str([2, 2], times=' x ')
1547  assert primes_str([2, 3], times=' x ') == '2 x 3', primes_str([2, 3], times=' x ')
1548  print('ok'); sys.stdout.flush()
1549 
1550 
1551  print('primes_to_n()...', end=''); sys.stdout.flush()
1552  assert primes_to_n([]) == 1, primaries([])
1553  assert primes_to_n([2]) == 2, primaries([0])
1554  assert primes_to_n([2, 2]) == 4, primaries([1])
1555  assert primes_to_n([2, 2, 2, 5]) == 40, primaries([2])
1556 
1557  for n in range(1, 1000):
1558  assert primes_to_n(primes(n)) == n, (n, primes_to_n(primes(n)), primes(n))
1559  if debug.assertspeed >= debug.ASSERT_NORMAL:
1560  for n in range(2**16 - 50, 2**16 + 100):
1561  assert primes_to_n(primes(n)) == n, (n, primes_to_n(primes(n)), primes(n))
1562  for n in range(2**24 - 5, 2**24 + 10):
1563  assert primes_to_n(primes(n)) == n, (n, primes_to_n(primes(n)), primes(n))
1564  else:
1565  print(debug.assertspeed_str(), end='')
1566  print('ok'); sys.stdout.flush()
1567 
1568 
1569  print('primes_to_primaries()...', end=''); sys.stdout.flush()
1570  assert primes_to_primaries([]) == [], primes_to_primaries([])
1571  assert primes_to_primaries([2]) == [(2, 1)], primes_to_primaries([2])
1572  assert primes_to_primaries([2, 2]) == [(2, 2)], primes_to_primaries([2, 2])
1573  assert primes_to_primaries([2, 2, 3]) == [(2, 2), (3, 1)], primes_to_primaries([2, 2, 3])
1574  for n in range(1, 1000):
1575  assert primes_to_primaries(primes(n)) == primaries(n), \
1576  (n, primes_to_primaries(primes(n)), primaries(n), primes(n))
1577  print('ok'); sys.stdout.flush()
1578 
1579 
1580  print('properdivisors_sum()...', end=''); sys.stdout.flush()
1581  assert properdivisors_sum([]) == 0, properdivisors_sum([])
1582  assert properdivisors_sum([(2, 1)]) == 1, properdivisors_sum([(2, 1)])
1583  assert properdivisors_sum([(2, 2)]) == 3, properdivisors_sum([(2, 2)])
1584  assert properdivisors_sum([(2, 1), (3, 1)]) == 6, \
1585  properdivisors_sum([(2, 1), (3, 1)])
1586  assert properdivisors_sum([(2, 5), (3, 2)]) == 63*13 - 32*9, \
1587  properdivisors_sum([(2, 5), (3, 2)])
1588 
1589  for n in range(1, 1000):
1590  f = primaries(n)
1591  d = natural.divisors(n)[:-1]
1592  assert properdivisors_sum(f) == natseq.sum(d), \
1593  (n, properdivisors_sum(f), natseq.sum(d), f, d)
1594  print('ok'); sys.stdout.flush()
1595 
1596 
1597  print('rad()...', end=''); sys.stdout.flush()
1598  assert rad([]) == [], rad([])
1599  for p in nat32.PRIME_16_ARRAY[:100]:
1600  for k in range(1, 10):
1601  assert rad([(p, k)]) == [(p, 1)], rad([(p, k)])
1602  for n in range(1, 100):
1603  f = primaries(n)
1604  r = rad(f)
1605  assert primaries_is(r), (n, r)
1606  assert len(f) == len(r), (n, len(f), len(r), f, r)
1607  assert squarefree_is(r), (n, r, f)
1608  print('ok'); sys.stdout.flush()
1609 
1610 
1611  print('squarefree_is()...', end=''); sys.stdout.flush()
1612  assert squarefree_is([])
1613  assert squarefree_is([(2, 1)])
1614  assert squarefree_is([(3, 1)])
1615  assert squarefree_is([(5, 1)])
1616  assert not squarefree_is([(2, 2)])
1617  assert not squarefree_is([(2, 33)])
1618  assert not squarefree_is([(5, 2)])
1619  assert not squarefree_is([(5, 33)])
1620  for n in range(100):
1621  assert squarefree_is(feb_primaries(n)), (n, feb_primaries(n))
1622  assert squarefree_is([(nat32.PRIME_16_ARRAY[n], 1)]), (n, nat32.PRIME_16_ARRAY[n])
1623  for k in range(2, 10):
1624  assert not squarefree_is([(nat32.PRIME_16_ARRAY[n], k)]), \
1625  (n, k, nat32.PRIME_16_ARRAY[n])
1626  print('ok'); sys.stdout.flush()
1627 
1628 
1629  print('totient()...', end=''); sys.stdout.flush()
1630  assert totient([]) == 1, totient([])
1631  assert totient([(2, 1)]) == 1, totient([(2, 1)])
1632  assert totient([(3, 1)]) == 2, totient([(3, 1)])
1633  assert totient([(2, 2)]) == 2, totient([(2, 2)])
1634  assert totient([(5, 1)]) == 4, totient([(5, 1)])
1635  assert totient([(2, 1), (3, 1)]) == 2, totient([(2, 1), (3, 1)])
1636  assert totient([(2, 3), (3, 1)]) == 8, totient([(2, 3), (3, 1)])
1637  for n in range(1, 500 if debug.assertspeed >= debug.ASSERT_NORMAL else 100):
1638  f = primaries(n)
1639  t = totient(f)
1640  assert t== len(natural.totatives(n)), \
1641  (n, t, len(natural.totatives(n)), natural.totatives(n), f)
1642  assert t + nontotient(f) == n, (n, t, nontotient(n), f)
1643  for p in nat32.PRIME_16_ARRAY:
1644  for k in range(1, 10):
1645  assert totient([(p, k)]) == (p - 1) * p**(k - 1), (p, k, totient([(p, k)]))
1646  if debug.assertspeed < debug.ASSERT_NORMAL:
1647  print(debug.assertspeed_str(), end='')
1648  print('ok'); sys.stdout.flush()
1649  debug.test_end()
1650 
1651  main_test()
1652 ##\endcond MAINTEST