MPO n°71 : Wilson au carré

[gege]

Bonjour,
D'après un gars nommé Wilson, pour tout nombre premier p, ((p-1)! + 1) est divisible par p.
(p-1)! étant la factorielle de (p-1).
Ok, mais il se trouve qu'on connaît trois nombres premiers inférieurs à 200 000 pour lesquels ((p-1)! + 1) est divisible par p² (oui, p au carré).

Saurez-vous trouver ces trois nombres ? Deux indices :
- Ils sont inférieurs à 1000
- Les nombres non premiers ne satisferont jamais à la relation (inutile de tester seulement les nombres premiers)

Bonne chasse !!
G.E.

EDIT : voici le sommaire des MPO : Sommaire des MPO

[emond]

Les seuls nombres premiers de Wilson connus sont 5, 13, et 563.

Un article intéressant sur le sujet ici

Et un article sur ce blog

[zpalm]

Sur une machine classique les deux premières valeurs se trouvent facilement, pour la troisième c'est plus dur, même sur la WP 34S en double précision ... j'ai du sortir la HP Prime:

Code : Tout sélectionner

EXPORT MPO70()
BEGIN
  LOCAL j:=2, s;
  WHILE j<=1000 DO
    s:="irem(("+j+"-1)!+1,"+j+"^2)"; 
    IF NOT CAS(s) THEN PRINT(j) END;
    j:=nextprime(j);
  END;
END;

[gege]

Les seuls nombres premiers de Wilson connus sont 5, 13, et 563.

Euh oui mais non, là il y a un carré...
La liste est exacte cependant !

Sur une machine classique les deux premières valeurs se trouvent facilement, pour la troisième c'est plus dur, même sur la WP 34S en double précision ... j'ai du sortir la HP Prime

Ben justement, il faudrait le faire sur une machine "normale"...
G.E.

[dprtl]

Est-ce qu'une tablette, ou un téléphone, avec un interpréteur Python embarqué pourrait battre une calculatrice "normale" ? Peut-être pas avec le code ci-dessous :

Code : Tout sélectionner

from math import factorial
import time

def rwh_primes1(n):
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

print time.asctime(time.localtime(time.time()))
max = 10000
tableau = rwh_primes1(max)
for a in tableau:
  reste = (factorial(a-1)+1) % (a*a)
  if (reste == 0):
    print a
print time.asctime(time.localtime(time.time()))

