>

```Warning, the protected names norm and trace have been redefined and unprotected
```

Everything.txt is part of Gibbons code. Lets you use notation of gasket.

#P_energy calculates p-energy on the m-th level of the gasket. fnum is the function #number, it is just an arbitrary number greater than 4. Numbers less than four. #represent the neighbors of vertices.p is obviously the general power one is #working with.

>
p_energy:=proc(T,m,fnum,p)
local vm,ver,i,j,k,tot;
tot:=0;
vm:=vertices(m):
make_neighbors(T,m,vm):
for j from 0 to 2 do
for k from 1 to 2 do
tot:=tot+abs(to_power((T[[j],fnum]-T[T[[j],k],fnum]),p,-1));
end do;
end do;
vm:=subsop(1=NULL,2=NULL,3=NULL,vm);
for ver in vm do
for i from 1 to 4 do
tot:=tot+abs(to_power((T[ver,fnum]-T[T[ver,i],fnum]),p,-1));
end do;
end do;
print(tot);
end:

#This program gives the "crude" initial p-energy and its vertices.How it does this is it sets initial vertices to 0 except for the boundary values, and then unassigns one value and minimizes the energy with respect to it, and this is the vertex's new value.This process is then cycled for all the vertices (non-boundary) k-times. m is again the level of the gasket one is working on. fnum is the same as above.p is the power one is looking at. k is the number of times the optimization is cycled. a,b,c are the boundary values.

> optimize:=proc(T,m,fnum,p,k,a,b,c,pnn,frst,secd,thrd)
local ans11,subT,vm,big,ver,i,j,ans1,ans2;
Digits:=15;
vm:=vertices(m):
make_neighbors(T,m,vm):
vm:=subsop(1=NULL,2=NULL,3=NULL,vm);
for ver in vm do
T[ver,fnum]:=0;
end do;
T[[0],fnum]:=a;
T[[1],fnum]:=b;
T[[2],fnum]:=c;
for i from 1 to k do
for ver in vm do
unassign('T[ver,fnum]');
subT:=T[ver,fnum];
ans1:=proc(subT) ans1(subT):=to_power((subT-T[T[ver,1],fnum]),(p-1),-1)+to_power((subT-T[T[ver,2],fnum]),(p-1),-1)+to_power((subT-T[T[ver,3],fnum]),(p-1),-1)+to_power((subT-T[T[ver,4],fnum]),(p-1),-1):
end;
ans11:=op(bisect_go(ans1,-10..10,.000001));

assign(T[ver,fnum], ans11);

end do;
end do;
#ans2:=p_energy(T,m,fnum,p);
#print(ans2);
for big from 0 to 2 do
unassign('T[[big],fnum]');
end do;
#for ver in vm do
#print(T[ver,fnum]);
#end do;
pnn:=ans2;
frst:=T[[0,1],fnum];
secd:=T[[1,2],fnum];
thrd:=T[[0,2],fnum];
end:

This thing below uses a least squares method and the program above to get an initial crude p-energy polynomial. To use it for the opt2 and polymake, you will probably have to copy and paste it, like NwN_30. This also gives intial p-harmonic extensions, B1 is for T[[0,1],fnum], B2 is for T[[1,2],fnum], and B3 is for T[[0,2],fnum].

> N30:=[]; N30_1:=[]; N30_2:=[]; N30_3:=[];
for i from 0 to 10 do
unassign('pnn,frst,secd,thrd');
optimize(T,1,8,12/5,50,-1,i/10,1,pnn,frst,secd,thrd);
N30:=[op(N30),pnn];
N30_1:=[op(N30_1),frst];
N30_2:=[op(N30_2),secd];
N30_3:=[op(N30_3),thrd];
end do;
times:=10;
f30:=proc(x)
f30(x):=x/N30[11];
end;
N30:=map(f30,N30);
print(N30);
An30:=proc(x)
An30(x):=N30[1]+b*x^2+c*x^4+d*x^6+e*x^8+f*x^10+h*x^12+m*x^14+n*x^16;
end:

