/*
test7.c
L'objectif est d'avoir une soubroutine en C
qui peut être accédé depuis le langage Python.

Pour compiler depuis un Terminal :
gcc -shared -o test7.so -fPIC test7.c

Il faut indiquer une fois pour toute où se trouve la librairie "Python.h"
Pour ajouter des répertoires de "include" (.h) :
C_INCLUDE_PATH="/usr/include/python3.5"
export C_INCLUDE_PATH

Depuis Python, exécuter le script test7.py qui est donné dans un autre fichier.

Pour les types que Python reconnait et peut convertir en C :
https://www.tutorialspoint.com/python/python_further_extensions.htm
***********************************************************************/

#include "Python.h"

double local_sommeInv(long lNbMin, long lNbMax) {
//===============================================
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
long lNbr;
double vRep;

// Le résultat de la somme est nettement plus précis si on
// commence par additionner les plus petits nombres
vRep = 0.0;
for (lNbr=lNbMax-1; lNbr>=lNbMin; lNbr--) vRep = vRep + 1.0 / lNbr;

return vRep;
} // local_sommeInv

static PyObject *test7_sommeInv(PyObject *self, PyObject *args) {
//===============================================================
// Définit la fonction en C qui sera appelée depuis Python
// PyObject *x_obj;  // autre possibilité.
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
// lNbMin et lNbMax sont les deux arguments.
long lNbMin;
long lNbMax;

if (!PyArg_ParseTuple(args, "ll", &lNbMin, &lNbMax)) return NULL; 
return PyFloat_FromDouble(local_sommeInv(lNbMin, lNbMax));
// https://docs.python.org/3.3/c-api/float.html
} // test7_sommeInv

double local_sommeInv9(long lNbMin, long lNbMax) {
//================================================
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
// Tous les nombres ayant le chiffre 9 dans leur écriture décimale
// sont enlevés de la somme.
long lNbr;
long lDixPow;
long lDixPowMax;
double vRep;

// Calcule la plus grande puissance de 10 qui est
// de justesse plus grande ou égale à lNbMax
lDixPowMax = 1;
while (lDixPowMax < lNbMax) lDixPowMax *=10;

// Somme en commençant par les plus petits nombres.
vRep = 0.0;
lNbr=lNbMax;
while (lNbr > lNbMin) {
  lNbr--;
  
  // Test si lNbr contient un chiffre "9".
  // Si "oui", passe au nombre précédent qui ne contient par de "9".
  lDixPow = lDixPowMax;
  while (lDixPow >= 10) {
    lDixPow /= 10; 
    if (lNbr/lDixPow % 10 == 9) lNbr -= lDixPow;
    } // while

  vRep += 1.0 / lNbr;
  } // while

// Cas exeptionnelle où  lNbMin est un nombre qui contient le chiffre "9"
if (lNbMin < lNbMax) { // ça c'est habituelle.
  lDixPow = lDixPowMax;
  while (lDixPow > 1) {
    lDixPow /= 10; 
    if (lNbMin/lDixPow % 10 == 9) {
      vRep -= 1.0 / lNbMin;
      break;
      }
    } // while
  } // if
  
// Pour des tests, cas où lNbMax < lNbMin
// Somme commençant pas les plus grands nombres. C'est moins précis.
lNbr=lNbMax;
while (lNbr<lNbMin) {
  vRep += 1.0 / lNbr;
  lNbr++;

  // Test si lNbr contient un chiffre "9".
  // Si "oui", passe au nombre suivant qui ne contient par de "9".
  if (lNbr % 10 == 9) { lNbr ++;
    if (lNbr/10 % 10 == 9) { lNbr += 10;
      if (lNbr/100 % 10 == 9) { lNbr += 100;
        // J'ai laissé quelques "if" pour comprendre la structures des tests.
        lDixPow = 1000;
        while (lNbr/lDixPow % 10 == 9) {
          lNbr += lDixPow;
          lDixPow *=10;
          }
        }      
      }
    }
  } // while

return vRep;
} // local_sommeInv9

static PyObject *test7_sommeInv9(PyObject *self, PyObject *args) {
//===============================================================
// Définit la fonction en C qui sera appelée depuis Python
// PyObject *x_obj;  // autre possibilité.
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
// Tous les nombres ayant le chiffre 9 dans leur écriture décimale
// sont enlevés de la somme.
// Il y a deux paramètrs en argument : lNbMin, lNbMax
long lNbMin;
long lNbMax;

if (!PyArg_ParseTuple(args, "ll", &lNbMin, &lNbMax)) return NULL; 
return PyFloat_FromDouble(local_sommeInv9(lNbMin, lNbMax));
// https://docs.python.org/3.3/c-api/float.html
} // test7_sommeInv9

