Misez p'tit Optimisez n°79 : Carré Parfait Interdit !

[C.Ret]

Bonjour à toutes et à tous...



L'idée de ce nouvel M.P.O. m'est venue alors que je divaguais sur la Toile à la recherche de quelques récréations mathématiques.




Je cherchais à résoudre quelques problèmes et à répondre à quelques questions que je pensais faciles à résoudre et escomptais rapidement gagner en notoriété et gloire. Peut-être même glaner ici ou là quelques récompenses en espèces sonnantes et trébuchantes, accompagnée de quelques bouteilles de champagne.

Je compte donc sur votre aimable et désintéressée contribution pour m'aider à programmer mes antiques assistants numériques boutonneux afin que ces derniers puissent me servir à calculer les termes de la suite f(n) ainsi définie.

On désigne par f(n) le n-ième  terme de la suite des entiers naturels qui n'est pas un carré parfait.

Ainsi f(1) = 2, f(2) = 3, f(3) = 5, f(4) = 6,etc...

Ce qui me permettra de vérifier que 2062 est le 2017-ième entier non carré parfait.

Comme à notre habitude en ce lieu, toutes vos propositions de code seront les bienvenues et seront passées à la loupe, analysées, discutées et commentées selon notre tradition maintenant bien établie. Je compte sur vous pour optimiser tout cela, miser petit et efficace.

En effet, une fois armé de ce premier code, le challenge qui nous est proposé est d'emboiter les termes de cette suite f :


On désigne par f^k(n), avec k entier quelconque > 1, la fonction f composée avec elle-même k fois.

Ainsi f^2(2) = f(f(2)) = f(3) = 5
et f^3(17) = f(f^2(17)) = f(f(f(17))) = f(f( 21 )) = f( 26 ) = 31
etc.


Tout cela me permettra de répondre aux questions du problème A2976 proposé sur le site Diophante.fr

A savoir :

Q₁ Calculer f^15(1406).
Q₂ Existe-t-il un entier n tel que f^20(n) = 2017?
Q₃ Trouver l'entier n tel que f^50(n) = 4867.
Q₄ Trouver l'entier n tel que f^2017(n) = 4072325

Pour rappel, les carrés parfait sont interdits, ainsi que toute calculatrice carrée ou parfaite :

[Xerxes]

Salut,

If I've understand the problem correcly, n+rnd(sqrt(n)) gives the n-th element of the sequence and n-int(sqrt(n)) gives the position of a non square number in the sequence.

[C.Ret]

Hi,

You are right this is a way to do it. There is other solutions, which may not be formulae but algorithms. What are in interest in a M.P.O. is how you implement this into code and how to optimize a given code on a specific calculator. The challenge is to present the shortest code.

A good point is you propose the function f(n) and the reciprocal function as well.

[Xerxes]

Ok, I've also tried an assembly language friendly way, using increment and decrement only. Maybe there is a shorter way,
but here the algorithm in PASCAL for f(n):

Code : Tout sélectionner

  readln(n);
  s:=1;
  i:=1;
  repeat
    inc(n);
    inc(i,2);
    inc(s,i);
  until s>n;
  writeln(n);

And the inverse of f(n):

Code : Tout sélectionner

  readln(n);
  s:=pred(n);
  i:=3;
  repeat
    dec(n);
    dec(s,i);
    inc(i,2);
  until s<0;
  writeln(n);

[Xerxes]

The algorithm on the FX-180P for f(n):

Code : Tout sélectionner

 01  1
 02  M+
 03  2
 04  Kin+1
 05  Kout1
 06  Kin+2
 07  Kout2
 08  x<=M
 09  MR

Usage example:

2017 Min 1 Kin1 Kin2 P1

And the inverse of f(n):

Code : Tout sélectionner

 01  1
 02  Kin-1
 03  Kout3
 04  Kin-2
 05  2
 06  Kin+3
 07  Kout2
 08  x>0
 09  Kout1

Usage example:

2062 Kin1 Kin2 1 Kin-2 3 Kin3 P1

[dprtl]

Voici ma première proposition pour Casio PB1000 (et compatibles), basée sur les formules de Xerxes :

Code : Tout sélectionner