Résultat (copie d'écran) :

[Gilles59]

Hp49g+. Mode exact. C'est assez rapide sur la vraie calc.

Code : Tout sélectionner

1
 1 3 START
   DO NEXTPRIME DUPDUP 1 - ! 1 + SWAP SQ UNTIL MOD NOT END
   DUP 
 NEXT

[/size] 63 octets

Retourne les 3 solutions sur la pile

[Gilles59]

Une variante qui n'utilise pas les fonctions "nombres premiers" (donc teste tout) mais évite de recalculer à chaque fois la factorielle. C'est 3 fois plus rapide (60 s.) et nul besoin de calculer la suite des nombres premiers.

Mode Exact

Code : Tout sélectionner

1 1 -> p f
<<
 1 3 START
  DO p 'f' STO* f 1 + 'p' INCR SQ UNTIL MOD NOT END
  p
 NEXT
>>

[/size]

Pour optimiser la vitesse on peut aussi éviter de recalculer le carré à chaque boucle sachant que
(a+1)^2=a^2+a+a+1
Est ce plus rapide ? a voir...

[gege]

Bonjour,
Superbe !
Et on doit pouvoir éviter de tester les nombres pairs en déroulant deux fois la boucle.
J'essaye.
G.E.

[dprtl]

Le premier MPO 70 était ici :

viewtopic.php?f=46&t=40139

Celui-ci sera donc le MPO 71 dans le sommaire des MPOs !..

[gege]

Bonjour,
Merci Daniel, l'erreur provenait du sommaire auquel manquait une ligne.
Tout ira bien dès que l'ami Ledudu aura modifié ce dernier.

zpalm, pourquoi n'as-tu pas fait un programme CAS ?

Code : Tout sélectionner

#cas
MPO71g():=
BEGIN
LOCAL j;2►j;
WHILE j<=1000 DO
 IF 0==irem((j-1)!+1,j^2) THEN PRINT(j) END;
 nextprime(j)►j;
END;
END;
#end

<Mode ronchon=ON>
Personnellement je déteste la notation := pour l'affectation. Soit = est l'affectation, soit c'est la comparaison d'égalité.
Mais avoir := pour l'affectation et == pour le test, c'est envoyer à l'utilisateur un gros "je t'emm*** !!".
Vu la fréquence relative d'utilisation de l'affectation et du test, la bonne décision de design est totalement évidente.
Il y a cependant encore des imbéciles en notre siècle pour penser autrement.
Va comprendre Charles !!
<Mode ronchon=OFF>

G.E.

[tyann]

Bonjour
Sur la Prime, en mode Home = fonctionne pour le test
d'égalité.
Essaye ceci :

Code : Tout sélectionner

EXPORT Ess(n)
BEGIN
 PRINT();
 IF n=3 THEN
  PRINT("n=3");
  END;
 PRINT(n=3);
END;

Perso, je m'accommoderai fort bien d'un signe unique pour l'égalité et l'affectation.
Je crois me souvenir qu'en GFA basic sur ST == était un test 'presque égal' qui arrondissait
les nombres réels à 5 ou 6 décimales.

[zpalm]

zpalm, pourquoi n'as-tu pas fait un programme CAS ?

Tout simplement parce que je n'ai pas l'habitude de faire des programmes CAS ... j'ai essayé rapidement, j'ai eu une erreur et je suis passé sur un programme non-CAS. En ce moment je travaille sur un programme BASIC sur Sharp ...

Je suis d'accord avec toi sur le :=, je préfèrerai largement le = tout seul. Par contre je n'arrive pas à me faire au ► et à l'inversion des arguments que cela implique.

[Gilles59]

<Mode ronchon=ON>
Personnellement je déteste la notation := pour l'affectation. Soit = est l'affectation, soit c'est la comparaison d'égalité.
Mais avoir := pour l'affectation et == pour le test, c'est envoyer à l'utilisateur un gros "je t'emm*** !!".
Vu la fréquence relative d'utilisation de l'affectation et du test, la bonne décision de design est totalement évidente.
Il y a cependant encore des imbéciles en notre siècle pour penser autrement.
Va comprendre Charles !!
<Mode ronchon=OFF>
G.E.

:= pour l'affection
= pour les équations
== pour les tests

pour moi les 3 choses sont différentes.

Mais je suis d'accord avec toi sur un point : sur une calculatrice, le nombre de touches à utiliser devraient être réduit au minimum et de ce point de vue, la Prime malgré toutes ses qualités n'est pas optimisée, son langage est en effet très verbeux

[gege]

Bonjour,
Il y a des choses à améliorer sur les fonctions CAS de la Prime.
D'abord, la syntaxe est bizarre :
#cas
toto(a,b):=
BEGIN
...
END;
#end

Ensuite, ces fonctions ne figurent pas dans le catalogue et sont donc difficiles à appeler, il faut taper le nom alphabétiquement sans se tromper.

Enfin, pour un rien ça part en erreur et la Prime est très capricieuse.
J'ai souvent effacé une fonction pour la retaper à force de ne pas comprendre.
On a des bugs bizarres, si vous essayez la fonction MPO71g ci-dessus et que vous l'interrompez, ensuite la fonction irem() renverra systématiquement une erreur "Programme interrompu" lorsqu'appellée depuis l'écran Home, et la fonction MPO71g ell-même plantera pour cette raison.
Pour débloquer la situation, il suffit d'appeler irem() une fois depuis l'écran CAS.
BIZARRE !!!!

Sans parler de l'utilité de taper restart de temps en temps, ça m'a ramené le menu Mem qui autrement plantait la Prime (extinction)... et certains gros programmes ont perdu une bonne partie de la fin...
Zarbi ce truc.
Comme on le sait, on ne perd rien mais c'est quand même seulement 99% stable.
On aurait aimé 100% !
G.E.

[C.Ret]

Les seuls nombres premiers de Wilson connus sont 5, 13, et 563.

Euh oui mais non, là il y a un carré...
La liste est exacte cependant !

Sur une machine classique les deux premières valeurs se trouvent facilement, pour la troisième c'est plus dur, même sur la WP 34S en double précision ... j'ai du sortir la HP Prime

Ben justement, il faudrait le faire sur une machine "normale"...
G.E.

Effectivement, un retranscrivant l'énoncé directement sur mon SHARP PC-1211, j'obtins naïvement le code suivant :

Code : Tout sélectionner

1 : F=1,P=1                                                 10 o
2 : P=P+1,Q=(F+1)/PP,F=FP:IF Q<>INT Q GOTO 2                32 o
3 : PRINT F,P:IF P<Ễ3 GOTO 2                                15 o
                                                            57 octets (et 3 registres)

L'idée est de parcourir tous les nombres p, y compris ceux non premier, afin de calculer petit à petit les factorielles. Il n'y a pas de fonction factorielle sur ce SHARP.

Et les résultats aberrants suivants:

Code : Tout sélectionner

Affichage                Saisie             Commentaires
3:PRINT F,P:IF P<Ễ3GOTO  [M][E][M][ENTER]   Vérifie emprise mémoire du programme
1367STEPS 170MEMORIES    [MODE][MODE]       Active le mode RUN 
>                        [R][U][N] [ENTER]  Lance le programme en mode RUN
        120.          5. [ENTER]            OK trouve le premier nombre (et affiche sa factorielle)
 6227020800.         13. [ENTER]            OK trouve le second nombre
 1.30767Ễ 12         15. [ENTER]            Aïe là ça coince !!!
 2.09228Ễ 13         16.  ...                      ...aïe aïe aïe et ainsi de suite

J'ai donc bien trouvé les deux premières valeurs et je me suis rendu compte de mon erreur. Pour trouver la troisième, il faut que je passe mon SHARP PC-1211 en mode 'calcul exact'. Car il me faut plus de précision dans le calcul des factorielles.
Il faut donc deux lignes supplémentaires:

Code : Tout sélectionner

 1:B=1,E=6,F=1                                                                                14 o
 2:B=B+1,D=0:FOR A=E TO 6 STEP -1:C=Ễ6D+A(A)+(A=6),D=C,C=C/BB,D=D-BB*INT C:NEXT A:D=0         67 o
 3:FOR A=6 TO E:A(A)=BA(A)+D,D=INT Ễ-6A(A),A(A)=A(A)-Ễ6D:NEXT A:IF D LET E=E+1,A(E)=D         66 o
 4:IF C<>INT C GOTO 2                                                                         10 o
 5:PRINT B:IF B<Ễ3 GOTO 2                                                                     13 o
                                                                                             170 octets

Qui revient à faire comme pour la première version, mais en mémorisant tous les chiffres de la factorielle.

Code : Tout sélectionner

 MEM : Variable   Commentaires
-----  ---------- ------------
    A: i          pour itération
    B: p          
    C: q 
    D: r          reste des divisions successives / retenue pour calcul factoriel
    E:  ^f        pointeur du dernier registre contenant le détail de la factorielle
006  :
      ... F(p-1)  factorielle mémorisée en 'full precision'
429  :

Cela marche mieux, je n'obtiens comme résultats que 5. et 13. et j'attends toujours le suivant.
Les temps de calculs sont bien trop longs et de plus, je ne suis pas sûr d'aboutir, il n'y a pas 429 registres disponibles !!!

Je cherche donc toujours une solution ad'hoc ...
J'ai pris rendez-vous avec Lagrange, Euler, Gauss et Stirling...

Pour le moment, toujours aucune réponse, ni de ces derniers, ni sur mon pocket