double local_sommeInvx(long lNbMin, long lNbMax, long lc) {
//=========================================================
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
long lNbr;
long lDixPow;
long lDixPowMax;
double vRep;

//if (lc == 9) return local_sommeInv9(lNbMin, lNbMax);

// Calcule la plus grande puissance de 10 qui est
// de justesse plus grande ou égale à lNbMax
lDixPowMax = 1;
while (lDixPowMax < lNbMax) lDixPowMax *=10;

// Somme en commençant par les plus petits nombres.
vRep = 0.0;
lNbr=lNbMax;
while (lNbr > lNbMin) {
  lNbr--;
  
  // Test si lNbr contient un chiffre "lc".
  // Si "oui", passe au nombre précédent qui ne contient par de "lc".
  lDixPow = lDixPowMax;
  while (lDixPow >= 10) {
    lDixPow /= 10; 
    if (lNbr/lDixPow % 10 == lc) lNbr -= lDixPow;
    } // while

  vRep += 1.0 / lNbr;
  } // while

// Cas exeptionnelle où  lNbMin est un nombre qui contient le chiffre "lc"
if (lNbMin < lNbMax) { // ça c'est habituelle.
  lDixPow = lDixPowMax;
  while (lDixPow > 1) {
    lDixPow /= 10; 
    if (lNbMin/lDixPow % 10 == lc) {
      vRep -= 1.0 / lNbMin;
      break;
      }
    } // while
  } // if

// Pour des tests, cas où lNbMax < lNbMin
// Somme commençant pas les plus grands nombres. C'est moins précis.
lNbr=lNbMax;
while (lNbr<lNbMin) {
  vRep = vRep + 1.0 / lNbr;
  lNbr++;

  // Test si lNbr contient un chiffre "lc".
  // Si "oui", passe au nombre suivant qui ne contient par de "lc".
  lDixPow = 1;
  while (lNbr/lDixPow >= 1) {
    if (lNbr/lDixPow % 10 == lc) {
      lNbr += lDixPow;
      if (lc != 9) break; // sort de la boucle while
      }
    lDixPow *=10;
    } // while
  } // while

return vRep;
} // local_sommeInvx

static PyObject *test7_sommeInvx(PyObject *self, PyObject *args) {
//===============================================================
// Définit la fonction en C qui sera appelée depuis Python
// PyObject *x_obj;  // autre possibilité.
// Somme des inverses des nombres entiers de lNbMin à lNbMax-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Il y a trois paramètrs en argument : lNbMin, lNbMax et lc
long lNbMin;
long lNbMax;
long lc;

if (!PyArg_ParseTuple(args, "lll", &lNbMin, &lNbMax, &lc)) return NULL; 
return PyFloat_FromDouble(local_sommeInvx(lNbMin, lNbMax, lc));
// https://docs.python.org/3.3/c-api/float.html
} // test7_sommeInvx


double local_pow(double vv, int rr) {
//===================================
// mise à la puissance, pas efficace.
// Retourne vv^rr
double vRep;

vRep = 1.0;
while (rr > 0) {vRep *= vv; rr--;}
return vRep;
} // local_pow

double local_Sc_pow(long li, long lj, long lc) {
//==============================================
// Somme des inverses des puissances lj des nombres entiers de 10^li à (10*10^li)-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Ici, le calcul est fait explicitement.
// Dans la fonction suivante, qui utilise celle-ci, le calcul
// est fait plus efficacement, si li > 4.
long lNbMin;
long lNbMax;
long lNbr;
long lDixPow;
double vRep;

// Calcule : lNbMin = 10^li
lNbMin = 1;
while (li > 0) { lNbMin = lNbMin*10, li--; }
lNbMax = 10*lNbMin;

// Somme en commençant par les plus petits nombres.
vRep = 0.0;

if (lc >= 0) { // Cas normal
  lNbr=lNbMax;
  while (lNbr>lNbMin) {
    lNbr--;
    
    // Test si lNbr contient un chiffre "lc".
    // Si "oui", passe au nombre précédent qui ne contient par de "lc".
    lDixPow = lNbMax; // lNbMax = 10^(li+1=
    while (lDixPow >= 10) {
      lDixPow /= 10; 
      if (lNbr/lDixPow % 10 == lc) lNbr -= lDixPow;
      } // while

    vRep += local_pow(1.0/lNbr, lj);
    } // while
  } // if
else { // Cas anormal, juste pour des tests.
  // Pour des tests, la somme commence par les plus grands nombres.
  lc = -lc;
  lNbr=lNbMin;
  while (lNbr<lNbMax) {
    vRep += local_pow(1.0/lNbr, lj);
    lNbr++;

    // Test si lNbr contient un chiffre "lc".
    // Si "oui", passe au nombre suivant qui ne contient par de "lc".
    lDixPow = 1;
    while (lNbr/lDixPow >= 1) {
      if (lNbr/lDixPow % 10 == lc) {
        lNbr += lDixPow;
        if (lc < 9) break; // sort de la boucle while
        }
      lDixPow *=10;
      } // while
    } // while
}

return vRep;
} // local_Sc_pow