> with(linalg):
tot:=0;
unassign('a,b,c,d,e,f,g,h,m,n'):
for i from 1 to times do
tot:=tot+(N30[i+1]-An30(i/(times)))^2:
end do:
l:=[b,c,d,e,f,h,m,n]:
DegreeAn30app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0},{b,c,d,e,f,h,m,n}):
NwAn30:=proc(x)
NwAn30(x):=eval(An30(x),DegreeAn30app):
end:
print(NwAn30(x));
Er:=[];
print(Errors);
for s from 0 to 10 do
Er:=[op(Er),abs(NwAn30(s/(times))-N30[s+1])]:
end do;
unassign('b,c,d,e,f,h,j,k,m');
E_3:=proc(x)
E_3(x):=N30_1[1]+b*x+c*x^2+d*x^3+e*x^4+f*x^5+h*x^6+j*x^7+k*x^8+m*x^9;
end:

l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 1 to times do
tot:=tot+(N30_1[i+1]-E_3(i/(times)))^2:
end do:
DegreeN_320app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m});
NwN_320:=proc(x)
NwN_320(x):=eval(E_3(x),DegreeN_320app);
end:
print(NwN_320(x));
unassign('b,c,d,e,f,h,j,k,m');
E_1:=proc(x)
E_1(x):=N30_2[1]+b*x+c*x^2+d*x^3+e*x^4+f*x^5+h*x^6+j*x^7+k*x^8+m*x^9;
end:
l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 1 to times do
tot:=tot+(N30_2[i+1]-E_1(i/(times)))^2:
end do:
DegreeN_120app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m}):
NwN_120:=proc(x)
NwN_120(x):=eval(E_1(x),DegreeN_120app);
end:
print(NwN_120(x));
unassign('b,c,d,e,f,h,j,k,m'):
E_2:=proc(x)
E_2(x):=b*x+c*x^3+d*x^5+e*x^7+f*x^9+h*x^11+j*x^13+k*x^15+m*x^17;
end:
l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 0 to times do
tot:=tot+(N30_3[i+1]-E_2(i/(times)))^2:
end do:
DegreeE_320app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m});
NwN_220:=proc(x)
NwN_220(x):=eval(E_2(x),DegreeE_320app);
end:
print(NwN_220(x));

>

A_p_non-sort is a program inside opt2.Its arguments are the boundary values, the p-power, and the initial polynomial. opt2 given a starting f1(I call NwN_30) p-energy polynomial will give a new polynomial closer to approximating the old p-energy polynomial. This should just be executed the rest of opt2 and polmake will take care of the rest.

> A_p_non_sort:=proc(a_1,a_2,a_3,p,f1)
local A,a1,a2,a3;
A:=[]: A:=sort([a_1,a_2,a_3]): a1:=A[1]: a2:=A[2]: a3:=A[3]:
A_p_non_sort(a1,a2,a3,p):=to_power((a3-a1)/2,p,-1)*f1(((2*a2-a1-a3)/(a3-a1)));
end:

opt2 is the important function for creating new p-energy polynomials. T,m,p are defined as above. est1 is used in decimate, which should also be executed before running any of these programs. It is the final radius in the decimate search. You can read about decimate in the decimate.mws code. est is absolute value of the difference between an initial vertex and its new value. itr is like k above it tells opt to cycle that many times. a,b,c are boundary values. pen,frst,secd,thrd, help in spitting values into the lists I have setup. f1 is the initial polynomial probably from the least squares method above.

> opt2:=proc(T,m,fnum,p,est1,est,itr,a,b,c,pen,frst,secd,thrd,f1)
local ans2,Enpoly,Enpoly_der,tp,Enpoly_non,cou,vm,z,Ave,EVLI,difT,MX,l,t,h,ver,tempor,tempie,G,P,E,tot,i,j,k,En,Myt,temp,pn1,fst1,scd1,trd1,T1,T2,T3,T_list,nx1,ny1,nz1;
tot:=0;
Digits:=11;
vm:=vertices(m):
make_neighbors(T,m,vm):
for ver in vm do
unassign('T[ver,fnum]');
end do;
for h from 0 to 2 do
unassign('T[[h],fnum]');
end do;
#vm:=subsop(1=NULL,2=NULL,3=NULL,vm);
#T[[0],fnum]:=a;
#T[[1],fnum]:=b;
#T[[2],fnum]:=c;
#T[[0,1],fnum]:=x1
#T[[1,2],fnum]:=y1
#T[[0,2],fnum}:=z1
if (b>.35) and (b <> 1) then
Enpoly_non:=proc(x1,y1,z1) Enpoly_non(x1,y1,z1):=A_p_non_sort(-1,x1,z1,p,f1)+A_p_non_sort(x1,y1,b,p,f1)+A_p_non_sort(z1,y1,1,p,f1);
end;
else
Enpoly_non:=proc(x1,y1,z1)
Enpoly_non(x1,y1,z1):=A_p_non_sort(-1,x1,z1,p,f1)+A_p_non_sort(x1,b,y1,p,f1)+A_p_non_sort(z1,y1,1,p,f1);
end;
end if;
optimize(T,m,fnum,p,itr,a,b,c,pn1,fst1,scd1,trd1):
unassign('nx1','ny1','nz1');
decimate_solve(Enpoly_non,T[[0,1],fnum],T[[1,2],fnum],T[[0,2],fnum],.1,est1,nx1,ny1,nz1):