10 INPUT"k";K
20 INPUT"n";N
30 R=N
40 IF K<0 THEN 70
50 FOR I=1 TO K:GOSUB 90:NEXT
60 GOTO 80
70 FOR I=1 TO(-K):GOSUB 100:NEXT
80 PRINT"f^";K;"(";N;") =";R:GOTO 10
90 R=INT(ROUND(R+SQR(R),-1)):RETURN
100 R=R-INT(SQR(R)):RETURN

Exemple d'utilisation :

Code : Tout sélectionner

k?-2017
n?4072325
f ^-2017 ( 4072325 ) = 1019091

[C.Ret]

Code : Tout sélectionner

  readln(n);
  s:=1;
  i:=1;
  repeat
    inc(n);
    inc(i,2);
    inc(s,i);
  until s>n;
  writeln(n);

I like this algorithm ! Computing squares without multiplication.
At first look, I was afraid a mistake.
But no, it is ok, the trick is to use variable n for entry and for the result !

J'aime cet algorithme! Calculer les carrés sans multiplication.
À première vue, j'avais peur d'une erreur.
Mais non, c'est bon, l'astuce consiste à utiliser la variable n pour l'entrée et pour le résultat!



En effet, les carrés n² et (n+1)² sont séparés de 2n-1 valeurs, la variable i est incrémentée de 2 en 2 et permet donc par simple addition de calculer le carré suivant à chaque tour de roue.

Si les carrés parfaits n'existaient pas, F(n) serait simplement la suite des entiers naturels.
Comme il y a de temps en temps un carré parfait, f(n) prend du "retard" par rapport à la suite des entiers. On perd une position à chaque fois que l'on dépasse un carré.

Voici en BASIC pour SHARP PC-1211 l'algorithme de Xerxes :

Code : Tout sélectionner

1 "F" AREAD N    
  : S=1            REM *** S pour 'square'  c'est à dire les carrrés 1² 2² 3² ... 
  , I=1            REM *** I  pour 'increment' c'est en fait la valeur 2n-1 qui permet de passer d'un carré au suivant
2 N=N+1            REM *** On incrémente f(n) de 1 car on vient de dépasser le carré  s (initialement 1²)
  ,I=I+2           REM *** On incrémente l"increment' de deux pour calculer le carré suivant
  ,S=S+I           REM *** On calcule le carré suivant
  :IF N>S GOTO 2   REM *** On boucle tant que n<s
3 PRINT N          REM *** On affiche le résultat
  :END             
 

et l'équivalent pour HP-41C:

Code : Tout sélectionner

001 LBL "MPO79"  1  ENTER^  ENTER^  R^
003 LBL 00  R^  ST+ T  ST+ T   ST+ Y   R^  ST+ T  RDN  RDN   x>y?  GTO 00
017 RTN

[dprtl]

Comme c'est un MPO, voici une autre version pour PB-1000, très proche de la précédente, mais dans laquelle j'ai supprimé vite fait tous les caractères et espaces inutiles :

Code : Tout sélectionner

10 INPUT"k";K:INPUT"n";N:R=N:IFK<0 THEN30
20 FORI=1 TOK:R=ROUND(R+SQR(R),-1):NEXT:GOTO40
30 FORI=1 TO-K:R=R-INTSQR(R):NEXT
40 PRINT"f^";K;"(";N;") =";R:GOTO10

Le résultat ci-dessous est calculé en 1 min 7 s environ :

Code : Tout sélectionner

f^-2017 ( 4072325 ) = 1019091

[Xerxes]

Shorter version for the FX-180P using the square function:

Code : Tout sélectionner

 01  1
 02  M+
 03  Kin+1
 04  Kout1
 05  x^2
 06  x<=M
 07  MR

Usage example:

2017 Min 1 Kin1 P1

[C.Ret]

Code : Tout sélectionner

10 INPUT"k";K:INPUT"n";N:R=N:IFK<0 THEN30
20 FORI=1 TOK:R=ROUND(R+SQR(R),-1):NEXT:GOTO40
30 FORI=1 TO-K:R=R-INTSQR(R):NEXT
40 PRINT"f^";K;"(";N;") =";R:GOTO10

Le résultat ci-dessous est calculé en 1 min 7 s environ :

Code : Tout sélectionner

f^-2017 ( 4072325 ) = 1019091

J'aime bien cette approche, où l'on entre une valeur positive ou négative de k en fonction du sens de la transformation et du nombre d'appels composés souhaité. Du coup, j'ai modifié mes codes initiaux car j'avais envisagé la chose différemment et cette dernière façon de faire est tellement plus pratique.

