La question (très simple) du dimanche après-midi.

Ici, on fait dans le petit, le LCD qui déchire sa race, on y cause même calculatrices quand on est en manque !

Modérateur : Politburo

Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

Bon, comme nous sommes vendredi soir, je vais donner quelques informations supplémentaires.

Les nombres que j'ai choisi sont volontairement les cas les plus défavorables pour la méthode par recherche des facteurs premiers. Ce qui explique certainement pourquoi la commande ISPRIME? du NewRPL est moins rapide que la méthode Rabin-Miller et autres méthode basée sur le petit théorème de Fermat.


Comme se sont de très grands nombres (respectivement 9 et 10 chiffres) c'est aussi un potentiel souci, car certaines méthodes probabilistes nécessitent de calculer une "Exponentiation Modulaire". Sur une machine dont le nombres de chiffres est limité, cela pose aussi un problème.

Lorsque a et b sont du même ordre de grandeur que n le nombre à tester , l'étape d'exponentiation a^b modulo n --- c'est à dire « a b * n MOD » ou (a*b MOD n) ou A*B-N*INT(A*B/N) --- posant un problème d'arrondi.

Problème d'arrondi qui peut être surmonté en décomposant le calcul comme cela fait proposée dans la Gazette n°2.
Mais la décomposition peut ralentir sérieusement ces méthodes et pourrait par là même les rendre aussi peu efficaces que la recherche systématique des éventuels facteurs premiers de n.
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Gilles59
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1602
Enregistré le : 27 oct. 2010 20:46

Re: La question (très simple) du dimanche après-midi.

Message par Gilles59 »

Bon je n'aurai pas le temps de faire le prog 602P.

Par contre je me suis 'amusé' avec la méthode probabiliste.Avec juste 3 nombres aléatoires je trouve les 10.000 premiers nombres premiers supérieur à 1.000.000 avec un taux de réussite de 100% … Idem avec seulement 2 nombres alétoires. Avec un seul, j'ai seulement 6 "faux positifs"...
(edit 100% de succès pour les nombres premiers < 10 millions avec 7 vérifications. Pourquoi 7 ? vu sur le net …)

