Page 1 of 3

quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 15:22
by loupiot
Voici un petit script bien pratique pour l'arithmétique si tu es au lycée :)
il contient une fonction pgcd, ppcm, diviseurs (renvoyant les diviseurs d'un nombre), quotient (liste de tous les quotients de l'algorithme d'Euclide), un test de primalité, une fonction renvoyant un couple (u,v) résolvant au+bv=1 et une fonction résolvant toutes les équations diophantiennes linéaires résolvables (celles vu au lycée), donc bien pratique pour la terminale S ;)
il s'appelle arithmetique.py
Code: Select all
from math import*

def isprime(n):
    if type(n)!=int or n<2: #test pour vérifier si n correspond bien à ce qu'on attend de lui. Il y en a beaucoup dans le script, souvent inutile mais au cas où :p
        return False
    for i in range(2,int(n**0.5)+1):
        if n%i==0:
           return False
    return True

def diviseurs(n):
    if type(n)!=int or n<1:
        return None
    if n==1:
        return [1]
    liste=[1]
    for i in range(2,int(n**0.5)+1):
        if n%i==0:
            liste.append(i)
            if n//i!=i: #pour tester si i²!=n, sans ce test pour un carré on aurait deux fois la racine
                liste.append(n//i)
    return liste+[n]

def pgcd(a,b):
    if type(a+b)!=int:
        return None
    while b!=0:
        a,b=b,a%b
    return abs(a)

def ppcm(a,b):
    if type(a+b)!=int:
        return None
    return a*b//pgcd(a,b)

def quotient(a,b):
    if type(a+b)!=int:
        return None
    liste=[]
    while b!=0:
        liste.append(a//b)
        a,b=b,a%b
    return liste

def bezout(a,b):
    if pgcd(a,b)!=1:
        return None
    if abs(a)<abs(b):
        a,b=b,a #u multiplie le max(a,b), ceci m'évite donc de faire cas par cas...
    liste_q=quotient(abs(a),abs(b))
    liste_u=[0]*len(liste_q)
    del liste_q[len(liste_q)-1] #dans le prochain algorithme pour trouver u, on ne veut pas le quotient du reste nul
    liste_u[1]=1
    u=0
    c=1
    for i in liste_q[1:]:
        liste_u[c+1]=-i*liste_u[c]+liste_u[c-1] #je ne sias pas pourquoi ça marche mais ça marche
    c+=1
    u=liste_u[-1]
    if a<0:
        u*=-1
    v=(-a*u+1)//b
    print(a,"*",u,"+",b,"*",v,"=1")
    return u,v

def diophant(a,b,c): #pour résoudre ax+by=c
    if type(a+b)!=int:
        return None
    d=pgcd(a,b)
    if pgcd(d,c)==d:
        if abs(a)<abs(b):
            a,b=b,a
        e,f,g=a//d,b//d,c//d
    else:
        return None
    u,v=bezout(e,f)
    print("(",u*g,"+",f,"k)*",a,"+","(",v*g,"-",e,"k)*",b,"=",c)

la fonction diophant n'est pas parfaite, pour diophant(17,-23,5), elle affiche ( -15 + 17 k)* -23 + ( -20 - - 23 k)* 17 = 5
le résultat est juste si on remplace - - par +
Ce script permet de combler l’absence de fonctions très pratiques absentes dans la version 1.3

Re: quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 16:32
by Bisam
Plusieurs remarques sur tes fonctions :
  • Il vaut mieux utiliser la commande assert test, texte_à_afficher_si_erreur pour faire les vérifications que tu fais au début de tes fonctions. Elle arrête l'exécution proprement et affiche une erreur plutôt que de renvoyer un résultat potentiellement aberrant dans l'exécution d'une autre fonction.
  • Il serait bon de dire que la fonction isprime n'est pas du tout efficace pour les nombres plus grands que quelques milliards.
  • Pour la liste des diviseurs, il y a un moyen simple de l'écrire dans l'ordre croissant, sans changer grand chose à la vitesse d'exécution.
  • Pour les deux dernières fonctions, tu utilises bien plus de variables que nécessaires dans l'algorithme d'Euclide étendu.

Si ça intéresse quelqu'un, je peux fournir ma propre version de ces fonctions classiques... et de bien d'autres.

Re: quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 18:39
by loupiot
oui je veux bien en voir certaines :p
notamment si t'as un meilleur test de primalité

Re: quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 18:46
by parisse

Re: quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 18:51
by Bisam
Je joins un fichier contenant mes fonctions "usuelles".

Ci-dessous, je recopie quelques-unes des fonctions que l'on peut y trouver :
Elles sont commentées et typées (puisque la syntaxe actuelle de Python le permet).
Code: Select all
def diviseurs(n: int) -> list:
    """Renvoie la liste des diviseurs de n"""

    assert n > 0, "n doit être un entier strictement positif"

    petits = []
    grands = []
    r = intsqrt(n)
    for i in range(1, r):
        if n % i == 0:
            petits.append(i)
            grands.insert(0, n // i)
    if n == r*r:
        petits.append(r)
    return petits + grands


def isprime(n: int) -> bool:
    """Teste si n est premier ou non par le test de Miller - Rabin.
    Fournit une réponse sûre lorsque n < 341,550,071,728,321."""

    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] # etc.
    if n < 2:
        return False
    for p in small_primes:
        if n < p * p:
            return True
        if n % p == 0:
            return False
    r = 0
    s = n - 1
    while s % 2 == 0:
        r += 1
        s //= 2
    if n < 1373653:
        needed_tests = 2
    elif n < 2152302898747:
        needed_tests = 5
    else:
        needed_tests = 7
    for a in small_primes[:needed_tests]:
        x = pow(a,  s,  n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

def eratho(N: int) -> list:
    """Renvoie une liste de True ou False, pour le crible d'Érathosthène"""

    t = [False, False, True] + [i % 2 == 1 for i in range(3, N+1)] # on initialise tous les pairs à False
    r = intsqrt(N)
    i = 3
    while i <= r:
        if t[i]: #si i est premier
            for j in range(i**2, N, i): #les multiples de i plus petits que i^2 sont déjà éliminés
                t[j] = False #on met à False tous les multiples de i
        i += 2
    return t


def primes_list(N: int) -> 'generator':
    """Renvoie (un générateur de) la liste des nombres premiers inférieurs ou égaux à N"""

    isprime = eratho(N) #création de l'oracle d'Ératosthène
    for a in range(N+1):
        if isprime[a]:
            yield a


def gcd(a: int, b: int) -> int:
    """Renvoie le pgcd de a et b"""

    while b:
       a, b = b, a % b
    return a


def bezout(a: int, b: int) -> (int, int):
    """Renvoie un couple (u,v) tel que u*a+v*b=pgcd(a,b)"""

    u0, v0, u1, v1 = 1, 0, 0, 1
    while b:
        a, (q, b) = b, divmod(a, b)
        u0, v0, u1, v1 = u1, v1, u0 - q * u1, v0 - q * v1
    return u0, v0


def intsolve(a: int, b: int, c: int) -> int:
    """Renvoie la plus petite solution entière positive de l'équation a*x=b[c]"""

    c0 = c
    u0, v0, u1, v1 = 1, 0, 0, 1
    while c:
        a, (q, c) = c, divmod(a, c)
        u0, v0, u1, v1 = u1, v1, u0 - q * u1, v0 - q * v1
    if b % a:
        return False
    return ((b // a) * u0) % (c0 // a)


def intsqrt(n: int) -> int:
    """Calcule la racine carrée entière de n,
    c'est-à-dire l'unique entier r tel que r²<=n<(r+1)²
    à l'aide uniquement d'additions et divisions."""

    q = n
    l = []
    while q != 0:
        q, r = divmod(q, 100)
        l.insert(0, r)
    res = 0
    j = 1
    r = 0
    for a in l:
        r = 100 * r + a
        u = 0
        while r >= j:
            r -= j
            j += 2
            u += 1
        res = 10 * res + u
        j = 10 * j - 9
    return res

Re: quelques fonctions d'arithmétique

Unread postPosted: 20 Feb 2018, 20:16
by parisse
Bisam, je ne vois pas bien l'interet d'utiliser un generateur pour renvoyer des nombres premiers, puisqu'on va verrouiller en memoire un tableau de N booleens. Pourquoi ne pas creer directement la liste des premiers plus petits que N a partir du crible? Sauf erreur, on stocke environ N/ln(N) premiers de taille <= ln(N), ca ne devrait pas occuper plus de memoire sauf pendant la conversion. Et ca permet plus d'usages directs, par exemple ithprime.
Ce sera interessant de voir jusqu'a quelle valeur de N on peut faire tourner le crible sur la Numworks avant le memory full. Comme la liste n'est pas typee, chaque booleen occupera au moins 1 mot 32 bits (en micropython, en Python il me semble que c'est 3 mots, soit 24 octets en 64 bits), on ne devrait pas pouvoir depasser 50 milles mais ca risque d'etre trop long de toutes facons. Pour ces valeurs de N, le stocakge utilisera 1 mot par entier, donc le tableau de premiers occupera moins de place que le crible.

Re: quelques fonctions d'arithmétique

Unread postPosted: 21 Feb 2018, 09:25
by Bisam
@parisse : Tu pointes effectivement un problème de ma fonction eratho.
C'est une des toutes premières fonctions que j'aie écrites en Python, et je ne me suis pas reposé la question de la mémoire.
Sur mon ordi, on peut stocker une liste de booléens comportant environ 10^8 valeurs. Ce serait déjà sans doute bien plus efficace si on utilisait un vrai tableau typé comme ceux de numpy.

Ensuite, tu as parfaitement raison sur le fait que renvoyer les premiers 1 par 1 n'est pas forcément le plus intéressant. Mon idée était d'avoir une réponse à la question "est-ce un nombre premier ?" en temps constant pour des valeurs raisonnables.

Ensuite, j'ai écrit l'algorithme de Miller-Rabin pour ce test... et il est très efficace, même sur un très grand nombre de valeurs.

[Edit] Je viens de faire un test avec les tableaux numpy... La différence en temps d'exécution n'est pas significative en dessous de 100000 valeurs, mais la taille en mémoire, en revanche, est drastiquement diminuée ! Par exemple, la liste eratho(300000) de mon premier jet occupait 2 400 048 octets, soit près de 2,4Mio, tandis que la version numpy n'occupe plus que 80 octets !!

D'ailleurs, c'est bizarre... Comment un tableau de longueur N rempli de 0 et de 1 peut-il prendre moins de N bits en mémoire ? Pire encore, comment se fait-il qu'un tableau numpy de type booléen prenne toujours 80 octets, quelque soit sa taille ?

Re: quelques fonctions d'arithmétique

Unread postPosted: 21 Feb 2018, 11:44
by parisse
La memoire d'un tableau numpy est certainement allouee hors de l'espace de visibilite de Python, numpy etant une librairie ecrite en C, mais elle n'est pas portee sur la Numworks et a mon avis elle a peu de chances de l'etre vu la taille de la flash. Concernant eratho, il me semble normal, visiblement un tableau d'objets Python occupe au minimum 8 octets par objet (64 bits), un tableau optimal (vector<bool> en C++) utiliserait 64 fois moins de memoire (pour toute valeur de n). C'est sur prime_list que je portais une reserve: je pense qu'il vaut mieux parcourir eratho une seule fois en creant la liste des premiers et liberer le tableau de booleens ensuite, la liste de premiers occupera alors moins de place et a l'avantage de pouvoir etre adressee en temps O(1) si on veut connaitre le j-ieme nombre premier.
Sinon, merci pour le code, je suis en train de modifier le source de giac pour qu'il passe dans mon interpreteur presque sans modifications. Pour l'avoir sans modification, il faudrait:
1/ utiliser des noms qui ne sont pas des noms de commandes giac, par exemple en mettant une majuscule a isprime, ou en utilisant des noms en francais (pgcd au lieu de gcd)
2/ ne pas utiliser + pour concatener deux listes (ni * pour repliquer une liste) car chez moi les listes sont des vecteurs et + et * ont ce sens. Mais je ne sais pas s'il existe un equivalent de concat en Python...

Re: quelques fonctions d'arithmétique

Unread postPosted: 21 Feb 2018, 11:56
by Bisam
@parisse : Pour concaténer deux listes l1 et l2, tu peux utiliser l1.extend(l2). Cela modifie l1, cependant... mais c'est souvent ce que l'on veut faire. Si on ne veut pas modifier l1, il faut d'abord faire une copie l3 de l1 puis étendre l3.
Tu peux aussi tricher en utilisant les méthodes cachées : l1.__add__(l2) renvoie une nouvelle liste qui est la concaténation de l1 et l2 (c'est exactement la même chose que l1 + l2 puisque c'est en fait la méthode qui est appelée lorsque le parseur reconnait cette écriture).

Re: quelques fonctions d'arithmétique

Unread postPosted: 21 Feb 2018, 14:13
by parisse
extend ira tres bien, il reste la duplication avec *, c'est utilise souvent?