Pour SHARP PC-1211, cela donne le code suivant:

Code : Tout sélectionner

1:" " AREAD N:F=N,K=1:INPUT "^";K      // Saisie de n et facultativement de k
2:FOR E=1 TO ABS K                     // Boucle pour chaque composée
    R=INT √F                           //   Selon le sens du caclul (càd le signe de k)
  , F=F+R*SGN K+(√F>.5+R)*(K>0)        //   ajoute ROUND(SQRT(n)) ou retire FLOOR(SQRT(n))      
 :NEXT E                               //   avec ROUND effectué en deux opérations car non disponible   
3:"=" PRINT F;" = F^";K;"(";N;")":END  // Affichage peut être simplifié

Mais, c'est toujours une résolution à l'aide d'une boucle. Et cela rend les déterminations d'autant plus longues que les composées sont nombreuses.


En étudiant de plus près certains cas de composées, je me rends compte qu'il y a une évolution régulière des valeurs f^k(n) en fonction de k.

D'où le code ci-dessus 'vectorialisé', c'est à dire sans boucle de calculs:

Code : Tout sélectionner

1:" " AREAD N:K=1:INPUT "^";K                                         // Saisie de n et facultativement de k
2:R=INT √N,R=R+(√N>.5+R),E=INT .5K,F=E-R+INT √N+(K>2E),F=N+NK+EF      // Calculs de R,E puis F
3:"=" PRINT F;" =^";K;N:END                                           // Affichage peut être affiné

Variables :

Code : Tout sélectionner

N  Entier anturel n saisi par l'utilisateur à l'aide de [shft][spc] en mode DEF

K  Nombre et sens des composé : k positif détermine les composées f^k(n) 
                                k négatif détermine les composées réciproques
   Par défaut K=1, le code détermine donc x = f(n) le n-ième entier non carré parfait.

R  Entier le plus proche de √n. Le SHARP PC1211 n'ayant pas de fonction ROUND, calcul en deux étapes   
E  Partie entière de K/2 qui permettra de calculer FLOOR(k/2) et CEIL(k/2) sur ce SHARP dépourvu
F  Résultat du calcul  F := n + k.ROUND(√n) + FLOOR(k/2).CEIL(k/2) 

Exemples d'applications:

Code : Tout sélectionner

Calcul de f(1e6):
 Affichage               Saisie             Remarques
  3:"="PRINT F;"=^";K;N  [mode][mode]            Après sasiie du code, passer en mode DEF
  >                      [Exp][ 6 ][shft][spc]   Saisie de n = 1 000 000
  ^_                     [ENTER]                 Pas de composé, taper ENTER directement
  1001000.=^1.1000000.                           Le millionième entier non carré parfait est 1001000.
                                                 il y a en effet 1000 carrés parfaits entre 1 et 10^6 !  

25/05/2017: EDIT : Attention ce qui suit est erroné, il faut lire 1629 et non pas 1627 - voir posts suivants :

Code : Tout sélectionner

Trouver le rang de 1669:
 Affichage               Saisie             Remarques
                         [cl]                    Efface affichage.
  >                      [1][6][6][9][shft][spc] Sasir 1669
  ^_                     [-][1][ENTER]           Demander la réciproque c.à.d. la composé (-1) fois
  1627.=^-1.1669.                                1669 est le 1627-ième entier non carré parfait
                         [shft][spc][ENTER]      Vérifions cette allégation
  1669.=^1.1627.                                 Le 1627-ième entier non carré est donc bien 1669

25/05/2017: EDIT : Attention ce qui suit est erroné, la réponse à Q2 est qu'il n'y a pas de solution - voir posts suivants :

Code : Tout sélectionner

Q1:1406[shft][spc]15[ENTER]        2017.=^15.1046            d'où f^15(1406) = 2017.
Q2:2017[shft][spc]-20[ENTER]       1207.=^-20.2017.          n = 1207
Q3:4867[shft][spc]-50[ENTER]       1967.=^-50.4867.          n = 1967
Q4:4072325[shft][spc]-2017[ENTER]  1019091.=^-2017.4072325.  n = 1019091

Faire attention cependant, ce code peu donner avec les composées inverses (k négatifs) des résultats aberrants.
En particulier dans le cas de détermination du rang d'un entier carré parfait:

Code : Tout sélectionner