Avec le simulateur, ces 10.000 premiers nombres premiers supérieur à 1.000.000 sont trouvés en … 1,8sec (420 sec sur la HP50g, soit un ratio nettement plus mauvais que d'habitude) . Le newRPL gérant nativement la multiprécision on peut aller jusque 2000 chiffres. Le prog fait 548 octets.

je me suis basé sur cet algorithme en ADA .

Code : Tout sélectionner

with Ada.Numerics.Discrete_Random;
 
package body Miller_Rabin is
 
   function Is_Prime (N : Number; K : Positive := 10)
                     return Result_Type
   is
      subtype Number_Range is Number range 2 .. N - 1;
      package Random is new Ada.Numerics.Discrete_Random (Number_Range);
 
      function Mod_Exp (Base, Exponent, Modulus : Number) return Number is
         Result : Number := 1;
      begin
         for E in 1 .. Exponent loop
            Result := Result * Base mod Modulus;
         end loop;
         return Result;
      end Mod_Exp;
 
      Generator : Random.Generator;
      D : Number := N - 1;
      S : Natural := 0;
      X : Number;
   begin
      -- exclude 2 and even numbers
      if N = 2 then
         return Probably_Prime;
      elsif N mod 2 = 0 then
         return Composite;
      end if;
 
      -- write N-1 as 2**S * D, with D mod 2 /= 0
      while D mod 2 = 0 loop
         D := D / 2;
         S := S + 1;
      end loop;
 
      -- initialize RNG
      Random.Reset (Generator);
      for Loops in 1 .. K loop
         X := Mod_Exp(Random.Random (Generator), D, N);
         if X /= 1 and X /= N - 1 then
        Inner : for R in 1 .. S - 1 loop
               X := Mod_Exp (X, 2, N);
               if X = 1 then return Composite; end if;
               exit Inner when X = N - 1;
            end loop Inner;
            if X /= N - 1 then return Composite; end if;
         end if;
      end loop;
 
      return Probably_Prime;
   end Is_Prime;
 
end Miller_Rabin;
EDIT : la "présomption de primalité' d'un nombre premier de 100 chiffres prends 0,73s sur la HP50, 0,006 sur l'émulateur. Commande ISPRIME? cherche encore … Bon connecté sur la prise USB je laisse chercher sans bouffer les piles ;D (edit, bon ca tourne depuis plus de 10 mn je stoppe)
Casio FX-502P /602P / 603P / FX180P+ / FX4000P / TI57 / TI66 / TI74 Basicalc / TI95 Procalc / HP12C / HP15C LE / DM41L / HP 30B / HP39GII / HP 48SX USA / 49G / 49g+ / 50G / 50G NewRPL / HP Prime / Oric 1 / Amstrad CPC 6128+ CM14 et MM12 / Alice 32
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

Intéressant, ton code en ADA semble utiliser le même algorithme type Miller Rabin que j'utilise sur mon HP-28S

Les nombres composés sont détecté en quelques secondes (le plus long semble d'imprimer la trace) et les nombres premiers après une ou deux minutes en itérant 9 fois le test (neuf nombres aléatoires).

C'est un progrès certains par rapport aux codes qui cherchent par décomposition en facteurs premiers qui mettent sur la même HP-28S quelques dizaines de minutes ou quelques heures.

Il faut dire que j'ai construit les nombres afin que cette méthode soit moins efficace, car pour les nombres composés, le premier facteur se trouve assez élevé, c'est à dire assez proche de la racine carrée du nombre.

Code : Tout sélectionner

  9831013529 = 2113 x 2153 x 2161 
  9831013559 est premier
  
 99875772649 = 4637 x 4639 x 4643
 99875772671 est premier
Avec la méthode déterministe, le plus long est de détecter les deux nombre premier car il faut tester les éventuel diviseurs jusqu'à respectivement 99151 et 316031. Affiner la méthode pour parcourir de façon les facteurs possibles jusqu'à ces limites vaut le coup.

Alors que pour les nombres composés, une fois le premier facteur trouvé, l'effort à fournir se réduit très rapidement. En effet, l'algorithme divise le nombre testé par le premier facteur ce qui réduit d'autant la limite supérieure (respectivement 2156 et 4640). Le dernier facteur étant le 'résidu' (respectivement 981013529/2113/2153 = 2161 et 99875772649/4637/4639 = 4643)


Je donne ci-dessous une version RPL de la décomposition en facteurs premiers avec un parcourt des valeurs testées peu évolué.

Code : Tout sélectionner

FCTRS:
« 'n' SWAP 2 WHILE DUP2 SQ >=
             REPEAT 0  ROT ROT
               WHILE DUP2 MOD NOT
               REPEAT ROT 1 + ROT 3 PICK / ROT
               END 4 ROLL 4 ROLL -> p « IF p THEN OVER IF p 1 > THEN 'x' ^ 3 p EXSUB END * END »
               ROT ROT DUP 2 > + 1 +
            END DROP 1 SWAP EXSUB »
le code est basé sur deux boucles enchevêtrées:
  • La boucle externe: les facteurs f = 2,3,5,7,9,11,13,... sont parcourus jusqu'à atteindre la racine carrée de n (cf. test n>=f² )
  • La boucle interne: qui compte le nombre p de fois où n est divisible par f tout en réduisant n par f
La dernière partie du code permet de mettre en forme la décomposition sous la forme d'une expression algébrique:
Par exemple 212340480 FCTRS donne '457*2^8*3*5*11^2'

Code : Tout sélectionner

   4:    3:       2:       1:  FCTRS:                           // Décompose en facteurs premiers
                            n  «                                //   
        'n'        n        2   'n' SWAP 2                      // Initialise produit et facteur f
      'eX*'        n        f   WHILE                           // ================================
'eX*'     n        f     n≥f²     DUP2 SQ >=                    //   Tant que n >= f²
      'eX*'        n        f   REPEAT                          //   ==============================
'eX*'     n        f        0     0  ROT ROT                    //   Initialise p et facteur f
'eX*'     p        n        f     WHILE                         //   ------------------------------
'eX*'     p        n        f       DUP2 MOD NOT                //     Tant que n est multiple de f
'eX*'     p        n        f     REPEAT                        //     ----------------------------
'eX*'     n        f      p+1       ROT 1 +                     //       - incrémente puissance p 
'eX*'     f      p+1      n/f       ROT 3 PICK /                //       - divise n par facteur f
'eX*'   p+1      n/f        f       ROT                         // 
'eX*'     p        n        f     END                           //   ------------------------------
    n     f    'eX*'        p     4 ROLL 4 ROLL -> p            //   réarrange et localise p
          n        f    'eX*'     « IF p THEN                   //     Si p est non nul 
    n     f    'eX*'        f         OVER                      //       copie le facteur
    n     f    'eX*'        f         IF p 1 > THEN             //     Si p>1 alors 
    n     f    'eX*'    'f^p'           'x' ^ 3 p EXSUB END     //       construit 'f^p'
          n        f 'eX*f^p'         * END »                   //     ajoute f ou f^p au produit 
      'eX*'        n        f     ROT ROT                       //
      'eX*'        n     f+df     DUP 2 > + 1 +                 //   augmente f de 1 (puis 2) 
      'eX*'        n        f   END                             // ================================
               'eX*'        n   DROP                            // supprime le dernier f
               'eX*'        n   1 SWAP EXSUB »                  // ajoute résidu n au produit final
                    'n*e*f^p'                                   // renvoi forme algébrique  
En effet, les facteurs testés sont 2,3 5 , 7 ... donc tous les nombres impairs après 2. Un nombre certain de multiples sont donc utilisés inutilement. L'efficacité de cet algorithme peut être sensiblement amélioré en sélectionnant mieux les facteurs parcourus. Par exemple parmi les nombre impairs présents entre 3 et 316031, 52671 sont des multiples de 3, 31602 sont des multiples de 5, 22573 sont multiples de 7 , 14364 multiples de 11 etc..
Donc, simplement en retirant de la liste parcourue ces multiples l'économie de tests est de 64845 soit une économie de 41%.
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Gilles59
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1602
Enregistré le : 27 oct. 2010 20:46

Re: La question (très simple) du dimanche après-midi.

Message par Gilles59 »

Version newRPL sur HP50. Algorithme de Miller-Rabin

Code : Tout sélectionner

« → n k 
 « 
  IF n 2 == THEN 1 EXIT END
  IF n 2 MOD NOT THEN 0 EXIT END 
  n 1 - 'd' LSTO 
  -1 's' LSTO 
  
  WHILE d 2 MOD NOT REPEAT 
   d 2 IQUOT 'd' STO 
   1 's' STO+ 
  END 

  1 SF 
  n MODSTO 

  1 k START 
   2 n 3 - RAND * IP + d POWMOD 'x' LSTO 
   IF x 1 ≠ x n 1 - ≠ AND THEN 
    1 s START  
     x 2 POWMOD 'x' LSTO 
     IF x 1 == THEN 1 CF EXIT END  
     IF x n 1 - == THEN EXIT END
    NEXT   
    IF x n 1 - ≠  THEN 1 CF EXIT END
    IF 1 FC? THEN EXIT END
   END
  NEXT
  IF 1 FS? THEN 1 ELSE 0 END 
 »
»
'Miller' STO
Ce programme ne fonctionne pas en RPL classique pour 2 raisons :
LSTO qui crée des variables locales (peut être remplacé par STO)
EXIT qui permet de sortir d'un niveau de boucle (ou du programme si pas dans une boucle) et n'existe pas en RPL HP. Pas d'équivalent. Il faudrait faire une boucle FOR NEXT et forcer la borne de sortie pas exemple. ou DO REPEAT … De toute façon ça sera 100x plus lent au moins et probablement beaucoup plus.

Si on veut tester des nombres de plus de 32 chiffres, faire xxx SETPREC (SET PRECision) avec xxx comme nombre de chiffres significatifs souhaités.
L'algo (et le NewRPL) est redoutablement efficace sur les très grands nombres. La primauté (?) probable d'un nombre de 100 chiffres est détectée presque immédiatement.


Usage :n k Miller ou 'Miller(n,k)'

n: nombre à tester
k: au moins 1, plus ce nombre est grand plus la fiabilité sera forte. k=7 donne une très grande fiabilité mais un simple k=3 donne déjà une bonne fiabilité.

A noter que le newRPL gère l'unicode. Il suffit de copier-Coller le programme dans l'émulateur et le transférer sur la calculatrice avec la commande USBSEND
Casio FX-502P /602P / 603P / FX180P+ / FX4000P / TI57 / TI66 / TI74 Basicalc / TI95 Procalc / HP12C / HP15C LE / DM41L / HP 30B / HP39GII / HP 48SX USA / 49G / 49g+ / 50G / 50G NewRPL / HP Prime / Oric 1 / Amstrad CPC 6128+ CM14 et MM12 / Alice 32
Gilles59
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1602
Enregistré le : 27 oct. 2010 20:46

Re: La question (très simple) du dimanche après-midi.

Message par Gilles59 »

Puisqu'il faut misez petit, un algo très simple pour Casio FX501..603P :

Code : Tout sélectionner

   
 *** P0  30 Steps...
  Min01 / 2 MinF = FRAC x=0 GOTO9
  MR01 √ Min1F 1 MinF
  LBL0
   2 M+F
   MR01 / MRF = FRAC x=0 GOTO9
   MR1F x>=F GOTO0
   MR01 MinF
 LBL9
   MRF
Nombre P0
retourne le même nombre si premier, le premier diviseur autrement

par ex :
9831013529 P0
retourne 2113 en 42s avec la 603P soit 25 boucles/secondes. On peut estimer qu'elle mettrait ~1heure pour prouver que 9831013559 est premier avec cet algo basique. On pourrait gagner ~20% de temps en supprimant la boucle des 3, un peu moins avec la boucle des 5, puis des 7 etc... Au lieu d'ajouter 2 à chaque boucle il faudrait ajouter 2 puis 4 puis … et reboucler à 2. Intégrer un crible d'érastotène en fait...
Modifié en dernier par Gilles59 le 10 juin 2019 17:10, modifié 11 fois.
Casio FX-502P /602P / 603P / FX180P+ / FX4000P / TI57 / TI66 / TI74 Basicalc / TI95 Procalc / HP12C / HP15C LE / DM41L / HP 30B / HP39GII / HP 48SX USA / 49G / 49g+ / 50G / 50G NewRPL / HP Prime / Oric 1 / Amstrad CPC 6128+ CM14 et MM12 / Alice 32
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

Les utilisateurs d'HP-28S vont m'en vouloir car il ne pourront pas transférer directement les codes ci-dessous.

Par contre ceux-ci seront utilisables sur tous les RPL sans bibliothèques spéciales.

La méthode de Miller-Rabin est implémentée en trois parties, le programme principal 'MilRab' utilisant deux sous-programmes.

Code : Tout sélectionner

MilRab:                     Usage:  n k MilRab --> { n 'Cmp.' OR 'Prm?' } 
« → n k « n 1 - 0 WHILE OVER 2 MOD NOT 
                  REPEAT 1 + SWAP 2 / SWAP END
                  → d s « 1 k START
                               1 SF n 2 - RAND * 2 + IP d n BEM 
                               IF DUP 1 == OVER n 1 - == OR
                               THEN 1 CF 
                               ELSE 1 s START 2 n BEM
                                         IF DUP n 1 - == THEN 1 CF END
                                    NEXT END DROP
                               IF 1 FS? THEN n {"Cmp."} + KILL END
                         NEXT  n { "Prm?" } + » » »
- 'BEM' qui implémente l'exponentiation modulaire a e ^m MOD selon la méthode efficace décrite dans WikiPédia à savoir la décomposition binaire de l'exposant.

Code : Tout sélectionner

BEM:                                    Usage:  a e m BEM --> a^e (MOD m) 
« → m
  « 1 WHILE OVER
      REPEAT
        IF OVER 2 MOD THEN 3 PICK m MMOD END
        ROT DUP m MMOD ROT 2 / IP ROT
      END ROT ROT DROP2 » »
- ''MMOD' qui réalise les multiplications modulaires a b * m MOD chiffre à chiffre (comme dans le sous-programme de la Gazetten°2) afin d'éviter tout débordement de capacité.

Code : Tout sélectionner

MMOD:                                  Usage: a b m MMOD -->  a*b (MOD m)
« → a b m 
  « a b MAX a b MIN → a b
    « 0 LOG IP 0 FOR k
        10 * m MOD a b k ALOG / IP 10 MOD * + m MOD
     -1 STEP » » »
Voici les traces des résolutions selon la méthode probabiliste de Miller-Rabin pour les quatre nombres de la question : pour chaque trace, les Témoins de Miller utilisés (ceux-ci sont tirés pseudo-aléatoirement à l'aide de l'instruction RAND) et le détail des exponentiations modulaires sont imprimés.

Pour 9831013529 j'ai lancé en tapant 9831013529 9 [MilRab]
9831013529.gif
9831013529.gif (42.31 Kio) Vu 12803 fois
9831013529 - 1 = 1228876691 * 2^3 donc d=1228876691 et s=3
Le seul témoin pseudo-aléatoire utilisé est 91951 dont l'exponentiation modulaire 91951^1228876691 mod 9831013529 est calculée ainsi que les trois carrés des restes modulaires successifs qui en découle (s=3).
Aucun de ces restes modulaires n'étant égal à 983101529-1, le programme est arrêté (KILL) en indiquant que 9831013529 est un nombre composé.
La décomposition FCTRS en facteurs premiers le confirme en quelques minutes.

De même, 9831013559 9 [MilRab] produit la trace suivante :
9831018559.gif
9831018559.gif (76.58 Kio) Vu 12803 fois
On a 98310113559 = 4915506779 * 2^1 .
Neufs témoins de Miller sont tirés au sort: 82367, 43492, 5403, 55552, 3503, 77537, 7107, 28700 et 386.
Mais aucun n'est un indicateur convenable, le reste modulaire de l'exponentiation par 4915506779 conduisant à chaque fois soit à un reste de 1 soit à 98310113558 = 98310113559-1
D'où la conclusion qui indique une très forte probabilité que 98310113559 soit premier.
La décomposition en facteur premier, qui met quelques heures, le confirme.


Je donne ci-dessou les traces pour 99875772649 et 99875772671 dont l'analyse ne prends que quelques secondes pour le premier et moins de deux minutes pour le second.
99875772649.gif
99875772649.gif (45.23 Kio) Vu 12803 fois
99875772671.gif
99875772671.gif (78.02 Kio) Vu 12803 fois
Là aussi, les décompositions en facteurs premiers confirment les diagnostiques de Miller-Rabin en moins de cinq minutes pour le premier et plus de 6 heures pour le second. Le gain de temps est plus qu'appréciable :)



P.S.: Je donne le listing des codes RPL en détailant le mouvement des arguments (c'est toujours délicat à voir directement dans le code) et donnant quelques commentaires pour repérer les étapes de l'algorithme:

Code : Tout sélectionner

    :     :        :        :                                  ::
 
   4:    3:       2:       1:  FCTRS:                           // Décompose en facteurs premiers
                            n  «                                //   
        'n'        n        2   'n' SWAP 2                      // Initialise produit et facteur f
      'eX*'        n        f   WHILE                           // ================================
'eX*'     n        f     n≥f²     DUP2 SQ >=                    //   Tant que n >= f²
      'eX*'        n        f   REPEAT                          //   ==============================
'eX*'     n        f        0     0  ROT ROT                    //   Initialise p et facteur f
'eX*'     p        n        f     WHILE                         //   ------------------------------
'eX*'     p        n        f       DUP2 MOD NOT                //     Tant que n est multiple de f
'eX*'     p        n        f     REPEAT                        //     ----------------------------
'eX*'     n        f      p+1       ROT 1 +                     //       - incrémente puissance p 
'eX*'     f      p+1      n/f       ROT 3 PICK /                //       - divise n par facteur f
'eX*'   p+1      n/f        f       ROT                         // 
'eX*'     p        n        f     END                           //   ------------------------------
    n     f    'eX*'        p     4 ROLL 4 ROLL -> p            //   réarrange et localise p
          n        f    'eX*'     « IF p THEN                   //     Si p est non nul 
    n     f    'eX*'        f         OVER                      //       copie le facteur
    n     f    'eX*'        f         IF p 1 > THEN             //     Si p>1 alors 
    n     f    'eX*'    'f^p'           'x' ^ 3 p EXSUB END     //       construit 'f^p'
          n        f 'eX*f^p'         * END »                   //     ajoute f ou f^p au produit 
      'eX*'        n        f     ROT ROT                       //
      'eX*'        n     f+df     DUP 2 > + 1 +                 //   augmente f de 1 (puis 2) 
      'eX*'        n        f   END                             // ================================
               'eX*'        n   DROP                            // supprime le dernier f
               'eX*'        n   1 SWAP EXSUB »                  // ajoute résidu n au produit final
                    'n*e*f^p'                                   // renvoi forme algébrique  

    :     :        :        :                                  ::
  
   4:    3:       2:       1:  MMOD:                            // P := a b * m MOD
         a         b        m  « -> a b m                       // mémorise a b et m
           MAX(a,b)  MIN(a,b)    « a b MAX a b MIN -> a b       // Trie  a > b
                            0      « 0                          // Initialise P
         p     lg b         0        LOG IP 0 FOR k             // ------ Pour k de longeur(b) à 0
                      10p % m          10 * m MOD               //    P := 10.P MOD m
         p        a       b_k          a b k ALOG / IP 10 MOD   //    extrait k-ième chiffre de b
                    p+a*b_k%m          * + m MOD                //    P := (P+a.b_k) MOD m 
                            P        -1 STEP » » »              // ------ boucle sur k

    :     :        :        :                                  ::
 
   4:    3:       2:       1:  BEM:                             //  R := b^e MOD m   
         b        e         m  « → m                            //  Mémorise m < 1.E12
         b        e         1    « 1                            //  Initialise R
         b        e         r      WHILE                        // ================================
   b     e        r      e>0?        OVER                       //  Tant que e>0 ? 
         b        e         r      REPEAT                       //   ==============================
   b     e        R  e impair        IF OVER 2 MOD              //     Si e est impair
         b        e r*b (% m)        THEN 3 PICK m MMOD END     //             alors R := R*b MOD m
         e        r   b²(% m)        ROT DUP m MMOD             //     b := b^2 MOD m 
         r       b²       e>>        ROT 2 / IP                 //     e := e DIV 2
         b        e         r        ROT                        //
                                   END                          // ================================
                            R      ROT ROT DROP2 » »            // vide pile pour ne laisser que R 

    :     :        :        :                                  ::

   4:    3:       2:       1:  MilRab:                          // Test Miller-rabin
                  n         k  « → n k                          // tester n en k itérations 
                n-1         0    « n 1 - 0                      // Initialise calcul de d et s
       n-1        s   n pair?      WHILE                        // ................................
                                     OVER 2 MOD NOT             // Tant que n-1 est pair
                n-1         s      REPEAT                       //   ..............................
                                     1 + SWAP 2 / SWAP          //   incrémente s et divise (n-1)
                                   END                          // ................................
                                   → d s                        // mémorise d et s
                                   « 1 k START                  // ================================
                            a          1 SF n 2 - RAND * 2 + IP //   a est témoin (drapeau 1 lévé)
         a        d         n          d n BEM                  //   x := a^d MOD n 
                    a^d MOD n          IF DUP 1 == OVER         //   si
                  x = 1? n-1?            n 1 - == OR            //      x=1 ou x=n-1
                                       THEN 1 CF                //   alors a n'est pas témoin
                            x          ELSE                     //   sinon boucler s fois
                                         1 s START              //      ---------------------------
                    x^2 MOD n              2 n BEM              //      x := x^2 MOD n
                   x^2 ≡ n-1 ?             IF DUP n 1 - ==      //      si x=n-1
                                           THEN 1 CF END        //      alors a n'est pas témoin
                            x            NEXT END DROP          //      ---------------------------
                                     IF 1 FS?                   //   si a est témoin
                   { n Cmp. }        THEN n {"Cmp."} + KILL END //   alors arrêter n est composé
                                   NEXT                         // ================================
                   { n Prm? }      n { "Prm?" } + » » »         // n est probablement premier 

    :     :        :        :                                  ::
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Gilles59
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1602
Enregistré le : 27 oct. 2010 20:46

Re: La question (très simple) du dimanche après-midi.

Message par Gilles59 »

C.Ret a écrit : 10 juin 2019 14:30 Les utilisateurs d'HP-28S vont m'en vouloir car il ne pourront pas transférer directement les codes ci-dessous.

Par contre ceux-ci seront utilisables sur tous les RPL sans bibliothèques spéciales.
- 'BEM' qui implémente l'exponentiation modulaire a e ^m MOD selon la méthode efficace décrite dans WikiPédia à savoir la décomposition binaire de l'exposant. (...) ::
Sauf erreur de ma part POWMOD est natif en RPL depuis la 48. Mais je me trompe peut-être. Sur depuis la 49G en tout cas
Casio FX-502P /602P / 603P / FX180P+ / FX4000P / TI57 / TI66 / TI74 Basicalc / TI95 Procalc / HP12C / HP15C LE / DM41L / HP 30B / HP39GII / HP 48SX USA / 49G / 49g+ / 50G / 50G NewRPL / HP Prime / Oric 1 / Amstrad CPC 6128+ CM14 et MM12 / Alice 32
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

A zut, si je l'avais su plus tôt, j'aurais appelé mon programme POWMOD au lieu de BEM !
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Gilles59
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1602
Enregistré le : 27 oct. 2010 20:46

Re: La question (très simple) du dimanche après-midi.

Message par Gilles59 »

C.Ret a écrit : 11 juin 2019 17:38 A zut, si je l'avais su plus tôt, j'aurais appelé mon programme POWMOD au lieu de BEM !
Je me suis un peu avancé. A priori ca existe dans les 48gII mais pas dans les "historiques" (48s et g). Probablement une fonction héritée du CAS.
Casio FX-502P /602P / 603P / FX180P+ / FX4000P / TI57 / TI66 / TI74 Basicalc / TI95 Procalc / HP12C / HP15C LE / DM41L / HP 30B / HP39GII / HP 48SX USA / 49G / 49g+ / 50G / 50G NewRPL / HP Prime / Oric 1 / Amstrad CPC 6128+ CM14 et MM12 / Alice 32
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

Gilles59 a écrit : 10 juin 2019 14:28Puisqu'il faut misez petit, un algo très simple pour Casio FX501..603P :
Sur fx-602p, le test de divisibilité simple à l'aide de FRAC x=0 pose problème et une erreur d'arrondi donne des résultats parfois incorrect pour nos grands nombres.

Je me suis donc inspiré du programme de décomposition qui se trouve au début du CASIO PROGRAM LIBRARY fx-601P / fx-602p (page 1 et 2)
Une note indique le souci du test avec 2512549139. Pour corriger cela, un second test refaisant le calcul détecte un faux positif dû à l'erreur d'arrondi. Cela alonge un peu le code, mais ne ralenti pas la recherche des facteurs car ce second test n'est effectué que dans les rares cas où FRAC x=0 est positif.

Mon code fait de même, il est donc plus long et plus complexe que celui de Gilles qui utilise le test simple. Pour les nombres petits ( < 10^8 ) on peut se passer de la vérification d'erreur. Dans ce cas, mon code peut être simplifié, il suffit de retirer les codes supplémentaires surlignés en jaune dans P5 :
fx-602p pFACTORS.gif
fx-602p pFACTORS.gif (41.81 Kio) Vu 12690 fois
Le programme principal P0 utilise une boucle qui permet de mieux choisir les facteurs à tester. En incrémentant f de +4+2+4+2+4+6+2+6 etc, j'évite les multiples de 2 3 et 5. Ce qui réduit significativement les facteurs testés. Tout ce passe comme si l'on avançait en moyenne de 3.75 au lieu d'aller de deux en deux. Celle accélération permet d'avoir la décomposition de 9831013529 en environ 2min sur une fx-602p
fx-602p pFACTORS ex.gif
fx-602p pFACTORS ex.gif (42.41 Kio) Vu 12637 fois
Modifié en dernier par C.Ret le 15 juin 2019 16:44, modifié 1 fois.
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Avatar du membre
dizzy33
Fonctionne à 2400 bauds
Fonctionne à 2400 bauds
Messages : 1463
Enregistré le : 13 févr. 2007 20:39
Localisation : Bordeaux

Re: La question (très simple) du dimanche après-midi.

Message par dizzy33 »

Vous me faites peur ... Vous me rappelez certains gars de ma classe, au collège, qui avaient parfois des conversations que je ne comprenais pas, et je me demandais "pourquoi se faire du mal ?"...
**** COMMODORE 64 BASIC V2 ****
64K RAM SYSTEM 38911 BASIC BYTES FREE
READY.
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

dizzy33 a écrit : 13 juin 2019 22:40[…] "pourquoi se faire du mal ?"...
Parceque ça fait tellement de bien quand on s'arrête … :D

Bon, je vois à sa signature que DIZZY33 a de quoi rejoindre notre club de j'me fais du bien en arrêtant de faire mal des maths :
DIZZI's CBM-C64ko color game console.gif
DIZZI's CBM-C64ko color game console.gif (67.61 Kio) Vu 12633 fois
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

Dans une vidéo récente de NUMBER PHILE, Matt Parker utilise la méthode de Miller-Rabin pour déterminer la primalité de 747.

Dès la première minute, il fait une erreur stratégique primordiale en allant chercher une CASIO fx-991 EX.
Witness Numbers (and the truthful 1,662,803) - Numberphile

Puis quelques dizaines de secondes plus loin dans son exposé, il utilise Wolfram Alpha pour calculer 23^373 MOD 747.

J'en déduis qu'il n'a pas d'HP-prime ni d'HP-28S et il n'a certainement jamais entendu parler de mon code MILRAB. Mais je ne comprends pas pourquoi il n'a pas utilisé sa CASIO pour calculer 23^373 MOD 747 fait 131.

Witness Numbers (and the truthful 1,662,803) - Numberphile

Sur une HP 28S cela se fait très facilement et immédiatement [alpha] 1 1 373 START 23 * 747 MOD NEXT [ENTER] affiche 131 en moins de 15 secondes ??
La CASIO fx 991 EX ne sait pas faire de même ??


Par contre, la suite de son exposé parle d'un sujet fort intéressant concernant le choix des 'Témoins de Primalité à utiliser. Il semblerai qu'avec un nombre très réduit de témoins bien choisis, on peut effectuer le test de Miller-Rabin tout en ayant une bonne présomption du résultat jusqu'à des nombres premiers de plusieurs milliards. Hummm!
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Avatar du membre
Over_score
Fonctionne à 300 bauds
Fonctionne à 300 bauds
Messages : 166
Enregistré le : 26 mars 2019 14:55
Localisation : Pas loin de Smartville

Re: La question (très simple) du dimanche après-midi.

Message par Over_score »

Il y a plusieurs ensembles de témoins de primalité ici : https://en.wikipedia.org/wiki/Miller%E2 ... s_of_bases
Avatar du membre
C.Ret
Fonctionne à 9600 bauds
Fonctionne à 9600 bauds
Messages : 3400
Enregistré le : 31 mai 2008 23:43
Localisation : N 49°22 E 6°10

Re: La question (très simple) du dimanche après-midi.

Message par C.Ret »

C'est très intéressant, car du fait de ces listes réduite de témoins bien choisi, il est facile d'écrire un code pour sa calculatrice qui teste rapidement la primalité d'un nombre à partir de quelques témoins.

Par exemple, mon SHARP PC-1360 ne calcule que sur 10 chiffres. Utiliser les témoins 2, 3, 5 et 7 devrais me permettre sans me tromper de tester la primalité des entiers impairs jusqu'à 1'428'571'427 inclus !!

Code : Tout sélectionner

 1 Q=X/M,R=X-M*INT Q: RETURN                          // Calcule le quotien Q=X/M et le reste R=X MOD M 
10 "M" AREAD N: WAIT 0: CLS :X=N-1,M=2,D=0                        // Input N et initialise
12 GOSUB 1:IF  R=0 LET D=D+1,X=X/2: GOTO 12                       // Réduit N = P * 2^D + 1 
14 P=X,M=N,T=2: GOSUB 20: FOR T=3 TO 7 STEP 2: GOSUB 20: NEXT T   // Témoins 2, 3, 5 et 7
16 PRINT N,"PRIME !?": END                            // N est probablement premier.   

20 S=D,R=1: FOR I=1 TO P:X=R*T : GOSUB 1:NEXT I:IF R=1 RETURN    // T suivant si T^P MOD N=1
22 IF R=N-1 RETURN                                               // ou T^P MOD N=-1 ou T^(2×…×2×P) MOD N=-1  
24 S=S-1:IF S>0 LET X=R*R : GOSUB 1: GOTO 22                     // Calcule X² MOD N et teste   
26 PRINT N,"COMPOSITE." : END                         // N est composite.
Utilisation, en mode RUN saisir le nombre à tester et presser DEF M pour lancer le test MILLER-RABIN.

747 DEF-M affiche: 747=373*2+1. puis 747.COMPOSITE car 2^373 MOD 747 = 731.
91 DEF-M affiche: 91=45*2+1. puis 91.COMPOSITE car 2^45 MOD 91 = 57.
1447 DEF-M affiche 1447=723*2^1+1. puis 1447.PRIME !? car on a successivement 2^45 MOD 1447=1, 3^45 MOD 1447=1446, 5^45 MOD 1447=1446 et 7^45 MOD 1447=1.
SHARP PC-1211 PC-1360 EL-5150 PC-E500 | Commodore C=128D | Texas Instruments Ti-57LCD Ti-74BASICalc Ti-92II Ti-58c Ti-95PROCalc Ti-30XPROMathPrint | Hewlett-Packard HP-28S HP-41C HP-15C HP-Prime HP-71B | CASIO fx-602p | NUMWORKS | Graphoplex Rietz Neperlog | PockEmul | Sommaire des M.P.O. | Ma...dov'il sapone.
Répondre

Retourner vers « Tous les Pockets »