/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        R E C H E R C H E   D E   P O I N T S   R A T I O N N E L S   S U R   U N E   C O U R B E   E L L I P T I Q U E  :         */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/CourbesElliptiques.11$c' :                                                                                 */
/*                                                                                                                                   */
/*                    Jean-Francois Colonna (LACTAMME, 20211123141952).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  "INCLUDES.01.I"

extern    double floor();

#include  "CourbesElliptiques.01.I"
                                        /* Definition du polynome :                                                                  */
                                        /*                                                                                           */
                                        /*                     3      2      1      0                                                */
                                        /*                  A.x  + B.x  + C.x  + D.x                                                 */
                                        /*                                                                                           */

#define   XPN                                                                                                                           \
                    (-1.0)
#define   XPD                                                                                                                           \
                    (+1.0)
#define   YPN                                                                                                                           \
                    (+1.0)
#define   YPD                                                                                                                           \
                    (+1.0)
                                        /* Point rationnel 'P' = {-1/1,+1/1} ('v $xtc/CourbesElliptiques.01$c').                     */
#define   XQN                                                                                                                           \
                    (+1.0)
#define   XQD                                                                                                                           \
                    (+4.0)
#define   YQN                                                                                                                           \
                    (+7.0)
#define   YQD                                                                                                                           \
                    (+8.0)
                                        /* Point rationnel 'Q' = {+1/4,+7/8} ('v $xtc/CourbesElliptiques.01$c').                     */

#define   EPSILON                                                                                                                       \
                    0.000001

#define   CHECK_POINT(x,y)                                                                                                              \
                    {                                                                                                                   \
                    double    Gdroite=y;                                                                                                \
                    double    Ddroite=PQa*x+PQb;                                                                                        \
                    double    Gcourbe=y*y;                                                                                              \
                    double    Dcourbe=A*(x*x*x)+B*(x*x)+C*(x)+D;                                                                        \
                                                                                                                                        \
                    printf("Le point {%+f,%+f} ",x,y);                                                                                  \
                                                                                                                                        \
                    if        (ABSO(Gdroite-Ddroite) < EPSILON)                                                                         \
                              {                                                                                                         \
                              printf("appartient a la droite");                                                                         \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              printf("n'appartient pas a la droite %f.x+%f (%f#%f)",PQa,PQb,Gdroite,Ddroite);                           \
                              }                                                                                                         \
                                                                                                                                        \
                    printf(" et ");                                                                                                     \
                                                                                                                                        \
                    if        (ABSO(Gcourbe-Dcourbe) < EPSILON)                                                                         \
                              {                                                                                                         \
                              printf("appartient a la courbe elliptique");                                                              \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              printf("n'appartient pas a la courbe elliptique (%f#%f)",Gcourbe,Dcourbe);                                \
                              }                                                                                                         \
                                                                                                                                        \
                    printf(".\n");                                                                                                      \
                    }

#define   PLUS_GRAND_DENOMINATEUR_POSSIBLE                                                                                              \
                    1000000

int       An,Ad;

ApproximationRationnelle(x)
double    x;
          {
          int       iterer=VRAI;

          double    a0=1,a1,a;
          double    b0=0,b1=1,b=1;

          double    en,q;

          a1=floor(x);
          en=x-a1;

          printf("\nApproximations rationnelles de %f :\n\n",x);

          while     (iterer == VRAI)
                    {
                    q=floor(1/en);
                    en=(1/en)-q;
                    a=(q*a1)+a0;
                    b=(q*b1)+b0;

                    if        (b <= PLUS_GRAND_DENOMINATEUR_POSSIBLE)
                              {
                              An=a;
                              Ad=b;

                              printf("%d/%d\n",An,Ad);

                              a0=a1;
                              b0=b1;

                              a1=a;
                              b1=b;
                              }
                    else
                              {
                              iterer=FAUX;

                              printf("\nsoit en flottant : %f\n",a/b);
                              }
                    }
          }

void      main()
          {
          double    PQa,PQb;
                                        /* Definition de la droite 'PQ'.                                                             */
          double    XR,YR;
          int       XRn,XRd,YRn,YRd;
                                        /* Coordonnees du troisieme point d'intersection 'R'.                                        */

          PQa=((YQN/YQD)-(YPN/YPD))/((XQN/XQD)-(XPN/XPD));
          PQb=(YPN/YPD)-(PQa*(XPN/XPD));

          XR=-((D-(PQb*PQb))/((XPN/XPD)*(XQN/XQD)));
          YR=PQa*XR+PQb;
                                        /* Equation de la droite 'PQ' :                                                              */
                                        /*                                                                                           */
                                        /*                  y = PQa.x + PQb                                                          */
                                        /*                                                                                           */
                                        /* Equation de la courbe elliptique :                                                        */
                                        /*                                                                                           */
                                        /*                   2      3      2      1      0                                           */
                                        /*                  y  = A.x  + B.x  + C.x  + D.x                                            */
                                        /*                                                                                           */
                                        /* On reporte l'equation de la droite 'PQ' dans celle de la courbe elliptique :              */
                                        /*                                                                                           */
                                        /*                               2      3      2      1      0                               */
                                        /*                  (PQa.x + PQb)  = A.x  + B.x  + C.x  + D.x                                */
                                        /*                                                                                           */
                                        /*                     3         2   2                  1         2   0                      */
                                        /*        [1]       A.x  + (B-PQa ).x  + (C-2.PQa.PQb).x  + (D-PQb ).x  = 0                  */
                                        /*                                                                                           */
                                        /* On cherche le troisieme point d'intersection 'R' en prolongement de 'P' et 'Q'. Les       */
                                        /* abscisses sont solutions de l'equation du troisieme degre :                               */
                                        /*                                                                                           */
                                        /*                  (x-XP).(x-XQ).(x-XR) = 0                                                 */
                                        /*                                                                                           */
                                        /* soit :                                                                                    */
                                        /*                                                                                           */
                                        /*                   3               2                              1               0        */
                                        /*        [2]       x  - (XP+XQ+XR).x  + [(XP.XQ)+(XQ.XR)+(XR.XP)].x  - (XP.XQ.XR).x  = 0    */
                                        /*                                                                                           */
                                        /* En identifiant [1] et [2], on obtient par exemple :                                       */
                                        /*                                                                                           */
                                        /*                                      2                                                    */
                                        /*                  -(XP.XQ.XR) = (D-PQb )                                                   */
                                        /*                                                                                           */
                                        /* d'ou :                                                                                    */
                                        /*                                                                                           */
                                        /*                               2                                                           */
                                        /*                        -(D-PQb )                                                          */
                                        /*                  XR = -----------                                                         */
                                        /*                          XP.XQ                                                            */
                                        /*                                                                                           */
                                        /* puis enfin :                                                                              */
                                        /*                                                                                           */
                                        /*                  YR = PQa.XR + PQb                                                        */
                                        /*                                                                                           */

          printf("XR=%f YR=%f\n",XR,YR);

          ApproximationRationnelle(XR);
          XRn=An;
          XRd=Ad;

          ApproximationRationnelle(YR);
          YRn=An;
          YRd=Ad;

          printf("\n");
          CHECK_POINT(XPN/XPD,YPN/YPD);
          CHECK_POINT(XQN/XQD,YQN/YQD);
          CHECK_POINT(((double)XRn)/((double)XRd),((double)YRn)/((double)YRd));
          }