Nullstellen von Polynomen - Bairstow-Verfahren

[p. specht]

Thema Kombinatorik
-----------------------
N<16 Bauern aus der schönen Steiermark*), alle mit typischem Steirerhut auf, gehen Sonntags in die Kirche. Der Pfarrer schaut böse, weil sie beim Betreten des Gotteshauses noch immer die Hüte aufhaben! Schnell werfen sie die Hüte mangels Ablage rechts in die Ecke. Am Ende, nach dem Amen und Hosianna-Frohlocken, nehmen sie sich rasch wieder jeder einen Hut und suchen das Weite (sprich: Sie eilen zum Frühschoppen).

Dummerweise sehen Steierhüte alle gleich aus. Aus hygienischen Gründen wäre es natürlich gut, wenn jeder seinen eigenen Hut erwischt hätte. Wie hoch stehen die Chancen dafür, daß K Hüte jeweils auf dem richtigen Kopf gelandet sind?

WindowTitle "RENCONTRES-Zahl, DERANGEMENT und SUBFAKULTÄT"
' (L)Copyleft 2011ff P.Specht für Paules PC Forum
' Versuch einer Umsetzung des Wikipedia-Artikels betr. Rencontres-Zahl
' in XProfan 11.2a, KEINE GEWÄHR - No warranty whatsoever!
Font 2
declare p!,n!,k!,Ren!,i&
Weiter:
cls
Print " In der KOMBINATORIK versteht man unter der RENCONTRES-Zahl   "
print " (französisch für 'Begegnungen') die mit D(n;k) bezeichnete   "
print " Anzahl jener PERMUTATIONEN einer Menge n unterscheidbarer    "
print " Elemente, bei der genau k Elemente ihren ursprünglichen bzw. "
print " einen bestimmten gewünschten Platz einnehmen (und n-k nicht)."
print " Ren=D(n;k)=n!/k!*SUM[i=0..(n-k)](-1)^i/i!=(n OVR k)*D(n-k;0) "
print "                                                              "
print " Für den Fall, dass KEINES der n Elemente seinen Platz ein-   "
print " nimmt bzw. 'wiederfindet', ergibt sich als Sonderfall die    "
print " Formel für die Zahl möglicher DERANGEMENTS oder 'Totalver-   "
print " setzungen' aller n Elemente zu  !n = 'SUBFAKULTÄT von n'     "
print " nach der Formel:  !n = D(n;0) = n! * SUM[i=0..n](-1)^i/i!    "  
' print " {Interessant: lim[n..+Inf](SUM[i=0..n](-1)^i/i!))= 1/exp(1)} "
Print "                                                              "
Print " Bsp: Anzahl n der zu permutierenden Elemente eingeben: ";:input n!
if n!>15
  print " Wegen oberer Integer-Grenze bitte nur Zahlen bis 15 - Sorry! "
  WaitInput
  goto "weiter"
endif
print " Prinzipiell gäbe es "; int(fakul(n!)); " Positions-Permutationen."
Print "                                                              "
Print " Wieviele Elemente sollen in Wunschposition stehen?: 0";:input k!
'print "                                                              "
print " Dann gibt es genau ";
  Ren!=Rencontres_D(n!,k!)
  set("decimals",0)
  print Ren!;" solche Permutationen."
set("decimals",3)
print " Die Wahrscheinlichkeit für so eine Stellung ist ";100*Ren!/fakul(int(n!));"%"
  set("decimals",0)
WaitInput
goto "Weiter"

Proc Rencontres_D : parameters n!,k!
var n&=int(n!)
var k&=int(k!)
var p!=1
whileLoop k&+1,n&
p!=p!*&Loop
EndWhile
var s!=0
var i&=0
while i&<=(n&-k&)
s! = s! + (1.0-2.0*(i& mod 2)) / fakul(i&)
inc i&
endwhile
'print "Vorfaktor: ";p!
'print "    Summe: ";s!
return p! * s!
EndProc

Proc fakul
parameters p&
var prd!=1
case p&<1 : p&=1
case p&>169 :prd! = -1
case prd!<0: goto "back"
whileloop p&,1,-1
prd!=prd!*&Loop
endwhile
back:
return prd!
EndProc

Alle kamen übrigens wohlbehalten heim, und nur drei machten Bekanntschaft mit dem Nudelholz der Frau Gemahlin... Bis nächsten Sonntag waren die Beulen aber dank der solide gearbeiteten Kopfbedeckung wieder abgeklungen. Dadurch wird auch der Begriff "Tracht Prügel" verständlich.
Gruss
____
*) ....wo übrigens auch Schwarzenegger herkommt!

[p. specht]

Thema Polynome bis ca. 12. Grades lösen - Das Bairstow-Verfahren
===========================================
Von einem wahren Könner abgekupfert und von mir in Profan11.2a übertragen.
Braucht man ja äh... fast täglich , sozusagen.. Gruss

'Var Koeff$="1,10"
'Var Koeff$="10,0,-1000"
'Var Koeff$="1,-10,0,100,0,10,0,1,0,-1000"
Var Koeff$="10,0,-20.2,1.8,-40,300,-1,0"

