Scientific Calculator

gio/Polinomi

By gio

// 07/03/2020
// i coefficienti dei polinomi sono ordinati dal termine noto alla potenza più alta come una lista
// in tutte le funzioni di polinomi.

//2 funzioni di prova, la seconda è ricorsiva
realp=function(x){real(x)};
divn=function(x,n){if (x mod n ==0)-> divn(x/n,n) else x};

//dati i polimomi "ff"==P(x) e "gg"== (x=t-a) trova i polinomio composto Y(t)=P((t-a))
polytrasl= function(ff,gg){var (k,j,r,s)=[0,0,[],0];
	while (k< len(ff)) ->(j=k;
		while (j<len(ff))->(s += nCr(j,k)*ff[j]*(gg[0])^(j-k);j+=1);
			r=concat(r,[s*(gg[1])^k]);s=0;k+=1;);r};

// trova i coefficienti del polinomio prodotto i cui fattori (numeri o polinomi)
// sono gli elementi di tipo lista della lista "lst"
polymul=function(lst){var (mul,ls,ww)=[[1],clone(lst),0];
	var MU=function(a,b){var (k,h,I,J,r,s)=[0,0,len(a)-1,len(b)-1,[],0];
		while (k<=I+J)->(h=max(0,k-J);
			while (h<=min(k,I))->(s+=a[h]*b[k-h];h+=1;);
			r=concat(r,[s]);s=0;k+=1);r};
    while (len(ls)>0) ->(ww=shift(ls);
        if (instring(system:type(ww),"integer-float-complex-rational"))[0]->(mul=mul*ww;)
        else (if (instring(system:type(ww),"list-matrix")[0])-> mul=MU(mul,ww))
             );mul};

// costruisce i coefficienti della derivata di un polinomio di "lst"
polyder=function(lst){var (k,d,ls)=[0,[],clone(lst)];d=shift(lst);
	d=map(x->(k=k+1;k*x),lst)};
// calcola il polinomio integrale indefinito del polinomio "lst".
// il termine noto è posto a zero
polyinteg=function(lst){var (k,ls)=[0,clone(lst)];
	concat([0],map(x->(k+=1;x/k),ls))};
// da il valore di un polinomio di coefficienti "cc" nel valore "x"
// c'è la funzione destra e quella sinistra in quanto funziona anche con polinomi
// di matrici
polydx=function(x,cc) {var (k,y,c)=[0,0,clone(cc)];
	while c != []->(y +=shift(c)*x^k;k+=1);
	y};
polysx=function(x,cc) {var (k,y,c)=[0,0,clone(cc)];
	while c != []->(y +=(x^k)*shift(c);k+=1);
	y};
	
// divide un polinomio (numeri o matrici) per (x-x0).
// in una lista da il resto e i coefficienti del polinomio ottenuto
polydivdx=function(x0,ll){var (l,q)=[clone(ll),[]];
	while(l != []) ->(q=concat(q,[polydx(x0,l)]);shift(l));
	[shift(q),q]};
polydivsx=function(x0,ll){var (l,q)=[clone(ll),[]];
	while(l != []) ->(q=concat(q,[polysx(x0,l)]);shift(l));
	[shift(q),q]};
	
// minimi quadrati
// lxx = lista valori x
// lyy = lista valori y
// gr = grado del polinomio interpolatore
// return=[lista coefficienti del polinomio,delta=deviazione standart degli scarti/deviazione standart delle y]
// se gr = -1 l' interpolazione viene eseguita sul logaritmo della funzione
// y=A*e^(alfa*x) [ln(y)=ln(A)+alfa*x]
// return [[A,alfa],delta=deviazione standart degli scarti/deviazione standart delle y] 
// se gr = -2 l' interpolazione viene eseguita sul logaritmo della funzione
// y=A*x^n*e^(alfa*x) [ln(y) =ln(A)+n*ln(x)+alfa*x]
// return [[A,n,alfa],delta=deviazione standart degli scarti/deviazione standart delle y] 
MINQUAD=function(lxx,lyy,gr){var (lx,ly)=[clone(lxx),clone(lyy)];
    var (dlx,dly)=[dim(lx)[0],dim(ly)[0]];var g=clone(gr);
    if (dlx<=g)-> g=dlx-1;
    media=function(lst,f){sum(lst,f)/dim(lst)[0]};
	POISSON=function(){
        var M=[ones(1,dim(lx)[0])[0],map(x->ln(x),lx),lx];
        var (A,b,c)= [M*M',M*map(x->ln(x),ly)',[ly,lx]'];
        b=((1/A)*b)'[0];b[0]= e^b[0];
        var sc=map(x->(x[0]-b[0]*(x[1]^b[1])*e^(b[2]*x[1])),[lyy,lxx]');
		var delta=sqrt(media(sc,x->x^2)/(media(ly,x->x^2)-media(ly)^2));
    [b,delta]};
    EXPO=function(){ly=map(x->ln(x),ly);g=1;var b=POLY()[0];b[0]=e^b[0];
        var sc=map(x->(x[0]-b[0]*e^(b[1]*x[1])),[lyy,lxx]');
		var delta=sqrt(media(sc,x->x^2)/(media(lyy,x->x^2)-media(lyy)^2));
        [b,delta]};
    POLY=function(){var (A,b,c)=[zeros(g+1,g+1),[],[ly,lx]'];var (k,j)=[0,0];
		while (k<=g)-> (
        j=k; while (j<=g) -> (
             A[k][j]=A[j][k]=sum(lx,x->x^(k+j));j+=1);
		b=push(b,sum(c,x->x[0]*x[1]^k));k+=1);
		b=((1/A)*b')'[0];
		var sc=map(x->(x[0]-polydx(x[1],b)),c);
		var delta=sqrt(media(sc,x->x^2)/(media(ly,x->x^2)-media(ly)^2));
    [b,delta]};
	CASE([dlx != dly,x->error("dati non congruenti");g== -1,EXPO;g== -2,POISSON;true,POLY])
    };

// soluzione di una equazione di secondo grado 
eq2=function(lst){
	var (c,b,a)=lst;[(-b-sqrt(b*b-4*a*c))/(2*a),(-b+sqrt(b*b-4*a*c))/(2*a)]};

// soluzione di una equazione di terzo grado col metodo di Cardano
eq3=function(lst,err){if dim(lst)[0] != 4 -> error("i coeff. devono essere 4.Da eq3");
	var er=(err==null ? 10 : err);
	var (f_i,ls,s,ww)=[true,clone(lst),[],e^(2/3*pi*1i)];	
	if ifzero(ls[3],er) ->(f_i=false;pop(ls);s=eq2(ls)) 
		else (ls=ls/ls[3];var b3 = -ls[2]/3;pt=polytrasl(ls,[b3,1]);var(q,p,o,a)=pt);		
	if (f_i and ifzero(p,er) and ifzero(q,er))->(f_i=false;s=[b3,b3,b3]);	
	if (f_i and ifzero(4*p^3+27*q^2,er))->(f_i=false;s=b3+map(x->3*q*x/(2*p),[2,-1,-1]));	
	if (f_i and ifzero(p,er))->(f_i=false;s=b3-q^(1/3)*[1,ww,1/ww])
		else (var z12=eq2([-(p^3)/27,q,a]););		
	if f_i -> (var z1=z12[0]^(1/3);var z2=-p/(3*z1);
	s=b3+[z1+z2,z1*ww+z2/ww,z1/ww+z2*ww]);
	s};

spam? | offensive?