AM Un petit defi: ecris un programme qui donne la valeur de zeta(3) a un nombre de decimales arbitraires (n'utilisant que Tcl, pas [exec])
GS Pour l'instant, j'arrive à 10 décimales :-(
set tcl_precision 17
proc zeta3 {} {
set z 0
set c1 1
set c2 5
set n 0
while {$n != 30} {
incr n
set c1 [expr {2*$c1*(2*$n-1)/$n}]
set n2 [expr {$n*$n}]
set w [expr {$c2/(2.*$c1*$n2*$n)}]
set z [expr {$z + $w}]
incr n
set c1 [expr {2*$c1*(2*$n-1)/$n}]
set n2 [expr {$n*$n}]
set w [expr {$c2/(2.*$c1*$n2*$n)}]
set z [expr {$z - $w}]
}
puts "$z"
}
zeta3KBK La valeur de zeta(3) se trouve dans une série curieuse:
On peut évaluer le coefficient binomial (2k,k) par la récurrence:
Il faudrait des procédures pour calculer dans précision arbitraire. Sans ces calculs, l'on peut obtenir 16 décimales par additioner 23 termes de la série:
set tcl_precision 17
proc zeta3 {} {
set z [expr wide(0)]
set c [expr wide(1)]; # Le nombre Catalan C_k; C_1 == 1
set n 5
for { set k [expr wide(1)] } { $k < 23 } { incr k } {
set ch [expr { ( $k + wide(1) ) * $c }]
# Le coefficient binomial (2k, k)
set d [expr { wide(2) * $ch * $k * $k * $k }]
set z [expr { $z + double($n) / double($d) }]
set c [expr { $c * wide(2) * ( wide(2) * $k + wide(1) )
/ ( $k + wide(2) ) }]
set n [expr { - $n }]
}
return $z
}
puts [zeta3]AM En effet, il y a plusieurs paquets d'arithmetique de precision arbitraire, par exemple MPA par Stephane Arnold.
KBK Et maintenant, je suis arrivé à 2.000 décimales - par ce script-ci:
# Addition des grands entiers:
proc add { a b } {
if { [llength $a] < [llength $b] } {
foreach { b a } [list $a $b] break
}
while { [llength $b] < [llength $a] } {
lappend b 0
}
set s {}
set c 0
foreach p $a q $b {
set r [expr { $p + $q + $c}]
if { $r > 99999999 } {
set c 1
incr r -100000000
} else {
set c 0
}
lappend s $r
}
if { $c } {
lappend s 1
}
return $s
}
# Soustraction des grands entiers:
proc sub { a b } {
if { [llength $a] < [llength $b] } {
return -code error "a < b"
}
while { [llength $b] < [llength $a] } {
lappend b 0
}
set d {}
set c 0
set zc 0
foreach p $a q $b {
set r [expr { $p - $q - $c }]
if { $r < 0 } {
set c 1
incr r 100000000
} else {
set c 0
}
if { $r != 0 } {
while { $zc > 0 } {
lappend d 0
incr zc -1
}
lappend d $r
} else {
incr zc
}
}
return $d
}
# Multiplication d'un grand entier x par un entier de précision single a
proc shmul { a x } {
set r {}
set c [expr {wide(0)}]
foreach p $x {
set prod [expr { wide($a) * wide($p) + $c }]
lappend r [expr { int( $prod % 100000000 ) }]
set c [expr { int( $prod / 100000000 ) }]
}
if { $c != 0 } {
lappend r $c
}
return $r
}
# Multiplication d'un grand entier a par un entier de précision single x.
# La valeur de la proc est une liste { a/x a%x }
proc shdivmod { a x } {
set f 0
set q $a
set r [expr {wide(0)}]
for { set i [expr { [llength $a] - 1 }] } { $i >= 0 } { incr i -1 } {
set p [lindex $a $i]
set y [expr { wide(100000000) * $r + $p }]
set qd [expr { int ( $y / $x ) }]
if { $qd != 0 && !$f } {
set l $i
set f 1
}
lset q $i [expr { int( $y / $x ) }]
set r [expr { int( $y % $x ) }]
}
if { !$f } {
return [list 0 $r]
} else {
return [list [lrange $q 0 $l] $r]
}
}
# Division a/x
proc shdiv { a x } {
return [lindex [shdivmod $a $x] 0]
}
# Multiplication d'un grand entier par une suite des entiers de
# précision single.
proc shmuln { a args } {
set p $a
foreach b $args {
set p [shmul $b $p]
}
return $p
}
# Division d'un grand entier par une suite des entiers de
# précision single.
proc shdivn { a args } {
set q $a
foreach d $args {
set q [shdiv $q $d]
}
return $q
}
# Comparaison des grands entiers.
proc >= { a b } {
if { [llength $a] > [llength $b] } {
return 1
} elseif { [llength $a] < [llength $b] } {
return 0
} else {
for { set i [expr { [llength $a] - 1 }] } { $i >= 0 } { incr i -1 } {
if { [lindex $a $i] > [lindex $b $i] } {
return 1
} elseif { [lindex $a $i] < [lindex $b $i] } {
return 0
}
}
return 1
}
}
# Le fonction zêta(3) de Riemann: on démande $nd décimales:
proc zeta3 { nd } {
set ch 20; # C(2k,k), k == 3
set ns [list 19]; # numérateurs des termes
set ds [list 96]; # Dénominateurs des termes
set nt [expr { int ( $nd * log( 10 ) / log( 4 ) + 1 ) }]
# Nombre des termes
set guard [expr { int ( log( $nt ) / log( 10 ) + 2 ) }]
# Compter $nt termes de la série
for { set k 3 } { $k < $nt } { incr k 2 } {
set twok [expr { $k + $k }]
set threek [expr { $twok + $k }]
set kp1 [expr { $k + 1 }]
set denom [shmuln $ch [expr { $twok + 1 }] $kp1 $kp1 \
$k $twok $twok]
set numer [shmul 5 \
[add 2 \
[shmuln [expr {$threek + 4}] [expr { $k+2 }] $k]]]
lappend ns $numer
lappend ds $denom
set ch [shdivn \
[shmuln $ch 4 [expr { 1 + $twok }] [expr { 3 + $twok}]] \
[expr { 1 + $k }] [expr { 2 + $k }]]
}
# Additioner la série et construire la réprésentation décimale:
set w 0
set line 1.
for { set i 0 } { $i < $nd + $guard + 16 } { incr i } {
# Construire la décimale $i:
set add 0
set newns $ns
set j 0
foreach n $ns d $ds {
set n [shmul 10 $n]
while { [>= $n $d] } {
set n [sub $n $d]
incr add
}
lset newns $j $n
incr j
}
# Additioner vers les décimales plus significantes:
set w [expr { wide(10) * $w + $add }]
if { $i > 16 && $i < $nd + 16} {
# ZUT! Si 9999999999999999 s'apparaît dans la réprésentation,
# ce code-ci ne marchera plus. (Il n'y a aucun 9999999999999999
# dans les 1.000.000 premières décimales.)
append line [expr { $w / 10000000000000000 }]
set w [expr { $w % 10000000000000000 }]
if { [string length $line] == 78 } {
puts $line
set line {}
}
}
set ns $newns
}
# Imprimer la dernière ligne incomplète:
if { [string length $line] != 0 } {
puts $line
}
}
zeta3 154AM Merveilleux! Non seulement un script qui peut faire la calculation en precision arbitraire mais aussi un petit bibliothèque pour des calculations a nombres entiers arbitraires. Ca vaut soumission au bibliothèque Tcllib, a mon avis.
David Cobac Pour montrer les effets des calculs itératifs décimaux, voici un petit exemple des termes d'une suite "classique" (exemple adapté d'un article de l'APMEP de JF Colonna) calculés de trois manières différentes, on comprend l'intérêt des bibliothèques de précision arbitraire :
set tcl_precision 17
proc suite1 { R u_n } {
return [expr { ($R+1)*$u_n - $R*pow($u_n,2) }]
}
proc suite2 { R u_n} {
return [expr { ($R+1)*$u_n - $R*$u_n*$u_n }]
}
proc suite3 { R u_n } {
return [expr { $u_n*($R+1-$R*$u_n) } ]
}
set n 0
set R 3
set u(1,$n) 0.5
set u(2,$n) $u(1,$n)
set u(3,$n) $u(1,$n)
foreach i {1 2 3} {puts "u($i,$n)=$u($i,$n)"}
puts ""
while {$n<=100} {
incr n
set prec [expr {$n-1}]
foreach i {1 2 3} {
set u($i,$n) [suite$i $R $u($i,$prec)]
if {$n % 10 == 0} {
puts "u($i,$n)=$u($i,$n)"
if {$i==3} {puts ""}
}
}
}AM Est-que tu as le reference pour cet article? (J'ai trouve le magasin bimensuel d'APMEP, mais pas l'article qui decrit un tel exemple).
David Cobac Je l'ai lu récemment mais je n'ai plus la brochure chez moi...Néanmoins, publimath permet de le retrouver et voilà : http://publimath.irem.univ-mrs.fr/biblio/AAA03063.htm Les tests décrits dans l'article sont écrits en C. Puis l'auteur donne quelques résultats sous différentes plateformes et avec différentes versions de gcc.
GS Un tout petit programme pour calculer 1000 décimales de e en utilisant un algorithme compte-goutte (spigot algorithm).
# 1000 décimales de e avec un algorithme compte-goutte
proc e1 {} {
set f {}
set n 1000
set b [expr {10 * $n / 4}]
for {set i 0} {$i <= $b} {incr i} {lappend f 1}
incr n -1
puts -nonewline "2."
for {set j 1} {$j <= $n} {incr j} {
set q 0
for {set i $b} {$i >= 1} {incr i -1} {
set k [expr {$i + 1}]
set x [expr {[lindex $f $i]*10 + $q}]
lset f $i [expr {$x % $k}]
set q [expr {$x / $k}]
}
puts -nonewline $q
flush stdout
}
}
e1Catégorie Mathématiques | Catégorie Concours