Quelques opérations sur les polynomes

Commenter

Contenu du snippet

Tout est dans le titre :p
Ce n'est pas un addon a proprement parler , notamment a cause de son inutilité totale pour le moment/pour certains , et le fait que les syntaxes soient un peu particulières. Les algorythmes utilisés sont simplistes : ils font les opérations que l'on ferait a la main ...
Enfin ce snippet s'inscrit dans une optique plus large d'addon de calcul formel qui n'est pas géré ici mais pas inhibé non plus par les fonctions : les simplifications d'expressions se font "en fraction" si les expressions sont totalement numériques ou ne se font pas pour l'instant si l'expression contient des lettres ...
Un parser permettant la lecture/simplification de telles expressions est en cours mais c la partie 'difficile' du travail :)

Source / Exemple : 


;### Gestion de polynomes
;division euclidienne (naturels)
Alias iquo {
  if (!$regex($1,^[0-9]+$)) || (!$regex($2,^[0-9]+$)) || (!$2) || (!$isid) { halt }
  if ($1 < $2) { return 0 $1 }
  var %len = $calc($len($1) + 1) , %res = $1 , %quo
  while (%len >= 0) {
    var %i = 9
    while (($calc(%res - $2 * %i * 10 ^ %len ) < 0) && (%i >= 1)) { dec %i }
    %quo = %quo $+ %i
    %res = $calc(%res - $2 * %i * 10^ %len )
    dec %len
  }
  return $calc(%quo) %res
}

;algorithme d'euclide pour le pgcd de deux entiers non nuls
Alias pgcd {
  if (!$regex($1,^-?[0-9]+$)) || (!$regex($2,^-?[0-9]+$)) || (!$1) || (!$2) || (!$isid) { halt }
  if ($abs($1) >= $abs($2)) { var %a = $abs($1) , %b = $abs($2) } | else { var %b = $abs($1) , %a = $abs($2) }
  var %r = $gettok($iquo(%a,%b),2,32)
  while (%r) {
    %reste = $gettok($iquo(%a,%b),2,32)
    %a = %b
    %b = %reste
    %r = $gettok($iquo(%a,%b),2,32)
  }
  return %b
}

;réduction d'un rationnel
Alias reduc {
  if (!$regex($1,^-?[0-9]+$)) || (!$regex($2,^-?[0-9]+$)) || (!$2) || (!$isid) { halt }
  if (!$abs($1)) { return 0 1 }
  return $calc($1 / $pgcd($1,$2)) $calc($2 / $pgcd($1,$2))
}

;calcul d'une expression dans Q
Alias Qcalc {
  if (!$isid) || ($count($1,$chr(40)) != $count($1,$chr(41))) { echo -ta error dans Qcalc | halt } 
  var %q = $remove($1,$chr(32))
  if (!$regex(%q,^[-0-9\+\*\/\(\)^]+$)) { return %q }
  if ($regex(%q,^[-0-9\+\*\(\)^]+$)) { return $calc(%q) }
  var %res = 1 , %i = 1 , %x = $regex($1,/([0-9]+)/g)
  while ($regml(%i) != $null) { if ($ifmatch != 0) { %res = %res * $ifmatch } | inc %i }
  var %res2 = $calc(%res * ( %q ))
  var %num = $gettok($reduc(%res2,%res),1,32) , %den = $gettok($reduc(%res2,%res),2,32)
  return %num $+ $iif(%den != 1 && %den,/ $+ %den)
}