double local_Sc(long li, long lj, long lc, int nFlag) {
//=====================================================
// Si nFlag == 0 :
// Somme des inverses des puissances lj des nombres entiers de 10^li à (10*10^li)-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Si nFlag == 1 :
// Somme des inverses des puissances lj des nombres entiers de 1 à (10*10^li)-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Le code n'est pas optimisé, mais il est suffisemment rapide.
int jMax = 20; // Valeur max pour l'indice "j"
int jj;
int ii;
int rr;
int pp;
double vC_r_j;
double avD[31]; // = dc(r) / 10^r
double avS[31];
double avS2[31];
double avTemp[31]; // Pour mémoriser des valeurs à sommer
double vSum;
double vSommeTot;  // Somme de 1 à (10*10^li)-1 des ..., c.f. nFlag == 1

if (li <= 3) {
  avS[lj] = local_Sc_pow(li, lj, lc);
  vSommeTot =  avS[lj];
  while (li > 0) { li--; vSommeTot += local_Sc_pow(li, lj, lc); }
  }
else {
  // C'est ici que l'on utilise la théorie qui permet d'optimiser les calculs.
  // Initialisation du vecteur : ( Sc(3,1) ; Sc(3,2) ; Sc(3,3) ; ... ; Sc(3,jMax))
  for (jj=1; jj<=jMax; jj++) avS[jj] = local_Sc_pow(3, jj, lc);
  vSommeTot = avS[lj] + local_Sc_pow(0, lj, lc) + local_Sc_pow(1, lj, lc) + local_Sc_pow(2, lj, lc);
  
  // Calculs de  dc(r)/10^r, pour r=0..jMax
  for (rr=0; rr<=jMax; rr++) {
    vSum = 0;
    for (pp=9; pp>=0; pp--) if (pp != lc) vSum += local_pow(0.1*pp, rr);
    avD[rr] = vSum;
    }

  // Par réccurence, calcule le vecteur : ( Sc(li,1) ; Sc(li,2) ; Sc(li,3) ; ... ; Sc(li,jMax))
  for (ii=4; ii<=li; ii++) {
    for (jj=1; jj<=jMax; jj++) {
      // Calcul de S(ii+1, jj)
      vSum = 0.0;
      vC_r_j = 1.0;
      for (rr=0; rr<=jMax-jj; rr++) {
        avTemp[rr] = vC_r_j * avD[rr] * avS[rr+jj];
        vC_r_j *= -1.0*(rr+jj) / (rr+1); // Calcul c(r+1, j)
        }

      // Pour sommer du plus petit au plus grand.
      for (rr=jMax-jj; rr>=0; rr--) vSum += avTemp[rr];
        
      avS2[jj] = vSum * local_pow(0.1, jj);
      }

    // Copie le vecteur  avS2  dans  avS
    for (jj=1; jj<=jMax; jj++) avS[jj] = avS2[jj];

    vSommeTot+=avS[lj];
    }
  } // else
  
if (nFlag == 0) return avS[lj]; // = Sc(li, lj)
else return vSommeTot;
} // local_Sc

static PyObject *test7_Sc(PyObject *self, PyObject *args) {
//=========================================================
// Définit la fonction en C qui sera appelée depuis Python
// PyObject *x_obj;  // autre possibilité.
// Somme des inverses puissances lj des nombres entiers de 10^li à (10*10^li)-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Il y a trois paramètrs en argument : li, lj et lc
long li;
long lj;
long lc;

if (!PyArg_ParseTuple(args, "lll", &li, &lj, &lc)) return NULL; 
return PyFloat_FromDouble(local_Sc(li, lj, lc, 0));
// https://docs.python.org/3.3/c-api/float.html
} // test7_Sc

