/*
Libreria per le radici dei polinomi
Molto commentato: http://www.elegio.it/calcolatrice/quarto-grado-1.html
Questo: http://www.elegio.it/calcolatrice/quarto_grado_esterno.js.html
Applicato in : http://www.elegio.it/calcolatrice/quarto_grado_esterno-.html
Anche: http://www.elegio.it/calcolatrice/quarto_grado_solo_esterno-1.html
e anche: http://www.elegio.it/omnia/js/meglioradici-interno.html
*/
//
var polsol={versione:"20170221",tartaglia:1};
//
polsol.polinomio=function(p,z){
//
// z rappresenta la variabile indipendente e
// puo' essere sia un numero che una Array,
// z[0] rappresenta la parte reale del numero complesso
// e z[1], se c'e' rappresenta la parte immaginaria.
//
// Il risultato e' sempre una Array di due elementi,
// con l'elemento 0 parte reale e l'elemento 1 parte
// immaginaria.
//
// Qui il polinomio ha come termine di grado 0 il
// primo elemento, e l'ultimo elemento e' quello
// di grado piu' elevato.
//
var x,y,rx,ry,rn,j;
if(Array.isArray(z)){x=z[0];
if(z.length>1)y=z[1];
else{y=0}}
else{x=z;y=0;}
rx=1;ry=0;
for(j=p.length-1;j>-1;j--){
if(Array.isArray(p[j])){rn=rx*x-ry*y+p[j][0];
if(p[j].length>1){ry=rx*y+ry*x+p[j][1];}
else{ry=rx*y+ry*x;}}
else{rn=rx*x-ry*y+p[j];ry=rx*y+ry*x;}
rx=rn;
}
//
// Un numero complesso e' indicato con [a,b] dove
// a e b sono numeri reali e dunque e' una Array di
// due elementi: l'elemento 0 e' la parte reale e
// l'elemento 1 e' la parte immaginaria.
//
return [rx,ry];
};
//
// Crea un polinomio di quarto grado di cui sa anche
// le sue quattro radici e fa in modo che il termine
// noto non possa mai essere esattamente nullo,
// anche se piccolissimo.
// Dati quattro argomenti reali calcola il polinomio
// di quarto grado con 1 sottinteso come fattore
// della quarta potenza.
//
polsol.perverifiche=function(a,bb,c,dd){
var b,d,pp,r1,r2,r3,r4,det;
if(Math.abs(bb)==0)b=1.e-5;
else b=bb;
if(Math.abs(dd)==0)d=1.e-5;
else d=dd;
pp=[[b*d,0],[2*(a*d+b*c),0],
[b+d+4*a*c,0],[2*(a+c),0]];
det=a*a-b;
if(0>det){r1=[-a,Math.sqrt(-det)];r2=[-a,-Math.sqrt(-det)]}
else{ r1=[-a+Math.sqrt(det),0];r2=[-a-Math.sqrt(det),0]}
det=c*c-d;
if(0>det){r3=[-c,Math.sqrt(-det)];r4=[-c,-Math.sqrt(-det)]}
else{r3=[-c+Math.sqrt(det),0];r4=[-c-Math.sqrt(det),0]}
return [pp,[r1,r2,r3,r4]]; }
//
// Equazione di secondo grado
//
polsol.secondogrado=function(r){
//
// Radici equaz. di secondo grado a coeff. complessi:
// (y + r[1] )*y + r[0] = 0
//
// r[0][0] rappresenta la parte reale del termine noto.
// r[0][1] e' la parte immaginaria del termine noto.
// Se r[0] non e' una Array, r[0] e' il termine noto reale.
// r[1][0] e' la parte reale del termine di primo grado.
// r[1][1] e' la parte immaginaria del termine di primo grado.
// Se r[1] non e' una Array r[1] e' il termine di primo grado reale.
// Il termine di secondo grado e' sottinteso e vale 1.
//
var p,q,s,t,x,a,b;
if(Array.isArray(r[0]))b=r[0];
else{b=[r[0],0]};
if(Array.isArray(r[1]))a=r[1];
else{a=[r[1],0]};
//
if(0>=b[0]*b[0]+b[1]*b[1]) return [ [0,0],[-a[0],-a[1]] ];
p=a[0]*a[0]-a[1]*a[1]-4*b[0];
q=2*a[0]*a[1]-4*b[1];
with(Math) {
if(p>0){
s=sqrt((sqrt(p*p+q*q)+p)*0.5);
t=0.5*q/s;
}
else {
t=sqrt((sqrt(p*p+q*q)-p)*0.5);
// radici coincidenti
if(0 >= t) return [[-a[0]/2,-a[1]/2],[-a[0]/2,-a[1]/2]];
s=0.5*q/t;
}
// Fa direttamente la radice di maggior modulo
if(a[0]*s+a[1]*t>0) x=[-0.5*(a[0]+s), -0.5*(a[1]+t)];
else x=[-0.5*(a[0]-s), -0.5*(a[1]-t)];
// Deduce l'altra dal termine noto.
p=1/( x[0]*x[0] + x[1]*x[1] );
}
return [ x, [p*(b[0]*x[0]+b[1]*x[1]), p*(-b[0]*x[1]+b[1]*x[0])] ];
};
//
// Equazione cubica.
//
polsol.cubica=function(poli){
//
// radici di polinomio di terzo grado con
// coefficienti reali.
// Se il polinomio ha coefficienti con parte immaginaria la
// parte immaginaria viene ignorata.
//
var p,q,r,a,b,d,aa,bb,rd,nr,psi;
var an,sn,qn,pn,segno;
if(3>poli.length) return [ [0,1],[0,1],[0,1] ];
if(Array.isArray(poli[2]))p=poli[2][0];
else{p=poli[2]};
if(Array.isArray(poli[1]))q=poli[1][0];
else{q=poli[1]};
if(Array.isArray(poli[0]))r=poli[0][0];
else{r=poli[0]};
a=(p*p-3*q)/9;
b=(p*(9*q-2*p*p)-27*r)/54;
d=b*b-a*a*a;
//
sn=p/3;
qn=r-sn*(q-sn*(p-sn));
pn=q-sn*(2*p-3*sn);
if(0>=pn*pn){
// cubo perfetto, tre radici coincidenti !
if(qn>0) segno=1;
else if(0>qn) segno=-1;
else segno=0;
an=segno*Math.pow(Math.abs(qn),1/3)-sn;
return [[an,0],[an,0],[an,0]];};
with(Math){
if(0>=d) { // Tre radici reali
rd=sqrt(-d);
nr=sqrt(-d+b*b);
if(abs(b)>=rd) {
psi=asin(rd/nr);
if(0>b) psi=PI-psi;
}
else psi=acos(b/nr);
aa=pow(nr,1/3);
bb=aa*sin(psi/3);
aa=aa*cos(psi/3);
return [ [ 2*aa-p/3,0],
[ bb*sqrt(3)-aa-p/3,0],
[-bb*sqrt(3)-aa-p/3,0] ];
}
else { // Una sola radice reale
if(b>0) {
aa = pow(sqrt(d)+b,1/3);
bb = a/aa;
}
else {
bb = -pow(sqrt(d)-b,1/3);
aa = a/bb;
}
return [ [aa+bb-p/3,0],
[-0.5*(aa+bb)-p/3,-0.5*sqrt(3)*(aa-bb)],
[-0.5*(aa+bb)-p/3, 0.5*sqrt(3)*(aa-bb)] ];
};
};
};
//
// Se la variabile tartaglia vale -1...
//
polsol.cubicabis=function(poli){
//
// radici di polinomio di terzo grado con
// coefficienti reali e calcolati con
// algoritmo basato su funzioni trigonometriche
// ed esponenziali
var pp0,pp1,pp2;
var a,p,q,s,am,bm,cm, dm, ca,sina,cosa,expa,segno,sqrt3;
sqrt3=Math.sqrt(3);
if(3>poli.length) return[ [0,1],[0,1],[0,1] ];
// Lo faccio uscire se il polinomio e' mal fatto.
if(Array.isArray(poli[2]))pp2=poli[2][0];
else{pp2=poli[2]};
if(Array.isArray(poli[1]))pp1=poli[1][0];
else{pp1=poli[1]};
if(Array.isArray(poli[0]))pp0=poli[0][0];
else{pp0=poli[0]};
s=pp2/3;
q=pp0-s*(pp1-s*(pp2-s));
p=pp1-s*(2*pp2-3*s);
if(0>=p*p){
// cubo perfetto, tre radici coincidenti !
if(q>0) segno=1;
else if(0>q) segno=-1;
else segno=0;
a=segno*Math.pow(Math.abs(q),1/3)-s;
return [[a,0],[a,0],[a,0]];};
if(p>0){
bm=Math.sqrt(4*p/3);
dm=-q*Math.sqrt(27/(4*p))/p;
a=Math.log(dm+Math.sqrt(dm*dm+1))/3;
expa=Math.exp(a);
sina=bm*(expa-1/expa)/2;
cosa=bm*(expa+1/expa)/2;
return [[sina-s,0],
[-sina/2-s,sqrt3*cosa/2],
[-sina/2-s,-sqrt3*cosa/2]];};
//
am=Math.sqrt(-4*p/3);
cm=-q*Math.sqrt(-27/(4*p))/p;
ca=Math.abs(cm);
if(ca>1){
if(cm>0) segno=1;
else if(0>cm) segno=-1;
else segno=0;
a=Math.log(ca+Math.sqrt(ca*ca-1))/3;
expa=Math.exp(a);
sina=segno*am*(expa-1/expa)/2;
cosa=segno*am*(expa+1/expa)/2;
return [[-cosa-s,0],
[cosa/2-s,-sqrt3*sina/2],
[cosa/2-s,sqrt3*sina/2]];};
// tre radici reali
a=Math.asin(cm)/3;
sina=am*Math.sin(a);
cosa=am*Math.cos(a);
return [[sina-s,0],
[-sina/2+sqrt3*cosa/2-s,0],
[-sina/2-sqrt3*cosa/2-s,0]];
};
//
polsol.cubicarobusta=function(poli){
// radici di polinomio di terzo grado con
// coefficienti reali. Calcolo preciso anche con
// termine noto molto piccolo.
var p,q,r,a,b,d,aa,bb,x1,x2,x3,rd,nr,psi;
var an,sn,qn,pn,segno;
if(3>poli.length) return [ [0,1],[0,1],[0,1] ];
if(Array.isArray(poli[2]))p=poli[2][0];
else{p=poli[2]};
if(Array.isArray(poli[1]))q=poli[1][0];
else{q=poli[1]};
if(Array.isArray(poli[0]))r=poli[0][0];
else{r=poli[0]};
// alert("In cubicarobusta"+" p= "+p+" q= "+q+" r= "+r);
a=(p*p-3*q)/9;
b=(p*(9*q-2*p*p)-27*r)/54;
d=b*b-a*a*a;
//
sn=p/3;
qn=r-sn*(q-sn*(p-sn));
pn=q-sn*(2*p-3*sn);
if(0>=pn*pn){
// cubo perfetto, tre radici coincidenti !
if(qn>0) segno=1;
else if(0>qn) segno=-1;
else segno=0;
an=segno*Math.pow(Math.abs(qn),1/3)-sn;
return [[an,0],[an,0],[an,0]];};
with(Math){
if(0>d) { // Tre radici reali
rd=sqrt(-d);
nr=sqrt(-d+b*b);
if(abs(b)>rd) {
psi=asin(rd/nr);
if(0>b) psi=PI-psi;
}
else psi=acos(b/nr);
aa=pow(nr,1/3);
bb=aa*sin(psi/3);
aa=aa*cos(psi/3);
x1=2*aa-p/3;
x2=bb*sqrt(3)-aa-p/3;
x3=-bb*sqrt(3)-aa-p/3;
//
if(!(!abs(x2)>abs(x1)
|| !abs(x3)>abs(x1))){
if(x2*x3!=0)x1=-r/(x2*x3);}
//
else if(!(!abs(x3)>abs(x2)
|| !abs(x1)>abs(x2))){
if(x1*x3!=0)x2=-r/(x1*x3);}
//
else if(!(!abs(x1)>abs(x3)
|| !abs(x2)>abs(x3))){
if(x1*x2!=0)x3=-r/(x1*x2);}
// alert("Tre radici reali:\n"+x1+"\n"+x2+"\n"+x3);
return [ [x1,0], [x2,0], [x3,0] ]
}
else { // Una sola radice reale
if(b>0) {
aa = pow(sqrt(d)+b,1/3);
bb = a/aa;
}
else {
bb = -pow(sqrt(d)-b,1/3);
aa = a/bb;
}
x1=aa+bb-p/3;
x2=-0.5*(aa+bb)-p/3;
x3=-0.5*sqrt(3)*(aa-bb);
if(x2*x2+x3*x3>x1*x1) x1=-r/(x2*x2+x3*x3);
// alert("La radice reale e' una sola: "+x1);
return [ [x1, 0], [x2, x3], [x2,-x3] ];
};
};
};
//
// Equazione di quarto grado.
//
polsol.tetra=function(p){
// radici di polinomio di quarto grado con
// coefficienti reali.
var a,b,c,d,q,u,s,t,v,w,e,f,g,h,m,n;
if(Array.isArray(p[3]))a=p[3][0];
else{a=p[3]};
if(Array.isArray(p[2]))b=p[2][0];
else{b=p[2]};
if(Array.isArray(p[1]))c=p[1][0];
else{c=p[1]};
if(Array.isArray(p[0]))d=p[0][0];
else{d=p[0]};
q=[[4*d*b-c*c-d*a*a,0],[c*a-4*d,0],[-b,0]];
if(polsol.tartaglia >0)u=polsol.cubica(q);
else u=polsol.cubicabis(q);
s=a*a/4+u[0][0]-b;
t=u[0][0]*u[0][0]/4-d;
with(Math){
v=abs(s); w=abs(t);
if(w>v) {
w=sqrt(w);
v=(a*u[0][0]-2*c)/(4*w);
}
else if(v>w) {
t=s;
v=sqrt(v);
w=(a*u[0][0]-2*c)/(4*v);
}
}
if(t >= 0) {
e=[u[0][0]/2-w,0];
f=[u[0][0]/2+w,0];
g=[a/2-v,0];
h=[a/2+v,0];
}
else {
e=[u[0][0]/2,w];
f=[u[0][0]/2,-w];
g=[a/2,-v];
h=[a/2,v];
}
m=polsol.secondogrado([e,g]);
n=polsol.secondogrado([f,h]);
return [n[0],n[1],m[0],m[1]];
};
//
polsol.tetrarobusta=function(p){
//
// radici di polinomio di quarto grado con
// coefficienti reali. Tratta bene il caso
// del termine noto quasi nullo.
//
var a,b,c,d,q,r,s,t,u,v,w,e,f,g,h,m,n;
if(Array.isArray(p[3]))a=p[3][0];
else{a=p[3]};
if(Array.isArray(p[2]))b=p[2][0];
else{b=p[2]};
if(Array.isArray(p[1]))c=p[1][0];
else{c=p[1]};
if(Array.isArray(p[0]))d=p[0][0];
else{d=p[0]};
//
// Correzione empirica per evitare situazioni critiche.
// Se d valesse esattamente 0 l'equazione sarebbe di
// terzo grado e una radice nulla e se d valesse
// esattamente 1 potrebbe essere una quarta potenza
// della equazione z=1 cioe' quattro radici coincidenti
// e problema simile se d valesse esattamente -1.
// COMUNQUE IN QUESTI CASI SI OTTENGONO PRECISIONI PESSIME...
// LE SOLUZIONI ANDREBBERO RAFFINATE ITERATIVAMENTE...
//
if(1.e-16>Math.abs(d))d=1.e-16;
if(5.e-16>Math.abs(d-1))d=1+5.e-16;
if(5.e-16>Math.abs(d+1))d=-1+5.e-16;
//
q=[[4*d*b-c*c-d*a*a,0],[c*a-4*d,0],[-b,0]];
if(polsol.tartaglia >0)u=polsol.cubicarobusta(q);
else u=polsol.cubicabis(q);
// alert("Trovo prima radice "+u[0]);
r=u[0][0];
with(Math){
if(!(!(0>=abs(u[1][1]))
|| !(abs(u[1][0]) > abs(r)) )) r=u[1][0];
if(!(!(0>=abs(u[2][1]))
|| !(abs(u[2][0]) > abs(r)) )) r=u[2][0];
//
s=a*a/4+r-b;
t=r*r/4-d;
v=abs(s); w=abs(t);
if(w>v) {
w=sqrt(w);
v=(a*r-2*c)/(4*w);
}
else if(v>w) {
v=sqrt(v);
w=(a*r-2*c)/(4*v);
}
if( abs(r/2-w) > abs(r/2+w) ) {
e=[r/2-w,0];
f=[d/e[0],0];
}
else if( abs(r/2+w) > abs(r/2-w) ) {
f=[r/2+w,0];
e=[d/f[0],0];
}
else {
e=[r/2-w,0];
f=[r/2+w,0];
}
g=[a/2-v,0];
h=[a/2+v,0];
}
m=polsol.secondogrado([e,g]);
n=polsol.secondogrado([f,h]);
return [n[0],n[1],m[0],m[1]];
};
//
// Se si usa questa funzione fa prove a caso delle varie funzioni della libreria polsol
//
polsol.provacaso=function(){
var j,n,z,xbr,poli=[],zeri=[];
var casi=[[2,polsol.secondogrado],
[3,polsol.cubica],[3,polsol.cubicabis],
[3,polsol.cubicarobusta],[4,polsol.tetra],
[4,polsol.tetrarobusta]];
xbr="\u003cbr/\u003e";
n=Math.round(100*Math.random())%casi.length;
for(j=0;casi[n][0]>j;j++)poli[j]=
[Math.round(2000*Math.random()-1000),0];
if(n==0)poli[0][1]=Math.round(2000*Math.random()-1000);
polsol.tartaglia=1-2*(Math.round(100*Math.random())%2);
z=casi[n][1](poli);
for(j=0;casi[n][0]>j;j++){
zeri[j]=polsol.polinomio(poli,z[j])
}
return (xbr+"Esempio "+n+xbr+
"Polinomio: "+poli.join(";"+xbr)+
xbr+"Radici: "+z.join("; "+xbr)+
xbr+"Primo Zero: "+
zeri.join(xbr+"Altro Zero: ") );
}
//
//