25[shft]-1[ENTER] affiche 20.=^-1.25.

Ce qui ne signifie pas que 25 est le 20-ième entier non carré parfait. En réalité 25 est un carré parfait qui se trouve entre f(20) et f(21) (soit entre 24 et 26 ) !!
Il doit être possible d'ajouter un test en ligne 1: pour détecter une telle situation ?


Voici une version équivalente pour HP41C.

Code : Tout sélectionner

001 LBL"MPO79"    011 2             021 ST- Y         031 ST+
002 FIX 0         012 R^            022 ST- Y         032 *
003 RCL Y         013 ST+ L         023 x<>y          033 +
004 SQRT          014 X<> L         024 ABS           034 END
005 INT           015 RDN           025 ST+ Z
006 LastX         016 /             026 X<> Y         Usage :
007 RND           017 ST* L         027 x<0?          1406 [ENTER^] 15 [XEQ][Alpha]MPO79[Alpha]
008 -             018 X<> L         028 DSE Y          affiche 2017.
009 x<>y          019 LastX         029 ABS 
010 ST* L         020 INT	    030 RDN 

La difficulté de retranscription, rencontrée en traduisant le code BASIC, provient du fait du comportement diffèrent de l'instruction INT sur HP-41C par rapport aux pockets SHARP

En effet, INT -7.5 donne avec le SHARP la valeur entière -8. L'instruction INT est donc bien une sorte de fonction FLOOR.
Alors que l'INT de l'HP-41C donne la partie entière de l'argument, ce qui signifie que 7.5 CHS INT donne -7

[Xerxes]

Faire attention cependant, ce code peu donner avec les composées inverses (k négatifs) des résultats aberrants.
En particulier dans le cas de détermination du rang d'un entier carré parfait:

Code : Tout sélectionner

25[shft]-1[ENTER] affiche 20.=^-1.25.

Ce qui ne signifie pas que 25 est le 20-ième entier non carré parfait. En réalité 25 est un carré parfait qui se trouve entre f(20) et f(21) (soit entre 24 et 26 ) !!
Il doit être possible d'ajouter un test en ligne 1: pour détecter une telle situation ?

It's obvious that a square number cannot be the argument of the inverse of f(n). But this problem can also happen, while repeatedly
calling the function. An example is to start with 30, that is the 25th element of the sequence. After the first iteration we have a
logical problem of how to continue now, without using an illegal argument.

f^50(2017)=4867

f^-50(4867)=2017 is correct, because there is no square number in the interim results.

f^-20(2017)=1225 is wrong, because one of the interim results is 37^2.

f^20(1225)=2015 and f^20(1226)=2026 so there cannot be a solution for f^-20(2017).

My conclusion is, that the occurrence of a square interim result, while calculating f^k(n) for negative k, means there is no solution.

[C.Ret]

Good catch !

So the answer to Q2 is NO, there is no integer n that verify f^20(n) = 2017.

I have to modify my codes in order to track for perfect square in interim results when using k<0 in sequence.

Par exemple en insérant un test :

Code : Tout sélectionner

1:" " AREAD N:F=N,K=1:INPUT "^";K    // Saisie de n et facultativement de k
2:FOR E=1 TO ABS K:                  // Boucle pour chaque composée
    R=INT √F:                        //   Selon le sens du caclul (càd le signe de k)
    IF K<0 IF RR=F BEEP 1:STOP       //   Vérifie validité des résultats intermédiaires
3:  F=F+R*SGN K+(√F>.5+R)*(K>0):     //   ajoute ROUND(SQRT(n)) ou retire FLOOR(SQRT(n))
 :NEXT E                             //   avec ROUND effectué en deux opérations car non disponible   
4:"=" PRINT F;" = F^";K;"(";N;")"    // Affichage peut être simplifié

Ce code s'arrête lors de la tentative de déterminer le rang d'un carré:
Ainsi lors de la détermination de f^-20(2017), il s'arrête après avoir calculé 17 valeurs intermédiaires aboutissant sur 1369 = 37*37.
2017<-1973<-1929<-1886<-1843<-1801<-1759<-1718<-1677<-1637<-1597<-1558<-1519<-1481<-1443<-1406<-1369 /!\

On vérifie et on corrige donc les réponses aux questions posées:

Code : Tout sélectionner

