Quelques calculs...

ça marche

Grâce à une discussion avec Derich (qui est à l'autre bout du monde, merci à mon chat de m'avoir réveillé, ce qui m'a fait profiter du décalage horaire ) j'ai une optimisation de l'algorithme qui simplifie les calculs, réduit les chances d'erreur et surtout réduit le nombre de cas à traiter.

Qu'est-ce qui est constant entre deux mises à jour ? La valeur de taxation et les seuils, entre autres.
Mais ce n'est pas tout ! Comme on aura une fonction de prévision apte à déterminer la date où la population [1] atteint une certaine valeur (pour gérer l'évolution des bonhommes et donc des salaires ou de la consommation, par exemple), je peux aussi considérer que le profil de la courbe est constant : il est croissant, stable ou décroissant. Par exemple, si la population est croissante, il suffira de prévoir quand elle atteindra le seuil [2] pour forcer une mise à jour à ce moment-là et donc avoir des profils constants entre mises à jour.

Voici la formule calculant la population d'une planète entre deux mises à jours t0 et t1, pour la richesse il suffit de remplacer a par b et les s par des r :

time = 30 jours
taux = log(time)

T = t1 - t0
if pop0 == s1
  pop(t) = s1 pour t entre 0 et T
if pop0 < s1
  a = exp(taux*pop0/s1)
  pop(t) = s1*log(a+t)/taux pour t entre 0 et T
if pop0 > s1
  a = exp(taux*(1 - pop0/s0))
  pop(t) = s0*(1 - log(a+t)/taux) pour t entre 0 et T

Donc pour l'algo il n'y a que 3*3 = 9 cas à traiter. Par soucis de clarté j'introduis également deux sous-programmes :

Code:

time = 30 jours
taux = log(time)

// intégrale de log(a+t)/taux pour t entre 0 et T
Ilog(a,T) =
  return ((a + T)*log(a + T) - a*log(a) - T)/taux

// approximation de l'intégrale de log(a+t)*log(b+t)/taux^2 pour t entre 0 et T
Ibilog(a,b,T) =
   return (((a + b)/2 + T)*log(a + T)*log(b + T) - (a + T)*log(a + T) - (b + T)*log(b + T) + 2*T - ((a + b)/2)*log(a)*log(b) + (a)*log(a) + (b)*log(b))/taux^2

imposition(planete,t1) =
   t0 = planete.lastUpdate
   T = t1 - t0
   pop = planete.population
   s0 = planete.dernierSeuilPop
   s1 = planete.seuilPop
   rich = planete.richesse
   r0 = planete.dernierSeuilRich
   r1 = planete.seuilRich
   if pop == s1
      if rich == r1
         impot = s1*r1*T
      if rich < r1
         b = exp(taux*rich/r1)
         impot = s1*r1*Ilog(b,T)
      if rich > r1
         b = exp(taux*(1 - rich/r0))
         impot = s1*r0*(T - Ilog(b,T))
   if pop < s1
      a = exp(taux*pop/s1)
      if rich == r1
         impot = s1*r1*Ilog(a,T)
      if rich < r1
         b = exp(taux*rich/r1)
         impot = s1*r1*Ibilog(a,b,T)
      if rich > r1
         b = exp(taux*(1 - rich/r0))
         impot = s1*r0*(Ilog(a,T) - Ibilog(a,b,T))
   if pop > s1
      a = exp(taux*(1 - pop/s0))
      if rich == r1
         impot = s0*r1*(T - Ilog(a,T))
      if rich < r1
         b = exp(taux*rich/r1)
         impot = s0*r1*(Ilog(b,T) - Ibilog(a,b,T)))
      if rich > r1
         b = exp(taux*(1 - rich/r0))
         impot = s0*r0*(T - Ilog(a,T) - Ilog(b,T) + Ibilog(a,b,T)))
   return impot

Il me reste encore à écrire un vrai programme et à faire des tests pour évaluer l'approximation du bilog et l'influence du nombre de mises à jours sur une période donnée. Il faut de plus que j'explicite la fonction de prévision.

