Initial commit
This commit is contained in:
parent
51c1149d60
commit
634e684081
20
README.md
20
README.md
@ -1,2 +1,22 @@
|
||||
# asymptote-config
|
||||
This is my personal config directory for asymptote, the descriptive vector graphics language.
|
||||
|
||||
# Contents
|
||||
My config directory contains packages from the following authors:
|
||||
- Christophe Grospellier for `cg_inequations.asy`.
|
||||
- G. Marris for `gm_probabilitytree.asy`, `gm_tableaux.asy` and `gm.asy`.
|
||||
- Jens Schwaiger for `js_polyhedron` and `js_tube.asy`.
|
||||
- Philippe Ivaldi for all files starting with `pi`.
|
||||
- Mine for `jl_lib.asy`, which only use some of the previous packages.
|
||||
- My global config file for asymptote which produce `.svg` files.
|
||||
|
||||
# Installation
|
||||
1. Clone this repository into your $HOME as the hidden config directory for asymptote :
|
||||
|
||||
`git clone https://github.com/jefflance/asymptote-config.git $HOME/.asy`
|
||||
|
||||
2. Enjoy !
|
||||
|
||||
# Troubleshootings
|
||||
This files aren't really up to date and may contains some bugs which are actually removed.
|
||||
Feel free to tell me.
|
||||
|
||||
163
cg_inequations.asy
Normal file
163
cg_inequations.asy
Normal file
@ -0,0 +1,163 @@
|
||||
// v4 29/10/2011, Christophe Grospellier.
|
||||
|
||||
// INSTALLATION:
|
||||
// copier ce fichier dans le sous-répertoire $HOME/.asy
|
||||
// Move this file in the sub-directory $HOME/.asy
|
||||
|
||||
//CODE:
|
||||
|
||||
import graph;
|
||||
import patterns;
|
||||
|
||||
int lsol=0, rsol=1;
|
||||
real hHatch=.2*cm; // Demi-hauteur des hachures
|
||||
real hBrack=.3*cm, wBrack=2*hBrack/3; // Demi-hauteur et largeur du crochet
|
||||
real extrem, bornebis;
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// a="écriture LaTeX de la solution" , borne=borne de l'inéquation //
|
||||
// position=position du label (N par défaut) //
|
||||
// dirsol=lsol (0) ou rsol : solutions à gauche ou à droite //
|
||||
// brack="]" ou "[" //
|
||||
// solcolor=couleur des solutions //
|
||||
// xMin et xMax=abscisses mini et maxi de l'axe //
|
||||
// ticks=Ticks("%",Step=...,step=...) ou Noticks() //
|
||||
// solutions=true ou false : écrire "solutions" //
|
||||
// hach=true ou false : dessiner les hachures //
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// Place un crochet "où on veut"
|
||||
void bracket(picture pic=currentpicture, Label L="", real a, pair pos=N, string sensc, pen p=currentpen){
|
||||
real sens;
|
||||
label(pic,Label(L,p),(a,0),3.5*pos);
|
||||
if (sensc=="]"){sens=-1;}
|
||||
else {sens=1;}
|
||||
pair hB=(0,hBrack), wB=(wBrack,0);
|
||||
pic.add(new void(frame f, transform t) {
|
||||
pair z=t*(a,0);
|
||||
draw(f,z-hB--z+hB,p);
|
||||
draw(f,z-hB--z-hB+sens*wB,p);
|
||||
draw(f,z+hB--z+hB+sens*wB,p);});
|
||||
}
|
||||
|
||||
|
||||
// Trace un intervalle
|
||||
void interval(picture pic=currentpicture, Label L1="", real a1, pair pos1=N, string sensc1,
|
||||
Label L2="", real a2, pair pos2=N, string sensc2,
|
||||
real xMin=min(a1,a2)-2.5, real xMax=max(a1,a2)+2.5,
|
||||
ticks ticks=Ticks(), pen p=currentpen){
|
||||
xlimits(pic,min=xMin,max=xMax);
|
||||
xaxis(pic,Label("$x$",align=N),ticks,Arrow);
|
||||
bracket(pic,L1,a1,pos1,sensc1,p);
|
||||
bracket(pic,L2,a2,pos2,sensc2,p);
|
||||
draw(pic,(a1,0)--(a2,0),p);
|
||||
}
|
||||
|
||||
// Trace des hachures "où on veut"
|
||||
void hatching(picture pic=currentpicture,
|
||||
real mini, real maxi, real hHatching=2mm,
|
||||
string myHatch){
|
||||
pair hH=(0,hHatch);
|
||||
pic.add(new void(frame f, transform t) {
|
||||
pair m=t*(mini,0), M=t*(maxi,0);
|
||||
fill(f,m-hH--m+hH--M+hH--M-hH--cycle,pattern(myHatch));});
|
||||
}
|
||||
|
||||
|
||||
// Solutions d'inéquation
|
||||
void solonaxis(picture pic=currentpicture,
|
||||
string a="",
|
||||
real borne,
|
||||
pair position=N,
|
||||
int dirsol,
|
||||
string brack,
|
||||
pen solcolor=currentpen,
|
||||
real xMin=borne-3, real xMax=borne+3,
|
||||
ticks ticks=Ticks(),
|
||||
bool solutions=true,
|
||||
bool hach=true)
|
||||
{
|
||||
// Tracé de l'axe // Possibilité de changer le repérage dans Ticks()
|
||||
xlimits(pic,min=xMin,max=xMax);
|
||||
xaxis(pic,Label("$x$",align=N),ticks,Arrow);
|
||||
// Crochet
|
||||
bracket(pic, a, borne, position, brack, 1.5bp+solcolor);
|
||||
// Solutions en couleur
|
||||
if (dirsol==0){extrem=xMin;
|
||||
bornebis=xMax;}
|
||||
else {extrem=xMax;
|
||||
bornebis=xMin;}
|
||||
draw(pic,(borne,0)--(extrem,0),1.5bp+solcolor);
|
||||
// Écrire "solutions"
|
||||
if (solutions){
|
||||
label(pic,Label("\scriptsize solutions",solcolor),((borne+extrem)/2,0),2*N);}
|
||||
// Définition de la hachure
|
||||
add("hach_cg",hatch(H=3mm,dir=NE,solcolor));
|
||||
if (hach){
|
||||
hatching(pic,borne,bornebis,"hach_cg");}
|
||||
}
|
||||
|
||||
//
|
||||
// solonaxis(pic,a,borne,position,dirsol,brack,solcolor=currentpen,
|
||||
// xMin=borne-3,xMax=borne+3,ticks=Ticks(),
|
||||
// solutions=true,hach=true)
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// Solutions d'inéquations
|
||||
void sol2onaxis(picture pic=currentpicture,
|
||||
string a1="",
|
||||
real borne1,
|
||||
pair pos1=N,
|
||||
int dirsol1,
|
||||
string bracket1,
|
||||
pen solcolor1=currentpen,
|
||||
string a2="",
|
||||
real borne2,
|
||||
pair pos2=N,
|
||||
int dirsol2,
|
||||
string bracket2,
|
||||
pen solcolor2=currentpen,
|
||||
real xMin=min(borne1,borne2)-2.5, real xMax=max(borne1,borne2)+2.5,
|
||||
ticks ticks=Ticks(),
|
||||
bool hach=true)
|
||||
{
|
||||
solonaxis(pic,a1,borne1,pos1,dirsol1,bracket1,solcolor1,xMin,xMax,ticks,solutions=false,hach);
|
||||
// Crochet
|
||||
bracket(pic,a2, borne2, pos2, bracket2, 1.5bp+solcolor2);
|
||||
// Solutions en couleur
|
||||
if (dirsol2==0){draw(pic,(borne2,0)--(xMin,0),1.5bp+solcolor2);}
|
||||
else {draw(pic,(borne2,0)--(xMax,0),1.5bp+solcolor2);}
|
||||
// Intervalle solution (solcolor=solcolor1+solcolor2)
|
||||
if (dirsol1==0 && dirsol2==0){
|
||||
draw(pic,(min(borne1,borne2),0)--(xMin,0),(1.5bp+solcolor1+solcolor2));
|
||||
}
|
||||
if (dirsol1!=0 && dirsol2!=0){
|
||||
draw(pic,(max(borne1,borne2),0)--(xMax,0),(1.5bp+solcolor1+solcolor2));
|
||||
}
|
||||
if (borne1<borne2){
|
||||
if (dirsol1!=0 && dirsol2==0){
|
||||
draw(pic,(borne1,0)--(borne2,0),(1.5bp+solcolor1+solcolor2));
|
||||
}
|
||||
}
|
||||
if (borne1>borne2){
|
||||
if (dirsol1==0 && dirsol2!=0){
|
||||
draw(pic,(borne2,0)--(borne1,0),(1.5bp+solcolor1+solcolor2));
|
||||
}
|
||||
}
|
||||
// Définition de la hachure
|
||||
add("hachback_cg",hatch(H=3mm,dir=NW,solcolor2));
|
||||
real bornebis2;
|
||||
if (hach){
|
||||
if (dirsol2==0){bornebis2=xMax;}
|
||||
else {bornebis2=xMin;}
|
||||
hatching(pic,borne2,bornebis2,"hachback_cg");}
|
||||
}
|
||||
|
||||
//
|
||||
// sol2onaxis(pic,a1,borne1,pos1,dirsol1,bracket1,solcolor1=currentpen,
|
||||
// a2,borne2,pos2,dirsol2,bracket2,solcolor2=currentpen,
|
||||
// xMin=min(borne1,borne2)-2.5,xMax=max(borne1,borne2)+2.5,
|
||||
// ticks=Ticks(),hach=true)
|
||||
|
||||
3
config.asy
Normal file
3
config.asy
Normal file
@ -0,0 +1,3 @@
|
||||
import settings;
|
||||
settings.outformat = "svg";
|
||||
//settings.tex = "pdflatex";
|
||||
156
dc_stat.asy
Normal file
156
dc_stat.asy
Normal file
@ -0,0 +1,156 @@
|
||||
//DCstat-1
|
||||
//02/03/2009
|
||||
|
||||
|
||||
pen[] couleur={mediumred,mediumcyan,mediumgreen,mediummagenta,mediumblue,yellow,
|
||||
orange,fuchsia,brown,heavycyan,heavygreen,heavymagenta,deepblue,white,black};
|
||||
|
||||
struct axes{pair min; pair max;real graduationx;real numerotationx;real taillex;real graduationy;real numerotationy;real tailley;}
|
||||
|
||||
axes calculaxes(real[] xi,real[] ni,pair ptO=(0,0)){
|
||||
axes ax;
|
||||
ax.min=ptO;
|
||||
ax.max=(max(xi)*1.05,max(ni)*1.05);
|
||||
|
||||
//axe des x
|
||||
if (ax.max.x-ax.min.x<1){ax.graduationx=0.025;ax.numerotationx=0.1;}
|
||||
else if (ax.max.x-ax.min.x<10){ax.graduationx=0.25;ax.numerotationx=1;}
|
||||
else if (ax.max.x-ax.min.x<100){ax.graduationx=2.5;ax.numerotationx=10;}
|
||||
else if (ax.max.x-ax.min.x<1000){ax.graduationx=25;ax.numerotationx=100;}
|
||||
else if (ax.max.x-ax.min.x<10000){ax.graduationx=250;ax.numerotationx=1000;}
|
||||
else {ax.graduationx=2500;ax.numerotationx=10000;}
|
||||
|
||||
//données axe y
|
||||
ax.tailley=0.0125*(ax.max.x-ax.min.x);
|
||||
if (ax.max.y-ax.min.y<1){ax.graduationy=0.025;ax.numerotationy=0.1;}
|
||||
else if (ax.max.y-ax.min.y<10){ax.graduationy=0.25;ax.numerotationy=1;}
|
||||
else if (ax.max.y-ax.min.y<100){ax.graduationy=2.5;ax.numerotationy=10;}
|
||||
else if (ax.max.y-ax.min.y<1000){ax.graduationy=25;ax.numerotationy=100;}
|
||||
else if (ax.max.y-ax.min.y<10000){ax.graduationy=250;ax.numerotationy=1000;}
|
||||
else {ax.graduationy=2500;ax.numerotationy=10000;}
|
||||
|
||||
//retour axe des x
|
||||
ax.taillex=0.0125*(ax.max.y-ax.min.y);
|
||||
return ax;}
|
||||
|
||||
void traceaxe(axes lesaxes,bool axedesy=true,bool origine=true,pair pointO=(0,0), string libellex,string libelley=""){
|
||||
//axe des abscisses
|
||||
draw((lesaxes.min.x,pointO.y)--(lesaxes.max.x,pointO.y),Arrow(1.5mm));
|
||||
|
||||
//tracé des graduations
|
||||
for (int n = round(lesaxes.min.x/lesaxes.graduationx) ; n <= (0.98*lesaxes.max.x)/lesaxes.graduationx; ++ n) {
|
||||
draw((n*lesaxes.graduationx,pointO.y-0.5*lesaxes.taillex)--(n*lesaxes.graduationx,pointO.y+0.5*lesaxes.taillex));}
|
||||
|
||||
//tracé des numérotations
|
||||
for (int n = round(lesaxes.min.x/lesaxes.numerotationx) ; n <= (0.98*lesaxes.max.x)/lesaxes.numerotationx; ++ n) {
|
||||
draw((n*lesaxes.numerotationx,pointO.y-lesaxes.taillex)--(n*lesaxes.numerotationx,pointO.y+lesaxes.taillex));
|
||||
if(origine){label(scale(0.6)*format(n*lesaxes.numerotationx),(n*lesaxes.numerotationx,pointO.y-lesaxes.taillex),S);}
|
||||
else {if(n!=0){label(scale(0.6)*format(n*lesaxes.numerotationx),(n*lesaxes.numerotationx,pointO.y-lesaxes.taillex),S);}}}
|
||||
|
||||
//Nom de l'axe
|
||||
label(scale(0.9)*libellex,(lesaxes.max.x,pointO.y-5*lesaxes.taillex),SW);
|
||||
|
||||
//axe des ordonnées
|
||||
if (axedesy) {draw((pointO.x,lesaxes.min.y)--(pointO.x,lesaxes.max.y),Arrow(1.5mm));
|
||||
|
||||
//tracé des graduations
|
||||
for (int n = round(lesaxes.min.y/lesaxes.graduationy) ; n <= (0.98*lesaxes.max.y)/lesaxes.graduationy; ++ n) {
|
||||
draw((pointO.x-0.5*lesaxes.tailley,n*lesaxes.graduationy)--(pointO.x+0.5*lesaxes.tailley,n*lesaxes.graduationy));}
|
||||
|
||||
//tracé des numérotations
|
||||
for (int n = round(lesaxes.min.y/lesaxes.numerotationy) ; n <= (0.98*lesaxes.max.y)/lesaxes.numerotationy; ++ n) {
|
||||
draw((pointO.x-lesaxes.tailley,n*lesaxes.numerotationy)--(pointO.x+lesaxes.tailley,n*lesaxes.numerotationy));
|
||||
label(scale(0.6)*format(n*lesaxes.numerotationy),(pointO.x-lesaxes.tailley,n*lesaxes.numerotationy),W);}
|
||||
|
||||
//Nom de l'axe
|
||||
label(scale(0.9)*rotate(90)*libelley,(pointO.x-(log10(lesaxes.max.y)+3)*1.5*lesaxes.tailley,lesaxes.max.y),SW);}
|
||||
}
|
||||
|
||||
void diagramme_baton(real[] xi,real[] ni,pen p1=0.5bp+blue,pair origine=(0,0),bool affichey=true,string nomdesx,string nomdesy="Effectif"){
|
||||
|
||||
real largeur=0.1;
|
||||
real hauteur=0;
|
||||
for(int i=0; i < ni.length;++i){hauteur=hauteur+xi[i];
|
||||
filldraw(((xi[i]-largeur),origine.y)--((xi[i]-largeur),ni[i])--((xi[i]+largeur),ni[i])--((xi[i]+largeur),origine.y)--cycle,p1);}
|
||||
|
||||
traceaxe(calculaxes(xi,ni,origine),origine,nomdesx,nomdesy);
|
||||
}
|
||||
|
||||
|
||||
void diagramme_barre(string[] modalite,real[] ni,pen p1=0.5bp+blue,pair origine=(0,0),bool affichey=true,string nomdesy="Effectif"){
|
||||
|
||||
//tracé des axes
|
||||
real[] tbxi;
|
||||
for (int i=0;i<modalite.length;++i){tbxi[i]=i+1;}
|
||||
axes ax=calculaxes(tbxi,ni,origine);
|
||||
ax.graduationx=(abs(ax.max.x)+abs(ax.min.x))*2;
|
||||
ax.numerotationx=(abs(ax.max.x)+abs(ax.min.x))*2;
|
||||
ax.max=(ax.max.x+0.5,ax.max.y);
|
||||
traceaxe(ax,true,false,origine,"",nomdesy);
|
||||
|
||||
//les rectangles
|
||||
real largeur=0.2;
|
||||
for(int i=0; i < ni.length;++i)
|
||||
{filldraw(((i+1-largeur),origine.y)--((i+1-largeur),ni[i])--
|
||||
((i+1+largeur),ni[i])--((i+1+largeur),origine.y)--cycle,p1);
|
||||
label(scale(0.8)*rotate(90)*modalite[i],(i+1,origine.y),S);}
|
||||
}
|
||||
|
||||
void histogramme(real[] tabxi,real[] tabni,bool axedesy=false,bool unitedaire =true,real uniteaire,pair origine=(0,0),string libellecaractere,string libelleunite="effectif \'egal \`a "+format(uniteaire),pen p1=.8lightgreen,pen p2=1bp+blue){
|
||||
//calculs
|
||||
// Une variable utile pour déterminer la plus petite amplitude :
|
||||
real largeurunite=abs(tabxi[1]-tabxi[0]);
|
||||
// ... et une autre pour le numéro de classe correspondant :
|
||||
int iclasse=0;
|
||||
|
||||
// Calcul des hauteurs (et de la plus petite amplitude de classe) :
|
||||
real[] tabhi;
|
||||
for(int i=0; i < tabni.length; ++i){
|
||||
tabhi[i]=tabni[i]/(tabxi[i+1]-tabxi[i]);
|
||||
if (largeurunite>abs(tabxi[i+1]-tabxi[i])) {largeurunite=abs(tabxi[i+1]-tabxi[i]);iclasse=i;}}
|
||||
|
||||
// Hauteur du rectangle le plus haut pour placer l'unité au dessus.
|
||||
real hauteurmaxi=max(tabhi);
|
||||
|
||||
// Calcul de la hauteur à donner au rectangle unité d'aire
|
||||
real hauteurunite=(uniteaire/tabni[iclasse])*tabhi[iclasse];
|
||||
|
||||
//tracé des axes
|
||||
axes ax=calculaxes(tabxi,tabhi,origine);
|
||||
traceaxe(ax,axedesy,true,origine,libellecaractere,"hauteurs");
|
||||
|
||||
// Hauteur du rectangle le plus haut pour placer l'unité au dessus.
|
||||
real hauteurmaxi=max(tabhi);
|
||||
|
||||
// Calcul de la hauteur à donner au rectangle unité d'aire
|
||||
real hauteurunite=(uniteaire/tabni[iclasse])*tabhi[iclasse];
|
||||
|
||||
// Tracé de l'histogramme
|
||||
for (int i=0; i < tabhi.length;++i){
|
||||
//tracé d'un rectangle
|
||||
filldraw(box((tabxi[i],origine.y),(tabxi[i+1],tabhi[i])),p1);
|
||||
draw((tabxi[i],origine.y)--(tabxi[i],tabhi[i])--(tabxi[i+1],tabhi[i])--(tabxi[i+1],origine.y),p2);}
|
||||
|
||||
// Tracé de l'unité d'aire et son étiquette
|
||||
if (unitedaire) {filldraw(box((0.05*ax.max.x,1.1*hauteurmaxi),(0.05*ax.max.x+largeurunite,1.1*hauteurmaxi+hauteurunite)),p1,p2);
|
||||
label(scale(0.8)*libelleunite,(0.05*ax.max.x+1.2*largeurunite,1.1*hauteurmaxi+0.5*hauteurunite),E);}
|
||||
|
||||
|
||||
}
|
||||
|
||||
void diagrammecirculaire(string[] tabmod,real[] tabeff,string nomdugraph=""){
|
||||
|
||||
real[] tabangle,tabanglecumule,tabanglelabel;
|
||||
tabanglecumule[0]=0;
|
||||
int n=tabeff.length;
|
||||
for(int i=0; i<n; ++i) {
|
||||
tabangle[i]=tabeff[i]*360/sum(tabeff);
|
||||
tabanglecumule[i+1]=tabanglecumule[i]+tabangle[i];
|
||||
tabanglelabel[i]=tabanglecumule[i]+tabangle[i]/2;
|
||||
path secteur=(0,0)--arc((0,0),1,tabanglecumule[i],tabanglecumule[i+1])--cycle;
|
||||
filldraw(secteur,couleur[i]);
|
||||
label(scale(0.9)*tabmod[i],(1.6,(i+1)/n),E);
|
||||
filldraw((1.4,(i+0.9)/n)--(1.4,(i+1.2)/n)--(1.6,(i+1.2)/n)--(1.6,(i+0.9)/n)--cycle,couleur[i]);}
|
||||
label(nomdugraph,(0.5,-1.2),S);}
|
||||
|
||||
|
||||
340
gm.asy
Normal file
340
gm.asy
Normal file
@ -0,0 +1,340 @@
|
||||
// gm.asy - Un fourre-tout de fonctions (récentes ou anciennes) pour Asymptote.
|
||||
//
|
||||
// Dernière modification : 04/07/15 pour l'ajout d'une ancienne fonction cote3D(...)
|
||||
// proposée sur le forum le 10/04/10.
|
||||
|
||||
import stats;
|
||||
import geometry;
|
||||
import graph;
|
||||
import three;
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
////////// STATISTIQUES
|
||||
///////////////////////////////////////////////////////////
|
||||
//
|
||||
// Remarque de février 2012 : ce qui suit est amené à disparaitre
|
||||
// Une extension gm_stats devrait voir le jour avant la fin 2012.
|
||||
// Tout sera réécrit depuis une feuille blanche...
|
||||
// ... mais si cela peut dépanner jusque là.
|
||||
//
|
||||
// Rappel de la histogramm défini dans stats.asy :
|
||||
// void histogram(picture pic=currentpicture, real[] bins, real[] count,
|
||||
// real low=-infinity,
|
||||
// pen fillpen=nullpen, pen drawpen=nullpen, bool bars=false,
|
||||
// Label legend="", real markersize=legendmarkersize)
|
||||
|
||||
// Une version perso (avec un e à la fin), plus adaptée pour faire
|
||||
// des histogrammes pour un de math de lycée :
|
||||
void histogramme(picture pic=currentpicture,
|
||||
real[] tabxi, // les xi définissants les classes
|
||||
real[] tabni, // les effectifs (ou fréquence) par classe.
|
||||
bool bars=true,
|
||||
pen p1=lightgray,
|
||||
pen p2=.8bp+blue,
|
||||
string libellecaractere="Valeurs du caract\`ere",
|
||||
real minaxe=min(tabxi),
|
||||
real maxaxe=max(tabxi),
|
||||
real uniteaxe=(maxaxe-minaxe)/4,
|
||||
string libelleunite="",
|
||||
bool afficherUniteAire=true,
|
||||
bool valeurUniteAire=true,
|
||||
bool frequence=false,
|
||||
bool pourcent=false,
|
||||
real uniteaire=sum(tabni)/100){
|
||||
real uniteairetempo;
|
||||
if(frequence==false) uniteairetempo=uniteaire;
|
||||
if(frequence==true && pourcent==false) {
|
||||
if(uniteaire>=1) uniteaire=uniteaire/100;
|
||||
uniteairetempo=uniteaire*sum(tabni);
|
||||
}
|
||||
if(frequence==true && pourcent==true) {
|
||||
uniteairetempo=uniteaire*sum(tabni)/100;
|
||||
}
|
||||
// Une variable utile pour déterminer la plus petite amplitude :
|
||||
real largeurunite=abs(tabxi[1]-tabxi[0]);
|
||||
// ... et une autre pour le numéro de classe correspondant :
|
||||
int iclasse=0;
|
||||
// Calcul des hauteurs (et de la plus petite amplitude de classe) :
|
||||
real[] tabhi;
|
||||
for(int i=0; i < tabni.length; ++i){
|
||||
tabhi[i]=tabni[i]/(tabxi[i+1]-tabxi[i]);
|
||||
if (largeurunite>abs(tabxi[i+1]-tabxi[i])) {
|
||||
largeurunite=abs(tabxi[i+1]-tabxi[i]);
|
||||
iclasse=i;
|
||||
}
|
||||
}
|
||||
// Hauteur du rectangle le plus haut pour placer l'unité au dessus.
|
||||
real hauteurmaxi=max(tabhi);
|
||||
// Calcul de la hauteur à donner au rectangle unité d'aire
|
||||
real hauteurunite=(uniteairetempo/tabni[iclasse])*tabhi[iclasse];
|
||||
|
||||
int nbrlignes = ceil(hauteurmaxi/hauteurunite);
|
||||
for(int k=0; k<=nbrlignes; ++k){
|
||||
draw((minaxe,k*hauteurunite)--(maxaxe,k*hauteurunite),.5bp+gray);
|
||||
}
|
||||
for(int k=0; k<=((maxaxe-minaxe)/largeurunite); ++k){
|
||||
draw((minaxe+k*largeurunite,0)--(minaxe+k*largeurunite,nbrlignes*hauteurunite),.5bp+gray);
|
||||
}
|
||||
// Tracé de l'histogramme
|
||||
histogram(tabxi,tabhi,low=0,bars=bars,p1,p2);
|
||||
// Tracé de l'unité d'aire et son étiquette
|
||||
if(afficherUniteAire){
|
||||
filldraw(shift(truepoint(N+W))*box((0,0),(largeurunite,hauteurunite)),p1,p2);
|
||||
if(valeurUniteAire)
|
||||
label(libelleunite+(string) uniteaire+(pourcent==true ? "\%": ""),truepoint(N+W)+(largeurunite,-hauteurunite/2),E);
|
||||
}
|
||||
// Ajout de l'axe gradué
|
||||
xaxis(libellecaractere, Bottom, minaxe,maxaxe,
|
||||
RightTicks(Label(currentpen+fontsize(6)),
|
||||
Step=uniteaxe),above=true);
|
||||
}
|
||||
|
||||
/////// bam et multibam pour des diagrammes en boite.
|
||||
|
||||
void bam(real[] xi, real h=1, real dh=0, bool lab=true, real Step=1){
|
||||
xi=sort(xi);
|
||||
real n=xi.length;
|
||||
real xmin=min(xi),
|
||||
q1= (floor(.25*n)==.25*n) ? xi[floor(.25*n)-1] : xi[floor(.25*n)],
|
||||
me= (n%2==0) ? (xi[floor(n/2)-1]+xi[floor(n/2)])/2 : xi[floor((n+1)/2)-1],
|
||||
q3= (floor(.75*n)==.75*n) ? xi[floor(.75*n)-1] : xi[floor(.75*n)],
|
||||
xmax=max(xi);
|
||||
xaxis(xmin,xmax,Ticks(Step=Step,Size=2));
|
||||
draw(box((q1,1+dh),(q3,1+dh+h)));
|
||||
draw((me,1+dh)--(me,1+dh+h));
|
||||
real hm=(1+dh+1+dh+h)/2;
|
||||
draw((xmin,hm)--(q1,hm)^^(q3,hm)--(xmax,hm));
|
||||
draw((xmin,hm-.1)--(xmin,hm+.1)^^(xmax,hm-.1)--(xmax,hm+.1));
|
||||
draw((xmin,0)--(xmin,hm)^^(q1,0)--(q1,hm)^^(me,0)--(me,hm)^^(q3,0)--(q3,hm)
|
||||
^^(xmax,0)--(xmax,hm),linetype("4 4"));
|
||||
if(lab){
|
||||
label("$X_{min}$",(xmin,hm),S,UnFill());
|
||||
label("$Q_1$",(q1,hm-h/2),S,UnFill());
|
||||
label("$M$",(me,hm-h/2),S,UnFill());
|
||||
label("$Q_3$",(q3,hm-h/2),S,UnFill());
|
||||
label("$X_{max}$",(xmax,hm),S,UnFill());
|
||||
}
|
||||
}
|
||||
|
||||
picture bam0(real[] xi, real larg=1,
|
||||
bool comment_b=false,
|
||||
bool titre_b =false,
|
||||
Label titre_l ="",
|
||||
pair titre_pos=max(xi),
|
||||
align titre_dir=NoAlign,
|
||||
pen stylo=currentpen, pen backgr=nullpen
|
||||
){
|
||||
picture pic;
|
||||
xi=sort(xi);
|
||||
real n=xi.length;
|
||||
real xmin=min(xi),
|
||||
q1= (floor(.25*n)==.25*n) ? xi[floor(.25*n)-1] : xi[floor(.25*n)],
|
||||
me= (n%2==0) ? (xi[floor(n/2)-1]+xi[floor(n/2)])/2 : xi[floor((n+1)/2)-1],
|
||||
q3= (floor(.75*n)==.75*n) ? xi[floor(.75*n)-1] : xi[floor(.75*n)],
|
||||
xmax=max(xi);
|
||||
filldraw(pic,box((q1,-larg/2),(q3,larg/2)),backgr,stylo);
|
||||
draw(pic,(me,-larg/2)--(me,larg/2),stylo);
|
||||
draw(pic,(xmin,0)--(q1,0)^^(q3,0)--(xmax,0),stylo);
|
||||
draw(pic,(xmin,-.1)--(xmin,.1)^^(xmax,-.1)--(xmax,.1),stylo);
|
||||
if(titre_b) label(pic,titre_l,titre_pos,titre_dir);
|
||||
if(comment_b){
|
||||
label(pic,"$X_{min}$",(xmin,0),S,UnFill());
|
||||
label(pic,"$Q_1$",(q1,-larg/2),S,UnFill());
|
||||
label(pic,"$M$",(me,-larg/2),S,UnFill());
|
||||
label(pic,"$Q_3$",(q3,-larg/2),S,UnFill());
|
||||
label(pic,"$X_{max}$",(xmax,0),S,UnFill());
|
||||
}
|
||||
return pic;
|
||||
}
|
||||
|
||||
void multibam(picture pic=currentpicture,
|
||||
real[][] ij,
|
||||
real[] larg =new real[],
|
||||
real[] esp =new real[],
|
||||
bool tit_b =false,
|
||||
Label[] tit =new Label[],
|
||||
align[] tit_d =new align[],
|
||||
pen[] sty =new pen[],
|
||||
pen[] bkg =new pen[],
|
||||
real Step=0, real step=0, real xmargin=0,
|
||||
bool vertical=false,
|
||||
bool twoaxis=false){
|
||||
picture pict;
|
||||
int n=ij.length;
|
||||
if(larg.length==0)larg.push(1);
|
||||
if(esp.length==0) esp.push(.5);
|
||||
while(larg.length<n) larg.push(larg[larg.length-1]);
|
||||
while(esp.length<n+1)esp.push(esp[esp.length-1]);
|
||||
while(tit.length<n) tit.push("");
|
||||
while(tit_d.length<n) tit_d.push((vertical)?plain.W:plain.E);
|
||||
while(sty.length<n) sty.push(currentpen);
|
||||
while(bkg.length<n) bkg.push(nullpen);
|
||||
real[] espa=new real[n+1];
|
||||
for(int k=0; k<n; ++k){
|
||||
espa[k]=((k!=0)?espa[k-1]+larg[k-1]/2:0)+esp[k]+larg[k]/2;
|
||||
add(pict,shift(0,espa[k])*bam0(
|
||||
ij[k],
|
||||
larg = larg[k],
|
||||
titre_b = false,
|
||||
titre_l = tit[k],
|
||||
titre_dir = tit_d[k],
|
||||
stylo = sty[k],
|
||||
backgr = bkg[k])
|
||||
);
|
||||
}
|
||||
espa[n]=espa[n-1]+larg[n-1]/2+esp[n];
|
||||
if(vertical){
|
||||
add(pic,reflect(Oy)*rotate(90)*pict);
|
||||
addMargins(pic,0,xmargin);
|
||||
xlimits(pic,0,espa[n]);
|
||||
ylimits(pic,min(ij)-xmargin,max(ij)+xmargin);
|
||||
yaxis(pic,BottomTop,
|
||||
Ticks("%",Step=Step,step=step,extend=true,.2bp+dashed));
|
||||
yaxis(pic,Left,Ticks(Step=Step,step=step,Size=2,size=1,pTick=.7bp+blue,ptick=black));
|
||||
if(twoaxis) yaxis(pic,Right,Ticks(Step=Step,step=step,Size=2,size=1,pTick=.7bp+blue,ptick=black));
|
||||
if(tit_b) for(int k=0; k<n; ++k)
|
||||
label(pic,tit[k],(espa[k],min(ij)-xmargin),plain.S);
|
||||
} else {
|
||||
add(pic,pict);
|
||||
addMargins(pic,xmargin);
|
||||
ylimits(pic,0,espa[n]);
|
||||
xlimits(pic,min(ij)-xmargin,max(ij)+xmargin);
|
||||
xaxis(pic,LeftRight,
|
||||
Ticks("%",Step=Step,step=step,extend=true,.2bp+dashed));
|
||||
xaxis(pic,Bottom,Ticks(Step=Step,step=step,Size=2,size=1,pTick=.7bp+blue,ptick=black));
|
||||
if(twoaxis) xaxis(pic,Top,Ticks(Step=Step,step=step,Size=2,size=1,pTick=.7bp+blue,ptick=black));
|
||||
if(tit_b) for(int k=0; k<n; ++k)
|
||||
label(pic,tit[k],(max(ij)+xmargin,espa[k]),plain.E);
|
||||
}
|
||||
}
|
||||
|
||||
void diabandes(real X, real Y,
|
||||
real ymin, real ymax, real ystep,
|
||||
real tickwidth, string yformat,
|
||||
Label LX, Label LY, Label[] LLX,
|
||||
real[] height,
|
||||
pen p=nullpen){
|
||||
draw((0,0)--(0,Y),EndArrow);
|
||||
draw((0,0)--(X,0),EndArrow);
|
||||
label(LX,(X,0),plain.SE,fontsize(9));
|
||||
label(LY,(0,Y),plain.N,fontsize(9));
|
||||
real yscale=Y/(ymax+ystep);
|
||||
for(real y=ymin; y<ymax; y+=ystep) {
|
||||
draw((-tickwidth,yscale*y)--(0,yscale*y));
|
||||
label(format(yformat,y),(-tickwidth,yscale*y),plain.W,fontsize(9));
|
||||
}
|
||||
int n=LLX.length;
|
||||
real xscale=X/(2*n+2);
|
||||
for(int i=0;i<n;++i) {
|
||||
real x=xscale*(2*i+1);
|
||||
path P=(x,0)--(x,height[i]*yscale)--(x+xscale,height[i]*yscale)--(x+xscale,0)--cycle;
|
||||
fill(P,p);
|
||||
draw(P);
|
||||
label(LLX[i],(x+xscale/2),plain.S,fontsize(10));
|
||||
}
|
||||
for(int i=0;i<n;++i)
|
||||
draw((0,height[i]*yscale)--(X,height[i]*yscale),dashed);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
/// TRACER UNE COURBE AVEC V.I.
|
||||
//////////////////////////////////////////////////////////////////
|
||||
|
||||
guide[] graphgm(picture pic=currentpicture, real f(real), real a, real b,
|
||||
int n=ngraph, real T(real)=identity,
|
||||
bool3 cond(real), interpolate join=operator --)
|
||||
{
|
||||
if(T == identity)
|
||||
return graph(join,cond)(new pair(real x) {
|
||||
return (x,pic.scale.y.T(f(pic.scale.x.Tinv(x))));},
|
||||
pic.scale.x.T(a),pic.scale.x.T(b),n);
|
||||
else
|
||||
return graph(join,cond)(new pair(real x) {
|
||||
return Scale(pic,(T(x),f(T(x))));},
|
||||
a,b,n);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
/// SIMILITUDES
|
||||
/// - simili_Ckphi(point C, real k, real phi)
|
||||
/// si sont connus : centre, rapport et angle
|
||||
/// - simili_ABCD(point A, point B, point C, point D)
|
||||
/// si sont connues les images C et D de deux points A et B.
|
||||
//////////////////////////////////////////////////////////////////
|
||||
|
||||
transform simili_Ckphi(point C, real k, real phi){
|
||||
return rotate(phi,C)*scale(k,C);
|
||||
}
|
||||
transform simili_ABCD(point A, point B, point C, point D){
|
||||
point a=(D-C)/(B-A), b=C-a*A, centre=b/(1-a);
|
||||
real rapport=abs(a), angle=degrees(a);
|
||||
return simili_Ckphi(centre,rapport,angle);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
///
|
||||
/// 3D : le nom des fonctions parlent d'elles-même.
|
||||
///
|
||||
//////////////////////////////////////////////////////////////////
|
||||
projection cavaliereX(real k=.5, real angle=45)
|
||||
{
|
||||
transform3 t={{-k*Cos(angle),1,0,0},
|
||||
{-k*Sin(angle),0,1,0},
|
||||
{1,0,0,-1},
|
||||
{0,0,0,1}};
|
||||
return projection((1,Cos(angle)^2,1-Cos(angle)^2),normal=(1,0,0),
|
||||
new transformation(triple,triple,triple) { return transformation(t);});
|
||||
}
|
||||
projection cavaliereX=cavaliereX();
|
||||
projection cavaliereYZ=cavaliereX();
|
||||
|
||||
projection cavaliereY(real k=.5, real angle=45)
|
||||
{
|
||||
transform3 t={{1,k*Cos(angle),0,0},
|
||||
{0,k*Sin(angle),1,0},
|
||||
{0,-1,0,-1},
|
||||
{0,0,0,1}};
|
||||
return projection((Cos(angle)^2,-1,1-Cos(angle)^2),normal=(0,-1,0),
|
||||
new transformation(triple,triple,triple) { return transformation(t);});
|
||||
}
|
||||
projection cavaliereY=cavaliereY();
|
||||
projection cavaliereXZ=cavaliereY();
|
||||
|
||||
projection cavaliereZ(real k=.5, real angle=45)
|
||||
{
|
||||
transform3 t={{1,0,-k*Cos(angle),0},
|
||||
{0,1,-k*Sin(angle),0},
|
||||
{0,0,1,-1},
|
||||
{0,0,0,1}};
|
||||
return projection((Cos(angle)^2,1-Cos(angle)^2,1),up=(0,1,0),normal=(0,0,1),
|
||||
new transformation(triple,triple,triple) { return transformation(t);});
|
||||
}
|
||||
projection cavaliereZ=cavaliereZ();
|
||||
projection cavaliereXY=cavaliereZ();
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
///////
|
||||
/////// Inspirée par la définition de la fonction distance(...)
|
||||
/////// que l'on trouve dans l'extension geometry,
|
||||
/////// une fonction qui permet d'indiquer une distance entre deux triples :
|
||||
///////
|
||||
//////////////////////////////////////////////////////////////////
|
||||
|
||||
void cote3D(picture pic=currentpicture, Label L="", triple A, triple B,
|
||||
triple v=O, real offset=1cm,
|
||||
pen p=currentpen, pen joinpen=dotted, arrowbar3 arrow=Arrows3){
|
||||
triple A=A, B=B, vAB=B-A;
|
||||
path3 g=A--B;
|
||||
if(v==O) v=(abs(vAB.x)>abs(vAB.y))? ((abs(vAB.y)>abs(vAB.z))?Z:Y) : ((abs(vAB.x)>abs(vAB.z))?Z:X);
|
||||
transform3 Tp=shift(offset*v);
|
||||
pic.add(new void(picture f, transform3 t) {
|
||||
picture opic;
|
||||
path3 G=Tp*t*g;
|
||||
Label L=L.copy();
|
||||
draw(opic,L,G,p,arrow);
|
||||
triple Ap=t*A, Bp=t*B;
|
||||
draw(opic,(Ap--Tp*Ap)^^(Bp--Tp*Bp), joinpen);
|
||||
add(f,opic);
|
||||
}, true);
|
||||
}
|
||||
147
gm_probabilitytree.asy
Normal file
147
gm_probabilitytree.asy
Normal file
@ -0,0 +1,147 @@
|
||||
/*
|
||||
AVERTISSEMENT du 05/02/11 :
|
||||
après près de deux ans d'utilisation, je songe à revopir
|
||||
en profondeur cette extension(ette) pour en faire
|
||||
une vraie extension avec davantage de possibilités.
|
||||
Tout est susceptible de changer !
|
||||
*/
|
||||
|
||||
// probabilitytree.asy - v. 1.1 - 19/04/09
|
||||
// Extension destinée à dessiner un arbre de probabilités
|
||||
// C'est une première version, non documentée, en phase de test.
|
||||
// G. Marris
|
||||
|
||||
real DistVertEntreNoeuds = .5cm;
|
||||
real DistHoriEntreNiveaux = 2cm;
|
||||
real HauteurNoeudMinimale = .5cm;
|
||||
pen StyleEvenParDefaut = 1bp+black;
|
||||
pen StyleBackGrParDefaut = invisible;
|
||||
pen StyleProbParDefaut = black;
|
||||
pen StyleBrParDefaut = black;
|
||||
|
||||
struct Noeud {Noeud parent;
|
||||
Noeud[] enfants;
|
||||
Label evenement;
|
||||
Label probabilite;
|
||||
pen coul_pr;
|
||||
pen coul_br;
|
||||
pair lieu;
|
||||
frame cadre;
|
||||
real decalage_y;
|
||||
}
|
||||
|
||||
Noeud addN( Noeud parent = null,
|
||||
Label evenement = "",
|
||||
pen[] evpen={StyleEvenParDefaut,StyleBackGrParDefaut},
|
||||
Label probabilite = null,
|
||||
pen prpen=StyleProbParDefaut,
|
||||
pen brpen=StyleBrParDefaut)
|
||||
{
|
||||
frame fr;
|
||||
Noeud nouvelenfant = new Noeud;
|
||||
if (evpen.length==1) evpen[1]=StyleBackGrParDefaut;
|
||||
label(fr,evenement,evpen[0]);
|
||||
roundbox(fr,
|
||||
xmargin=1mm,
|
||||
Fill(evpen[1]));
|
||||
nouvelenfant.cadre = fr;
|
||||
nouvelenfant.probabilite = probabilite;
|
||||
nouvelenfant.coul_pr = prpen;
|
||||
nouvelenfant.coul_br = brpen;
|
||||
if( parent != null ) {
|
||||
nouvelenfant.parent = parent;
|
||||
parent.enfants.push( nouvelenfant );
|
||||
}
|
||||
return nouvelenfant;
|
||||
}
|
||||
|
||||
real CalculHauteur( int niveau, Noeud noeud )
|
||||
{
|
||||
real hauteurcadrenoeud=(max(noeud.cadre)-min(noeud.cadre)).y;
|
||||
if( noeud.enfants.length > 0 ) {
|
||||
real hauteur[] = new real[noeud.enfants.length];
|
||||
real H = 0;
|
||||
|
||||
for( int i = 0; i < noeud.enfants.length; ++i ) {
|
||||
hauteur[i] = CalculHauteur(niveau+1, noeud.enfants[i] );
|
||||
noeud.enfants[i].lieu =(niveau*DistHoriEntreNiveaux, H-hauteur[i]/2);
|
||||
H -= hauteur[i]+DistVertEntreNoeuds;
|
||||
}
|
||||
real hauteurramifications = (sum(hauteur)
|
||||
+DistVertEntreNoeuds*(hauteur.length-1));
|
||||
for( int i = 0; i < noeud.enfants.length; ++i ) {
|
||||
noeud.enfants[i].decalage_y = hauteurramifications/2;
|
||||
}
|
||||
return max(hauteurcadrenoeud,sum(hauteur)
|
||||
+DistVertEntreNoeuds*(hauteur.length-1));
|
||||
}
|
||||
else {
|
||||
return max(hauteurcadrenoeud,HauteurNoeudMinimale);
|
||||
}
|
||||
}
|
||||
void TracerRecursivement( Noeud noeud, frame f )
|
||||
{
|
||||
pair pos;
|
||||
if( noeud.parent != null )
|
||||
pos = (0, noeud.parent.lieu.y + noeud.decalage_y);
|
||||
else
|
||||
pos = (0, noeud.decalage_y);
|
||||
noeud.lieu += pos;
|
||||
|
||||
noeud.cadre = shift(noeud.lieu)*noeud.cadre;
|
||||
|
||||
if( noeud.parent != null ) {
|
||||
pair dirtrait=noeud.parent.lieu-noeud.lieu;
|
||||
path p = point(noeud.cadre,dirtrait)--point(noeud.parent.cadre,E);
|
||||
draw(f,p,noeud.coul_br);
|
||||
if( noeud.probabilite != null ) {
|
||||
label(f,scale(.9)*noeud.probabilite,
|
||||
relpoint(p,.3),noeud.coul_pr,Fill(1,white));
|
||||
}
|
||||
}
|
||||
|
||||
add( f, noeud.cadre );
|
||||
|
||||
for( int i = 0; i < noeud.enfants.length; ++i )
|
||||
TracerRecursivement( noeud.enfants[i], f );
|
||||
}
|
||||
void TracerArbre( Noeud racine, pair pos=(0,0) )
|
||||
{
|
||||
frame f;
|
||||
racine.lieu = (0,0);
|
||||
CalculHauteur( 1, racine );
|
||||
TracerRecursivement( racine, f );
|
||||
add(f,pos);
|
||||
}
|
||||
void Bernouilli( Label Succes="$S$", Label probS="$p$",
|
||||
Label Echec="$\overline{S}$", Label probE="$q$",
|
||||
int repet=2, pair pos=(0,0) )
|
||||
{
|
||||
int n = 2^repet-1, j;
|
||||
Noeud[] N;
|
||||
N[0] = addN( );
|
||||
for (int i=1; i<=2*n; i+=2) {
|
||||
j=floor((i-1)/2);
|
||||
N[i] = addN( N[j], Succes, probS );
|
||||
N[i+1] = addN( N[j], Echec, probE );
|
||||
}
|
||||
TracerArbre( N[0], pos );
|
||||
}
|
||||
void GagnerAToutPrix( Label Succes="$G$", Label probS="$\frac{1}{2}$",
|
||||
Label Echec="$P$", Label probE="$\frac{1}{2}$",
|
||||
int repet=4, pair pos=(0,0) )
|
||||
{
|
||||
int n = 2*repet;
|
||||
Noeud[] N;
|
||||
N[0] = addN( );
|
||||
for (int i=0; i<=repet; ++i) {
|
||||
N[2*i+1] = addN( N[2*i], Succes, probS );
|
||||
N[2*i+2] = addN( N[2*i], Echec, probE );
|
||||
}
|
||||
TracerArbre( N[0], pos );
|
||||
}
|
||||
|
||||
pair operator cast(Noeud P)
|
||||
{
|
||||
return P.lieu;
|
||||
}
|
||||
157
gm_tableaux.asy
Normal file
157
gm_tableaux.asy
Normal file
@ -0,0 +1,157 @@
|
||||
// Un début de quelque chose qui pourrait être une extension, un jour,
|
||||
// pour tracer des tableaux de signes et/ou variations.
|
||||
|
||||
// MAIS, il y a un grand MAIS : si je devais me lancer vraiment
|
||||
// là dedans... tout serait à repenser pour donner la souplesse
|
||||
// nécessaire à l'ajout de nouvelles possibilités.
|
||||
|
||||
import fontsize;
|
||||
|
||||
restricted bool cr=true;
|
||||
restricted bool decr=false;
|
||||
|
||||
picture tabvar(string var="x",
|
||||
string fonct="f",
|
||||
string deriv="f'(x)",
|
||||
string[] lx,
|
||||
string[] ly,
|
||||
bool3 psv=default,
|
||||
real[] v={0,1},
|
||||
bool affderivee=false,
|
||||
string[] ld={" "},
|
||||
bool m=true,
|
||||
real x1=1.2,
|
||||
real x2=2.5,
|
||||
real y1=1,
|
||||
real y2=2,
|
||||
pen styletrait=1bp+black,
|
||||
pen stylefleche=styletrait) {
|
||||
picture pic;
|
||||
real fx(real x) {return 3/2*x1+x*x2;}
|
||||
real l=x1/20;
|
||||
string d=(m==true)?"$":"";
|
||||
int n=lx.length-1;
|
||||
while (ly.length<lx.length) ly.push(" ");
|
||||
if(psv==cr) {v.delete(); v.push(0);}
|
||||
if(psv==decr) {v.delete(); v.push(1);}
|
||||
if(v.length==1) (v[0]>0)?v.push(0):v.push(1);
|
||||
while (v.length<lx.length) v.push(v[v.length-2]);
|
||||
label(pic,d+var+d,(x1/2,-y1/2));
|
||||
real delta=0;
|
||||
if(affderivee==true) {
|
||||
label(pic,d+deriv+d,(x1/2,-3/2*y1));
|
||||
delta=-y1;
|
||||
draw(pic,(0,-2y1)--(2*x1+n*x2,-2y1),styletrait);
|
||||
for(int i=0; i<ld.length; ++i) {
|
||||
real xd=3/2*x1+i*x2/2;
|
||||
if(ld[i]=="O") {
|
||||
draw(pic,circle((xd,-3*y1/2),y1/7));
|
||||
draw(pic,(xd,-y1)--(xd,-2*y1),styletrait); }
|
||||
else if(ld[i]=="0") {
|
||||
label(pic,d+ld[i]+d,(xd,-3*y1/2));
|
||||
draw(pic,(xd,-y1)--(xd,-2*y1),styletrait); }
|
||||
else if(ld[i]=="VI")
|
||||
draw(pic,(xd-l,-y1)--(xd-l,-2*y1)^^(xd+l,-y1)--(xd+l,-2*y1),styletrait);
|
||||
else label(pic,d+ld[i]+d,(xd,-3*y1/2));
|
||||
}
|
||||
}
|
||||
label(pic,d+fonct+d,(x1/2,-3/2*y1+delta-y2/2));
|
||||
draw(pic,box((0,0),(2*x1+n*x2,-2*y1+delta-y2)),styletrait);
|
||||
draw(pic,(0,-y1)--(2*x1+n*x2,-y1)^^(x1,0)--(x1,-2*y1+delta-y2),styletrait);
|
||||
real ymin=min(v), ymax=max(v);
|
||||
real[] vn;
|
||||
object[] im;
|
||||
for(real y:v) vn.push(-3/2*y1+delta-y2+(y-ymin)/(ymax-ymin)*y2);
|
||||
int decal=0;
|
||||
bool[] relier;
|
||||
for(int i=0; i<n+1; ++i) {
|
||||
label(pic,d+((m==true)?replace(lx[i],"inf","\infty"):lx[i])+d,(fx(i),-y1/2));
|
||||
int pos=rfind(ly[i],"VI");
|
||||
if (pos>-1) {
|
||||
draw(pic,(fx(i)-l,-y1+delta)--(fx(i)-l,-2*y1+delta-y2)^^(fx(i)+l,-y1+delta)--(fx(i)+l,-2*y1+delta-y2),styletrait);
|
||||
string[] lim;
|
||||
lim=split(ly[i],"VI");
|
||||
int zu=1; if(i==0) zu=0;
|
||||
if(i!=0) {
|
||||
im.push(draw(pic,Label(d+((m==true)?replace(lim[0],"inf","\infty"):lim[0])+d,align=W),box,(fx(i),vn[i+decal]),invisible));
|
||||
relier.push(false);
|
||||
}
|
||||
if(i!=n) {
|
||||
im.push(draw(pic,Label(d+((m==true)?replace(lim[1],"inf","\infty"):lim[1])+d,align=E),box,(fx(i),vn[i+decal+zu]),invisible));
|
||||
relier.push(true);
|
||||
}
|
||||
if(i!=0&&i!=n) {
|
||||
++decal;
|
||||
vn.push(vn[vn.length-2]);
|
||||
}
|
||||
}
|
||||
else {
|
||||
im.push(draw(pic,d+((m==true)?replace(ly[i],"inf","\infty"):ly[i])+d,box,(fx(i),vn[i+decal]),invisible));
|
||||
relier.push(true);
|
||||
}
|
||||
}
|
||||
add(pic,new void(picture pic2, transform t) {
|
||||
for(int i=0; i<relier.length-1; ++i){
|
||||
if(relier[i]){
|
||||
pair d=(fx(i+1),vn[i+1])-(fx(i),vn[i]);
|
||||
draw(pic2,point(im[i],d,t)--point(im[i+1],-d,t),stylefleche,Arrow(SimpleHead,size=8*linewidth(stylefleche)));
|
||||
}
|
||||
}
|
||||
});
|
||||
return pic;
|
||||
}
|
||||
|
||||
picture tabsgn(string[][] tab, real ul=1, real kx=1.1, real ky=1.2 ,pen p=1bp+black, bool m=true){
|
||||
picture pic;
|
||||
string d=(m==true)?"$":"";
|
||||
pen stylo=fontsize(ul*10pt);
|
||||
int n=tab.length;
|
||||
object[] obj=new object[tab.length];
|
||||
real larg,haut;
|
||||
for(int k=0; k<n; ++k){
|
||||
obj[k]=object(Label("$"+tab[k][0]+"$",stylo));
|
||||
larg=max(larg,(max(obj[k])-min(obj[k])).x);
|
||||
haut=max(haut,(max(obj[k])-min(obj[k])).y);
|
||||
}
|
||||
string[] lx;
|
||||
lx=split(tab[0][1]," ");
|
||||
object[] lxo=new object[lx.length];
|
||||
real xlarg,xhaut,dxl=1.1;
|
||||
for(int i=0; i<lx.length; ++i){
|
||||
lxo[i]=object(Label(d+((m==true)?replace(lx[i],"inf","\infty"):lx[i])+d,stylo));
|
||||
xlarg=max(xlarg,(max(lxo[i])-min(lxo[i])).x);
|
||||
}
|
||||
larg*=kx/2; haut*=ky/2; xlarg*=kx*1.5;
|
||||
draw(pic,(-larg,haut-2n*haut)--(-larg,haut)--(larg,haut)--(larg,haut-2n*haut),p);
|
||||
for(int k=0; k<n; ++k){
|
||||
add(pic,shift(0,-2k*haut)*obj[k]);
|
||||
draw(pic,(-larg,haut-2(k+1)*haut)--(larg,haut-2(k+1)*haut),p);
|
||||
real L=(lxo.length-1/2)*xlarg;
|
||||
if(k==0) {
|
||||
draw(pic,(larg,haut)--(dxl^2*larg+L,haut)--(dxl^2*larg+L,-haut)--(larg,-haut),p);
|
||||
for(int i=0; i<lxo.length; ++i)
|
||||
add(pic,shift(dxl*larg+(1/4+i)*xlarg,0)*lxo[i]);
|
||||
} else {
|
||||
string[] ly;
|
||||
real l=xlarg/30;
|
||||
ly=split(tab[k][1]," ");
|
||||
for(int i=0; i<ly.length; ++i){
|
||||
real xla = dxl*larg+(1/4+i/2)*xlarg;
|
||||
if(ly[i]=="O") {
|
||||
draw(pic,circle((xla,-k*2haut),haut/2),linewidth(1bp));
|
||||
draw(pic,(xla,-k*2haut-haut)--(xla,-k*2haut+haut),p); }
|
||||
else if(ly[i]=="0") {
|
||||
label(pic,shift(xla,-k*2haut)*(d+ly[i]+d),stylo);
|
||||
draw(pic,(xla,-k*2haut-haut)--(xla,-k*2haut+haut),p); }
|
||||
else if(ly[i]=="|")
|
||||
draw(pic,(xla,-k*2haut-haut)--(xla,-k*2haut+haut),p);
|
||||
else if(ly[i]=="||")
|
||||
draw(pic,(xla-l,-k*2haut-haut)--(xla-l,-k*2haut+haut)
|
||||
^^(xla+l,-k*2haut-haut)--(xla+l,-k*2haut+haut),linewidth(1bp));
|
||||
else label(pic,shift(xla,-k*2haut)*(d+ly[i]+d),stylo);
|
||||
}
|
||||
draw(pic,(dxl^2*larg+L,haut-k*2haut)--(dxl^2*larg+L,-haut-k*2haut)--(larg,-haut-k*2haut),p);
|
||||
}
|
||||
}
|
||||
return pic;
|
||||
}
|
||||
41
jl_lib.asy
Normal file
41
jl_lib.asy
Normal file
@ -0,0 +1,41 @@
|
||||
import graph_pi;
|
||||
import gm_probabilitytree;
|
||||
import gm_tableaux;
|
||||
|
||||
marker croix(pen sColor)
|
||||
{
|
||||
return marker(scale(2)*rotate(45)*cross(4), sColor);
|
||||
}
|
||||
|
||||
|
||||
void loadaxis(real xmin, real xmax, real ymin, real ymax)
|
||||
{
|
||||
//xlimits(xmin,xmax,Crop);
|
||||
//ylimits(ymin,ymax,Crop);
|
||||
|
||||
xaxis(Label("$x$",position=EndPoint, align=.2*SE),
|
||||
xmin=xmin,xmax=xmax,
|
||||
Ticks(scale(.7)*Label(align=E),
|
||||
NoZero,
|
||||
//begin=false,
|
||||
end=false,
|
||||
//beginlabel=false,
|
||||
endlabel=false,
|
||||
Step=1,step=.25,
|
||||
Size=1mm, size=.5mm,
|
||||
pTick=black,ptick=gray),
|
||||
Arrow);
|
||||
|
||||
yaxis(Label("$y$",position=EndPoint, align=.2*NW),
|
||||
ymin=ymin,ymax=ymax,
|
||||
Ticks(scale(.7)*Label(),
|
||||
NoZero,
|
||||
//begin=false,
|
||||
end=false,
|
||||
//beginlabel=false,
|
||||
endlabel=false,
|
||||
Step=1,step=.25,
|
||||
Size=1mm, size=.5mm,
|
||||
pTick=black,ptick=gray),
|
||||
Arrow);
|
||||
}
|
||||
4431
js_polyhedron.asy
Normal file
4431
js_polyhedron.asy
Normal file
File diff suppressed because it is too large
Load Diff
182
js_tube.asy
Normal file
182
js_tube.asy
Normal file
@ -0,0 +1,182 @@
|
||||
// Copyright (c) 2007, author: Jens Schwaiger
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or (at
|
||||
// your option) any later version.
|
||||
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
// INSTALLATION:
|
||||
// Paste this file in the sub-directory $HOME/.asy
|
||||
|
||||
// Code:
|
||||
|
||||
import graph3;
|
||||
|
||||
// similar to roundedpath for 3D
|
||||
guide3 roundedguide(guide3 A, real r=0.2){
|
||||
guide3 rounded;
|
||||
triple before, after, indir, outdir;
|
||||
int len=length(A);
|
||||
bool guideclosed=cyclic(A);
|
||||
if(len<2){return A;};
|
||||
if(guideclosed){rounded=point(point(A,0)--point(A,1),r);}
|
||||
else {rounded=point(A,0);};
|
||||
for(int i=1;i<len;i=i+1){
|
||||
before=point(point(A,i)--point(A,i-1),r);
|
||||
after=point(point(A,i)--point(A,i+1),r);
|
||||
indir=dir(point(A,i-1)--point(A,i),1);
|
||||
outdir=dir(point(A,i)--point(A,i+1),1);
|
||||
rounded=rounded--before{indir}..{outdir}after;
|
||||
}
|
||||
if(guideclosed) {
|
||||
before=point(point(A,0)--point(A,len-1),r);
|
||||
indir=dir(point(A,len-1)--point(A,0),1);
|
||||
outdir=dir(point(A,0)--point(A,1),1);
|
||||
rounded=rounded--before{indir}..{outdir}cycle;}
|
||||
else rounded=rounded--point(A,len);
|
||||
|
||||
return rounded;
|
||||
};
|
||||
|
||||
// returns a triple orthogonal to the triple p
|
||||
triple orthv(triple p=(0,0,1)){
|
||||
if(abs((p.x,p.y))>0)
|
||||
{return unit((-p.y,p.x,0));} else {return (1,0,0);};
|
||||
};
|
||||
// used in constructin the array of Bishop frames
|
||||
triple nextnormal(triple p, triple q){
|
||||
triple nw=p-(dot(p,q)*q);
|
||||
if(abs(nw)<0.001){return p;} else {return unit(nw);}
|
||||
};
|
||||
|
||||
// Bishop frame itself; for closed curves a modification guarantees
|
||||
// smoothness also at the
|
||||
// "closing" position of the guide3 g
|
||||
// tw<>0 means some additional twist (measured in radians)
|
||||
// for closed g twist should be a multiple of 2pi
|
||||
/*
|
||||
See http://ada.math.uga.edu/research/software/tube/tube.html
|
||||
*/
|
||||
triple[][] bframe(guide3 g, int subdiv=20, real tw=0){
|
||||
triple[][] bf=new triple[subdiv+1][3];
|
||||
real lg=arclength(g);
|
||||
for(int i=0;i<subdiv+1;i=i+1){bf[i][0]=dir(g,arctime(g,(i/subdiv)*lg));}
|
||||
bf[0][1]=orthv(bf[0][0]);
|
||||
bf[0][2]=cross(bf[0][0],bf[0][1]);
|
||||
for(int i=1;i<subdiv+1;i=i+1){bf[i][1]=nextnormal(bf[i-1][1],bf[i][0]);
|
||||
bf[i][2]=cross(bf[i][0],bf[i][1]);
|
||||
|
||||
};
|
||||
|
||||
if(cyclic(g)){// Modify frame, such that surface closes smoothly
|
||||
triple[] startframe=new triple[3];
|
||||
triple[] endframe=new triple[3];
|
||||
startframe=bf[0]; endframe=bf[subdiv];
|
||||
|
||||
pair tmp=(dot(endframe[1],startframe[1]),-dot(endframe[2],startframe[1]));
|
||||
real alpha=angle(unit(tmp));
|
||||
for(int i=1;i<subdiv+1;i=i+1){
|
||||
bf[i][1]=rotate(-alpha*180/pi*i/subdiv,bf[i][0])*bf[i][1];
|
||||
bf[i][2]=rotate(-alpha*180/pi*i/subdiv,bf[i][0])*bf[i][2];
|
||||
};
|
||||
};
|
||||
for(int i=1;i<subdiv+1;i=i+1){
|
||||
bf[i][1]=rotate(tw*180/pi*i/subdiv,bf[i][0])*bf[i][1];
|
||||
bf[i][2]=rotate(tw*180/pi*i/subdiv,bf[i][0])*bf[i][2];
|
||||
};
|
||||
return bf;
|
||||
};
|
||||
|
||||
typedef guide crosssec(real);
|
||||
guide cs0(real s){return scale(0.3)*unitcircle;};
|
||||
|
||||
// produces a tubelike surface around g; sc ist the radius of the tube
|
||||
/*
|
||||
See http://ada.math.uga.edu/research/software/tube/tube.html
|
||||
*/
|
||||
surface spacetube(guide3 g, int nx=20, int ny=12,
|
||||
crosssec cs=cs0,
|
||||
real twist=0, bool cover=false)
|
||||
{
|
||||
triple[][] bf=bframe(g,nx,twist);
|
||||
triple[] pt=new triple[];
|
||||
real lg=arclength(g);
|
||||
for(int i=0;i<nx+1;i=i+1){
|
||||
pt[i]=relpoint(g,i/nx);}
|
||||
triple[][] surfc=new triple[nx+1][ny+1];
|
||||
for(int i=0;i<nx+1;i=i+1)
|
||||
for(int j=0;j<ny+1;j=j+1){
|
||||
guide rhox=cs((i/nx));
|
||||
if(cover){
|
||||
if((!cyclic(g))&&(i==0||i==nx)){rhox=(0,0);};};
|
||||
pair prhox=relpoint(rhox,j/ny);
|
||||
real scxx=prhox.x;
|
||||
real scyy=prhox.y;
|
||||
surfc[i][j]=pt[i]+scxx*bf[i][1]+
|
||||
scyy*bf[i][2];
|
||||
};
|
||||
return surface(surfc, new bool[][] {});
|
||||
}
|
||||
|
||||
surface spacetube(guide3 g, int nx=20,
|
||||
path cs,
|
||||
real twist=0, bool cover=false)
|
||||
{
|
||||
surface sf;
|
||||
// path3 sec=path3(cs,ZXplane);
|
||||
triple[][] bf=bframe(g,nx,twist);
|
||||
triple[] pt=new triple[];
|
||||
real lg=length(g), r=abs(max(cs)-min(cs))/2;
|
||||
int n=length(cs);
|
||||
path3 tmp1,tmp2;
|
||||
for(int i=0;i<nx+1;i=i+1) pt[i]=relpoint(g,i/nx);
|
||||
// triple pt1, pt2;
|
||||
// for(int i=0; i < n-1; ++i) {
|
||||
// real S=straightness(g,i);
|
||||
// if(S < epsilon*r) {
|
||||
// pt1=point(g,i);
|
||||
// pt2=point(g,i+1);
|
||||
// triple[][] bf=bframe(subpath(g,i,i+1),3,twist);
|
||||
// for (int k=0; k < (cyclic(cs) ? n : n-1); ++k) {
|
||||
// path sec=subpath(cs,k,k+1);
|
||||
// tmp1=path3(sec,new triple(pair z){return pt1+z.x*bf[1][1]+z.y*bf[1][2];});
|
||||
// tmp2=path3(sec,new triple(pair z){return pt2+z.x*bf[2][1]+z.y*bf[2][2];});
|
||||
// sf.append(surface(tmp1--reverse(tmp2)--cycle));
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// triple[][] surfc=new triple[nx+1][ny+1];
|
||||
path3 tmp1,tmp2;
|
||||
for(int i=0; i < nx; ++i) {
|
||||
// for(int j=0; j < ny-1; ++j) {
|
||||
// if(cover){
|
||||
// if((!cyclic(g))&&(i==0||i==nx)){rhox=(0,0);};};
|
||||
// pair prhox=relpoint(rhox,j/ny);
|
||||
// real scxx=prhox.x;
|
||||
// real scyy=prhox.y;
|
||||
// surfc[i][j]=pt[i]+scxx*bf[i][1]+scyy*bf[i][2];
|
||||
// surface tmp;
|
||||
// path3 tmp1=shift(pt[i])*align(cross(bf[i][1],bf[i][2]))*sec;
|
||||
// path3 tmp2=shift(pt[i+1])*align(cross(bf[i+1][1],bf[i+1][2]))*sec;
|
||||
for (int k=0; k < (cyclic(cs) ? n : n-1); ++k) {
|
||||
path sec=subpath(cs,k,k+1);
|
||||
tmp1=path3(sec,new triple(pair z){return pt[i]+z.x*bf[i][1]+z.y*bf[i][2];});
|
||||
tmp2=path3(sec,new triple(pair z){return pt[i+1]+z.x*bf[i+1][1]+z.y*bf[i+1][2];});
|
||||
// tmp=surface();
|
||||
// sf.append(surface(subpath(tmp1,k,k+1)--subpath(tmp2,k,k+1)--cycle));
|
||||
sf.append(surface(tmp1--reverse(tmp2)--cycle));
|
||||
}
|
||||
// draw(tmp1);
|
||||
// return surface(surfc, new bool[][] {});
|
||||
}
|
||||
return sf;
|
||||
}
|
||||
159
pi_arrows.asy
Normal file
159
pi_arrows.asy
Normal file
@ -0,0 +1,159 @@
|
||||
// Copyright (c) 2007, Philippe Ivaldi.
|
||||
// Last modified: Mon Dec 31 15:29:19 CET 2007
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or (at
|
||||
// your option) any later version.
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
|
||||
arrowhead EdgeHead()
|
||||
{
|
||||
arrowhead oa;
|
||||
oa.head=new path(path g, position position, pen p=currentpen, real size=0,
|
||||
real angle=arrowangle)
|
||||
{
|
||||
if(size == 0) size=arrowsize(p);
|
||||
bool relative=position.relative;
|
||||
real position=position.position.x;
|
||||
if(relative) position=reltime(g,position);
|
||||
path r=subpath(g,position,0.0);
|
||||
pair x=point(r,0);
|
||||
real t=arctime(r,size);
|
||||
pair y=point(r,t);
|
||||
path base=y+2*size*I*dir(r,t)--y-2*size*I*dir(r,t);
|
||||
path left=rotate(-angle,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,r);
|
||||
pair denom=point(right,T[1])-y;
|
||||
real factor=denom != 0 ? length((point(left,T[0])-y)/denom) : 1;
|
||||
path left=rotate(-angle,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,r);
|
||||
return subpath(left,0,T[0])--y--cycle;
|
||||
};
|
||||
return oa;
|
||||
}
|
||||
|
||||
arrowhead EdgeHead=EdgeHead();
|
||||
|
||||
arrowhead EdgeHookHead(real dir=arrowdir, real barb=arrowbarb)
|
||||
{
|
||||
arrowhead oa;
|
||||
oa.head=new path(path g, position position, pen p=currentpen, real size=0,
|
||||
real angle=arrowangle)
|
||||
{
|
||||
if(size == 0) size=arrowsize(p);
|
||||
angle *= arrowhookfactor;
|
||||
bool relative=position.relative;
|
||||
real position=position.position.x;
|
||||
if(relative) position=reltime(g,position);
|
||||
path r=subpath(g,position,0);
|
||||
pair x=point(r,0);
|
||||
real t=arctime(r,size);
|
||||
pair y=point(r,t);
|
||||
path base=y+2*size*I*dir(r,t)--y-2*size*I*dir(r,t);
|
||||
path left=rotate(-angle,x)*r;
|
||||
path right=rotate(angle,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,right);
|
||||
pair denom=point(right,T[1])-y;
|
||||
real factor=denom != 0 ? length((point(left,T[0])-y)/denom) : 1;
|
||||
path left=rotate(-angle,x)*r;
|
||||
path right=rotate(angle*factor,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,right);
|
||||
left=subpath(left,0,T[0]);
|
||||
right=subpath(right,T[1],0);
|
||||
pair pl0=point(left,0), pl1=relpoint(left,1);
|
||||
pair pr0=relpoint(right,0), pr1=relpoint(right,1);
|
||||
pair M=(pl1+pr0)/2;
|
||||
pair v=barb*unit(M-pl0);
|
||||
pl1=pl1+v; pr0=pr0+v;
|
||||
left=pl0{dir(-dir+degrees(M-pl0))}..pl1--M;
|
||||
right=M--pr0..pr1{dir(dir+degrees(pr1-M))};
|
||||
return left--y--cycle;
|
||||
};
|
||||
return oa;
|
||||
}
|
||||
arrowhead EdgeHookHead=EdgeHookHead();
|
||||
|
||||
|
||||
arrowhead EdgeSimpleHead(real dir=arrowdir, real barb=arrowbarb)
|
||||
{
|
||||
arrowhead oa;
|
||||
oa.head=new path(path g, position position, pen p=currentpen, real size=0,
|
||||
real angle=arrowangle)
|
||||
{
|
||||
if(size == 0) size=arrowsize(p);
|
||||
bool relative=position.relative;
|
||||
real position=position.position.x;
|
||||
if(relative) position=reltime(g,position);
|
||||
path r=subpath(g,position,0);
|
||||
pair x=point(r,0);
|
||||
real t=arctime(r,size);
|
||||
pair y=point(r,t);
|
||||
path base=y+2*size*I*dir(r,t)--y-2*size*I*dir(r,t);
|
||||
path left=rotate(-angle,x)*r;
|
||||
path right=rotate(angle,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,right);
|
||||
pair denom=point(right,T[1])-y;
|
||||
real factor=denom != 0 ? length((point(left,T[0])-y)/denom) : 1;
|
||||
path left=rotate(-angle,x)*r;
|
||||
path right=rotate(angle*factor,x)*r;
|
||||
real[] T=arrowbasepoints(base,left,right);
|
||||
return subpath(left,T[0],0);
|
||||
};
|
||||
return oa;
|
||||
}
|
||||
arrowhead EdgeSimpleHead=EdgeSimpleHead();
|
||||
|
||||
|
||||
private real position(position position, real size, path g, bool center)
|
||||
{
|
||||
bool relative=position.relative;
|
||||
real position=position.position.x;
|
||||
if(relative) {
|
||||
position *= arclength(g);
|
||||
if(center) position += 0.5*size;
|
||||
position=arctime(g,position);
|
||||
} else if(center)
|
||||
position=arctime(g,arclength(subpath(g,0,position))+0.5*size);
|
||||
return position;
|
||||
}
|
||||
|
||||
arrowbar EdgeArrows(arrowhead head=EdgeHead,real size=0, real angle=arrowangle,
|
||||
filltype filltype=FillDraw, position position=EndPoint,
|
||||
real space=infinity)
|
||||
{
|
||||
return new bool(picture pic, path g, pen p, margin margin) {
|
||||
pair direction;
|
||||
real sg=sgn(dot(N,space*I*dir(g,length(g)/2)));
|
||||
space = (space == infinity) ? 2*linewidth(p) : space/2;
|
||||
if (sg>=0)
|
||||
{
|
||||
direction=-space*I*dir(g,length(g)/2);
|
||||
sg=1;
|
||||
}
|
||||
else direction=space*I*dir(g,length(g)/2);
|
||||
picture tpic;
|
||||
tpic.add(new void (frame f, transform t) {
|
||||
drawarrow(f,head,t*shift(inverse(t)*(-direction))*g,p,
|
||||
size,sg*angle,filltype,position,true,margin,false);
|
||||
drawarrow(f,head,t*shift(inverse(t)*direction)*reverse(g),p,
|
||||
size,sg*angle,filltype,position,true,margin,false);
|
||||
});
|
||||
tpic.addPath(g,p);
|
||||
real sz=size;
|
||||
real gle=angle;
|
||||
filltype fl=filltype;
|
||||
addArrow(tpic,head,g,p,sz,gle,fl,position(position,size,g,false));
|
||||
add(pic,tpic);
|
||||
return false;
|
||||
};
|
||||
};
|
||||
|
||||
arrowbar EdgeArrows=EdgeArrows();
|
||||
482
pi_base.asy
Normal file
482
pi_base.asy
Normal file
@ -0,0 +1,482 @@
|
||||
// Copyright (c) 2007, Philippe Ivaldi.
|
||||
// Version: $Id: base_pi.asy,v 0.0 "2007/01/27 10:35:52" Philippe Ivaldi Exp $
|
||||
// Last modified: Sun Oct 14 22:36:12 CEST 2007
|
||||
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or (at
|
||||
// your option) any later version.
|
||||
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
// Commentary:
|
||||
|
||||
// THANKS:
|
||||
|
||||
// BUGS:
|
||||
|
||||
// INSTALLATION:
|
||||
|
||||
// Code:
|
||||
|
||||
import math;
|
||||
|
||||
pair userMin(picture pic=currentpicture){return (pic.userMin().x,pic.userMin().y);}
|
||||
pair userMax(picture pic=currentpicture){return (pic.userMax().x,pic.userMax().y);}
|
||||
|
||||
// A rotation in the direction dir limited to [-90,90]
|
||||
// This is useful for rotating text along a line in the direction dir.
|
||||
private transform rotate(explicit pair dir)
|
||||
{
|
||||
real angle=degrees(dir);
|
||||
if(angle > 90 && angle < 270) angle -= 180;
|
||||
return rotate(angle);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.......................Structures......................*
|
||||
// *=======================================================*
|
||||
/*<asyxml><variable type="guide" signature="Straight(... guide[])"><code></asyxml>*/
|
||||
guide Straight(... guide[])=operator --;/*<asyxml></code><documentation></documentation></variable></asyxml>*/
|
||||
/*<asyxml><variable type="guide" signature="Spline(... guide[])"><code></asyxml>*/
|
||||
guide Spline(... guide[])=operator ..;/*<asyxml></code><documentation></documentation></variable></asyxml>*/
|
||||
|
||||
typedef guide interpolate(... guide[]);
|
||||
|
||||
/*<asyxml><struct signature="rational"><code></asyxml>*/
|
||||
struct rational
|
||||
{/*<asyxml></code><documentation>
|
||||
'p' est le numérateur, 'q' est le dénominateur.
|
||||
'ep' est la précision avec laquelle le rationnel a été obtenu dans
|
||||
le cas où il y convertion à partir d'un irrationnel.
|
||||
..................................................
|
||||
'p' is the numerator, 'q' is the denominator.
|
||||
'ep' is the precision with which the rational was obtained in the case of a
|
||||
convertion from irrational.
|
||||
</documentation></asyxml>*/
|
||||
int p;
|
||||
int q;
|
||||
real ep;
|
||||
}/*<asyxml></struct></asyxml>*/
|
||||
|
||||
rational operator init() {return new rational;}
|
||||
|
||||
/*ANCrational(real,real)ANC*/
|
||||
/*<asyxml><function type="rational" signature="rational(real,real)"><code></asyxml>*/
|
||||
rational rational(real x, real ep=1/10^5)
|
||||
{/*<asyxml></code><documentation>Retourne le rationnel qui approxime 'x' tel que 'abs(p/q-x)<=ep'.
|
||||
..................................................
|
||||
Return the rational which approximates 'x' such as
|
||||
'abs(p/q-x)<=ep'.
|
||||
</documentation></function></asyxml>*/
|
||||
rational orat;
|
||||
int q=1;
|
||||
while (abs(round(q*x)-q*x)>ep)
|
||||
{
|
||||
++q;
|
||||
}
|
||||
orat.p=round(q*x);
|
||||
orat.q=q;
|
||||
orat.ep=ep;
|
||||
return orat;
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *...................Calculus routines...................*
|
||||
// *=======================================================*
|
||||
/*<asyxml><function type="int" signature="pgcd(int,int)"><code></asyxml>*/
|
||||
int pgcd(int a, int b)
|
||||
{/*<asyxml></code><documentation>Greatest common divisor.</documentation></function></asyxml>*/
|
||||
int a_=abs(a), b_=abs(b), r=a_;
|
||||
if (b_>a_) {a_=b_; b_=r; r=a_;}
|
||||
while (r>0)
|
||||
{
|
||||
r=a_%b_;
|
||||
a_=b_;
|
||||
b_=r;
|
||||
}
|
||||
return a_;
|
||||
}
|
||||
|
||||
/*<asyxml><function type="int" signature="gcd(int,int)"><code></asyxml>*/
|
||||
int gcd(int a, int b)
|
||||
{/*<asyxml></code><documentation>Greatest common divisor.</documentation></function></asyxml>*/
|
||||
return pgcd(a,b);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.................Extend point routine..................*
|
||||
// *=======================================================*
|
||||
/*<asyxml><function type="pair[]" signature="points(path g,real[])"><code></asyxml>*/
|
||||
pair[] points(path g, real[] t)
|
||||
{/*<asyxml></code><documentation>Extend 'point(path, real)' routine to array of 'real'.</documentation></function></asyxml>*/
|
||||
pair [] op;
|
||||
for (int i=0; i < t.length; ++i) {
|
||||
op.push(point(g,t[i]));
|
||||
}
|
||||
return op;
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="points(path,int[])"><code></asyxml>*/
|
||||
pair [] points(path g, int[] t)
|
||||
{/*<asyxml></code><documentation>Extend 'point(path, int)' routine to array of 'int'.</documentation></function></asyxml>*/
|
||||
pair [] op;
|
||||
for (int i=0; i < t.length; ++i) {
|
||||
op.push(point(g,t[i]));
|
||||
}
|
||||
return op;
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.........................join..........................*
|
||||
// *=======================================================*
|
||||
/*<asyxml><function type="guide" signature="join(pair[],interpolate)"><code></asyxml>*/
|
||||
guide join(pair[] a, interpolate join=operator --)
|
||||
{/*<asyxml></code><documentation></documentation></function></asyxml>*/
|
||||
guide og;
|
||||
for(int i=0; i < a.length; ++i) og=join(og,a[i]);
|
||||
return og;
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *...................Extend intersect....................*
|
||||
// *=======================================================*
|
||||
/*<asyxml><function type="real" signature="intersectp(path,pair,int,real)"><code></asyxml>*/
|
||||
real intersectp(path g, pair a, int n=1, real fuzz=0)
|
||||
{/*<asyxml></code><documentation>Retourne le "temps" par rapport à 'g' du premier point
|
||||
d'intersection de 'g' avec le plus petit cercle de centre 'a'
|
||||
coupant 'g'.
|
||||
La précision du découpage peut être augmentée en augmentant 'n'.
|
||||
..................................................
|
||||
Return the time along 'g' of the first intersection point of the path
|
||||
"g" with the smaller circle centered in 'a' and which is intersecting
|
||||
'g'.
|
||||
The cutting precision is increased by increasing n.</documentation></function></asyxml>*/
|
||||
real r=0;
|
||||
real [] ip=intersect(g,(path)a,fuzz);
|
||||
while (ip.length < 2)
|
||||
{
|
||||
r+=1/(50*n);
|
||||
ip=intersect(g,shift(a)*scale(r)*unitcircle,fuzz);
|
||||
}
|
||||
return ip[0];
|
||||
}
|
||||
|
||||
/*<asyxml><function type="real[]" signature="intersectsv(path,real)"><code></asyxml>*/
|
||||
real[] intersectsv(path p, real x)
|
||||
{/*<asyxml></code><documentation>Retourne les "temps" par rapport à 'g' de tous les points
|
||||
d'intersection de 'g' avec la droite verticale passant par (x,0).
|
||||
..................................................
|
||||
Return the times along "g" of all intersection points of the path
|
||||
"g" with the vertical line passing by (x,0).</documentation></function></asyxml>*/
|
||||
return intersections(p,(x,0),(x,1));
|
||||
}
|
||||
|
||||
/*<asyxml><function type="real[]" signature="intersectsh(path g, real y)"><code></asyxml>*/
|
||||
real[] intersectsh(path p, real y)
|
||||
{/*<asyxml></code><documentation>Retourne les "temps" par rapport à 'g' de tous les points
|
||||
d'intersection de 'g' avec la droite horizontale passant par (0,y).
|
||||
La précision du découpage peut être augmentée en augmentant 'n'.
|
||||
..................................................
|
||||
Return the times along 'g' of all intersection points of the path
|
||||
'g' with the horizontal line passing by (0,y).</documentation></function></asyxml>*/
|
||||
return intersections(p,(0,y),(1,y));
|
||||
}
|
||||
|
||||
/*<asyxml><function type="real[]" signature="intersectsd(path,pair,pair)"><code></asyxml>*/
|
||||
real[] intersectsd(path g, pair a, pair b)
|
||||
{/*<asyxml></code><documentation>Retourne les "temps" par rapport à 'g' de tous les points
|
||||
d'intersection de la demi-droite [ab) avec 'g'.
|
||||
..................................................
|
||||
Return the times along 'g' of all intersection points of the
|
||||
half-line from 'a' towards 'b' with the path 'g'.</documentation></function></asyxml>*/
|
||||
real[] ot, ott;
|
||||
ott=intersections(g,a,b);
|
||||
pair ab=b-a;
|
||||
for(real t:ott) if(dot(point(g,t)-a,ab) >= 0) ot.push(t);
|
||||
return ot;
|
||||
}
|
||||
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="intersectionpointsv(path,real)"><code></asyxml>*/
|
||||
pair[] intersectionpointsv(path g, real x)
|
||||
{/*<asyxml></code><documentation>Retourne tous les points d'intersection de 'g' avec la droite
|
||||
verticale passant par (x,0).
|
||||
..................................................
|
||||
Return all the intersection points of the path
|
||||
"g" with the vertical line passing by (x,0).</documentation></function></asyxml>*/
|
||||
return points(g,intersectsv(g,x));
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="intersectionpointsh(path,real)"><code></asyxml>*/
|
||||
pair[] intersectionpointsh(path g, real y)
|
||||
{/*<asyxml></code><documentation>
|
||||
Retourne tous les points d'intersection de 'g' avec la droite
|
||||
horizontale passant par (0,y).
|
||||
..................................................
|
||||
Return all the intersection points of the path
|
||||
"g" with the horizontal line passing by (0,y).</documentation></function></asyxml>*/
|
||||
return points(g,intersectsh(g,y));/*IDOCpoints(path,real[])IDOC*/
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="intersectionpointsd(path,pair,pair)"><code></asyxml>*/
|
||||
pair[] intersectionpointsd(path g, pair a, pair b)
|
||||
{/*<asyxml></code><documentation>
|
||||
Retourne tous les points d'intersection de la demi-droite [ab) avec
|
||||
'g'.
|
||||
..................................................
|
||||
Return all the intersection points of the
|
||||
half-line from 'a' towards 'b' with the path 'g'.</documentation></function></asyxml>*/
|
||||
return points(g,intersectsd(g,a,b));
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="intersectionpoints(path,pair,pair)"><code></asyxml>*/
|
||||
pair[] intersectionpoints(path g, pair a, pair b)
|
||||
{/*<asyxml></code><documentation>
|
||||
Retourne tous les points d'intersection de la droite (ab) avec
|
||||
'g'.
|
||||
..................................................
|
||||
Return all the intersection points of the line (ab) with the path
|
||||
'g'.</documentation></function></asyxml>*/
|
||||
return points(g,intersections(g,a,b));
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.......................Fractions.......................*
|
||||
// *=======================================================*
|
||||
|
||||
/*<asyxml><function type="string" signature="texfrac(int,int,string,bool,bool,bool,bool)"><code></asyxml>*/
|
||||
string texfrac(int p, int q,
|
||||
string factor="",
|
||||
bool signin=false, bool factorin=true,
|
||||
bool displaystyle=false,
|
||||
bool zero=true)
|
||||
{/*<asyxml></code><documentation> Retourne le code LaTeX pour écrire la fraction p/q*factor.
|
||||
Si 'signin' vaut 'true' le signe '-' est dans la fraction (au
|
||||
numérateur).
|
||||
Si 'displaystyle' vaut 'true' le code est en mode 'displaystyle'.
|
||||
Si 'zero' vaut 'false' et 'p' vaut 0, le code génère 0/p*factor; 0
|
||||
si 'zero' vaut 'true'.
|
||||
..................................................
|
||||
Return the LaTeX code to write the fraction p/q*factor.
|
||||
If 'signin' is 'true' the sign '-' is inside the fraction (within
|
||||
the numerator).
|
||||
If 'displaystyle' is 'true' the code is in mode 'displaystyle'.
|
||||
If 'zero' is 'false' and 'p' is 0, the code generates 0/p*factor; 0
|
||||
if 'zero' is 'true'.</documentation></function></asyxml>*/
|
||||
if (p==0) return (zero ? "$0$" : "");
|
||||
string disp= displaystyle ? "$\displaystyle " : "$";
|
||||
int pgcd=pgcd(p,q);
|
||||
int num= round(p/pgcd), den= round(q/pgcd);
|
||||
string nums;
|
||||
if (num==1)
|
||||
if (factor=="" || (!factorin && (den !=1))) nums="1"; else nums="";
|
||||
else
|
||||
if (num==-1)
|
||||
if (factor=="" || (!factorin && (den !=1))) nums="-1"; else nums="-";
|
||||
else nums= (string) num;
|
||||
if (den==1) return "$" + nums + factor + "$";
|
||||
else
|
||||
{
|
||||
string dens= (den==1) ? "" : (string) den;
|
||||
if (signin || num>0)
|
||||
if (factorin)
|
||||
return disp + "\frac{" + nums + factor + "}{" + (string) dens + "}$";
|
||||
else
|
||||
return disp + "\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
|
||||
else
|
||||
{
|
||||
if (num==-1)
|
||||
if (factor=="" || !factorin) nums="1"; else nums="";
|
||||
else nums=(string)(abs(num));
|
||||
if (factorin)
|
||||
return disp + "-\frac{" + nums + factor + "}{" + (string) dens + "}$";
|
||||
else
|
||||
return disp + "-\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*<asyxml><function type="string" signature="texfrac(rational,string,bool,bool,bool,bool)"><code></asyxml>*/
|
||||
string texfrac(rational x,
|
||||
string factor="",
|
||||
bool signin=false, bool factorin=true,
|
||||
bool displaystyle=false,
|
||||
bool zero=true)
|
||||
{/*<asyxml></code><documentation></documentation></function></asyxml>*/
|
||||
return texfrac(x.p, x.q, factor, signin, factorin, displaystyle, zero);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *......................About paths......................*
|
||||
// *=======================================================*
|
||||
|
||||
/*<asyxml><function type="void" signature="drawline(picture,Label,pair,bool,pair,bool,align,pen,arrowbar,arrowbar,margin,Label,marker)"><code></asyxml>*/
|
||||
void drawline(picture pic=currentpicture, Label L="",pair P, bool dirP, pair Q, bool dirQ,
|
||||
align align=NoAlign, pen p=currentpen,
|
||||
arrowbar arrow=None, arrowbar bar=None, margin margin=NoMargin,
|
||||
Label legend="", marker marker=nomarker)
|
||||
{/*<asyxml></code><documentation>Ajoute les deux paramètres 'dirP' et 'dirQ' à la routine native
|
||||
'drawline' du module 'math'.
|
||||
La segment [PQ] sera prolongé en direction de P si 'dirP=true',
|
||||
en direction de Q si 'dirQ=true'.
|
||||
Si 'dirP=dirQ=true', le comportement est celui du 'drawline' natif.
|
||||
Ajoute tous les autres paramètres de 'draw'.
|
||||
..................................................
|
||||
Add the two parameters 'dirP' and 'dirQ' to the native routine
|
||||
'drawline' of the module 'maths'.
|
||||
Segment [PQ] will be prolonged in direction of P if 'dirP=true', in
|
||||
direction of Q if 'dirQ=true'.
|
||||
If 'dirP=dirQ=true', the behavior is that of the native 'drawline'.
|
||||
Add all the other parameters of 'Draw'.</documentation></function></asyxml>*/
|
||||
pic.add(new void (frame f, transform t, transform, pair m, pair M) {
|
||||
picture opic;
|
||||
// Reduce the bounds by the size of the pen.
|
||||
m -= min(p); M -= max(p);
|
||||
|
||||
// Calculate the points and direction vector in the transformed space.
|
||||
pair z=t*P;
|
||||
pair q=t*Q;
|
||||
pair v=t*Q-z;
|
||||
path g;
|
||||
real cp = dirP ? 1:0;
|
||||
real cq = dirQ ? 1:0;
|
||||
// Handle horizontal and vertical lines.
|
||||
if(v.x == 0) {
|
||||
if(m.x <= z.x && z.x <= M.x)
|
||||
g= dot(v,(z.x,m.y))<0 ?
|
||||
(z.x,z.y+cp*(m.y-z.y))--(z.x,q.y+cq*(M.y-q.y)):
|
||||
(z.x,q.y+cq*(m.y-q.y))--(z.x,z.y+cp*(M.y-z.y));
|
||||
} else if(v.y == 0) {
|
||||
if(m.y <= z.y && z.y <= M.y)
|
||||
g=(m.x,z.y)--(M.x,z.y);
|
||||
g= dot(v,(m.x,z.y))<0 ?
|
||||
(z.x+cp*(m.x-z.x),z.y)--(q.x+cq*(M.x-q.x),z.y):
|
||||
(q.x+cq*(m.x-q.x),z.y)--(z.x+cp*(M.x-z.x),z.y);
|
||||
} else {
|
||||
// Calculate the maximum and minimum t values allowed for the
|
||||
// parametric equation z + t*v
|
||||
real mx=(m.x-z.x)/v.x, Mx=(M.x-z.x)/v.x;
|
||||
real my=(m.y-z.y)/v.y, My=(M.y-z.y)/v.y;
|
||||
real tmin=max(v.x > 0 ? mx : Mx, v.y > 0 ? my : My);
|
||||
real tmax=min(v.x > 0 ? Mx : mx, v.y > 0 ? My : my);
|
||||
pair pmin=z+tmin*v;
|
||||
pair pmax=z+tmax*v;
|
||||
if(tmin <= tmax)
|
||||
g= z+cp*tmin*v--z+(cq==0 ? v:tmax*v);
|
||||
}
|
||||
if (length(g)>0) draw(opic, L=L, g=g, align=align, p=p,
|
||||
arrow=arrow, bar=bar, margin=margin,
|
||||
legend=legend, marker=marker);
|
||||
add(f,opic.fit());
|
||||
});
|
||||
}
|
||||
|
||||
/*<asyxml><function type="void" signature="drawline(picture,Label,path,bool,bool,align,pen,arrowbar,arrowbar,margin,Label,marker)"><code></asyxml>*/
|
||||
void drawline(picture pic=currentpicture, Label L="",path g, bool begin=true, bool end=true,
|
||||
align align=NoAlign, pen p=currentpen,
|
||||
arrowbar arrow=None, arrowbar bar=None, margin margin=NoMargin,
|
||||
Label legend="", marker marker=nomarker)
|
||||
{/*<asyxml></code><documentation></documentation></function></asyxml>*/
|
||||
drawline(pic, L, point(g,0), begin, point(g,length(g)), end,
|
||||
align, p, arrow, bar, margin, legend, marker);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *....................Rotated labels.....................*
|
||||
struct rotatedLabel
|
||||
{
|
||||
Label L;
|
||||
};
|
||||
|
||||
rotatedLabel rotatedLabel(string s, string size="",
|
||||
align align=NoAlign,
|
||||
pen p=nullpen, filltype filltype=NoFill)
|
||||
{
|
||||
rotatedLabel OL;
|
||||
OL.L.init(s,size,align,p,Rotate,filltype);
|
||||
return OL;
|
||||
}
|
||||
|
||||
rotatedLabel rotatedLabel(Label L, explicit position position,
|
||||
align align=NoAlign,
|
||||
pen p=nullpen, filltype filltype=NoFill)
|
||||
{
|
||||
rotatedLabel OL;
|
||||
OL.L=Label(L,align,p,Rotate,filltype);
|
||||
OL.L.position(position);
|
||||
return OL;
|
||||
}
|
||||
|
||||
rotatedLabel rotatedLabel(Label L, pair position,
|
||||
align align=NoAlign,
|
||||
pen p=nullpen, filltype filltype=NoFill)
|
||||
{
|
||||
return rotatedLabel(L,(position) position,align,p,filltype);
|
||||
}
|
||||
|
||||
void draw(picture pic=currentpicture, rotatedLabel L, path g, align align=NoAlign,
|
||||
pen p=currentpen, arrowbar arrow=None, arrowbar bar=None,
|
||||
margin margin=NoMargin, Label legend="", marker marker=nomarker)
|
||||
{
|
||||
Label LL=L.L.copy();
|
||||
bool relative=LL.position.relative;
|
||||
real position=LL.position.position.x;
|
||||
if(LL.defaultposition) {relative=true; position=0.5;}
|
||||
if(relative) position=reltime(g,position);
|
||||
LL.embed=Rotate(rotate(dir(g,position))*(1,0));
|
||||
LL.align.dir=LL.embed(identity())*LL.align.dir;
|
||||
align lalign=align.copy();
|
||||
lalign.dir=dir(g,position)*align.dir;
|
||||
draw(pic, LL, g, lalign, p, arrow, bar, margin, legend, marker);
|
||||
}
|
||||
// *...................End rotatedLabel....................*
|
||||
// *=======================================================*
|
||||
|
||||
|
||||
/*<asyxml><function type="void" signature="finalbounds(picture,pen)"><code></asyxml>*/
|
||||
void finalbounds(picture pic=currentpicture,pen p=currentpen)
|
||||
{/*<asyxml></code><documentation>Write the final bounding box of picture 'pic'.
|
||||
This routine is useful to determine the right top and left bottom
|
||||
point for enlarging manually the bounding box.</documentation></function></asyxml>*/
|
||||
pic.add(new void (frame f, transform t, transform, pair m, pair M) {
|
||||
// Reduce the bounds by the size of the pen and the margins.
|
||||
m += min(p); M -= max(p);
|
||||
transform T=inverse(t);
|
||||
write("box("+(string)(T*m)+", "+(string)(T*M)+")");
|
||||
},true);
|
||||
}
|
||||
|
||||
/*<asyxml><function type="bool" signature="isPrime(int)"><code></asyxml>*/
|
||||
bool isPrime(int num)
|
||||
{/*<asyxml></code><documentation></documentation></function></asyxml>*/
|
||||
if (num == 2)
|
||||
return true;
|
||||
else if (num % 2 == 0)
|
||||
return false;
|
||||
else
|
||||
{
|
||||
bool prime = true;
|
||||
int divisor = 3;
|
||||
int upperLimit = ceil(sqrt(num) + 1);
|
||||
while (divisor <= upperLimit)
|
||||
{
|
||||
if (num % divisor == 0)
|
||||
prime = false;
|
||||
divisor += 2;
|
||||
}
|
||||
return prime;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
60
pi_edvenn.asy
Normal file
60
pi_edvenn.asy
Normal file
@ -0,0 +1,60 @@
|
||||
// Copyright (c) 2007, Philippe Ivaldi.
|
||||
// Version: $Id: edvenn_pi.asy,v 0.0 2007/02/19 17:49:31 Philippe Ivaldi Exp $
|
||||
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or (at
|
||||
// your option) any later version.
|
||||
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
// Commentary:
|
||||
|
||||
// THANKS:
|
||||
|
||||
// BUGS:
|
||||
|
||||
// INSTALLATION:
|
||||
|
||||
// Code:
|
||||
|
||||
|
||||
// Venn diagram // Diagramme de Venn
|
||||
// Edwards' construction // Construction d'Edwards
|
||||
import roundedpath;
|
||||
size(10cm,0);
|
||||
path [] EdVenn(int n)
|
||||
{
|
||||
path [] opath;
|
||||
if (n>=1)
|
||||
opath.push(shift(-1.4,-.9)*roundedpath(xscale(2.8)*yscale(.9)*unitsquare,.1));
|
||||
if (n>=2)
|
||||
opath.push(shift(0,-.9)*roundedpath(xscale(1.4)*yscale(1.8)*unitsquare,.1));
|
||||
if (n>=3)
|
||||
opath.push(scale(.5)*unitcircle);
|
||||
for (int i=1; i<=n-3; ++i)
|
||||
{
|
||||
pair pcle=point(opath[2],1/(2^i)),
|
||||
ccle=intersectionpoint(pcle--(pcle-dir(opath[2],1/(2^i))), (0,0)--(1,0));
|
||||
path cle=shift(ccle)*scale(abs(pcle-ccle))*unitcircle;
|
||||
real[] p1=intersect(cle, opath[2]);
|
||||
path ocle=subpath(cle,-p1[0],p1[0]);
|
||||
guide tpath;
|
||||
real step=360/(2^i), a=0;
|
||||
for (int j=0; j<2^i; ++j)
|
||||
{
|
||||
tpath=tpath..rotate(a)*ocle;
|
||||
a+=step;
|
||||
}
|
||||
opath.push(tpath..cycle);
|
||||
}
|
||||
return opath;
|
||||
}
|
||||
748
pi_graph.asy
Normal file
748
pi_graph.asy
Normal file
@ -0,0 +1,748 @@
|
||||
// Copyright (c) 2007, Philippe Ivaldi.
|
||||
// Version: $Id: graph_pi.asy,v 0.0 "2007/01/27 10:35:52" Philippe Ivaldi Exp $
|
||||
// Last modified: Fri Mar 28 15:57:10 CET 2008
|
||||
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or (at
|
||||
// your option) any later version.
|
||||
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
// Commentary:
|
||||
|
||||
// THANKS:
|
||||
|
||||
// BUGS:
|
||||
|
||||
// INSTALLATION:
|
||||
// copier ce fichier dans le sous-répertoire $HOME/.asy
|
||||
// Move this file in the sub-directory $HOME/.asy
|
||||
|
||||
//CODE:
|
||||
|
||||
import graph;
|
||||
import base_pi;
|
||||
import markers;
|
||||
usepackage("mathrsfs");
|
||||
|
||||
int None=0, onX=1, onY=2, onXY=3;
|
||||
// Real functions
|
||||
typedef real realfunction(real);
|
||||
|
||||
struct graphicrules
|
||||
{// used to comunicate graphicrules to cartesianaxis.
|
||||
real xmin,xmax;
|
||||
real ymin,ymax;
|
||||
bool xcrop,ycrop;
|
||||
void set(picture pic=currentpicture){
|
||||
xlimits(pic,xmin, xmax, xcrop);
|
||||
ylimits(pic,ymin, ymax, ycrop);
|
||||
}
|
||||
};
|
||||
graphicrules init() {return new graphicrules;}
|
||||
graphicrules graphicrules = new graphicrules;
|
||||
void graph_pi_exitfunction(){};
|
||||
|
||||
/*ANCgraphicrulesANC*/
|
||||
void graphicrules(picture pic=currentpicture, real unit=1cm,
|
||||
real xunit=unit != 0 ? unit : 0,
|
||||
real yunit=unit != 0 ? unit : 0,
|
||||
real xmin=-infinity, real xmax=infinity, real ymin=-infinity, real ymax=infinity,
|
||||
bool crop=NoCrop, bool xcrop=crop, bool ycrop=crop)
|
||||
{
|
||||
graphicrules.xmin=xmin;
|
||||
graphicrules.xmax=xmax;
|
||||
graphicrules.ymin=ymin;
|
||||
graphicrules.xmax=xmax;
|
||||
graphicrules.ymax=ymax;
|
||||
graphicrules.xcrop=xcrop;
|
||||
graphicrules.ycrop=ycrop;
|
||||
pic.unitsize(x=xunit,y=yunit);
|
||||
graphicrules.set(pic);
|
||||
graph_pi_exitfunction = new void() {
|
||||
graphicrules.set(pic);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
ticklabel NoZero(string s=defaultformat) {
|
||||
return new string(real x) {
|
||||
if (x!=0) return format(s,x);
|
||||
else return "";
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
ticklabel labelfrac(real ep=1/10^5, real factor=1,
|
||||
string symbol="",
|
||||
bool signin=false, bool symbolin=true,
|
||||
bool displaystyle=false,
|
||||
bool zero=true)
|
||||
{
|
||||
return new string(real x) {
|
||||
return texfrac(rational(x/factor), symbol, signin, symbolin, displaystyle, zero);
|
||||
};
|
||||
}
|
||||
|
||||
ticklabel labelfrac=labelfrac();
|
||||
|
||||
// *=======================================================*
|
||||
// *....................Graph paper....................*
|
||||
// *=======================================================*
|
||||
/*ANCmillimeterpaperANC*/
|
||||
picture millimeterpaper(picture pic=currentpicture, pair O=(0,0),
|
||||
real xmin=infinity, real xmax=infinity,
|
||||
real ymin=infinity, real ymax=infinity,
|
||||
pen p=.5bp+orange)
|
||||
{
|
||||
picture opic;
|
||||
real
|
||||
cofx = pic.xunitsize/cm,
|
||||
cofy = pic.yunitsize/cm;
|
||||
real
|
||||
xmin = (xmin == infinity) ? pic.userMin().x*cofx : xmin*cofx,
|
||||
xmax = (xmax == infinity) ? pic.userMax().x*cofx : xmax*cofx,
|
||||
ymin = (ymin == infinity) ? pic.userMin().y*cofy : ymin*cofy,
|
||||
ymax = (ymax == infinity) ? pic.userMax().y*cofy : ymax*cofy;
|
||||
path
|
||||
ph = (xmin*cm, 0)--(xmax*cm, 0),
|
||||
pv = (0, ymin*cm)--(0, ymax*cm);
|
||||
real [] step={5, 1, .5, .1};
|
||||
pen [] p_={ p, scale(.7)*p, scale(.4)*p, scale(.2)*p};
|
||||
|
||||
for (int j=0; j<4; ++j)
|
||||
{
|
||||
for (real i=O.y; i <= ymax; i += step[j]) {
|
||||
draw(opic, shift(0, i*cm) * ph, p_[j]);
|
||||
}
|
||||
for (real i=O.y; i >= ymin ; i -= step[j]) {
|
||||
draw(opic, shift(0, i*cm) * ph, p_[j]);
|
||||
}
|
||||
for (real i=O.x; i <= xmax; i += step[j]) {
|
||||
draw(opic, shift(i*cm, 0) * pv, p_[j]);
|
||||
}
|
||||
for (real i=O.x; i >= xmin; i -= step[j]) {
|
||||
draw(opic, shift(i*cm, 0) * pv, p_[j]);
|
||||
}
|
||||
}
|
||||
|
||||
return opic;
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.....................Axis and grid.....................*
|
||||
// *=======================================================*
|
||||
/*ANCgridANC*/
|
||||
void grid(picture pic=currentpicture,
|
||||
real xmin=pic.userMin().x, real xmax=pic.userMax().x,
|
||||
real ymin=pic.userMin().y, real ymax=pic.userMax().y,
|
||||
real xStep=1, real xstep=.5,
|
||||
real yStep=1, real ystep=.5,
|
||||
pen pTick=nullpen, pen ptick=grey, bool above=false)
|
||||
{
|
||||
draw(pic,box((xmin,ymin),(xmax,ymax)),invisible);
|
||||
xaxis(pic, BottomTop, xmin, xmax,
|
||||
Ticks("%",extend=true,Step=xStep,step=xstep,pTick=pTick,ptick=ptick),
|
||||
above=above, p=nullpen);
|
||||
yaxis(pic, LeftRight, ymin, ymax,
|
||||
Ticks("%",extend=true,Step=yStep,step=ystep,pTick=pTick,ptick=ptick),
|
||||
above=above, p=nullpen);
|
||||
}
|
||||
|
||||
/*ANCcartesianaxisANC*/
|
||||
void cartesianaxis(picture pic=currentpicture,
|
||||
Label Lx=Label("$x$",align=2S),
|
||||
Label Ly=Label("$y$",align=2W),
|
||||
real xmin=-infinity, real xmax=infinity,
|
||||
real ymin=-infinity, real ymax=infinity,
|
||||
real extrawidth=1, real extraheight=extrawidth,
|
||||
pen p=currentpen,
|
||||
ticks xticks=Ticks("%",pTick=nullpen, ptick=grey),
|
||||
ticks yticks=Ticks("%",pTick=nullpen, ptick=grey),
|
||||
bool viewxaxis=true,
|
||||
bool viewyaxis=true,
|
||||
bool above=true,
|
||||
arrowbar arrow=Arrow)
|
||||
{
|
||||
graphicrules.set(pic);
|
||||
xmin=(xmin == -infinity) ? pic.userMin().x : xmin;
|
||||
xmax=(xmax == infinity) ? pic.userMax().x : xmax;
|
||||
ymin=(ymin == -infinity) ? pic.userMin().y : ymin;
|
||||
ymax=(ymax == infinity) ? pic.userMax().y : ymax;
|
||||
extraheight= pic.yunitsize != 0 ? cm*extraheight/(2*pic.yunitsize) : 0;
|
||||
extrawidth = pic.xunitsize != 0 ? cm*extrawidth/(2*pic.xunitsize) : 0;
|
||||
if (viewxaxis)
|
||||
{
|
||||
yequals(pic, Lx, 0, xmin-extrawidth, xmax+extrawidth, p, above, arrow=arrow);
|
||||
yequals(pic, 0, xmin, xmax, p, xticks, above);
|
||||
}
|
||||
if (viewyaxis)
|
||||
{
|
||||
xequals(pic, Ly, 0, ymin-extraheight, ymax+extraheight, p, above, arrow=arrow);
|
||||
xequals(pic, 0, ymin, ymax, p, yticks, above);
|
||||
}
|
||||
}
|
||||
|
||||
real labelijmargin=1;
|
||||
|
||||
/*ANClabeloijANC*/
|
||||
void labeloij(picture pic=currentpicture,
|
||||
Label Lo=Label("$O$",NoFill),
|
||||
Label Li=Label("$\overrightarrow{\imath}$",NoFill),
|
||||
Label Lj=Label("$\overrightarrow{\jmath}$",NoFill),
|
||||
pen p=scale(2)*currentpen,
|
||||
pair diro=SW, pair diri=labelijmargin*S, pair dirj=labelijmargin*1.5*W,
|
||||
filltype filltype=NoFill, arrowbar arrow=Arrow(2mm),
|
||||
marker marker=dot)
|
||||
{
|
||||
if (Lo.filltype==NoFill) Lo.filltype=filltype;
|
||||
if (Li.filltype==NoFill) Li.filltype=filltype;
|
||||
if (Lj.filltype==NoFill) Lj.filltype=filltype;
|
||||
labelx(pic, Lo, 0, diro, p);
|
||||
draw(pic, Li, (0,0)--(1,0), diri, p, arrow);
|
||||
draw(pic, Lj, (0,0)--(0,1), dirj, p, arrow);
|
||||
if(marker != nomarker) draw(pic, (0,0), p, marker);
|
||||
}
|
||||
|
||||
real labelIJmargin=1;
|
||||
|
||||
/*ANClabeloIJANC*/
|
||||
void labeloIJ(picture pic=currentpicture,
|
||||
Label Lo=Label("$O$",NoFill),
|
||||
Label LI=Label("$I$",NoFill),
|
||||
Label LJ=Label("$J$",NoFill),
|
||||
pair diro=SW, pair dirI=labelIJmargin*S, pair dirJ=labelIJmargin*W,
|
||||
pen p=currentpen,
|
||||
filltype filltype=NoFill,
|
||||
marker marker=dot)
|
||||
{
|
||||
if (Lo.filltype==NoFill) Lo.filltype=filltype;
|
||||
if (LI.filltype==NoFill) LI.filltype=filltype;
|
||||
if (LJ.filltype==NoFill) LJ.filltype=filltype;
|
||||
labelx(pic, LI, 1, dirI, p);
|
||||
labely(pic, LJ, 1, dirJ, p);
|
||||
labelx(pic, Lo, 0, diro, p);
|
||||
if(marker != nomarker) draw(pic, (0,0), p, marker);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *....................Recursivegraph.....................*
|
||||
// *=======================================================*
|
||||
typedef void recursiveroutime(picture, real F(real), real, int, int,
|
||||
Label, align, pen, arrowbar, arrowbar,
|
||||
margin, Label, marker);
|
||||
|
||||
/*ANCrecursiveoptionANC*/
|
||||
recursiveroutime recursiveoption(Label L="u",
|
||||
bool labelbegin=true,
|
||||
bool labelend=true,
|
||||
bool labelinner=true,
|
||||
bool labelalternate=false,
|
||||
string format="",
|
||||
int labelplace=onX,
|
||||
pen px=nullpen,
|
||||
pen py=nullpen,
|
||||
bool startonyaxis=false,
|
||||
arrowbar circuitarrow=None,
|
||||
marker automarker=marker(cross(4)),
|
||||
marker xaxismarker=nomarker,
|
||||
marker yaxismarker=nomarker,
|
||||
marker xmarker=nomarker,
|
||||
marker fmarker=nomarker)
|
||||
{
|
||||
return new void(picture pic, real F(real), real u0, int n0, int n,
|
||||
Label L_, align align,
|
||||
pen p, arrowbar arrow, arrowbar bar,
|
||||
margin margin, Label legend, marker marker_)
|
||||
{
|
||||
real [] u;
|
||||
u[n0]=u0;
|
||||
for(int i=n0+1;i<n;++i) u[i]=F(u[i-1]);
|
||||
guide g= (labelplace==2 || labelplace==3 || startonyaxis) ? nullpath : (u[n0],0);
|
||||
bool addlabelautomark;
|
||||
for(int i=n0; i < n-1; ++i) {
|
||||
g=g--(u[i],u[i+1])--(u[i+1],u[i+1]);
|
||||
}
|
||||
g = ((labelplace==2 || labelplace==3 || startonyaxis) && u0<0) ? (0,u0)--(u0,u0)--(u0,0)--g : g;
|
||||
draw(pic, L_,g, align, p, arrow, bar, margin, legend, marker_);
|
||||
if (circuitarrow!=None)
|
||||
{
|
||||
if (labelplace==2 || labelplace==3 || startonyaxis)
|
||||
{
|
||||
draw(pic, (0,u[n0])--(u[n0],u[n0]), p, circuitarrow);
|
||||
draw(pic, (u[n0],u[n0])--(u[n0],u[n0+1]), p, circuitarrow);
|
||||
}
|
||||
else
|
||||
{
|
||||
draw(pic, (u[n0],0)--(u[n0],u[n0+1]), p, circuitarrow);
|
||||
}
|
||||
draw(pic, (u[n0],u[n0+1])--(u[n0+1],u[n0+1]), p, circuitarrow);
|
||||
for(int i=n0+1; i < n-1; ++i) {
|
||||
draw(pic, (u[i],u[i])--(u[i],u[i+1]), p, circuitarrow);
|
||||
draw(pic, (u[i],u[i+1])--(u[i+1],u[i+1]), p, circuitarrow);
|
||||
}
|
||||
}
|
||||
if (px==nullpen)
|
||||
if (labelplace==onX || labelplace==onXY) px=dotted; else px=invisible;
|
||||
if (py==nullpen)
|
||||
if (labelplace==onY || labelplace==onXY) py=dotted; else py=invisible;
|
||||
Label L=L.copy();
|
||||
|
||||
bool Le=(L.s=="");
|
||||
bool fe= (format=="");
|
||||
bool Lno=(L.s=="%");
|
||||
string Ls=L.s;
|
||||
for (int i=n0; i<n; ++i)
|
||||
{
|
||||
addlabelautomark=((i==n0 && labelbegin) || (i!=n0 && labelinner) || (i==n-1 && labelend));
|
||||
if (i>n0)
|
||||
{
|
||||
if (px!=invisible)
|
||||
{
|
||||
draw(pic,(u[i],u[i])--(u[i],0),px);
|
||||
if (addlabelautomark) add(pic,automarker.f, (u[i],u[i]));
|
||||
}
|
||||
}
|
||||
if (i>n0)
|
||||
{
|
||||
if (py!=invisible)
|
||||
{
|
||||
if (addlabelautomark) add(pic,automarker.f, (u[i-1],u[i]));
|
||||
draw(pic,(u[i-1],u[i])--(0,u[i]),py);
|
||||
}
|
||||
}
|
||||
if (Le)
|
||||
L.s=format(format == "" ? defaultformat : format,u[i]);
|
||||
else if (Lno) L.s=""; else
|
||||
if (fe) L.s="$" + Ls + "_{" + (string) i + "}$"; else
|
||||
L.s="$" + Ls + "_{" + (string) i + "}" + format(format,u[i]) + "$";
|
||||
if (labelplace==1 || labelplace==3)
|
||||
{//Label on xaxis
|
||||
L.position((u[i],0));
|
||||
L.align(L.align,S);
|
||||
if (labelalternate && i!=n0) {L.align(-L.align.dir);}
|
||||
L.s=baseline(L.s,"$p_{1234567890}$");
|
||||
if (addlabelautomark) label(pic,L);
|
||||
if (xaxismarker==nomarker && addlabelautomark) add(pic,automarker.f, (u[i],0));
|
||||
}
|
||||
if (labelplace==2 || labelplace==3)
|
||||
{//Label on yaxis
|
||||
L.position((0,u[i]));
|
||||
L.align(L.align,W);
|
||||
if (labelalternate && i!=n0) {L.align(-L.align.dir);}
|
||||
L.s=baseline(L.s,"$w_{1234567890}$");
|
||||
if (addlabelautomark) label(pic,L);
|
||||
if (yaxismarker==nomarker && addlabelautomark) add(pic,automarker.f, (0,u[i]));
|
||||
}
|
||||
add(pic,xaxismarker.f, (u[i],0));
|
||||
if (i>n0 || startonyaxis) add(pic,yaxismarker.f, (0,u[i]));
|
||||
if (i>n0 || startonyaxis) add(pic,xmarker.f, (u[i],u[i]));
|
||||
if (i>n0) add(pic,fmarker.f, (u[i-1],u[i]));
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
recursiveroutime DefaultRecursiveOption=recursiveoption();
|
||||
|
||||
struct recursivegraph
|
||||
{
|
||||
real f(real);
|
||||
real u0;
|
||||
int n0;
|
||||
int n;
|
||||
recursiveroutime recursiveroutime=DefaultRecursiveOption;
|
||||
void draw(picture pic=currentpicture,
|
||||
Label L, align align,
|
||||
pen p, arrowbar arrow, arrowbar bar,
|
||||
margin margin, Label legend, marker marker)
|
||||
{recursiveroutime(pic, f, u0, n0, n,
|
||||
L, align,
|
||||
p, arrow, bar,
|
||||
margin, legend, marker);};
|
||||
};
|
||||
recursivegraph operator init() {return new recursivegraph;}
|
||||
|
||||
recursivegraph recursivegraph(real F(real), real u0, int n0=0, int n)
|
||||
{
|
||||
recursivegraph orec= new recursivegraph;
|
||||
orec.f=F;
|
||||
orec.u0=u0;
|
||||
orec.n0=n0;
|
||||
orec.n=n;
|
||||
return orec;
|
||||
}
|
||||
|
||||
|
||||
void draw(picture pic=currentpicture, Label L="", recursivegraph g, recursiveroutime lr=DefaultRecursiveOption,align align=NoAlign,
|
||||
pen p=currentpen, arrowbar arrow=None, arrowbar bar=None,
|
||||
margin margin=NoMargin, Label legend="", marker marker=nomarker)
|
||||
{
|
||||
g.recursiveroutime=lr;
|
||||
g.draw(pic, L, align, p, arrow, bar, margin, legend, marker);
|
||||
}
|
||||
|
||||
guide graph(picture pic=currentpicture, real f(real),
|
||||
int n=ngraph, interpolate join=operator --)
|
||||
{
|
||||
return graph(pic, f, a=pic.userMin().x, b=pic.userMax().x, n, join);
|
||||
}
|
||||
|
||||
|
||||
/*ANCgraphpoint(...)ANC*/
|
||||
void graphpoint(picture pic=currentpicture,
|
||||
Label L="",
|
||||
real f(real), real xCoordinate,
|
||||
real xmin=0, real ymin=0,
|
||||
int draw=onXY,
|
||||
pen px=nullpen, pen py=px,
|
||||
arrowbar arrow=None, arrowbar bar=None,
|
||||
margin marginy=NoMargin, margin marginx=NoMargin,
|
||||
bool extend=false, bool extendx=extend, bool extendy=extend,
|
||||
Label legendx="", Label legendy="",
|
||||
marker markerx=nomarker, marker markery=nomarker)
|
||||
{/*DOC Mark a point on a curve defined by real f(real). DOC*/
|
||||
/*EXAgraphpoint(...)EXA*/
|
||||
real xmax,ymax;
|
||||
if (extendy) {
|
||||
xmin=pic.userMin().x;
|
||||
xmax=pic.userMax().x;
|
||||
} else xmax=xCoordinate;
|
||||
if (extendx) {
|
||||
ymin=pic.userMin().y;
|
||||
ymax=pic.userMax().y;
|
||||
} else ymax=f(xCoordinate);
|
||||
px = (px==nullpen) ? currentpen+linetype("6 6") : px;
|
||||
py = (py==nullpen) ? currentpen+linetype("6 6") : py;
|
||||
L.align(L.align,NE);
|
||||
label(pic, L, (xCoordinate,f(xCoordinate)));
|
||||
if (draw==onX || draw==onXY)
|
||||
draw(pic, (xCoordinate,ymin)--(xCoordinate,ymax),
|
||||
p=px, arrow=arrow, bar=bar, margin=marginx, legend=legendx, marker=markerx);
|
||||
if (draw==onY || draw==onXY)
|
||||
draw(pic, (xmin,ymax)--(xmax,ymax),
|
||||
p=py, arrow=arrow, bar=bar, margin=marginy, legend=legendy, marker=markery);
|
||||
}
|
||||
|
||||
// *=======================================================*
|
||||
// *.....................About tangent.....................*
|
||||
// *=======================================================*
|
||||
/*ANCtangentANC*/
|
||||
path tangent(path g, real x, path b=box(userMin(currentpicture),userMax(currentpicture)))
|
||||
{//Return the tangent with the maximun size allowed by b (cyclic path)
|
||||
if (!cyclic(b)) abort("tangent: path b is not a cyclic path...");
|
||||
pair pt=point(g,intersectsv(g,x)[0]);
|
||||
real t=intersectp(g,pt);
|
||||
real rt=intersectp(reverse(g),pt);
|
||||
pair dirr=dir(g,t);
|
||||
pair dll=intersectionpointsd(b,pt,shift(pt)*dirr)[0];
|
||||
pair dlr=intersectionpointsd(b,pt,shift(pt)*(-dirr))[0];
|
||||
return dll--dlr;
|
||||
}
|
||||
|
||||
/*ANCaddtangentANC*/
|
||||
void addtangent(picture pic=currentpicture,
|
||||
path g,
|
||||
pair pt,//Point on the path g
|
||||
real size=infinity,//ABSOLUTE size of the tangent line (infinity=maximun size according the size of the pic)
|
||||
bool drawright=true,//Draw the tangent at the right
|
||||
bool drawleft=true,//... left
|
||||
pair v=(infinity,infinity),//A finite value forces the value of the derivative
|
||||
pair vr=v,//A finite value forces the value of the derivative at right
|
||||
pair vl=v,//A finite value forces the value of the derivative at left
|
||||
arrowbar arrow=null,//null=automatic determination
|
||||
margin margin=NoMargin,//Useful with size=infinity
|
||||
Label legend="",
|
||||
pen p=currentpen,
|
||||
real dt=2,//Increase this number can help to discern tangent at the right and tgt at the left.
|
||||
bool differentiable=true)//Set it "true" maybe useful if you are sure that "g is differentiable" at this point.
|
||||
{
|
||||
arrowbar arrow_=arrow;
|
||||
pair dir_r,dir_l;
|
||||
if (intersect(g,pt).length<2) abort("addtangent: the point is not on the path.");
|
||||
real t=intersectp(g,pt);
|
||||
if (!differentiable) {
|
||||
path subpa,subpb;
|
||||
subpa=subpath(g,0,t-dt/2);
|
||||
subpb=subpath(g,t+dt/2,length(g));
|
||||
dir_l=(vl.x<infinity || vl.y<infinity) ? dir((0,0)--vl,.5) : dir(subpa,length(subpa));
|
||||
dir_r=(vr.x<infinity || vr.y<infinity) ? dir((0,0)--vr,.5) : dir(subpb,0);
|
||||
} else {
|
||||
if (v.x<infinity || v.y<infinity) dir_r=dir((0,0)--v,.5);
|
||||
else if (vr.x<infinity || vr.y<infinity) dir_r=dir((0,0)--vr,.5);
|
||||
else if (vl.x<infinity || vl.y<infinity) dir_r=dir((0,0)--vl,.5);
|
||||
else dir_r=dir(g,t);
|
||||
dir_l=dir_r;
|
||||
}
|
||||
pair dr_a,dl_a;
|
||||
pair dr_b,dl_b;
|
||||
dl_a=shift(pt)*(-dir_l);
|
||||
dl_b=shift(pt)*dir_l;
|
||||
dr_a=shift(pt)*dir_r;
|
||||
dr_b=shift(pt)*(-dir_r);
|
||||
if (size==infinity) {
|
||||
draw(pic,g,invisible);
|
||||
dl_a=intersectionpointsd(box(userMin(currentpicture),userMax(currentpicture)),pt,dl_a)[0];
|
||||
dl_b=intersectionpointsd(box(userMin(currentpicture),userMax(currentpicture)),pt,dl_b)[0];
|
||||
dr_a=intersectionpointsd(box(userMin(currentpicture),userMax(currentpicture)),pt,dr_a)[0];
|
||||
dr_b=intersectionpointsd(box(userMin(currentpicture),userMax(currentpicture)),pt,dr_b)[0];
|
||||
if (arrow_==null) arrow_=None;
|
||||
if (drawright && drawleft) {
|
||||
draw(dl_a--dl_b,p=p,arrow=arrow_,margin=margin, legend=legend);
|
||||
if (!differentiable) draw(dr_a--dr_b,p=p,arrow=arrow_);
|
||||
} else if (drawright) draw(pt--dr_a,p=p,arrow=arrow_,margin=margin, legend=legend);
|
||||
else if (drawleft) draw(pt--dl_a,p=p,arrow=arrow_,margin=margin, legend=legend);
|
||||
} else {//Fixed size
|
||||
dl_a=shift((0,0))*(-size*unit(pic.calculateTransform()*dir_l));
|
||||
dl_b=shift((0,0))*(size*unit(pic.calculateTransform()*dir_l));
|
||||
dr_a=shift((0,0))*(size*unit(pic.calculateTransform()*dir_r));
|
||||
dr_b=shift((0,0))*(-size*unit(pic.calculateTransform()*dir_r));
|
||||
picture pict;
|
||||
if (drawright && drawleft) {
|
||||
if (differentiable) {
|
||||
if (arrow_==null) arrow_=Arrows;
|
||||
draw(pict,dl_a--dl_b,p,arrow=arrow_,margin=margin, legend=legend);
|
||||
} else {
|
||||
if (arrow_==null) arrow_=Arrow;
|
||||
draw(pict,(0,0)--dl_a,p,arrow=arrow_,margin=margin, legend=legend);
|
||||
draw(pict,(0,0)--dr_a,p,arrow=arrow_,margin=margin, legend=legend);
|
||||
}
|
||||
} else {
|
||||
if (arrow_==null) arrow_=Arrow;
|
||||
if (drawleft) draw(pict,(0,0)--dl_a,p,arrow=arrow_,margin=margin, legend=legend);
|
||||
else if (drawright) draw(pict,(0,0)--dr_a,p,arrow=arrow_,margin=margin, legend=legend);
|
||||
}
|
||||
add(pic,pict,pt);
|
||||
}
|
||||
}
|
||||
|
||||
void addtangent(picture pic=currentpicture,
|
||||
path g,
|
||||
real x,//x-Coodinate
|
||||
real size=infinity,//ABSOLUTE size of the tangent line (infinity=maximun size according the size of the pic)
|
||||
bool drawright=true,//Draw the tangent at the right
|
||||
bool drawleft=true,//... left
|
||||
pair v=(infinity,infinity),//A finite (x OR y) value forces the value of the derivative
|
||||
pair vr=v,//A finite value forces the value of the derivative at right
|
||||
pair vl=v,//A finite value forces the value of the derivative at left
|
||||
arrowbar arrow=null,//null=automatic determination
|
||||
margin margin=NoMargin,//Useful with size=infinity
|
||||
Label legend="",
|
||||
pen p=currentpen,
|
||||
real dt=2,//Increase this number can help to discern tangent at the right and tgt at the left.
|
||||
bool differentiable=true)//Set it "true" maybe useful if you are sure that "g is differentiable" at this point.
|
||||
{
|
||||
addtangent(pic,g,point(g,intersectsv(g,x)[0]),size,drawright,drawleft,v,vr,vl,arrow,margin,legend,p,dt,differentiable);
|
||||
}
|
||||
|
||||
|
||||
// *=======================================================*
|
||||
// *.....................Special marks.....................*
|
||||
// *=======================================================*
|
||||
// On picture pic, add to path g the frame f rotated by the direction of path g
|
||||
// at the begin or the end of the path g.
|
||||
markroutine dirmarkextremroutine(bool begin=true, bool end=true) {
|
||||
return new void(picture pic=currentpicture, frame f, path g) {
|
||||
if(!begin && !end) return;
|
||||
else {
|
||||
real [] pos;
|
||||
if (begin){
|
||||
add(pic, rotate(degrees(pic.calculateTransform()*dir(g,arctime(g,0))))*f, point(g,0));
|
||||
}
|
||||
if (end) {
|
||||
add(pic, rotate(180+degrees(pic.calculateTransform()*dir(g,length(g))))*f, relpoint(g,1));
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
// A new marker constructor which uses the markroutine dirmarkendroutine.
|
||||
marker markerextrem(frame f,bool begin=true, bool end=true,bool above=true)
|
||||
{
|
||||
return marker(f=f,markroutine=dirmarkextremroutine(begin=begin,end=end),above=above);
|
||||
}
|
||||
|
||||
real graphmarksize=sqrt(2)*dotsize(currentpen);
|
||||
real graphmarksize(){return graphmarksize;}
|
||||
// *=======================================================*
|
||||
// *................How define a new marker................*
|
||||
// *=======================================================*
|
||||
// 1. Definition of a mark as frame.
|
||||
frame arcpicture(real radius=graphmarksize(), real angle=90, pen p=currentpen)
|
||||
{
|
||||
frame ofr;
|
||||
draw(ofr, shift((-radius,0))*arc((0,0),radius,-angle/2,angle/2),p);
|
||||
return ofr;
|
||||
}
|
||||
|
||||
// 2. Definition of the marker itself
|
||||
marker ArcMarkerExtrem(real radius=graphmarksize(), real angle=180,
|
||||
bool begin=true, bool end=true,
|
||||
pen p=currentpen, bool above=true)
|
||||
{
|
||||
return markerextrem(f=arcpicture(radius=radius,angle=angle,p=p),
|
||||
begin=begin,end=end,above=above);
|
||||
}
|
||||
|
||||
// 3. Definition of an alias to use default values.
|
||||
marker ArcMarkerExtrem=ArcMarkerExtrem();
|
||||
//End of section 'How define a new marker'
|
||||
|
||||
|
||||
// HookMarkerExtrem
|
||||
frame hookpicture(real height=graphmarksize(), real width=height/sqrt(2), real angle=90, pen p=currentpen)
|
||||
{
|
||||
// if (!(width<infinity)) width=graphmarksize()/2;
|
||||
frame ofr;
|
||||
draw(ofr, (0,height)--(0,-height),p);
|
||||
draw(ofr, (0,height)--(-width,height),p);
|
||||
draw(ofr, (0,-height)--(-width,-height),p);
|
||||
return ofr;
|
||||
}
|
||||
|
||||
marker HookMarkerExtrem(real height=graphmarksize(), real width=height/2,
|
||||
bool begin=true, bool end=true,
|
||||
pen p=currentpen, bool above=true)
|
||||
{
|
||||
return markerextrem(f=hookpicture(height=height,width=width,p=p),
|
||||
begin=begin,end=end,above=above);
|
||||
}
|
||||
|
||||
marker HookMarkerExtrem=HookMarkerExtrem();
|
||||
//// End HookMarkerExtrem
|
||||
|
||||
//CircleMarkerExtrem
|
||||
frame circlepicture(real radius=graphmarksize(), filltype filltype=NoFill, pen p=currentpen)
|
||||
{
|
||||
frame ofr;
|
||||
filltype filltype_=filltype;
|
||||
path cle=shift((-radius,0))*scale(radius)*unitcircle;
|
||||
filltype.fill(ofr,cle,p);
|
||||
if (filltype_==NoFill) draw(ofr, cle,p);
|
||||
return ofr;
|
||||
}
|
||||
|
||||
marker CircleMarkerExtrem(real radius=graphmarksize(), real angle=90,
|
||||
bool begin=true, bool end=true,
|
||||
pen p=currentpen, filltype filltype=NoFill,
|
||||
bool above=true)
|
||||
{
|
||||
return markerextrem(f=circlepicture(radius=radius,p=p,filltype=filltype),
|
||||
begin=begin,end=end,above=above);
|
||||
}
|
||||
|
||||
marker CircleMarkerExtrem=CircleMarkerExtrem();
|
||||
// End CircleMarkerExtrem
|
||||
|
||||
|
||||
exitfcn currentexitfunction=atexit();
|
||||
|
||||
atexit(new void() {
|
||||
graph_pi_exitfunction();
|
||||
if(currentexitfunction != null) currentexitfunction();
|
||||
});
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
145
pi_hull.asy
Normal file
145
pi_hull.asy
Normal file
@ -0,0 +1,145 @@
|
||||
// Copyright (c) 2008, Philippe Ivaldi.
|
||||
// Version: : hull_pi.asy,v 0.1 2008/18/11 Philippe Ivaldi Exp $
|
||||
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 3 of the License, or (at
|
||||
// your option) any later version.
|
||||
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
||||
// 02110-1301, USA.
|
||||
|
||||
// Commentary:
|
||||
// Graham scan method of computing a hull nodes of a given set of pairs.
|
||||
|
||||
// THANKS:
|
||||
|
||||
// BUGS:
|
||||
|
||||
// INSTALLATION:
|
||||
// Paste this file in the sub-directory ~/.asy
|
||||
|
||||
// Code:
|
||||
|
||||
/*<asyxml><function type="path" signature="polygon(pair[])"><code></asyxml>*/
|
||||
path polygon(pair[] c)
|
||||
{/*<asyxml></code><documentation>Join the nodes c with segments.</documentation></function></asyxml>*/
|
||||
guide g;
|
||||
for (int k=0; k < c.length; ++k)
|
||||
g=g--c[k];
|
||||
return g--cycle;
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair" signature="pivot(pair[])"><code></asyxml>*/
|
||||
pair pivot(pair[] c) {
|
||||
/*<asyxml></code><documentation>Return the point with the lowest y-coordinate.
|
||||
If there is a tie, the point with the lowest x-coordinate out of the tie breaking candidates is returned.</documentation></function></asyxml>*/
|
||||
real[][] coords;
|
||||
for (int i=0; i < c.length; ++i) coords.push(new real[] {c[i].y,c[i].x,i});
|
||||
return c[round(sort(coords)[0][2])];
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature=""><code></asyxml>*/
|
||||
pair[] polarSort(pair[] c, int pivot=-1)
|
||||
{/*<asyxml></code><documentation>Sort points by the polar angles in ascending order.
|
||||
If pivot < 0, use the pair returned by pivot(c) as origin else c[pivot].</documentation></function></asyxml>*/
|
||||
int n=c.length;
|
||||
if(pivot >= n) abort("polarSort: pivot to large.");
|
||||
pair O;
|
||||
real[][] polar;
|
||||
if(pivot < 0) {
|
||||
O=pivot(c);
|
||||
for (int i=0; i < n; ++i)
|
||||
polar.push(new real[] {degrees(c[i]-O,false),abs(c[i]-O),i});
|
||||
polar=sort(polar);
|
||||
} else {
|
||||
O=c[pivot]; // Origine of the polar coordinates system
|
||||
real[][] polp, polm;
|
||||
for (int i=0; i < n; ++i) {
|
||||
real d;
|
||||
if(i != pivot) {
|
||||
d=degrees(c[i]-O,false);
|
||||
if(d > 180) d=d-360;
|
||||
if(d > 0) // 0 <= angle <= 180
|
||||
polp.push(new real[] {d,abs(c[i]-O),i});
|
||||
else // 180 < angle < 0
|
||||
polm.push(new real[] {d,abs(c[i]-O),i});
|
||||
}
|
||||
}
|
||||
|
||||
// Sort the angles in ascending order;
|
||||
polp=sort(polp);
|
||||
polm=sort(polm);
|
||||
|
||||
void append(real[][] a, real[][] b, real[][] c)
|
||||
{
|
||||
polar=copy(a);
|
||||
polar.append(b);
|
||||
polar.append(c);
|
||||
}
|
||||
|
||||
if(polp.length > 0 && polm.length > 0) {
|
||||
pair p1=c[round(polp[0][2])],
|
||||
p2=c[round(polm[polm.length-1][2])],
|
||||
p3=1/3*(O+p1+p2);
|
||||
// We must be careful in rotation to join the paths
|
||||
if((p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) > 0)
|
||||
append(new real[][]{{0,0,pivot}}, polp, polm);
|
||||
else
|
||||
append(new real[][]{{0,0,pivot}}, polm, polp);
|
||||
} else
|
||||
append(new real[][]{{0,0,pivot}}, polp, polm);
|
||||
}
|
||||
return sequence(new pair(int i){return c[round(polar[i][2])];}, n);
|
||||
}
|
||||
|
||||
/*<asyxml><function type="pair[]" signature="hull(pair[],real,real,real,real,int)"><code></asyxml>*/
|
||||
pair[] hull(pair[] c, real depthMin=infinity, real depthMax=0,
|
||||
real angleMin=360, real angleMax=0, int pivot=-1)
|
||||
{/*<asyxml></code><documentation>Graham scan method of computing a hull nodes of a given set of points.
|
||||
With default parameter, return the convex hull.
|
||||
depthMin and depthMax control the minimum and the maximum depth of cracks from the bounding box of c when it's possible.
|
||||
angleMin and angleMax control the minimum and the maximum angle (in degrees) defined by three consecutive points when it's possible.
|
||||
The origin for sorting polar coordinates is the point returned by <html><a href="#pivot(pair[])">pivot(c)</a></html> if pivot < 0 else c[pivot].</documentation></function></asyxml>*/
|
||||
pair minb, maxb, center;
|
||||
if(depthMax > 0) {
|
||||
minb=minbound(c);
|
||||
maxb=maxbound(c);
|
||||
center=(minb+maxb)/2;
|
||||
}
|
||||
real dbound(pair M, pair dir)
|
||||
{
|
||||
return abs(M-minb-realmult(rectify(dir-center),maxb-minb));
|
||||
}
|
||||
|
||||
pair[] nodes;
|
||||
int n=c.length;
|
||||
|
||||
nodes=polarSort(c,pivot);
|
||||
nodes.cyclic=true;
|
||||
|
||||
bool modified;
|
||||
do {
|
||||
modified=false;
|
||||
for (int i=0; i < n; ++i) {
|
||||
pair p1=nodes[i], p2=nodes[i+1], p3=nodes[i+2];
|
||||
if((p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) < 0)
|
||||
if((depthMax <= 0 || dbound(p2,0.5*(p1+p3)) > depthMax) ||
|
||||
(depthMin == infinity || dbound(p2,0.5*(p1+p3)) < depthMin) ||
|
||||
(angleMin >=360 || (degrees(p3-p2)-degrees(p1-p2) < angleMin)) ||
|
||||
(angleMax <=0 || (degrees(p3-p2)-degrees(p1-p2) > angleMax))) {
|
||||
nodes.delete(i+1);
|
||||
modified=true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
} while(modified);
|
||||
return nodes;
|
||||
}
|
||||
313
pi_numeric.asy
Normal file
313
pi_numeric.asy
Normal file
@ -0,0 +1,313 @@
|
||||
pair sqrt(explicit pair z)
|
||||
{
|
||||
return sqrt(abs(z))*expi(angle(z,false)/2);
|
||||
}
|
||||
|
||||
struct laguer {
|
||||
pair x;
|
||||
int its;
|
||||
}
|
||||
|
||||
real EPS=10*realEpsilon;
|
||||
laguer laguer(pair[] a, int m=a.length-1, pair x)
|
||||
{/*
|
||||
Given the degree 'm' and 'm+1' complex coefficients a[0...m] of polynomial sum(a[i]*x^i),
|
||||
and given a complex value "x", this routine improves 'x' by Laguerre's method until it converges to a root of the given polynomial. The value of x found and the number of iterations taken is returned as 'laguer' structure.
|
||||
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
||||
*/
|
||||
|
||||
static real MR=8, MT=10, MAXIT=MT*MR;
|
||||
int its;
|
||||
laguer ol;
|
||||
real[] ret;
|
||||
int iter,j;
|
||||
real abx,abp,abm,err;
|
||||
pair dx,x1,b,d,f,g,h,sq,gp,gm,g2,ox=x;
|
||||
static real[] frac={0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1};
|
||||
|
||||
laguer ret(){
|
||||
ol.x=x;
|
||||
ol.its=its;
|
||||
return ol;
|
||||
}
|
||||
for (iter=1; iter <= MAXIT ; ++iter) {
|
||||
its=iter;
|
||||
b=a[m];
|
||||
err=abs(b);
|
||||
d=(0,0); f=d;
|
||||
abx=abs(x);
|
||||
for (int j=m-1; j >=0 ; --j) {
|
||||
f=x*f+d;
|
||||
d=x*d+b;
|
||||
b=x*b+a[j];
|
||||
err=abs(b)+abx*err;
|
||||
}
|
||||
err *= EPS;
|
||||
if(abs(b) <= err) return ret();
|
||||
g=d/b;
|
||||
g2=g^2;
|
||||
h=g2-2*f/b;
|
||||
sq=sqrt((m-1)*(m*h-g2));
|
||||
gp=g+sq;
|
||||
gm=g-sq;
|
||||
abp=abs(gp);
|
||||
abm=abs(gm);
|
||||
if(abp < abm) gp=gm;
|
||||
dx=(max(abp,abm) > 0.0) ? (m,0)/gp :
|
||||
(1+abx)*(cos(iter),sin(iter));
|
||||
x1=x-dx;
|
||||
if(x.x == x1.x && x.y == x1.y) return ret();
|
||||
if(iter%MT == 0) x=x1;
|
||||
else x=x-frac[floor(iter/MT)]*dx;
|
||||
}
|
||||
abort("Too many iteration in laguer.");
|
||||
return ret();
|
||||
}
|
||||
|
||||
pair[] zroots(pair[] a, bool polish=true)
|
||||
{/*
|
||||
Given the complex coefficients a[0...m] of the polynomial sum(a[i]*x^i),
|
||||
this routine returns all the complex roots by Laguer's method.
|
||||
If 'polish' is true, the roots are polished (also by Laguer's method');
|
||||
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
||||
*/
|
||||
laguer L;
|
||||
int m=a.length-1,its=L.its;
|
||||
pair[] ad=copy(a);
|
||||
pair x=L.x,b,c;
|
||||
pair[] roots;
|
||||
for (int j=m; j >= 1 ; --j) {
|
||||
x=(0,0);
|
||||
L=laguer(ad,j,L.x);
|
||||
x=L.x; its=L.its;
|
||||
if(abs(x.y) <= 2*EPS*abs(x.x)) x=(x.x,0);
|
||||
roots[j-1]=x;
|
||||
b=ad[j];
|
||||
for (int k=j-1; k >= 0 ; --k) {
|
||||
c=ad[k];
|
||||
ad[k]=b;
|
||||
b=x*b+c;
|
||||
}
|
||||
}
|
||||
|
||||
if(polish) {
|
||||
for (int j=0; j < m ; ++j) {
|
||||
L.x=roots[j];
|
||||
L.its=its;
|
||||
L=laguer(a,L.x);
|
||||
roots[j]=L.x;
|
||||
its=L.its;
|
||||
}
|
||||
}
|
||||
for (int j=0; j < m ; ++j) {
|
||||
pair x=roots[j];
|
||||
int i;
|
||||
for (i=j-1; i >= 0 ; --i) {
|
||||
if(roots[i].y == 0) break;
|
||||
roots[i+1]=roots[i];
|
||||
}
|
||||
roots[i+1]=x;
|
||||
}
|
||||
return roots;
|
||||
}
|
||||
|
||||
real sgn(real a, real b)
|
||||
{
|
||||
return (b == 0) ? abs(a): sgn(b)*abs(a);
|
||||
}
|
||||
|
||||
void hqr(real[][] a, int n, real[] wr=new real[]{}, real[] wi= new real[]{})
|
||||
{/*
|
||||
Finds all eigenvalues of an upper Hessenberg matrix a[1...n][1..n].
|
||||
The real and imaginary parts of eigenvalues are returned in wr[1...n] and wi[1...n], respectively.
|
||||
Adapted from http://www.nrbook.com/a/bookcpdf/c11-6.pdf
|
||||
*/
|
||||
int nn,m,l,k,j,its,i,mmin;
|
||||
real z,y,x,w,v,u,t,s,r,q,p,anorm=0;
|
||||
for (i=1; i <= n ; ++i)
|
||||
for (j=max(i-1,1); j <= n ; ++j) anorm += abs(a[i][j]);
|
||||
nn=n;t=0;
|
||||
while(nn >= 1) {
|
||||
its=0;
|
||||
do {
|
||||
for (l=nn; l >= 2 ; --l) {
|
||||
s=abs(a[l-1][l-1])+abs(a[l][l]);
|
||||
if(s == 0) s=anorm;
|
||||
if(abs(a[l][l-1])+s == s) {
|
||||
a[l][l-1]=0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
x=a[nn][nn];
|
||||
if(l == nn) {
|
||||
wr[nn]=x+t;
|
||||
wi[--nn]=0;
|
||||
} else {
|
||||
y=a[nn-1][nn-1];
|
||||
w=a[nn][nn-1]*a[nn-1][nn];
|
||||
if(l == (nn-1)) {
|
||||
p=0.5*(y-x);
|
||||
q=p^2+w;
|
||||
z=sqrt(abs(q));
|
||||
x += t;
|
||||
if(q >= 0) {
|
||||
z=p+sgn(z,p);
|
||||
wr[nn-1]=wr[nn]=x+z;
|
||||
if(z != 0) wr[nn]=x-w/z;
|
||||
wi[nn-1]=wi[nn]=0;
|
||||
} else {
|
||||
wr[nn-1]=wr[nn]=x+p;
|
||||
wi[nn-1]=-(wi[nn]=z);
|
||||
}
|
||||
nn -= 2;
|
||||
} else {
|
||||
if(its == 30) abort("Too many iterations in hqr...");
|
||||
if(its == 10 || its == 20) {
|
||||
t += x;
|
||||
for(i=1; i <= nn; ++i) a[i][i] -= x;
|
||||
s=abs(a[nn][nn-1])+abs(a[nn-1][nn-2]);
|
||||
y=x=0.75*s;
|
||||
w= -0.4375*s^2;
|
||||
}
|
||||
++its;
|
||||
for (m=nn-2; m >= l; --m ) {
|
||||
z=a[m][m];
|
||||
r=x-z;
|
||||
s=y-z;
|
||||
p=(r*s-w)/a[m+1][m]+a[m][m+1];
|
||||
q=a[m+1][m+1]-z-r-s;
|
||||
r=a[m+2][m+1];
|
||||
s=abs(p)+abs(q)+abs(r);
|
||||
p /= s; q /= s; r /= s;
|
||||
if(m == 1) break;
|
||||
u=abs(a[m][m-1])*(abs(q)+abs(r));
|
||||
v=abs(p)*(abs(a[m-1][m-1])+abs(z)+abs(a[m+1][m+1]));
|
||||
if(u+v == v) break;
|
||||
}
|
||||
for (i=m+2; i <= nn; ++i) {
|
||||
a[i][i-2]=0;
|
||||
if(i != m+2) a[i][i-3]=0;
|
||||
}
|
||||
for (k=m; k <= nn-1; ++k) {
|
||||
if(k != m) {
|
||||
p=a[k][k-1];
|
||||
q=a[k+1][k-1];
|
||||
r=0;
|
||||
if(k != nn-1) r=a[k+2][k-1];
|
||||
if((x=abs(p)+abs(q)+abs(r)) != 0) {
|
||||
p /= x;
|
||||
q /= x;
|
||||
r /=x;
|
||||
}
|
||||
}
|
||||
if((s=sgn(sqrt(p^2+q^2+r^2),p)) != 0) {
|
||||
if(k == m) {
|
||||
if(l != m) a[k][k-1]=-a[k][k-1];
|
||||
} else a[k][k-1]=-s*x;
|
||||
p += s;
|
||||
x=p/s; y=q/s; z=r/s; q /= p; r /= p;
|
||||
for (j=k; j <= nn; ++j) {
|
||||
p=a[k][j]+q*a[k+1][j];
|
||||
if(k != nn-1) {
|
||||
p += r*a[k+2][j];
|
||||
a[k+2][j] -= p*z;
|
||||
}
|
||||
a[k+1][j] -= p*y;
|
||||
a[k][j] -= p*x;
|
||||
}
|
||||
mmin=nn < k+3 ? nn : k+3;
|
||||
for (i=l; i <= mmin; ++i) {
|
||||
p=x*a[i][k]+y*a[i][k+1];
|
||||
if(k != nn-1) {
|
||||
p += z*a[i][k+2];
|
||||
a[i][k+2] -= p*r;
|
||||
}
|
||||
a[i][k+1] -= p*q;
|
||||
a[i][k] -= p;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} while(1 < nn-1);
|
||||
}
|
||||
}
|
||||
|
||||
real RADIX=8;
|
||||
void balanc(real[][] a, int n=a[0].length)
|
||||
{/* Given a matrix 'a[1...n][1..n]', this routine replaces it by a balanced matrix with identical eigenvalues.
|
||||
A symmetric matrix is already balanced and is unaffected by this procedure.
|
||||
The variable RADIX should be the machine's floating-point radix.
|
||||
Adapted from http://www.nrbook.com/a/bookcpdf/c11-5.pdf
|
||||
*/
|
||||
int last=0, j, i;
|
||||
real s, r, g, f, c, sqrdx=RADIX^2;
|
||||
while(last == 0) {
|
||||
last=1;
|
||||
for (i=1; i <= n; ++i) {
|
||||
r=c=0;
|
||||
for (j=1; j <= n; ++j)
|
||||
if(j != i) {
|
||||
c += abs(a[j][i]);
|
||||
r += abs(a[i][j]);
|
||||
}
|
||||
if(c != 0 && r != 0) {
|
||||
g=r/RADIX;
|
||||
f=1;
|
||||
s=c+r;
|
||||
while(c < g) {
|
||||
f *= RADIX;
|
||||
c *= sqrdx;
|
||||
}
|
||||
g=r*RADIX;
|
||||
while(c > g) {
|
||||
f /= RADIX;
|
||||
c /= sqrdx;
|
||||
}
|
||||
if((c+r)/f < 0.95*s) {
|
||||
last=0;
|
||||
g=1/f;
|
||||
for (j=1; j <= n; ++j) a[i][j] *= g;
|
||||
for (j=1; j <= n; ++j) a[j][i] *= j;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pair[] zroots(real[] a, bool polish=true)
|
||||
{/* Return all the roots of a polynomial with real coefficients, sum(a[i]x^i).
|
||||
The method is to construct an upper Hessenberg matrix whose eigenvalues are
|
||||
the desired roots, and the use the routine 'balanc' and 'hqr'.
|
||||
If polish is true, root-polishing by Newton-Raphson's method is applied.
|
||||
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
||||
*/
|
||||
int m=a.length-1;
|
||||
int j,k;
|
||||
real[][] hess=new real[m+1][m+1];
|
||||
real[] rtr=sequence(new real(int n){return 0;},m+1),
|
||||
rti=sequence(new real(int n){return 0;},m+1);
|
||||
real xr,xi;
|
||||
if(a[m] == 0) abort("Bad args in zroots.");
|
||||
for (k=1; k <= m; ++k) {
|
||||
hess[1][k]=-a[m-k]/a[m];
|
||||
for(j=2; j <= m; ++j) hess[j][k]=0;
|
||||
if(k != m) hess[k+1][k]=1;
|
||||
}
|
||||
balanc(hess,m);
|
||||
hqr(hess,m,rtr,rti);
|
||||
for (j=2; j <= m; ++j) {
|
||||
xr=rtr[j];
|
||||
xi=rti[j];
|
||||
for (k=j-1; k >= 1; --k) {
|
||||
if(rtr[k] <= xr) break;
|
||||
rtr[k+1]=rtr[k];
|
||||
rti[k+1]=rti[k];
|
||||
}
|
||||
rtr[k+1]=xr;
|
||||
rti[k+1]=xi;
|
||||
}
|
||||
pair[] roots=sequence(new pair(int n){return (rtr[n],rti[n]);},m+1);
|
||||
roots.delete(0);
|
||||
return roots;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user