MPO n°71 : Wilson au carré

[C.Ret]

Cette version est définitivement trop lente.

En effet, les temps de test de la divisibilité et du calcul des factorielles augmentent exponentiellement.
Voici une version bien plus rapide qui donne les trois résultats en moins de 5 heures

L'accélération est due au fait que le test de divisibilité n'est effectué que pour les nombres premiers.
Le test prend un peu de temps, mais ce temps est rapidement amorti par en évitant une palanquée de tests de divisibilité inutiles.

Par contre la boucle de recherche parcourt touts les entiers afin de construire toutes les factorielles et de mémoriser tous les chiffres.

Pour éviter les erreurs d'arrondi, je ne peux mémorisé dans chaque registre que 7 chiffres (sur les 10 d'un registre SHARP PC-1211). En effet, il faut pouvoir calculer le produit de N (au plus trois chiffres) avec les chiffres de la factorielle et ne pas dépasser 10 pour un résultat correct.

Ce qui pose un problème pour le test de divisibilité. Je ne peux l'effectuer sur les 7 chiffres de chaque registre de la factorielle, car le carré peut lui aussi faire jusqu'à 6 chiffres ce qui va dépasser la capacité de calcul limité à 10 chiffres. Je suis donc obligé de faire le test de divisibilité en deux étape pour chaque registre de la factorielle. En prenant 3 ou 4 chiffres j'obtiens un calcul du module correct.
C'est évidemment plus lent, mais comme le test n'est effectué que pour les nombres premiers, je gagne du temps, surtout vers la fin de la recherche où les nombres premiers sont de plus en plus rares.


Donc, voici le code pour SHARP PC-1211 qui lui fonctionne et donne les trois nombres de Wilson inférieurs à 1000 dont le carré est multiple de la factorielle précédente.

Code : Tout sélectionner

PROGRAMM  
1:F=7,G=2:FOR B=3 TO 997:C=B,D=2                                                                     
2:IF D<C LET C=B/D,D=D+1+(D>2):GOTO 2+2*(C=INT C)                                                    
3:C=0:FOR A=F TO 7 STEP -1:E=Ḛ3,D=A(A)/10E:GOSUB 9:E=Ḛ4,D=DE+(A=7):GOSUB 9:D=C:NEXT A:IF D=0 PRINT B 
4:C=0,E=Ḛ-7:FOR A=7 TO F:A(A)=BA(A)+C,C=INT EA(A),A(A)=A(A)-C/E:NEXT A:IF C LET F=F+1,A(F)=C         
5:NEXT B:END                                                                                         
9:C=CE+INT D,D=D-INT D,C=C-BB*INT (C/BB):RETURN                                                      

Code : Tout sélectionner

MEMORIES
1:  oo27    Initie  F=2!
2:  oo37    Test primalité de n
3:  oo68    Test divisabilité en deux étapes - calcul F(…) MOD p² 
4:  oo74    Calcule  F(n)=n*F(n-1)
5:  ooo7    Boucle pour tous les n
9:  oo28    SS-prg calcul modulo
TOT o241 octets

Code : Tout sélectionner

REGISTERS : VARIABLES
Ligne  2:Test Primalité    3:Test Divisibilité      4:Calcul Factorielle  9:Ss-Prg:Calcul Module   
A:      -                  i=indice d'F(i)          i=indice d'F(i)         -             -
B:     n=nombre recherché  p=nombre premier         m=multiplicateur      IN:p            -
C:     q=ratio n/d         r=retenue modulo n²      r=F(…) DIV 1E7        IN:retenue   Out:module n²  
D:     d=diviseurs testés  d= <- FFF.ffff(…)         -                    IN:donnée    Out:partie fractionaire    
E:      -                  e= 1e4 ou 1e3            e=1e7                 IN:puissance    -
F:      -                  f= ^n° dernier F(…)      f= ^n° dernier F(…)     -             - 
G:      -               F(7)=                  
…                       F(…)=  Chiffres de la factorielle
A(F):                   F(F)=

Code : Tout sélectionner

EXECUTION
   4'25                                                                24=!(   5.)    
  16'23                                                         479001600=!(  13.)    
4h32'23"  …925300632803498973864795497173769781248000000000000000000000…0=!( 563.)  

[C.Ret]

Pouff ! 5h sur SHARP PC 1211
Bon algorithme qui donne le même résultat en moins d'une seconde sur HP'

Code : Tout sélectionner

 5 :                                                                                                                   24
13 :                                                                                                            479001600
563:                       1128062102650147335604931877600826356732599785899156109948799220678512168507765234982483167383
7539058091447061350404912288224752676166232504987238597734405218500241472929420765736368417076000786819232752956209532030
1752682551980654899150765947042244162080245635559296669405015282601776448053508845595778363500203282317374589290155381842
1944325990177606926800286097414627054051864536041870228285721514685385860292377988279901751941225277578878135356899314842
0813759218092818613683428979748535035060534255480019959564308157354162160693833915878636825625939973652741227366357593618
8904743783303502232789030651288843466602506358615343735588012541284262325375959472028210304185657819879248669560192564140
0233914291862551103903663877005347420624372478858572341037671247912700942422558687330698793138312803725477763383475708100
5076422592059604736922016788198489589238019374441997355890966064057477037099668260531617003866321836270934196409597472684
3275715521979250368348435703471063132673555492450549572930312250812511626370105274275251208167297520077950422716994917370
4232465321579395032843200164749811794626780197808167144622306371292530063280349897386479549717376978124800000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Avec le code ci-dessous (les n° de ligne indiqués correspondent aux lignes du programme en BASIC correspondant)

Code : Tout sélectionner

      #cas
      WILSON():=
  1   BEGIN  LOCAL n,f:=2; FOR n FROM 3 TO 997 DO
  2     IF isprime(n) THEN 
  3        IF ((1+f) MOD (n^2))=0 THEN PRINT(n&":"&f); END; END;
  4     f*n►f;
  5   END; END; 
      #end

J'ai sorti mon HP-41C, mais je m'attends à ce que les DM 42s trouvent quelque chose bien avant moi.