WindowTitle "Nullstellen von Polynomen - Bairstow-Verfahren"
' Quelle: http://www.rhirte.de/vb/home2.htm
' Umsetzung des VB-Programms von Prof.em. Dr.Rolf Hirte, Technische Fachhochschule Wildau
' nach XProfan 11.2a; Nur zur Demonstration gedacht; Keine Garantie! Demoware only!
' Begleittext: Die Lösung von Polynomausdrücken (Nullstellensuche) wird unbequem, wenn das
' Polynom etwa vom Grade 5 oder z,B. 25 ist. Bairstow verfährt dabei folgendermaßen:
' Er spaltet in einer Iteration laufend die quadratischen Faktoren ab, die dann in bekannter Weise 
' (Satz von Vieta) gelöst werden - das so lange, bis das Restpolynom vom Grade 0 oder 1 ist.
' Der Anwender muß LEDIGLICH die einzelnen Koeffizienten des auf Normalform "Polynom = 0" gebrachten Polynoms  
' in fallender Potenz bereitstellen (Text in starker Anlehnung an die genannte Quelle).
'
' WIR HIER packen die Koeffizienten sowie das Absolutglied nach fallendem Exponenten in einen 
' Komma-separierten String, z.B. Var Koeff$="10,20,-2.2,1.8,-40,300,1.1111111,0", siehe 1. Zeile.
' Danach gehts los:

Cls
Font 2
AppendMenuBar 253,"Run Bairstow"
AppendMenuBar 252,"Finish"
print
Declare Grad&,A$[],A![],W!
A$[]=Explode(Koeff$,","):Grad&=SizeOf(A$[])-1
SetSize A![],Grad&
WhileLoop 0,Grad&
  A![&Loop]=val(A$[Grad&-&Loop])
EndWhile
WhileLoop Grad&,0,-1
  Set("NumWidth",1) :set("Decimals",0)
  if &Loop<Grad&:print " +";:else:print "  ";:endif
  Print "X^";&Loop,
  Set("NumWidth",26):set("Decimals",15):Print " mal ";A![&Loop]
EndWhile
Print:Print " Für Berechnungsstart bitte Taste drücken!"
WaitInput
CLS
Set("NumWidth",1) :set("Decimals",0)
Print "Das gegebene Polynom vom Grad ";Grad&;" hat folgende Nullstellen:":Print
    Bairstow(A![])
WaitInput
End 'Main


Proc Bairstow
Parameters A![]
Var Grad& = SizeOf(A![])-1
Declare i&,R!,P!,P1!,Q!,Q1!,Ce!,s!,t!
Declare B![Grad&],C![Grad&]
set("NumWidth",20):set("Decimals",15)
While Grad& > 2
  R! = 0
  P1! = 1
  Q1! = -1
  B![Grad&] = A![Grad&]
  C![Grad&] = A![Grad&]
  Repeat
    P! = P1!
    Q! = Q1!
    B![Grad&-1] = B![Grad&] * P! + A![Grad& - 1]
    C![Grad&-1] = B![Grad&-1] + B![Grad&] * P!

    Whileloop Grad& - 2 , 0 , -1 : i& = &Loop  'For
      B![i&] = B![i& + 2] * Q! + B![i& + 1] * P! + A![i&]
      C![i&] = C![i& + 2] * Q! + C![i& + 1] * P! + B![i&]
    EndWhile ' Next
    Ce! = C![2] * C![2] - C![1] * C![3]
    case Ce! = 0 : Print " Andere Startwerte nötig!"
    P1! = P! - (B![1] * C![2] - B![0] * C![3]) / Ce!
    Q1! = Q! - (B![0] * C![2] - B![1] * C![1]) / Ce!
    R! = R! + 1
    If R! > 4000
      Print:Print:Print "   Sorry, nach 4000 Runden keine (weitere) Konvergenz!"
      WaitInput
      End
    EndIf
  Until Abs(B![0]) + Abs(B![1]) < 10^-12 ' Innere loop

' Nullstelle des quad. Faktors
  s! = P1! / 2
  t! = P1!*P1! + 4 * Q1!
  If t! < 0
      Print s!;"   + ";0.5 * Sqrt(-t!);"*i "; : comment
      Print s!;"   - ";0.5 * Sqrt(-t!);"*i "; : comment
  Else
      Print s! + 0.5 * Sqrt(t!)
      Print s! - 0.5 * Sqrt(t!)
  EndIf
  whileloop 2,Grad&
    i& = &Loop
    A![i& - 2] = B![i&]
  endwhile
  Grad& = Grad& - 2

EndWhile ' Outer Loop


If Grad& = 1
  Print -A![0]/A![1]
Else
  s! = -0.5 * A![1] / A![2]
  t! = A![1] * A![1] - 4 * A![2] * A![0]
    If t! < 0
        Print s!;"   + ";0.5 * Sqrt(-t!) / A![2];"*i "; : comment
        Print s!;"   - ";0.5 * Sqrt(-t!) / A![2];"*i "; : comment
    Else
        Print s! + 0.5 * Sqrt(t!) / A![2]
        Print s! - 0.5 * Sqrt(t!) / A![2]
    EndIf
EndIf
EndProc


Proc comment
if nearly(s!,0,9)
  print "(Imaginär)"
else
  print "(Komplex)"
endif
endproc