-----------------
[1] ou la richesse, mutatis mutandis et cætera...
[2] en passant, pour la population il suffit de prévoir l'évolution vers le prochain milliard de population, puisque les seuils seront entiers (ou on se débrouillera pour qu'ils le soient, ça m'horripilerait que ma planète reste bloquée à 4 milliards ET DEMI de population)

Commençons par la fonction de prévision.

Il s'agit, après avoir fait une mise à jour en t0, de déterminer la date de la prochaine mise à jour t1. Pour cela on calcule t tel que t1 = t0 + t. Le t en question sera le minimum parmi différentes valeurs possibles, ici je vais calculer celles pour la richesse et la population. Prenons l'exemple de la population :

On souhaite qu'une mise à jour soit faite quand la population atteint une valeur v.
Si pop = s1 alors l'évolution de la population est constante et il n'y a rien à faire.
Si pop < v <= s1 alors a = exp(taux*pop/s1) et v = s1*log(a+t)/taux, donc t = exp(taux*v/s1) - exp(taux*pop/s1).
Si s1 <= v < pop alors a = exp(taux*(1 - pop/s0)) et v = s0*(1 - log(a+t)/taux), donc t = exp(taux*(1 - v/s0)) - exp(taux*(1 - pop/s0)).
Dans les cas où v ne respecte pas sa contrainte avec pop, t est alors négatif.
On peut aussi poser que quand t = 0 aucune mise à jour n'est faite.

Cette fonction ne sert pas que pour mettre à jour quand la planète atteint un seuil de population ou de richesse, mais peut servir pour n'importe quelle valeur à atteindre, par exemple pour déterminer la date où un travailleur disparaît :

Code:

// t > 0 si v est atteignable, sinon t = 0
nextUpdPop(planete,v) =
   t = 0
   pop = planete.population
   s0 = planete.dernierSeuilPop
   s1 = planete.seuilPop
   if pop < v <= s1
      t = exp(taux*v/s1) - exp(taux*pop/s1)
   if s1 <= v < pop
      t = exp(taux*(1 - v/s0)) - exp(taux*(1 - pop/s0))
   return t

// t > 0 si v est atteignable, sinon t = 0
nextUpdRich(planete,v) =
   t = 0
   rich = planete.richesse
   r0 = planete.dernierSeuilRich
   r1 = planete.seuilRich
   if rich < v <= r1
      t = exp(taux*v/r1) - exp(taux*rich/r1)
   if r1 <= v < rich
      t = exp(taux*(1 - v/r0)) - exp(taux*(1 - rich/r0))
   return t

//il nous faut ignorer les 0 lorsque le min est calculé entre deux valeurs positives
minz(a,b) =
   if a == 0
      return b
   elsif b == 0
      return a
   else
      return min(a,b)

nextUpdate(planete) =
   t = 0
   // si la population atteint son seuil
   s1 = planete.seuilPop
   tPop = nextUpdPop(planete,s1)
   t = minz(t,tPop)
   // si la richesse atteint son seuil
   r1 = planete.seuilRich
   tRich = nextUpdRich(planete,r1)
   t = minz(t,tRich)
   // si un travailleur disparaît
   popActive = ...
   tAct = nextUpdPop(planete,popActive)
   t = minz(t,tAct)
   // on peut continuer à ajouter des t sur le même modèle
   t0 = planete.lastUpdate
   if t > 0
      return t0 + t
   else
      //pas de nouvelle mise à jour

Voilà qui est fait, et j'ai aussi édité l'algo d'imposition au post précédent pour le rendre plus lisible.
Du coup, Nimé, tu pourras implémenter tout ça pendant mon départ.
Si j'ai le temps je ferai des tests pour traquer les erreurs d'inattention ou évaluer l'impact du nombre de mises à jour sur la valeur de l'imposition, mais ce n'est pas gagné.

PS : Pas la peine de te mettre martel en tête pour ma remarque sur les seuils entiers, ce n'est pas nécessaire.