T_list:=[];
T1:=abs(nx1-T[[0,1],fnum]);
T2:=abs(ny1-T[[1,2],fnum]);
T3:=abs(nz1-T[[0,2],fnum]);
T_list:=[T1,T2,T3];
tempor:=max(op(T_list));
T[[0,1],fnum]:=nx1;
T[[1,2],fnum]:=ny1;
T[[0,2],fnum]:=nz1;
while (tempor > est) do
unassign('nx1,ny1,nz1');
decimate_solve(Enpoly_non,T[[0,1],fnum],T[[1,2],fnum],T[[0,2],fnum],est1,(est1)/4,nx1,ny1,nz1):
#print(Enpoly_non(nx1,ny1,nz1));
T_list:=[];
T1:=abs(nx1-T[[0,1],fnum]);
T2:=abs(ny1-T[[1,2],fnum]);
T3:=abs(nz1-T[[0,2],fnum]);
T_list:=[T1,T2,T3];
tempor:=max(op(T_list));
T[[0,1],fnum]:=nx1;
T[[1,2],fnum]:=ny1;
T[[0,2],fnum]:=nz1;
print(inloop);
end do;
print(energy_is);
ans2:=Enpoly_non(T[[0,1],fnum],T[[1,2],fnum],T[[0,2],fnum]);
print(ans2);
pen:=ans2;
frst:=T[[0,1],fnum];
secd:=T[[1,2],fnum];
thrd:=T[[0,2],fnum];
print(vertices_are);
print(T[[0,1],fnum]); print(T[[1,2],fnum]); print(T[[0,2],fnum]);
end:

This cycles opt2 for different boundary values. This makes the final polynomial, and its p-harmonic extensions. p is the power. counte is the number of cycles of evolving the p-energy polynomial, and times is the number of points that correspond to different boundary values. I have some examples below.

> polymake:=proc(p,counte,times,f1)
local l,N30,N30_1,N30_2,Er,s,N30_3,f30,tot,w,Gr100,NwN_both,E_1,E_2,E_3,a,b,c,d,e,f,g,h,j,k,i,An30,DegreeAn30app,NwAn30,DegreeN_320app,NwN_320,DegreeN_120app,NwN_120,DegreeE_320app,NwN_220,f2,as,bs,cs,ds,es,fs,hs,js,ks,ls,ms,ns;
Digits:=11;
f2:=f1;
for w from 0 to counte do
unassign('a,b,c,d,e,f,g,h,k,l,m,n');
N30:=[]: N30_1:=[]: N30_2:=[]: N30_3:=[]:
for i from 0 to times do
unassign('pn');
unassign('fst');
unassign('scd');
unassign('trd');
opt2(T,1,8,p,.000001,.000001,35,-1,i/(times),1,pn,fst,scd,trd,f2);
N30:=[op(N30),pn];
print(N30);
N30_1:=[op(N30_1),fst];
print(N30_1);
N30_2:=[op(N30_2),scd];
print(N30_2);
N30_3:=[op(N30_3),trd];
print(N30_3);
end do;
f30:=proc(x)
f30(x):=x/N30[11];
end;
N30:=map(f30,N30);
print(N30);
unassign('a,b,c,d,e,f,h,m,n');
as:=N30[1]:
An30:=proc(x)
An30(x):=as+b*x^2+c*x^4+d*x^6+e*x^8+f*x^10+h*x^12+m*x^14+n*x^16;
end;