static PyObject *test7_ScTot(PyObject *self, PyObject *args) {
//============================================================
// Définit la fonction en C qui sera appelée depuis Python
// PyObject *x_obj;  // autre possibilité.
// Somme des inverses puissances lj des nombres entiers de 1 à (10*10^li)-1
// Tous les nombres ayant le chiffre lc dans leur écriture décimale
// sont enlevés de la somme.
// Il y a trois paramètrs en argument : li, lj et lc
long li;
long lj;
long lc;

if (!PyArg_ParseTuple(args, "lll", &li, &lj, &lc)) return NULL; 
return PyFloat_FromDouble(local_Sc(li, lj, lc, 1));
// https://docs.python.org/3.3/c-api/float.html
} // test7_ScTot

// Définit l'aide associée à la fonction.
#define AIDE "sommeInv prend deux nombres entiers Nmin et Nmax en entrée\n\
et retourne la somme des inverses des entiers de Nmin à (Nmax-1)."

// Définit l'aide associée à la fonction.
#define AIDE9 "sommeInv9 prend deux nombres entiers Nmin et Nmax en entrée\n\
et retourne la somme des inverses des entiers de Nmin à (Nmax-1)\n\
dans laquelle tous les nombres contenant le chiffre 9 ont été éliminés.\n\
Voici des nombres éliminés :\n\
9 ; 19 ; 29 ; 9xxx ; 19xxx ; 29xxx ; 999999.\n\
Nmin ne doit pas contenir de chiffre 9."

// Définit l'aide associée à la fonction.
#define AIDEx "sommeInvx prend trois nombres entiers Nmin, Nmax et Nc en entrée\n\
et retourne la somme des inverses des entiers de Nmin à (Nmax-1)\n\
dans laquelle tous les nombres contenant le chiffre Nc ont été enliminés\n\
c.f. simmeInv9 pour l'analogie.\n\
Nmin ne doit pas contenir de chiffre Nc."

#define AIDESc "Sc prend trois nombres entiers li, lj et lc en entrée\n\
et retourne la somme des inverses des entiers puissance lj de 10^li à (10*10^li)-1.\n\
dans laquelle tous les nombres contenant le chiffre lc."

#define AIDEScTot "Sc prend trois nombres entiers li, lj et lc en entrée\n\
et retourne la somme des inverses des entiers puissance lj de 1 à (10*10^li)-1.\n\
dans laquelle tous les nombres contenant le chiffre lc."

static PyMethodDef test7Methods[] = {
//===================================
{"sommeInv",  test7_sommeInv, METH_VARARGS, AIDE}, // Est un texte d'aide. Obtenu avec : print(test7.sommeInv.__doc__)
{"sommeInv9",  test7_sommeInv9, METH_VARARGS, AIDE9}, // Est un texte d'aide. Obtenu avec : print(test7.sommeInv9.__doc__)
{"sommeInvx",  test7_sommeInvx, METH_VARARGS, AIDEx}, // Est un texte d'aide. Obtenu avec : print(test7.sommeInvx.__doc__)
{"Sc",  test7_Sc, METH_VARARGS, AIDESc}, // Est un texte d'aide. Obtenu avec : print(test7.Sc.__doc__)
{"ScTot",  test7_ScTot, METH_VARARGS, AIDEScTot}, // Est un texte d'aide. Obtenu avec : print(test7.ScTot.__doc__)
{NULL, NULL, 0, NULL}  // Sentinel 
};

static struct PyModuleDef test7module = {
//=======================================
// Est une structure nécessaire pour Python ???
  PyModuleDef_HEAD_INIT,
  "test7",   /* name of module */
  NULL,  //test7_doc, /* module documentation, may be NULL */
  -1,       /* size of per-interpreter state of the module,
             or -1 if the module keeps state in global variables. */
  test7Methods
};

PyMODINIT_FUNC PyInit_test7(void) {
//=================================
// Fonction de création du module
return PyModule_Create(&test7module);
} // PyInit_test7

int main(int argc, char *argv[]) {
//================================
// Fonction principale, qui sera exécutée au chargement du module
wchar_t *program = Py_DecodeLocale(argv[0], NULL);
if (program == NULL) {
  fprintf(stderr, "Fatal error: cannot decode argv[0]\n");
  exit(1);
  }

/* Add a built-in module, before Py_Initialize */
PyImport_AppendInittab("test7", PyInit_test7);

/* Pass argv[0] to the Python interpreter */
Py_SetProgramName(program);

/* Initialize the Python interpreter.  Required. */
Py_Initialize();

/* Optionally import the module; alternatively,
   import can be deferred until the embedded script
   imports it. */
PyImport_ImportModule("test7");

PyMem_RawFree(program);
return 0;
} // main