;Polynômes
;la syntaxe employée est la liste des coefs par ordre DECROISSANT > 3*X^2-5*X+4 s'écrit 3;-5;4
;volontairement , les trois lois de composition rendent un résultat évalué de facon très simple (pas de factorisation et simplifications compliquées) , afin de permettre le calcul formel et éviter les erreur d'arrondis sur les fractions ($calc() rend 0 lorsque l'argument n'est pas une opération sur des nombres...)
;cf le parser de sax sur scriptsdb pour la construction d'un binaire respectant les priorités mais je n'en suis pas la :
;N-ième coef
alias coef {
  if (!$regex($2,^[0-9]+$)) || (!$isid) || (;; isin $1) { halt }
  if ($2 > $count($1,;)) { return 0 }
  if ($count($1,;) == 0) { return $1 }
  if ($prop == detail) { echo -ta Ici on cherche simplement l'emplacement entre points virgules n° $calc($2 +1) en partant de la fin. }
  if ($2 == $count($1,;)) { return $left($1,$calc($pos($1,;,1) -1)) }
  var %mask = ;([^;]+)(;[^;]+) $+ $chr(123) $+ $2 $+ $chr(125) $+ $
  if ($regex($1,%mask)) { return $regml(1) }
  return 0
}

;degré d'un polynome (sur un polynome numérique sinon l'évaluation ne se fait pas)
alias degree {
  if (!$isid) || (;; isin $1) { halt }
  var %x = 1 , %p = $simplifypoly($1) , %deg = $count(%p,;)
  if ($prop == detail) { echo -ta Simplification du polynôme de départ : $1 -> %p qui comporte $calc(%deg +1) coefficients et que l'on réduit par la fonction $ $+ reducpoly en $reducpoly(%p) de premier coefficient non nul d'ou le degré. }  
  %deg = $count($reducpoly(%p),;)
  if (%deg) { return %deg }
  return $iif($coef(%p,0),0,-Infinity)
}

;Multiplication d'un polynome par le monome X^N : ex: produit de 3*X^2+1 par X^4 : $uppoly(3;0;1,4) retourne 3*X^6+X^4 donc : 3;0;1;0;0;0;0
Alias uppoly { 
  if ($prop == detail) { echo -ta Simple ajout à la fin (on dit aussi concaténation) de $str(;0,$2) ce qui correspond a la multiplication par X^ $+ $2 }
  return $1 $+ $str(;0,$2) 
}