> with(linalg):
tot:=0;
unassign('a,b,c,d,e,f,g,h,m,n'):
for i from 1 to times do
tot:=tot+(N30[i+1]-An30(i/(times)))^2:
end do:
l:=[b,c,d,e,f,h,m,n]:
DegreeAn30app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0},{b,c,d,e,f,h,m,n}):
NwAn30:=proc(x)
NwAn30(x):=eval(An30(x),DegreeAn30app):
end:
f2:=NwAn30:
print(NwAn30(x));
Er:=[];
print(Errors);
for s from 0 to 10 do
Er:=[op(Er),abs(NwAn30(s/(times))-N30[s+1])]:
end do;
print(Er);
unassign('b,c,d,e,f,h,j,k,m');
E_3:=proc(x)
E_3(x):=N30_1[1]+b*x+c*x^2+d*x^3+e*x^4+f*x^5+h*x^6+j*x^7+k*x^8+m*x^9;
end;

l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 1 to times do
tot:=tot+(N30_1[i+1]-E_3(i/(times)))^2:
end do:
DegreeN_320app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m});
NwN_320:=proc(x)
NwN_320(x):=eval(E_3(x),DegreeN_320app);
end;
print(NwN_320(x));
unassign('b,c,d,e,f,h,j,k,m');
E_1:=proc(x)
E_1(x):=N30_2[1]+b*x+c*x^2+d*x^3+e*x^4+f*x^5+h*x^6+j*x^7+k*x^8+m*x^9;
end;
l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 1 to times do
tot:=tot+(N30_2[i+1]-E_1(i/(times)))^2:
end do:
DegreeN_120app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m}):
NwN_120:=proc(x)
NwN_120(x):=eval(E_1(x),DegreeN_120app);
end;
print(NwN_120(x));
unassign('b,c,d,e,f,h,j,k,m'):
E_2:=proc(x)
E_2(x):=b*x+c*x^3+d*x^5+e*x^7+f*x^9+h*x^11+j*x^13+k*x^15+m*x^17;
end;
l:=[b,c,d,e,f,h,j,k,m];
tot:=0;
for i from 0 to times do
tot:=tot+(N30_3[i+1]-E_2(i/(times)))^2:
end do:
DegreeE_320app:=fsolve({Gr100[1]=0,Gr100[2]=0,Gr100[3]=0,Gr100[4]=0,Gr100[5]=0,Gr100[6]=0,Gr100[7]=0,Gr100[8]=0,Gr100[9]=0},{b,c,d,e,f,h,j,k,m});
NwN_220:=proc(x)
NwN_220(x):=eval(E_2(x),DegreeE_320app);
end;
print(NwN_220(x));
end do;
NwN_both:=proc(x)
if (x>=0) then
NwN_both(x):=NwN_120(x):
else
NwN_both(x):=-NwN_320(-x):
end if:
end:
plot(NwN_both, -1..1, style=line,thickness=3);
plot(NwN_220, -1..1,-(1/5)...(1/5) ,style=line ,thickness=3);
end:

these are some examples of the program.

>
NwAn30:=proc(x)
NwAn30(x):=.77806771929+.23894976631*x^2-.77030730821e-1*x^4+.32033503861*x^6-.87405331971*x^8+1.3104349877*x^10-.99943270749*x^12+.30676026747*x^14-.403100300e-2*x^16:
end;
polymake(12/7,2,10,NwAn30);
print(next1);
NwAn30:=proc(x)
NwAn30(x):=.71763521754+.2320738035*x^2+.2487171443*x^4-1.0443644595*x^6+3.0082640570*x^8-5.2150093545*x^10+5.2546360752*x^12-2.8279587105*x^14+.6260061898*x^16:
polymake(12/5,5,10,NwAn30);
print(next1);
NwAn30:=proc(x)
NwAn30(x):=.7811453332+.208273466*x^2-.305952879*x^4+1.991207265*x^6-3.742578664*x^8-.787376029*x^10+10.78661085*x^12-12.44164914*x^14+4.51031982*x^16:
end:
polymake(16/9,5,10,NwAn30);
print(next1);
NwAn30:=proc(x)
NwAn30(x):=.7153884018+.296313196*x^2+.214162952*x^4-1.536024583*x^6+3.096039197*x^8-.57364606e-1*x^10-7.304440561*x^12+8.819441661*x^14-3.243515661*x^16:
end:
polymake(16/7,5,10,NwAn30);
print(endy);