Q1:1406[shft][spc]15[ENTER]        2017.=^15.1046            d'où f^15(1406) = 2017.
Q2:2017[shft][spc]-20[ENTER]       BREAK AT 2:                 Il n'existe pas d'entier n tel que f^20(n)=2017.   
Q3:4867[shft][spc]-50[ENTER]       1967.=^-50.4867.          n = 1967
Q4:4072325[shft][spc]-2017[ENTER]  1019091.=^-2017.4072325.  n = 1019091

[Xerxes]

I have to modify my codes in order to track for perfect square in interim results when using k<0 in sequence.

If it's possible to calculate f^k(n) without iterations, it would be more practical to check the result for k<0 with
f^-k(result) for plausibility. I wrote "if possible", because I noticed a wrong result for Q3 that should be 2017.

[Xerxes]

Calculating f^k(n) for k>0:

d=floor(sqrt(n)*2)-1
if frac(sqrt(n))=0 then d=d-1
n=n+round((k^2+2dk+2k)/4)


Calculating f^k(n) for k<0:

d=floor(sqrt(n)*2)-1
n=n-round((-k^2+2dk+2k)/4)

[C.Ret]

Correct, or both k<0 and k>O as I have done:

D'où le code ci-dessus 'vectorialisé', c'est à dire sans boucle de calculs:

Code : Tout sélectionner

1:" " AREAD N:K=1:INPUT "^";K                                         // Saisie de n et facultativement de k
2:R=INT √N,R=R+(√N>.5+R),E=INT .5K,F=E-R+INT √N+(K>2E),F=N+NK+EF      // Calculs de R,E puis F
3:"=" PRINT F;" =^";K;N:END                                           // Affichage peut être affiné

Variables :

Code : Tout sélectionner

N  Entier naturel n saisi par l'utilisateur à l'aide de [shft][spc] en mode DEF

K  Nombre et sens des composé : k positif détermine les composées f^k(n) 
                                k négatif détermine les composées réciproques
   Par défaut K=1, le code détermine donc x = f(n) le n-ième entier non carré parfait.

R  Entier le plus proche de √n. Le SHARP PC1211 n'ayant pas de fonction ROUND, calcul en deux étapes   
E  Partie entière de K/2 qui permettra de calculer FLOOR(k/2) et CEIL(k/2) sur ce SHARP dépourvu
F  Résultat du calcul  F := n + k.ROUND(√n) + FLOOR(k/2).CEIL(k/2) 

[...]
Voici une version équivalente pour HP41C.

Code : Tout sélectionner

001 LBL"MPO79"    011 2             021 ST- Y         031 ST+
002 FIX 0         012 R^            022 ST- Y         032 *
003 RCL Y         013 ST+ L         023 x<>y          033 +
004 SQRT          014 X<> L         024 ABS           034 END
005 INT           015 RDN           025 ST+ Z
006 LastX         016 /             026 X<> Y         Usage :
007 RND           017 ST* L         027 x<0?          1406 [ENTER^] 15 [XEQ][Alpha]MPO79[Alpha]
008 -             018 X<> L         028 DSE Y          affiche 2017.
009 x<>y          019 LastX         029 ABS 
010 ST* L         020 INT	    030 RDN 

Code : Tout sélectionner

  n   k     √n      r   R  i  m   e    n  +  k*R  + e*( e-i+m)          
1406, 14  37.49666  37 37  0  0   7  1406 + 14*37 + 7*  7      = 1924  + 49 = 1973
1406, 15  37.49666  37 37  0  1   7  1406 + 15*37 + 7*( 7  +1) = 1961  + 56 = 2017
1407, 14  37.50999  37 38  1  0   7  1407 + 14*38 + 7*( 7-1  ) = 1939  + 42 = 1981
1407, 15  37.50999  37 38  1  1   7  1407 + 15*38 + 7*( 7-1+1) = 1977  + 49 = 2026
1973,-14  44.41846  44 44  0  0  -7  1973 - 14*44 - 7*  7      = 1357  + 49 = 1406
2026,-15  45.01111  45 45  0  1  -8  2026 - 15*45 - 8*(-8  +1) = 1351  + 56 = 1407
1981,-14  44.50842  44 45  1  0  -7  1981 - 14*45 - 7*(-7-1  ) = 1351  + 56 = 1407
2017,-15  44.91102  44 45  1  1  -8  2017 - 15*45 - 8*(-8-1+1) = 1342  + 64 = 1406

i = R-r 
m= k MOD 2