;lci addition : somme de deux listes $sumpoly(poly1,poly2,[+|-]) (+ par défaut ...)
Alias sumpoly {
  if ($count($1,;) >= $count($2,;)) { var %p1 = $1 , %p2 = $str(0;,$calc($count($1,;) - $count($2,;))) $+ $2 , %p3 }
  if ($count($1,;) < $count($2,;)) { var %p2 = $2 , %p1 = $str(0;,$calc($count($2,;) - $count($1,;))) $+ $1 , %p3 }
  var %i = $count(%p1,;) , %op = $iif($3 == -,-,+) , %res
  if ($prop == detail) { echo -ta Les deux polynomes sont ramenés a un nombre de termes identiques si ce n'est pas le cas: %p1 et %p2 . Ensuite les termes sont ajoutés en degré décroissants }
  while (%i) {
    if ($prop == detail) { echo -ta Coefficient en X^ $+ %i : $+($chr(40),$coef(%p1,%i),$chr(41),%op,$chr(40),$coef(%p2,%i),$chr(41)) }
    %res = $+(%res,$chr(40),$coef(%p1,%i),$chr(41),%op,$chr(40),$coef(%p2,%i),$chr(41),;) 
    dec %i
  }
  if ($prop == detail) { echo -ta Coefficient en X^0 : $+($chr(40),$coef(%p1,0),$chr(41),%op,$chr(40),$coef(%p2,0),$chr(41)) }
  %res = $+(%res,$chr(40),$coef(%p1,0),$chr(41),%op,$chr(40),$coef(%p2,0),$chr(41))
  if ($prop == detail) { echo -ta D'où finalement le résultat non simplifié %res puis simplifié : } 
  return $simplifypoly(%res)
}

;Multiplication externe : $scalpoly(3;0;2,5) <-> 5*(3*X^2+2) = 15*X^2+10 <-> 5*3;5*0;5*2 qui évaluée donne : 15;0;10
;cette fonction est inutile dans l'absolu car par abus de notation le 1 scalaire se confond avec le 1 de Q[X] mais elle est utilisée pour plus de clarté dans les aliases ... (savoir si on manipule un scalaire ou un vecteur de Q[X])
Alias scalpoly {
  var %i = $count($1,;) , %res
  while (%i) {
    if ($prop == detail) { echo -ta Coefficient en X^ $+ %i : $+($chr(40),$2,$chr(41),*,$chr(40),$coef($1,%i),$chr(41)) }
    %res = $+(%res,$chr(40),$2,$chr(41),*,$chr(40),$coef($1,%i),$chr(41),;) 
    dec %i
  }
  if ($prop == detail) { echo -ta Coefficient en X^0 : $+($chr(40),$2,$chr(41),*,$chr(40),$coef($1,0),$chr(41)) }
  %res = $+(%res,$chr(40),$2,$chr(41),*,$chr(40),$coef($1,0),$chr(41))
  if ($prop == detail) { echo -ta D'où finalement le résultat non simplifié %res puis simplifié : }
  return $simplifypoly(%res)
}
;Multiplication interne : ((Q[X],.,+,*) est une Q-algèbre associative unitaire... la belle affaire :p... lol) : par abus de notation : réalise l'opération : (3*X^2+1)*(X+1) = 3X^3+3X^2+X+1 : $prodpoly(3;0;1,1;1) > 3;3;1;1)
Alias prodpoly {
  var %i = $count($1,;) , %res = 0 , %res2
  if ($prop == detail) { echo -ta On utilise ici la distributivité (à droite mais c'est ici arbitraire) de la multiplication de polynômes. On multiplie donc successivement le polynome $1 par les monomes composant $2 , a savoir $coef($2,0) * X^0 , $coef($2,1) * X^1 , etc grace aux fonctions $ $+ scalpoly  et $ $+ uppoly }
  while (%i >= 0) {
    if ($prop == detail) { echo -ta Multiplication de $1 par par $coef($1,%i) * X^ $+ %i soit $uppoly($scalpoly($2,$coef($1,%i)),%i) }
    %res2 = $uppoly($scalpoly($2,$coef($1,%i)),%i)
    %res = $sumpoly(%res,%Res2)
    dec %i
  }
  if ($prop == detail) { echo -ta D'où finalement le résultat non simplifié %res puis simplifié : }
  return $simplifypoly(%Res)
}
;réduction d'une liste jusqu'a un 'coefficient dominant non nul... (ne fonctionne que sur des polynomes numériques simplifiés. > 5-5 est un coeficient non nul pour cet alias ...)
Alias reducpoly {
  if (!$regex($1,^([^;]+;)*[^;]+$)) || (!$isid) { halt }
  if ($prop == detail) { echo -ta On se content ici d'éliminer le plus possible de 0; en début de chaine ce qui revient a chercher le premier coefficient non nul . (marche sur un polynome déja simplifié) }
  var %x = $regsub($1,^(0+;)+(.*)$,\2,%v) , %x2 = %v
  unset %v
  return %x2
}

;simplification d'écritures
alias simplifypoly {
  var %i = $count($1,;) , %res = $qcalc($coef($1,%i))
  dec %i
  if ($prop == detail) { echo -ta Pour chaque coefficient , on utilise la fonction $ $+ qcalc qui renvoie : le même résultat que $ $+ calc si on a que des nombres sans fractions. Un résultat fractionnaire si il y a des fractions (ce résultat peut etre erroné si des puissances interviennent) . Un résultat non évalué pour des expressions littérales. }
  while (%i >= 0) {
    %res = %res $+ ; $+ $qcalc($coef($1,%i))
    dec %i
  }
  return $reducpoly(%res)
}
;division de polynomes $iquopoly(3;3;1;2,1;1) retourne la partie entière : 3;0;1 et le reste 1 . (de facon a ce que 3;3;1;2 = 3;0;1 * 3;0;1 + 1 schématiquement)
alias iquopoly {
  if (!$regex($1,^([^;]+;)*[^;]+$)) || (!$isid) { halt }
  var %p1 = $reducpoly($simplifypoly($1)) , %p2 = $reducpoly($simplifypoly($2)) , %deg1 = $degree(%p1) , %deg2 = $degree(%p2) , %i = %deg1 , %quo , %res = %p1 , %poly
  if (%deg1 < %deg2) { if ($prop == detail) { echo -ta $1 est de degré inférieur a $2 donc la division des deux polynomes ne fournit pas de partie entière : } | return 0 $1 }
  if (%deg2 > 0) {
    var %p22 = $scalpoly(%p2,$+(1/,$chr(40),$coef(%p2,%deg2),$chr(41)))
    %p22 = 1; $+ $gettok(%p22,2-,59)
    if ($prop == detail) { echo -ta Le diviseur est rendu unitaire pour éviter de gérer des fractions dans la boucle : $2 > %p22 .L'algorithme qui suit est pratiquement le même qu' "à la main". }
    while (%i >= %deg2) {
      %poly = $uppoly($scalpoly(%p22,$coef(%res,%i)),$calc(%i - %deg2))
      ;mise a jour du quotient
      %quo = %quo $+ $coef(%res,%i) $+ $iif(%i > %deg2,;)
      ;différence 
      if ($prop == detail) { echo -ta Etape X^ $+ %i : Dans %res on a $coef(%res,%i)) fois X^ %i (on regarde les coefs dominants) on mets donc a jour le quotient en rajoutant $coef(%res,%i) pour X^ $+ %i et on retranche %poly au reste provisoire %res : $simplifypoly($sumpoly(%res,%poly,-)) }
      %res = $simplifypoly($sumpoly(%res,%poly,-))
      dec %i
    }
    if ($prop == detail) { echo -ta Le reste est maintenant de degré < au degré de %p2 . On a effectué la division par un polynôme unitaire de quotient %quo et de reste %res , on  se ramène maintenant a la division de départ : }
    %quo = $scalpoly(%quo,$+(1/,$chr(40),$coef(%p2,%deg2),$chr(41)))
    return $simplifypoly(%quo) $simplifypoly(%res)
  }
  if (%deg2 == 0) {
    if ($prop == detail) { echo -ta Simple division par un réel $2 , qui revient a un $ $+ scalpoly de $1 par $+(1/,$chr(40),%p2,$chr(41)) }
    return $simplifypoly($scalpoly(%p1,$+(1/,$chr(40),%p2,$chr(41)))) 0 
  }
}

;évaluation
alias evalpoly {
  if (!$regex($1,^([^;]+;)*[^;]+$)) || (!$regex($2,^-?[^;]+$)) || (!$isid) { halt }
  var %i = $count($1,;) , %res
  while (%i >= 0) {
    if ($prop == detail) { echo -ta > Ajout du coefficient en x^ $+ %i > $+(%res,$iif(%res != $null,+),$chr(40),$coef($1,%i),*,$2,^,%i,$chr(41)) > $+(%res,$iif(%res != $null,+),$qcalc($+($chr(40),$coef($1,%i),*,$qcalc($2 ^ %i),$chr(41)))) = $qcalc($+(%res,$iif(%res != $null,+),$chr(40),$coef($1,%i),*,$qcalc($2 ^ %i),$chr(41))) }
    %res = $qcalc($+(%res,$iif(%res != $null,+),$chr(40),$coef($1,%i),*,$qcalc($2 ^ %i),$chr(41)))
    dec %i
  }
  return $qcalc(%res)
}

;primitive
alias intpoly {
  if (!$regex($1,^([^;]+;)*[^;]+$)) || (($2 != $null) && (!$regex($2,^-?[^;]+$))) || (!$isid) { halt }
  var %p = $simplifypoly($1) , %i = $degree(%p) , %res , %x0 = $iif($2 != $null,$2,0)
  if ($prop == detail) { echo -ta > On intègre terme à terme les monômes composant le polynôme. }
  while (%i >= 0) {
    if ($prop == detail) { echo -ta > Terme en X^ $+ $calc(%i +1) : $+($chr(40),$coef(%p,%i),$chr(41),*1/,$chr(40),$calc(%i +1),$chr(41)) }
    %res = $+(%res,$chr(40),$coef(%p,%i),$chr(41),*1/,$chr(40),$calc(%i +1),$chr(41),;)
    dec %i
  }
  %Res = %res $+ 0
  if ($prop == detail) { echo -ta > on retranche maintenant à %res (s'annulant en 0) son évaluation en %x0 de facon a avoir une expression s'annulant pour cette valeur d'où après simplification : }
  %res = $sumpoly(%res,$evalpoly(%res,%x0),-)
  return $simplifypoly(%Res)
}

;dérivation
alias diffpoly {
  if (!$regex($1,^([^;]+;)*[^;]+$)) || (!$isid) { halt }
  var %p = $simplifypoly($1) , %i = $degree(%p) , %res
  if ($prop == detail) { echo -ta > On dérive terme à terme les monômes composant le polynôme et de degré supérieur à 1. }
  if (%i == -Infinity) || (%i == 0) { return 0 }
  while (%i) {
    if ($prop == detail) { echo -ta > Terme en X^ $+ $calc(%i - 1) : $+($chr(40),$coef(%p,%i),$chr(41),*,%i) }
    %res = $+(%res,$chr(40),$coef(%p,%i),$chr(41),*,%i,$iif(%i != 1,;))
    dec %i
  }
  if ($prop == detail) { echo -ta > D'où le résultat non simplifié %Res puis : }
  return $simplifypoly(%Res)
}

;aide
alias Expoly {
  echo -ta 
  echo -ta Opérations sur les nombres
  echo -ta Alias: 4 $ $+ iquo(A,B) : division euclidienne de A par B. Retourne le quotient et le reste :4 $ $+ iquo(7,2) , division de 7 par 2 a pour quotient 3 et pour reste 1 de telle sorte que 7=2*3+1 d'où :4 3 1 (A et B doivent être positifs)
  echo -ta Alias: 4 $ $+ pgcd(N,M) : retourne le plus grand diviseur commun de M et N (M et N peuvent être relatifs)
  echo -ta Alias: 4 $ $+ reduc(N,M) : réduction de la fraction rationnelle N/M. 4 $ $+ reduc(4,8) met la fraction 4/8 sous forme irréductible soit 1/2 et retourne donc >4 1 2 (M et N peuvent être relatifs)
  echo -ta Alias: 4 $ $+ qcalc(expression) : calcul d'une expression numérique en restant dans Q (rationnels)
  echo -ta 
  echo -ta Polynômes : un polynôme est dans ce remote une liste de coefficients donnée par ordre décroissants séparés par des points virgule: le polynôme 4X^2+8X-1 s'écrira 4;8;-1
  echo -ta Pour certaines les commandes ci dessous , vous pouvez ajouter l'attribut personnalisé "detail" qui fournit le détail de ce que fait la fonction. 4$Commande(arguments).detail
  echo -ta Alias: 4 $ $+ coef(poly,N) : retourne le coefficient en X^N du polynôme :4 $ $+ coef(4;3;5;6;9;8;7,3) récupère le coef en X^3 , soit4 6
  echo -ta Alias: 4 $ $+ degree(poly) : degré d'un polynome à coefficients numériques.
  echo -ta Alias: 4 $ $+ uppoly(poly,N) : multiplication d'un polynome par X^N
  echo -ta Alias: 4 $ $+ sumpoly(poly1,poly2,[+/-]) : somme de deux polynomes . le "+" n'est pas obligatoire , c'est le signe utilisé par défaut , le "-" exécute la différence.
  echo -ta Alias: 4 $ $+ scalpoly(poly,N) : produit polynôme/scalaire :4 $ $+ scalpoly(3;0;2,5) : produit 5*(3*X^2+2) = 15*X^2+10 retourne donc4 15;0;10
  echo -ta Alias: 4 $ $+ prodpoly(poly1,poly2) : produit de polynômes :4 $ $+ prodpoly(3;0;1,1;1) : produit de 3X^2+1 par X+1 :retourne donc 3X^3+3X^2+X+1 soit4 3;3;1;1
  echo -ta Alias: 4 $ $+ reducpoly(poly) : réduction d'un polynôme sous forme simplifiée jusqu'à un coefficient de plus haut degré non nul. (/!\ 5-5 n'est pas nul dans cette fonction) 
  echo -ta Alias: 4 $ $+ simplifypoly(poly) : simplification de polynomes (uniquement les expressions numériques pour l'instant)
  echo -ta Alias: 4 $ $+ iquopoly(poly1,poly2) : division euclidienne de polynomes :4 $ $+ iquopoly(3;4;1;2,1;1) : division de 3X^3+4X^2+X+2 par X+1 : retourne le quotient 3X^2+X suivi du reste 2 >4 3;1;0 2
  echo -ta Alias: 4 $ $+ evalpoly(poly,x0) : évaluation du polynôme en x0.
  echo -ta Alias: 4 $ $+ intpoly(poly,[x0]) : primitive du polynôme s'annulant en x0 . Si x0 n'est pas spécifié , primitive s'annulant en 0.
  echo -ta Alias: 4 $ $+ diffpoly(poly) : dérivation.
  echo -ta 
}

Conclusion :


c/c du code dans un de vos remote puis //expoly
Je poste ici surtout pour connaitre les bugs éventuels (un seul bug connu pour le moment , le calcul d'expression avec des puissances au dénominateur ne reste pas dans les rationels parfois ...) et toutes les améliorations/critiquesArgumentées que vous pourrez formuler :)

A voir également

S'inscrire maintenant

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.

Du même auteur (cs_PaDa)