text:
# This file sets up the letters of the alphabet as one-dimensional simplicial
# complexes.

letters := table():

alphabet_list := [["A","B","C","D","E"],
                  ["F","G","H","I","J"],
                  ["K","L","M","N","O"],
                  ["P","Q","R","S","T"],
                  ["U","V","W","X","Y"]]:

alphabet_type_list := [
 ["C","G","L","M","N","S","U","V","W"],
 ["E","F","J","T","Y"],
 ["H","I"],
 ["K","X"],
 ["O","D","A","Q"],
 ["B"]
]:

set_map := proc(T)
 local P,e,s,t;
 T["map"] := table():
 P := T["embedding"]:
 for e in T["max_simplices"] do
  T["map"][e] := unapply((1-t) *~ P[e[1]] +~ t *~ P[e[2]],t);
 od:
end:

set_alt_map := proc(T)
 local P,e,s,t;
 T["alt_map"] := table():
 P := T["alt_embedding"]:
 for e in T["max_simplices"] do
  T["alt_map"][e] := unapply((1-t) *~ P[e[1]] +~ t *~ P[e[2]],t);
 od:
end:

set_isotopy := proc(T)
 local P,Q,e,s,t;
 T["isotopy"] := table():
 for e in T["max_simplices"] do
  T["isotopy"][e] :=
   unapply((1-s) *~ T["map"][e](t) +~ s *~ T["alt_map"][e](t),s,t);
 od:
end:

rotate_isotopy := proc(T,a)
 local P0,P1,P2,e,v,s,t;
 P0 := eval(T["alt_embedding"]):
 P1 := eval(T["alt_map"]):
 P2 := eval(T["isotopy"]):
 T["alt_embedding"] := table():
 T["alt_map"] := table():
 T["isotopy"] := table():
 for v in T["vertices"] do
  T["alt_embedding"][v] := rot(a,P0[v]):
 od:
 for e in T["max_simplices"] do
  T["alt_map"][e] := unapply(rot(a,P1[e](t)),t):
  T["isotopy"][e] := unapply(rot(a*s,P2[e](s,t)),s,t):
 od:
end:

translate_letter := proc(T,a)
 local P0,P1,P2,P3,P4,e,v,s,t;
 P0 := eval(T["embedding"]):
 P1 := eval(T["map"]):
 P2 := eval(T["alt_embedding"]):
 P3 := eval(T["alt_map"]):
 P4 := eval(T["isotopy"]):
 T["embedding"] := table():
 T["map"] := table():
 T["alt_embedding"] := table():
 T["alt_map"] := table():
 T["isotopy"] := table():
 for v in T["vertices"] do
  T["embedding"][v] := P0[v] +~ a:
  T["alt_embedding"][v] := P3[v] +~ a:
 od:
 for e in T["max_simplices"] do
  T["map"][e] := unapply(P1[e](t) +~ a,t):
  T["alt_map"][e] := unapply(P3[e](t) +~ a,t):
  T["isotopy"][e] := unapply(P2[e](s,t) +~ a,s,t):
 od:
end:

pol := (r,t) -> evalf([r * cos(Pi*t),r * sin(Pi*t)]):
rot := (t,x) -> evalf([cos(Pi*t)*x[1]-sin(Pi*t)*x[2],sin(Pi*t)*x[1] + cos(Pi*t)*x[2]]):

show_map := (T,s) ->
 display(seq(plot([op(T["map"][e](t)),t=0..1]),e in T["max_simplices"]),axes=none,scaling=constrained):

show_alt_map := (T,s) ->
 display(seq(plot([op(T["alt_map"][e](t)),t=0..1]),e in T["max_simplices"]),axes=none,scaling=constrained):

show_isotopy := (T,s) ->
 display(seq(plot([op(T["isotopy"][e](s,t)),t=0..1]),e in T["max_simplices"]),axes=none,scaling=constrained):

show_movie := (T) ->
 display([seq(show_isotopy(T,s),s=0..1,0.02)],insequence=true,view=[-1.5..1.5,-1.5..1.5]):

show_list_map := proc(L)
 display(seq(seq(
   translate(show_map(letters[L[i][j]]),2*j,-3*i),
  j=1..nops(L[i])),i=1..nops(L)),scaling=constrained);
end:

show_list_alt_map := proc(L)
 display(seq(seq(
   translate(show_alt_map(letters[L[i][j]]),3*j,-3*i),
  j=1..nops(L[i])),i=1..nops(L)),scaling=constrained);
end:

show_list_isotopy := proc(L,s)
 display(seq(seq(
   translate(show_isotopy(letters[L[i][j]],s),3*j,-3*i),
  j=1..nops(L[i])),i=1..nops(L)),scaling=constrained);
end:

show_list_movie := proc(L)
 display([seq(show_list_isotopy(L,s),s=0..1,0.02)],insequence=true):
end:

tikz_map := proc(T)
 local s,e,u,n;
 s := "";
 for e in T["max_simplices"] do
  u := T["map"][e];
  if type(u(t),list(polynom(realcons,t))) and
     max(map(degree,u,t)) <= 1 then
   n := 1;
  else
   n := 16;
  fi;
  s := cat(s," \\draw ",
    StringTools[Join]([seq(sprintf("(%.3f,%.3f)",op(u(i/n))),i=0..n)]," -- "),
    ";\n");
 od:
 return s;
end:

tikz_alt_map := proc(T)
 local s,e,u,n;
 s := "";
 for e in T["max_simplices"] do
  u := T["alt_map"][e];
  if type(u(t),list(polynom(realcons,t))) and
     max(map(degree,u,t)) <= 1 then
   n := 1;
  else
   n := 16;
  fi;
  s := cat(s," \\draw ",
    StringTools[Join]([seq(sprintf("(%.3f,%.3f)",op(u(i/n))),i=0..n)]," -- "),
    ";\n");
 od:
 return s;
end:

tikz_list_map := proc(L,{width::numeric := 2},{height::numeric := 3})
 cat(seq(seq(
  cat(
   sprintf("\\begin{scope}[shift={(%f,%f)}]\n",width*j,-height*i),
   tikz_map(letters[L[i][j]]),
   "\\end{scope}\n"),
   j=1..nops(L[i])),i=1..nops(L)));
end:

tikz_list_alt_map := proc(L,{width::numeric := 3.2},{height::numeric := 3})
 cat(seq(seq(
  cat(
   sprintf("\\begin{scope}[shift={(%f,%f)}]\n",width*j,-height*i),
   tikz_alt_map(letters[L[i][j]]),
   "\\end{scope}\n"),
   j=1..nops(L[i])),i=1..nops(L)));
end:

######################################################################
# A

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[4,5],[2,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([
 1 = [-1,-1], 2 = [-0.5,0], 3 = [0,1], 4 = [0.5,0], 5 = [1,-1]
]):
T["alt_embedding"] := table([
 1 = [-1.5,0], 2 = [-1,0], 3 = [0,1], 4 = [1,0], 5 = [1.5,0]
]):
set_map(T):
set_alt_map(T):
T["alt_map"][[2,3]] := (t) -> [-cos(Pi/2*t),sin(Pi/2*t)]:
T["alt_map"][[3,4]] := (t) -> [ sin(Pi/2*t),cos(Pi/2*t)]:
T["alt_map"][[2,4]] := (t) -> [-cos(Pi*t),-sin(Pi*t)]:
set_isotopy(T):
letters["A"] := eval(T):

######################################################################
# B

T := table():
T["vertices"] := [1,2,3,4,5,6]:
T["max_simplices"] := [[1,2],[1,4],[2,3],[2,5],[3,6],[4,5],[5,6]]:
T["embedding_dim"] := 2:
T["embedding"] := table([
 1 = [-0.5,-1],2=[-0.5,0],3=[-0.5,1],4=[0,-1],5=[0,0],6=[0,1]
]):
T["alt_embedding"] := table([
 1 = pol(1,-2/3), 2 = pol(1,1), 3 = pol(1,2/3), 4 = pol(1,-1/3), 5 = pol(1,0), 6 = pol(1,1/3)
]):
set_map(T):
set_alt_map(T):
T["map"][[4,5]] := (t) -> [0,-0.5] +~ pol(0.5,t-0.5):
T["map"][[5,6]] := (t) -> [0, 0.5] +~ pol(0.5,t-0.5):
T["alt_map"][[1,2]] := (t) -> pol(1,(-2-t)/3):
T["alt_map"][[1,4]] := (t) -> pol(1,(-2+t)/3):
T["alt_map"][[2,3]] := (t) -> pol(1,(-3-t)/3):
T["alt_map"][[3,6]] := (t) -> pol(1,( 2-t)/3):
T["alt_map"][[4,5]] := (t) -> pol(1,(-1+t)/3):
T["alt_map"][[5,6]] := (t) -> pol(1,(   t)/3):
set_isotopy(T):
letters["B"] := eval(T):

######################################################################
# C

T := table():
T["vertices"] := [1,2]:
T["max_simplices"] := [[1,2]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1 = pol(1,0.4) *~ [0.8,1], 2 = pol(1,-0.4) *~ [0.8,1]]):
T["alt_embedding"] := table([1 = [0,1], 2=[0,-1]]):
T["map"] := table([[1,2] = ((t) -> pol(1,0.4 + 1.2*t) *~ [0.8,1])]):
set_alt_map(T):
set_isotopy(T):
letters["C"] := eval(T):

######################################################################
# D

T := table():
T["vertices"] := [1,2,3]:
T["max_simplices"] := [[1,2],[1,3],[2,3]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1 = [-0.5,-1], 2=[-0.5,1], 3 = [0.5,0]]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,1], 3 = [1,0]]):
set_map(T):
T["map"][[1,3]] := (t) -> [-0.5,0] +~ pol(1,-0.5 + 0.5*t):
T["map"][[2,3]] := (t) -> [-0.5,0] +~ pol(1, 0.5 - 0.5*t):
T["alt_map"] := table():
T["alt_map"][[1,2]] := (t) -> pol(1,-0.5 - t):
T["alt_map"][[1,3]] := (t) -> pol(1,-0.5 + 0.5*t):
T["alt_map"][[2,3]] := (t) -> pol(1, 0.5 - 0.5*t):
set_isotopy(T):
letters["D"] := eval(T):

######################################################################
# E

T := table():
T["vertices"] := [1,2,3,4,5,6]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[4,5],[3,6]]:
T["embedding_dim"] := 2:
T["embedding"] := table([
 1 = [0.5,-1], 2 = [-0.5,-1], 3 = [-0.5,0], 4 = [-0.5,1], 5 = [0.5,1], 6 = [0,0]
]):
T["alt_embedding"] := table([
 1 = [-1,-1], 2 = [-1,-0.5], 3 = [-1,0], 4 = [-1,0.5], 5 = [-1,1], 6 = [1,0]
]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
T["isotopy"][[1,2]] := (s,t) -> (1-s) *~ [-0.5,-1] +~ s *~ [-1,-0.5] +~ pol((1-0.5*s)*(1-t),-s/2):
T["isotopy"][[4,5]] := (s,t) -> (1-s) *~ [-0.5, 1] +~ s *~ [-1, 0.5] +~ pol((1-0.5*s)*(1-t), s/2):
rotate_isotopy(T,-0.5):
letters["E"] := eval(T):

######################################################################
# F

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[2,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([
 1 = [-0.5,-1], 2 = [-0.5,0], 3 = [-0.5,1], 4 = [0.5,1], 5 = [0,0]
]):
T["alt_embedding"] := table([
 1 = [-1,-1], 2 = [-1,0], 3 = [-1,0.5], 4 = [-1,1], 5 = [1,0]
]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
T["isotopy"][[3,4]] := (s,t) -> (1-s) *~ [-0.5, 1] +~ s *~ [-1, 0.5] +~ pol((1-0.5*s)*(1-t), s/2):
rotate_isotopy(T,-0.5):
letters["F"] := eval(T):

######################################################################
# G

T := table():
T["vertices"] := [1,2,3]:
T["max_simplices"] := [[1,2],[2,3]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1 = pol(1,0.4) *~ [0.8,1], 2 = pol(1,-0.4) *~ [0.8,1],
                         3 = pol(1,-0.4) *~ [0.8,1] +~ [0,-0.4]
                         ]):
T["alt_embedding"] := table([1 = [0,1], 2=[0,0], 3=[0,-1]]):
set_map(T):
T["map"][[1,2]] := ((t) -> pol(1,0.4 + 1.2*t) *~ [0.8,1]):
set_alt_map(T):
set_isotopy(T):
letters["G"] := eval(T):

######################################################################
# H

T := table():
T["vertices"] := [1,2,3,4,5,6]:
T["max_simplices"] := [[1,2],[2,3],[2,5],[4,5],[5,6]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.5,-1],2=[-0.5,0],3=[-0.5,1],
                         4=[ 0.5,-1],5=[ 0.5,0],6=[ 0.5,1]
                         ]):
T["alt_embedding"] := table([1 = [-1,-1], 2=[-1,0], 3=[-1,1],
                             4=[1,-1],5=[1,0],6=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["H"] := eval(T):

######################################################################
# I

T := table():
T["vertices"] := [1,2,3,4,5,6]:
T["max_simplices"] := [[1,2],[2,3],[2,5],[4,5],[5,6]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.2,-1],2=[0,-1],3=[0.2,-1],
                         4=[-0.2,1],5=[0,1],6=[0.2,1]
                         ]):
T["alt_embedding"] := table([1 = [-1,-1], 2=[0,-1], 3=[1,-1],
                             4=[-1,1],5=[0,1],6=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["I"] := eval(T):

######################################################################
# J

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,4],[3,4],[4,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.6,-0.7],2=[0,-0.7],3=[-0.5,1],
                         4=[0,1],5=[0.5,1]
                         ]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,0], 3=[-1,1],
                             4=[0,1],5=[1,1]]):
set_map(T):
T["map"][[1,2]] := (t) -> [-0.3,-0.7] +~ pol(0.3,-1+t):
set_alt_map(T):
set_isotopy(T):
letters["J"] := eval(T):

######################################################################
# K

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[2,4],[2,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.5,-1],2=[-0.5,0],3=[-0.5,1],
                         4=[0.5,-1],5=[0.5,1]
                         ]):
T["alt_embedding"] := table([1 = [-1,-1], 2=[0,0], 3=[-1,1],
                             4=[1,-1],5=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["K"] := eval(T):

######################################################################
# L

T := table():
T["vertices"] := [1,2,3]:
T["max_simplices"] := [[1,2],[2,3]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[0.5,-1],2=[-0.5,-1],3=[-0.5,1]]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,0], 3=[0,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["L"] := eval(T):

######################################################################
# M

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[4,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.6,-1],2=[-0.3,1],3=[0,-1],4=[0.3,1],5=[0.6,-1]]):
T["alt_embedding"] := table([1 = [-1,0], 2=[-0.5,0], 3=[0,0],4=[0.4,0],5=[1,0]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["M"] := eval(T):

######################################################################
# N

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[3,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.5,-1],2=[-0.5,1],3=[0.5,-1],4=[0.5,1]]):
T["alt_embedding"] := table([1 = [-1,0], 2=[-0.33,0], 3=[0.33,0],4=[1,0]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["N"] := eval(T):

######################################################################
# O

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[1,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=pol(1,0) *~ [0.6,1], 2=pol(1,0.5) *~ [0.6,1],
                         3=pol(1,1) *~ [0.6,1], 4=pol(1,1.5) *~ [0.6,1]]):
T["alt_embedding"] := table([1 = pol(1,0),2=pol(1,0.5),3=pol(1,1),4=pol(1,1.5)]):
set_map(T):
set_alt_map(T):
T["map"][[1,2]] := (t) -> pol(1,0.5* t   ) *~ [0.6,1]:
T["map"][[2,3]] := (t) -> pol(1,0.5*(t+1)) *~ [0.6,1]:
T["map"][[3,4]] := (t) -> pol(1,0.5*(t+2)) *~ [0.6,1]:
T["map"][[1,4]] := (t) -> pol(1,0.5*(4-t)) *~ [0.6,1]:
T["alt_map"][[1,2]] := (t) -> pol(1,0.5* t   ):
T["alt_map"][[2,3]] := (t) -> pol(1,0.5*(t+1)):
T["alt_map"][[3,4]] := (t) -> pol(1,0.5*(t+2)):
T["alt_map"][[1,4]] := (t) -> pol(1,0.5*(4-t)):
set_isotopy(T):
letters["O"] := eval(T):

######################################################################
# P

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[2,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.25,-1],2=[-0.25,0],3=[-0.25,1],4=[0.5,0.5]]):
T["alt_embedding"] := table([1=[0,-1.5],2=[0,-1],3=pol(1,2/3),4=pol(1,1/3)]):
set_map(T):
set_alt_map(T):
T["map"][[2,4]] := (t) -> pol(0.5,0.5*(t-1)) *~ [1.5,1] +~ [-0.25,0.5]:
T["map"][[3,4]] := (t) -> pol(0.5,0.5*(1-t)) *~ [1.5,1] +~ [-0.25,0.5]:
T["alt_map"][[2,3]] := (t) -> pol(1,-1/2-2/3*t):
T["alt_map"][[2,4]] := (t) -> pol(1,-1/2+2/3*t):
T["alt_map"][[3,4]] := (t) -> pol(1,5/6-2/3*t):
set_isotopy(T):
rotate_isotopy(T,0.5):
letters["P"] := eval(T):

######################################################################
# Q

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[1,3],[1,4],[1,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([
 1= pol(1,-1/4)   *~ [0.6,1],
 2= pol(1,5/12)   *~ [0.6,1],
 3= pol(1,13/12)  *~ [0.6,1],
 4= pol(1.3,-1/4) *~ [0.6,1],
 5= pol(0.7,-1/4) *~ [0.6,1]
]):
T["alt_embedding"] :=
table([
 1= pol(1,-1/4),
 2= pol(1,5/12),
 3= pol(1,13/12),
 4= pol(1.5,-1/4),
 5= pol(0.5,-1/4)
]):
set_map(T):
set_alt_map(T):
T["map"][[1,2]] := (t) -> pol(1,-1/4+2/3*t) *~ [0.6,1]:
T["map"][[1,3]] := (t) -> pol(1,-1/4-2/3*t) *~ [0.6,1]:
T["map"][[2,3]] := (t) -> pol(1,5/12+2/3*t) *~ [0.6,1]:
T["alt_map"][[1,2]] := (t) -> pol(1,-1/4+2/3*t):
T["alt_map"][[1,3]] := (t) -> pol(1,-1/4-2/3*t):
T["alt_map"][[2,3]] := (t) -> pol(1,5/12+2/3*t):
set_isotopy(T):
rotate_isotopy(T,0.25):
letters["Q"] := eval(T):

######################################################################
# R

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[2,4],[2,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.25,-1],2=[-0.25,0],3=[-0.25,1],4=[0.5,0.5],5=[0.5,-1]]):
T["alt_embedding"] := table([1=[-0.5,-1.5],2=[0,-1],3=pol(1,2/3),4=pol(1,1/3),5=[0.5,-1.5]]):
set_map(T):
set_alt_map(T):
T["map"][[2,4]] := (t) -> pol(0.5,0.5*(t-1)) *~ [1.5,1] +~ [-0.25,0.5]:
T["map"][[3,4]] := (t) -> pol(0.5,0.5*(1-t)) *~ [1.5,1] +~ [-0.25,0.5]:
T["alt_map"][[2,3]] := (t) -> pol(1,-1/2-2/3*t):
T["alt_map"][[2,4]] := (t) -> pol(1,-1/2+2/3*t):
T["alt_map"][[3,4]] := (t) -> pol(1,5/6-2/3*t):
set_isotopy(T):
rotate_isotopy(T,0.5):
letters["R"] := eval(T):

######################################################################
# S

T := table():
T["vertices"] := [1,2,3]:
T["max_simplices"] := [[1,2],[2,3]]:
T["embedding_dim"] := 2:
theta := -0.05;
T["embedding"] := table([
 1=rot(theta,[0,-0.5] +~ pol(0.5,-0.7)),
 2=rot(theta,[0,0],3=[0,0.5] +~ pol(0.5,0.3))
 ]):
T["alt_embedding"] := table([1=[0,-1],2=[0,0],3=[0,1]]):
set_map(T):
set_alt_map(T):
T["map"][[1,2]] := (t) -> rot(theta,[0,-0.5] +~ pol(0.5,-0.7+1.2*t)):
T["map"][[2,3]] := (t) -> rot(theta,[0, 0.5] +~ pol(0.5,-0.5-1.2*t)):
set_isotopy(T):
letters["S"] := eval(T):

######################################################################
# T

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[2,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[0,-1],2=[0,1],3=[-0.7,1],4=[0.7,1]]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,1], 3=[-1,1],4=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["T"] := eval(T):

######################################################################
# U

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[3,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.6,1],2=[-0.6,-0.4],3=[0.6,-0.4],4=[0.6,1]]):
T["alt_embedding"] := table([1 = [-1,0], 2=[-0.33,0], 3=[0.33,0],4=[1,0]]):
set_map(T):
set_alt_map(T):
T["map"][[2,3]] := (t) -> [0,-0.4] +~ pol(0.6,t-1):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["U"] := eval(T):

######################################################################
# V

T := table():
T["vertices"] := [1,2,3]:
T["max_simplices"] := [[1,2],[2,3]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.6,1],2=[0,-1],3=[0.6,1]]):
T["alt_embedding"] := table([1 = [-1,0], 2=[0,0], 3=[1,0]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["V"] := eval(T):

######################################################################
# W

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[3,4],[4,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.6,1],2=[-0.3,-1],3=[0,1],4=[0.3,-1],5=[0.6,1]]):
T["alt_embedding"] := table([1 = [-1,0], 2=[-0.5,0], 3=[0,0],4=[0.4,0],5=[1,0]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
rotate_isotopy(T,-0.5):
letters["W"] := eval(T):

######################################################################
# X

T := table():
T["vertices"] := [1,2,3,4,5]:
T["max_simplices"] := [[1,2],[2,3],[2,4],[2,5]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[-0.5,-1],2=[0,0],3=[-0.5,1],
                         4=[0.5,-1],5=[0.5,1]
                         ]):
T["alt_embedding"] := table([1 = [-1,-1], 2=[0,0], 3=[-1,1],
                             4=[1,-1],5=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["X"] := eval(T):

######################################################################
# Y

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[2,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[0,-1],2=[0,0],3=[-0.7,1],4=[0.7,1]]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,1], 3=[-1,1],4=[1,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["Y"] := eval(T):

######################################################################
# Z

T := table():
T["vertices"] := [1,2,3,4]:
T["max_simplices"] := [[1,2],[2,3],[3,4]]:
T["embedding_dim"] := 2:
T["embedding"] := table([1=[0.5,-1],2=[-0.5,-1],3=[0.5,1],4=[-0.5,1]]):
T["alt_embedding"] := table([1 = [0,-1], 2=[0,-0.33], 3=[0,0.33],4=[0,1]]):
set_map(T):
set_alt_map(T):
set_isotopy(T):
letters["Z"] := eval(T):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_mgga_x *)

n_mbeef := 5:
coefs_mbeef := [
    [ 1.17114923e+00,  1.15594371e-01, -5.32167416e-02, -2.01131648e-02,  1.41417107e-03],
    [-6.76157938e-02,  4.53837246e-02, -2.22650139e-02,  1.92374554e-02,  9.19317034e-07],
    [ 1.48659502e-02,  3.18024096e-02, -5.21818079e-03,  1.33707403e-07, -5.00749348e-07],
    [ 1.40794142e-03, -6.08338264e-03, -6.57949254e-07, -5.49909413e-08,  5.74317889e-08],
    [ 1.41530486e-04, -1.00478906e-07,  2.01895739e-07,  3.97324768e-09, -3.40722258e-09]
]:

$include "mbeef.mpl"

f := (rs, x, t, u) ->
  mgga_mbeef(x, t):

ext: mpl

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
input := FileTools:-Text:-ReadFile("AoC-2021-3-input.txt" ):
bins := map(b->[seq(ifelse(b[i]="0",0,1), i=1..length(b))],
            StringTools:-Split(StringTools:-Trim(input)) ):

mostcommon := proc(L)
    map(b->ifelse(b>=nops(L)/2, 1, 0),
        [seq(add(L[j][i], j=1..nops(L)), i=1..nops(L[1]))] );
end proc:
leastcommon := proc(L) map(b->ifelse(b=0, 1, 0), mostcommon(L)); end proc:
bin2dec := proc(L) add(L[-i]*2^(i-1), i=1..nops(L)); end proc:

bin2dec(mostcommon(bins)) * bin2dec(leastcommon(bins)); # part 1 answer

o2gen_bins := bins:
for i to nops(bins[1]) while nops(o2gen_bins)<>1 do
   mc := mostcommon(o2gen_bins);
   o2gen_bins := select(r->r[i]=mc[i], o2gen_bins);
end do:

co2scrub_bins := bins:
for i to nops(bins[1]) while nops(co2scrub_bins)<>1 do
   mc := leastcommon(co2scrub_bins);
   co2scrub_bins := select(r->r[i]=mc[i], co2scrub_bins);
end do:

bin2dec( o2gen_bins[] ) * bin2dec( co2scrub_bins[] ); # part 2 answer


text:
RG := proc(A::algebraic,var::symbol,n::integer)
   local B1,B2,B3,i;
   local B;
   B:=A;
   for i from 1 to n do
   B1:=simplify(B*var^(n-i+1));
   B2:=simplify(eval(B1,var=0));
   if (B2 <> 0) then
   B3:=eval(B2,0=var);
   lprint(B3/var^(n-i+1));
   B:=simplify(B-B3/var^(n-i+1));
   end if;
   end do;
   return B;
end proc;


text:
flat_torus_complex := proc(n::posint,m::posint)
 local f,g,T,V,E,F,R,r;

 f := (i) -> modp(i+1,n);
 g := (j) -> modp(j+1,m);

 V := [seq(seq([i,j],j=0..m-1),i=0..n-1)];
 E := [seq(seq([[i,j],[f(i),j]],j=0..m-1),i=0..n-1),
       seq(seq([[i,j],[i,g(j)]],j=0..m-1),i=0..n-1),
       seq(seq([[i,j],[f(i),g(j)]],j=0..m-1),i=0..n-1)];
 F := [seq(seq([[i,j],[f(i),j],[f(i),g(j)]],j=0..m-1),i=0..n-1),
       seq(seq([[i,j],[i,g(j)],[f(i),g(j)]],j=0..m-1),i=0..n-1)];

 T := table([]);

 T["vertices"] := V;
 T["edges"] := E;
 T["faces"] := F;
 T["max_simplices"] := F;
 T["simplices"] := [map(v -> [v],V),op(E),op(F)];

 R := 3;
 r := 1;

 T["embedding_dim"] := 3;
 T["embedding"] :=
  table([seq(seq([i,j] = evalf(torus_embedding(R,r)(i/n,j/m)),j=0..m-1),i=0..n-1)]);

 T["plot"] := `plot/simplicial_complex`(T["vertices"])(T["faces"],3,T["embedding"]);

 return eval(T);

end:


text:
[plugin]
pluginname = PrometheusExporter
description = Export MAD metrics compatible with Prometheus
author = Artanicus
url = https://github.com/Artanicus/mad_prometheus_exporter
version = 0.2.1


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

op_a1 := 1.5214:
op_a2 := 0.5764:
op_b1 := 1.1284:
op_b2 := 0.3183:

(* This wrapper is to avoid overflows in the OP functionals. The
   energy is not affected, since the value is only changed for
   densities that are screened away.  *)
op_b88_zab := (f_x, rs, z, xs0, xs1) ->
  my_piecewise3(
    b88_zab(1, op_enhancement, rs, z, xs0, xs1) = 0,
    DBL_EPSILON,
    b88_zab(1, op_enhancement, rs, z, xs0, xs1)
  ):

op_beta := (rs, z, xs0, xs1) ->
  op_qab/op_b88_zab(op_enhancement, rs, z, xs0, xs1):

op_f_s := (rs, z, xt, xs0, xs1) ->
  - (1 - z^2)*n_total(rs)/4.0
  * (op_a1*op_beta(rs, z, xs0, xs1) + op_a2)
  / (op_beta(rs, z, xs0, xs1)^4 + op_b1*op_beta(rs, z, xs0, xs1)^3 + op_b2*op_beta(rs, z, xs0, xs1)^2)
:

op_f := (rs, z, xt, xs0, xs1) ->
   my_piecewise3(1 - abs(z) <= p_a_zeta_threshold or (screen_dens(rs,z) and screen_dens(rs,-z)), 0,
                 op_f_s(rs, z_thr(z), xt, xs0, xs1)):

f := (rs, z, xt, xs0, xs1) ->
  op_f(rs, z, xt, xs0, xs1):


text:
 # a function returning (i + j) * -998
 func $foo (
# var %i xxx,
  var %i i32,
   var %j i32,
   var %k i32) i32 {
   return (
     select i32 (
      dread i32 %i,
      dread i32 %j,
      dread i32 %k))}

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
repeat process ACV_teleportation {
/scoreboard players reset @e ACV_TP_FROM
/scoreboard players reset @e ACV_TP_UNEVEN
/scoreboard players reset @e ACV_TP_TO

/execute @e[type=area_effect_cloud,name=ACV_PortalIn] ~ ~ ~ scoreboard players operation @e[tag=!ACV_TpProof,dy=0] ACV_TP_FROM = @e[type=area_effect_cloud,name=ACV_PortalIn,r=0,c=1] ACV_COLOR

// Calulate the target portal
/execute @e[score_ACV_TP_FROM_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_TO = @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_FROM
/execute @e[score_ACV_TP_FROM_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_UNEVEN = @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_FROM
/scoreboard players operation @e[score_ACV_TP_FROM_min=0] ACV_TP_UNEVEN %= 2 ACV_CONSTANT
/scoreboard players add @e[score_ACV_TP_UNEVEN_min=0,score_ACV_TP_UNEVEN=0] ACV_TP_TO 1
/scoreboard players remove @e[score_ACV_TP_UNEVEN_min=1,score_ACV_TP_UNEVEN=1] ACV_TP_TO 1

// Teleport the entity to the other portal
/tp @e[score_ACV_TP_TO_min=0,score_ACV_TP_TO=0] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=0,score_ACV_COLOR=0]
/tp @e[score_ACV_TP_TO_min=1,score_ACV_TP_TO=1] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=1,score_ACV_COLOR=1]
/tp @e[score_ACV_TP_TO_min=2,score_ACV_TP_TO=2] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=2,score_ACV_COLOR=2]
/tp @e[score_ACV_TP_TO_min=3,score_ACV_TP_TO=3] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=3,score_ACV_COLOR=3]
}


text:
`is_element/raw_matroids` := proc(M)
 global reason;
 local n,s,t,s1,t1,u,a,b;

 if not(type(M,set(set))) then
  reason := ["not a set of sets"];
  return false;
 fi;

 n := nops(M);
 for s in M do
  for t in M do
   if s <> t then
    s1 := s minus t;
    t1 := t minus s;
    if s1 = {} or t1 = {} then
     reason := ["s and t are nested",s,t];
     return false;
    fi;
    for a in s1 do
     u := select(b -> member((s minus {a}) union {b},M),t1);
     if u = {} then
      reason := ["exchange fails",s,t];
      return false;
     fi;
    od;
   fi;
  od;
 od;

 return true;
end:

`vertices/raw_matroid` := (M) -> map(op,{op(M)});

`is_basis/raw_matroid` := (M,s) -> member(s,M);

`is_independent/raw_simplicial_complex` := proc(M,s)
 local s0,u;

 if not(type(s,{list,set})) then return false; fi;

 s0 := {op(s)};

 for u in M do
  if s0 minus {op(u)} = {} then
   return true;
  fi;
 od:
 return false;
end:

`is_dependent/raw_simplicial_complex` := proc(M,s)
 local s0,u;

 if not(type(s,{list,set})) then return false; fi;

 s0 := {op(s)};

 for u in M do
  if s0 minus {op(u)} = {} then
   return false;
  fi;
 od:
 return true;
end:

`is_element/matroids` := proc(T)
 global reason;

 if not(type(T,table)) then
  return false;
 fi;

 if not (`is_element/raw_matroids`(T["bases"])) then
  return false;
 fi;

 if `vertices/raw_matroid`(T["bases"]) minus T["vertices"] <> {} then
  reason := ["vertex mismatch",T["vertices"],`vertices/raw_matroid`(T["bases"])];
  return false;
 fi;

 return true;
end:

`cook/raw_matroid` := proc(M,V_)
 local T;

 T := table():
 if nargs > 1 then
  T["vertices"] := V_;
 else
  T["vertices"] := `vertices/raw_matroid`(M);
 fi;
 T["bases"] := M;

end:

`set_rank/matroid` := proc(T)
 local m,r,n,P,A;

 m := nops(T["vertices"]);
 r := table():
 n := table():
 P := combinat[powerset](T["vertices"]);
 for A in P do
  r[A] := max(op(map(s -> nops(s intersect A),T["bases"])));
  n[A] := m - r[A];
 od:

 T["rank"] := eval(r);
 T["nullity"] := eval(r);
end:

`set_closure/matroid` := proc(T)
 local r,V,P,c,A,m;

 if not(type(T["rank"],table)) then `set_rank/matroid`(T); fi;

 r := eval(T["rank"]);
 V := T["vertices"];
 P := combinat[powerset](T["vertices"]);
 c := table():
 for A in P do
  m := r[A];
  c[A] := A union select(b -> r[A union {b}] = m + 1,V minus A);
 od:

 T["closure"] := eval(c):
 T["flats"] := select(A -> nops(c[A]) = nops(A),P);
end:

`res/matroid` := proc(T,A)
 local U,B,m;

 if not(type(A,set) and A minus T["vertices"] = {}) then
  error("A is not a subset of the vertices of T");
 fi;

 U := table():
 U["vertices"] := A;
 B := map(s -> s intersect A,T["bases"]);
 m := max(map(nops,B));
 U["bases"] := select(s -> nops(s) = m,T["bases"]);
end:

uniform_matroid := proc(A_::{nonnegint,set},k)
 local T,A;

 if type(A_,nonnegint) then
  A := {seq(i,i=1..A_)};
 else
  A := A_;
 fi;

 T := table():
 T["vertices"] := A;
 T["bases"] := combinat[choose](A,k);

 return eval(T);
end:

linear_matroid := proc(p,d)
 local T,V,i,j,B,f;

 V := [[]];

 for i from 1 to d do
  V := [seq(seq([op(v),j],j=0..p-1),v in V)];
 od:

 B := [seq([j],j=2..p^d)];
 for i from 2 to d do
  B := [seq(seq([op(b),j],j=b[-1]+1..p^d),b in B)];
 od:

 f := b -> modp(Determinant(Matrix(map(i -> V[i],b))),p);
 B := select(b -> f(b) <> 0,B);

 T := table():
 T["vectors"] := V;
 T["vertices"] := {seq(i,i=1..p^d)};
 T["bases"] := map(b -> {op(b)},{op(B)});

 return eval(T);
end:


text:
#@ Not autoload

# This gives a two-parameter family of embeddings of the trefoil in R^3

trefoil0 := proc(t,a,b)
 [(1 + a * cos(3 * t)) * cos(2 * t),
  (1 + a * cos(3 * t)) * sin(2 * t),
  b * sin(3 * t)]
end:

# This is a particular embedding from the above family for which the
# speed is just sqrt(8)*(1 + 2/3 * cos(3*t)).  For a typical element
# of the above family, the speed is a much more complicated function.

trefoil := (t) -> trefoil0(t,2/3,2/9*sqrt(7));

Trefoil := table():

Trefoil["x"] := trefoil(t,2/3,2*sqrt(7)/9);
Trefoil["t"] := t;
Trefoil["v"] := combine(simplify(expand(map(diff,Trefoil["x"],t)))) assuming t > 0;
Trefoil["speed"] := combine(simplify(factor(expand(`norm_2/R`(3)(Trefoil["v"]))))) assuming t > 0;
Trefoil["s"] := int(Trefoil["speed"],t);
Trefoil["v_hat"] := combine(factor(expand(Trefoil["v"] /~ Trefoil["speed"]))) assuming t > 0;
Trefoil["T"] := Trefoil["v_hat"];
Trefoil["a"] := map(diff,Trefoil["v"],t);
Trefoil["v_hat_dot"] := combine(factor(expand(map(diff,Trefoil["v_hat"],t)))) assuming t > 0;
Trefoil["curvature"] := combine(simplify(factor(expand(`norm_2/R`(3)(Trefoil["v_hat_dot"]))))) assuming t > 0;
Trefoil["radius_of_curvature"] := simplify(1/Trefoil["curvature"]);
Trefoil["N"] := factor(combine(factor(expand(rationalize(Trefoil["v_hat_dot"] /~ Trefoil["curvature"]))))) assuming t > 0;
Trefoil["B"] := factor(combine(factor(expand(`cross_product`(Trefoil["T"],Trefoil["N"]))))) assuming t > 0;


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_x_optx_params *params;

  assert(p->params != NULL);
  params = (gga_x_optx_params * )(p->params);
*)

optx_f := x-> params_a_a + params_a_b*(params_a_gamma*x^2/(1 + params_a_gamma*x^2))^2:

f := (rs, z, xt, xs0, xs1) -> gga_exchange(optx_f, rs, z, xs0, xs1):


text:
#indirect access are associated with pionters and arrays
#iassign<type><field-id>(<addr-expr>, <rhs-expr>)
# type gives the high level type of <addr-expr> and must be a pointer type
# <addr-expr> is computed to return an address
# <rhs-expr> return a value
#iread<prim-type><type><field-id>(<addr-expr>)
#
func $new ( var %n i32 )a64
func $strcpy ( var %dest <*i8>, var %src <*i8> )<*i8>
type $People < struct{
     @name < * i8 >,
     @number < * i32> }>
func $addition ( var %a i32, var %b i32 )i32 {
     return(
	add i32( dread i32 %a, dread i32 %b ))}
var $name < [6] i8 > = [ 104, 117, 97, 119, 101, 105 ]
var $temp<[6] i8> = [ 104, 117, 97, 119, 101, 105 ]
func $main ( ) i32 {
     var %cPtr < * i8 >
     dassign %cPtr ( constval a64 0)
     dassign %cPtr (
     	     addrof a64 $name)

     var %cPtr1 i8
     dassign %cPtr1(
     	     iread i32 < * i8 >( dread a64 %cPtr ))
     dassign %cPtr(
     	     add u64(
	     	 dread a64 %cPtr,
		 constval u64 1))
     var %cPtr2 <* i8>
     dassign %cPtr2(
     	     add u64(
	     	 dread a64 %cPtr,
		 constval u64 1))

     var %iPtr1 <* i32 >
     var %iPtr2 <* i32 >
     iassign < * i32>(
     	     dread a64 %iPtr1,
	     constval i32 10)
     iassign < * i32 >(
     	     dread a64 %iPtr2,
	     iread i32 < * i32> (
	     	   dread a64 %iPtr1))
     var %sPtr <* $People >
     #call $new(
     #	  iread <>)
     call &new (constval u32 128)
     dassign %sPtr(
         regread a64 %%retval
         )
     call &strcpy ( dread a64 %sPtr 0 ,
     	  	    addrof a64 $temp)
     iassign <* i32>(
     	     dread a64 %sPtr 0,
	     constval i32 5)
     var %i i32
     var %j i32
     var %aPtr1 < [10] i32 >
     var %aPtr2	< [10][10] i32 >
     var %aPtr3	< [10][10][10] i32 >
     dassign %i ( constval i32 0 )
     while(
       lt i32 i32 ( dread i32 %i, constval i32 10 )){
       	  iassign < *[10] i32 >(
	  	  array 1 a64 < *[10] i32 >(addrof a64 %aPtr1,
		  		          dread i32 %i),
                  constval i32 10 )

          iassign< *[10] i32 >(
	  	  array 1 a64 < *[10] i32 >( addrof a64 %aPtr1,
		  	                   dread i32 %i),
		  dread i32 %i )

          iassign < *[10] i32 >(
	  	  array 1 a64 < *[10] i32 >( addrof a64 %aPtr1,
		  	      	    	   dread i32 %i ),
		  iread i32 <* i32 >(
		  	addrof a64 %aPtr1 ))
       dassign %i(
	 	 add i32(
		     dread i32 %i,
		     constval i32 1))}

     dassign %i ( constval i32 0 )
     while(
	lt i32 i32 ( dread i32 %i, constval i32 10 )){
	dassign %j ( constval i32 0 )
	while(
		lt i32 i32 ( dread i32 %j, constval i32 10 )){
		   iassign < *[10][10] i32 >(
		   	   array 1 a64 < *[10][10] i32 >( addrof a64 %aPtr2,
		   	   	 dread i32 %i,
				 dread i32 %j),
		           constval i32 10 )
	           iassign < *[10][10] i32 >(
		   	  array 1 a64< *[10][10] i32 >( addrof a64 %aPtr2,
			  	dread i32 %i,
				dread i32 %j),
			  mul i32 (
			      dread i32 %i,
			      dread i32 %j))
		   iassign< *[10][10] i32 >(
		   	    array 1 a64 < *[10][10] i32 >( addrof a64 %aPtr2,
			    	  dread i32 %i,
				  dread i32 %j),
		            iread i32 < *[10] i32 >(
			    	  array 1 a64< *[10] i32 >(
				  	addrof a64 %aPtr1,
				  	rem i32 (
					    mul i32 (
					    	dread i32 %i,
						dread i32 %j),
					    constval i32 10))))
                dassign %j(
			add i32(
		 	    dread i32 %j,
		     	    constval i32 1))}
     dassign %i(
     	     add i32(
	     	 dread i32 %i,
		 constval i32 1))}
     var %k i32
     dassign %i ( constval i32 0)
     while(
	lt i32 i32 ( dread i32 %i, constval i32 10 )){
	dassign %j ( constval i32 0)
	while(
	lt i32 i32 ( dread i32 %j, constval i32 10)){
	   while(
		lt i32 i32 ( dread i32 %k, constval i32 10)){
		iassign < *[10][10][10] i32 >(
			array 1 a64 < *[10][10][10] i32 >(addrof a64 %aPtr3,
			dread i32 %i,
			dread i32 %j,
			dread i32 %k),
			constval i32 10)
		iassign < *[10][10][10] i32 >(
			array 1 a64 < *[10][10][10] i32>(addrof a64 %aPtr3,
			dread i32 %i,
			dread i32 %j,
			dread i32 %k),
			mul i32 (
			    mul i32 (
			    	dread i32 %i,
				dread i32 %j),
			    dread i32 %k))
	        iassign < *[10][10][10] i32 >(
			array 1 a64 < *[10][10][10] i32 >(addrof a64 %aPtr3,
			dread i32 %i,
			dread i32 %j,
			dread i32 %k),
			mul i32 (
			    iread i32 < *[10][10] i32>(
			    	  array 1 a64 <*[10][10] i32>( addrof a64 %aPtr1,
				  dread i32 %i,
				  dread i32 %j)),
				  constval i32 10))
 		dassign %k(
			add i32(
			    dread i32 %k,
			    constval i32 1))}
         dassign %j(
	 	 add i32(
		     dread i32 %j,
		     constval i32 1))}
    dassign %i(
    	    add i32(
	    	dread i32 %i,
		constval i32 1))}

    var %fPtr1 <func( i32, i32) i32 >
    var %fPtr2 <func( i32, i32) i32 >

    dassign %fPtr1(
          addrof a64 $addition)
    dassign %fPtr2(
          addrof a64 $addition)
    var %result i32
    icall ( dread a64 %fPtr1, constval i32 2, constval i32 3)
    dassign %result(
    	    regread i32 %%retval)
    return ( dread i32 %result)}

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
################################################################
# Utility function for irregular term extraction
# Usage:
#
#    gencoeffs(a/x+c+b*x+a*x^3,x);
#
#     Returns ->>  [-1, a], [0, c], [1, b], [2, 0], [3, a]
#
#  (i.e. it returns a list of powers of x and their coefficients
################################################################

gencoeffs := proc(expr::algebraic, var::name)
   local i;
   if (expr <> 0 ) then
   return seq([i,frontend(coeff,[expr,var,i])],
       i=frontend(ldegree,[expr,var]) .. frontend(degree,[expr,var]));
   else
    return [0,0];
   end if;
end proc;

#######################################################
# List of grid functions symbols defined in the problem
# See BSSN_Spher.mw or BSSN_Spher.mpl
#######################################################

grid_functions:={alpha,K,beta,phi,Axx,Athth,A,B,metxx,metthetatheta,DDLxx,DDLthetatheta,rho,TS,Lamx,em4phi,RRxx,RRthetatheta,JSxx,JSthth,Uxx,Pxx,Uthetatheta,Pthetatheta,C5s,C1s,psi,Sx,PI2,PHI2,U,rPI,iPI,rPHI,iPHI,a1,rpsi,ipsi,rPIb2,iPIb2,b1,a2,rpsidot,ipsidot,C5psi,divbeta};

################################################################
# List of tensors, will be used for getting the RHS of equations
# and looping over all of the equations
################################################################

all_tensors:=[
      [EMPH,0,dn],
   [DPH,0,dn] , [dg,2,dn] ,
         [met,2,dn], [DDL,2,dn],
         [DB,0,dn],
   [DTK,0,dn],
         [RR,2,dn], [U,2,dn] , [P,2,dn],
   [DAA,2,dn],
   [DLAM,1,up],
       [C1,0,dn],
       [C5,0,dn],
   [HS,0,dn],
   [HS3,0,dn],
   [MK,1,up]
             ];

# List of coordinates used: (x is r in spherical coordinate here)
coord:=[x,theta,phi];

######################################################################
# The coordiante coefficients that are defined
# in the definition of the generic 3-metric in curvilinear coordinate.
######################################################################

coordcoef:=[1,x^2,x^2*sin(theta)^2];

####################################################################
# Initiating a table for the diff expression of tensors
# tt[A][all][x] will store the x compononent of tensor A
# tt[A][sep][x] will store the x component of tensor A in seperated
# powers of 'r'
# tt[A][reg][x] will store the regular part of above
# tt[A][dis][x] will become the finite difference approximation of
# the above. All in all, 'tt' is our working table of everything!
####################################################################

tt:=table([]):

#########################################
# All of the following atomic variables
# will be used for indexing in table 'tt'
# and shall not be assigned anythin
#########################################

unassign('all');
unassign('sep');
unassign('tres');
unassign('reg');
unassign('dis');
unassign('disreg');
unassign('disinf');
unassign('inf');

#######  IMPORTANT NOTE: ###############################
# Never use i,j,k as indexing as when working with FD as
# they are protected by default to be used as indexing
# of coordinates x, y , z
########################################################

for jj from 1 to nops(all_tensors) do

  # Tensor name
  nm:= all_tensors[jj][1]:

  # Tensor rank
  rank := all_tensors[jj][2]:

  # Tensor type
  tp  := all_tensors[jj][3]:

  # rank-1 tensors
  if rank = 1 then
  tt[nm][all] := table( [ seq( coord[ii] = expand(grcomponent(nm(tp),[coord[ii]])/coordcoef[ii])  ,ii=1..nops(coord)) ] ):
  end if;

  #Assuming rank-2 tensors are all diagnoal (it is the case in spherical symmetry)
  if rank = 2 then
  tt[nm][all] := table( [ seq( coord[ii] = expand(grcomponent(nm(tp,tp),[coord[ii],coord[ii]])/coordcoef[ii])  ,ii=1..nops(coord))  ] ):
  end if;

  if rank = 0 then
  tt[nm][all] := grcomponent(nm):
  end if;

  #Breaking down the expressions to regular and irregular parts

  if rank <> 0 then
  tt[nm][sep] := table (  [  seq( coord[ii] =  [ gencoeffs(tt[nm][all][coord[ii]],'x')   ]  , ii=1..nops(coord) )   ]  ) ;

  end if;

  if rank = 0 then
    tt[nm][sep] := [gencoeffs(tt[nm][all],'x')];
  end if;
end do:

# Reading EOM for matter and adding it to tt table
read("eq_evol_matter.mpl");

evol_matters := [ [EQRPIb2,0,dn] , [EQIPIb2,0,dn] , [EQRPSI,0,dn]  ,[EQIPSI,0,dn], [DRPSI,0,dn], [DIPSI,0,dn] ]:

tt[EQRPIb2][all] := rhs(  eq_evol_rPIb2  ):
tt[EQIPIb2][all] := rhs( eq_evol_iPIb2   ):
tt[EQRPSI][all] := rhs(eq_evol_rpsi):
tt[EQIPSI][all] := rhs(eq_evol_ipsi):
tt[DRPSI][all] := diff(rpsi(t,x),x)+myzero*x:
tt[DIPSI][all] := diff(ipsi(t,x),x)+myzero*x:

for jj from 1 to nops(evol_matters) do
  nm:= evol_matters[jj][1]:
  #Breaking down the expressions to regular and irregular parts
  tt[nm][sep] := [gencoeffs(tt[nm][all],'x')];
end do:

# Updating all_tensors with newly added tensors:

all_tensors:=[
      [EMPH,0,dn],
   [DPH,0,dn] , [dg,2,dn] ,
         [met,2,dn], [DDL,2,dn],
         [DB,0,dn],
   [DTK,0,dn],
         [RR,2,dn], [U,2,dn] , [P,2,dn],
   [DAA,2,dn],
   [DLAM,1,up],
       [C1,0,dn],
       [C5,0,dn],
   [HS,0,dn],
   [HS3,0,dn],
   [MK,1,up],
   [EQRPIb2,0,dn],
   [EQIPIb2,0,dn],
   [EQIPSI,0,dn],
   [EQRPSI,0,dn],
   [DRPSI,0,dn],
   [DIPSI,0,dn]
             ];

#Checking that gencoeffs generated coefficients correctly:
check_res := true;
if (check_res) then

		for jj from 1 to nops(all_tensors) do
		  nm:= all_tensors[jj][1]:
		  rank := all_tensors[jj][2]:
		  tp  := all_tensors[jj][3]:
		  if rank <> 0 then
			tt[nm][tres] := table ( [ seq (  coord[ii] = simplify(expand( tt[nm][all][coord[ii]] - sum(x^tt[nm][sep][coord[ii]][k][1]*tt[nm][sep][coord[ii]][k][2],k=1..nops(tt[nm][sep][coord[ii]])   ) ))      ,ii=1..nops(coord) ) ] );
		  else
			tt[nm][tres] := simplify(expand( tt[nm][all] - sum(x^tt[nm][sep][k][1]*tt[nm][sep][k][2],k=1..nops(tt[nm][sep])   ) ))
		  end if;
		end do:

		for jj from 1 to nops(all_tensors) do
		  nm:= all_tensors[jj][1]:
		  rank := all_tensors[jj][2]:
		  tp  := all_tensors[jj][3]:
		  if rank <> 0 then
			 for ii from 1 to nops(coord) do
			  print(tt[nm][tres][coord[ii]]);
			 end do
		  else
			  print(tt[nm][tres]);
		  end if;

		end do:
end if;

for jj from 1 to nops(all_tensors) do
  nm:= all_tensors[jj][1]:
  rank := all_tensors[jj][2]:
  tp  := all_tensors[jj][3]:
  if rank <> 0 then
    for ii from 1 to nops(coord) do
        if tt[nm][sep][coord[ii]][1][1] >= 0 then
               print(nm + comp + coord[ii] + reg );
        else
                print( nm + comp + coord[ii] + ireg);
                kk := 1;
                while ( tt[nm][sep][coord[ii]][kk][1] < 0 ) do
                   print( tt[nm][sep][coord[ii]][kk][2]*x^tt[nm][sep][coord[ii]][kk][1] );
                    kk := kk + 1;
                end do
        end if;
    end do;
  else
       if tt[nm][sep][1][1] >= 0 then
          print(nm + comp + reg);
       else
          print(nm + comp + ireg);
          kk := 1;
          while ( tt[nm][sep][kk][1] < 0 ) do
                   print( tt[nm][sep][kk][2]*x^tt[nm][sep][kk][1] );
                    kk := kk + 1;
          end do
        end if;
  end if;

end do:

# End of residual check

#### NOTE: ############################################
# Below was effectively disabled by setting all to 0
# i.e. we decided to leave the expressions as they are
#######################################################

# Table to regularize the negative powers of see the following
reg_tbl:= table( [
 [dg,theta,-1] = 0 ,
 [dg,phi,-1]   = 0 ,
 [DDL,theta,-1] = 0 ,
 [DDL,phi,-1]  = 0 ,
 [RR,x,-2] = 0 ,
 [RR,x,-1] = 0,
 [RR,theta,-2] = 0,
 [RR,theta,-1] = 0,
 [RR,phi,-2] = 0,
 [RR,phi,-1] = 0,
 [P,theta,-1] = 0,
 [P,phi,-1]   = 0,
 [DLAM,x,-2] = 0,
 [DLAM,x,-1] = 0,
 [C1,-1] = 0,
 [HS,-1] = 0,
 [HS3,-1] = 0,
 [EQRPIb2,-1] = 0 ,
 [EQIPIb2,-1] = 0,
 [MK,x,-1] = 0,
 [DB,-1] = 0
] );

#################################################
# Building list of regular terms in expressions:
# Irregular term replacement was disabled
#################################################
for ii from 1 to nops(all_tensors) do
  nm:= all_tensors[ii][1]:
  rank := all_tensors[ii][2]:
  tp  := all_tensors[ii][3]:
  if rank <> 0 then
    for jj from 1 to nops(coord) do
        if tt[nm][sep][coord[jj]][1][1] >= 0 then
           tt[nm][reg][coord[jj]] := tt[nm][all][coord[jj]];
        else
           acm := 0;
           for kk from 1 to nops( tt[nm][sep][coord[jj]] ) do
             ord :=  tt[nm][sep][coord[jj]][kk][1];
             trm :=  tt[nm][sep][coord[jj]][kk][2];
             if ord < 0 then
               if not(assigned(reg_tbl[ [nm,coord[jj],ord] ] ) ) then
                  print("uncomplete reg table for:");
                  print(nm + coord[jj] + ord);
                  error("quiting...");
               else
                  acm := acm + reg_tbl[ [nm,coord[jj],ord] ];
               end if;
             else
               acm := acm + trm*x^ord;
             end if;
           end do;
           tt[nm][reg][coord[jj]] := acm;
        end if
     end do;
  else

        if tt[nm][sep][1][1] >= 0 then
           tt[nm][reg] := tt[nm][all];
        else
           acm := 0;
           for kk from 1 to nops( tt[nm][sep] ) do
             ord :=  tt[nm][sep][kk][1];
             trm :=  tt[nm][sep][kk][2];
             if ord < 0 then
               if not(assigned(reg_tbl[ [nm,ord] ] ) ) then
                  print("uncomplete reg table for:");
                  print(nm + ord);
                  error("quitting...");
               else
                  acm := acm + reg_tbl[ [nm,ord] ];
               end if;
             else
               acm := acm + trm*x^ord;
             end if;
           end do;
           tt[nm][reg]:= acm;
        end if

  end if;
end do:

# Treating 3-Reimmenn tensor seperately

not_all_tensors:=[
         [RR,2,dn]
             ];

R_rest := table([]);
R_ireg := table([]);

for ii from 1 to nops(not_all_tensors) do
  nm:= not_all_tensors[ii][1]:
  rank := not_all_tensors[ii][2]:
  tp  := not_all_tensors[ii][3]:
  if rank <> 0 then
    for jj from 1 to nops(coord) do
        rr_ireg := 0;
        rr_rest := 0;
        for kk from 1 to nops( tt[nm][sep][coord[jj]] ) do
             ord :=  tt[nm][sep][coord[jj]][kk][1];
             trm :=  tt[nm][sep][coord[jj]][kk][2];
             if ord < 0 then
                  rr_ireg:= rr_ireg + trm*x^ord;
             else
                  rr_rest := rr_rest + trm*x^ord;
             end if;
           end do;
        R_ireg[coord[jj]] := rr_ireg;
        R_rest[coord[jj]] := rr_rest;
     end do;
  end if;
end do:

check := true;

#############################################
# Following will print the difference
# between all and reg aka the irregular terms
# in the tensor components
#############################################

if check then
 for ii from 1 to nops(all_tensors) do
  nm:= all_tensors[ii][1]:
  rank := all_tensors[ii][2]:
  tp  := all_tensors[ii][3]:
  if rank <> 0 then
     for jj from 1 to nops(coord) do
       print(nm + coord[jj] );
       ress:=simplify(expand(tt[nm][all][coord[jj]] - tt[nm][reg][coord[jj]]));
       print(expand(ress));
     end do:
  else
      print(nm);
      ress:=simplify(expand(tt[nm][all] - tt[nm][reg]));
      print(expand(ress));
  end if;

 end do;

end if;


text:
Digits := 100:
interface(quiet=true):
read "common-procedures.mpl":
with(orthopoly):

mkdir("TEMPTRIG"):

# - Evaluation scheme :
# case 1 : return x
# case 2 (or Fast) : compute a simple polynomial
# case 3 : do an argument reduction...

########################################################
#  Case 1 : Small arguments
# return x for sine and tan, return 1 for cos
########################################################

xmax_return_x_for_sin := 2^(-26):
xmax_return_1_for_cos_RN := sqrt(2^(-53)):
xmax_return_1_for_cos_RDIR:=2^(-26):
one_rounded_down := evalf(1-ulp(1/2)):

xmax_return_x_for_tan := 2^(-27):

########################################################
# Case 2 : simple polynomial approximation
########################################################

# We want to use the same polynomial in case 2 and 3.
# So see after arg red

#################################################
#   Case 3 : Argument reduction
#################################################

#################################################
# CODY and WAITE  Argument reduction

C := Pi/256:
invC:= nearest(1/C):
reminvC := evalf(1/C - invC):
expC:=ieeedouble(C)[2]:
epsinvC := abs(reminvC*C):

# There are three sets of constants :
#  - split redC into two constants, for small values when we are concerned with absolute error
#  - split redC into three constants, for larger values
#  - split redC into three doubles, for the cases when we need
#     good relative precision on the result and fear cancellation

# Fastest reduction using two-part Cody and Waite (up to |k|=2^22)

bitsCh_0:=34:

# 1/2 <= C/2^(expC+1) <1
Ch:= round(evalf(  C * 2^(bitsCh_0-expC-1))) / (2^(bitsCh_0-expC-1)):
# recompute bitsCh in case we are lucky (and we are for bitsCh_0=32)
bitsCh:=1+log2(op(2,ieeedouble(Ch)[3])) :  # this means the log of the denominator

Cl:=nearest(C - Ch):
# Cody and Waite argument reduction will work for |k|<kmax_cw2
kmax_cw2:=2^(53-bitsCh):

# The constants to move to the .h file
RR_CW2_CH := Ch:
RR_CW2_MCL := -Cl:
XMAX_CODY_WAITE_2 := nearest(kmax_cw2*C):

# The error in this case (we need absolute error)
delta_repr_C_cw2   := abs(C-Ch-Cl):
delta_round_cw2    := kmax_cw2* 1/2 * ulp(Cl) :
delta_cody_waite_2 := kmax_cw2 * delta_repr_C_cw2 + delta_round_cw2:
# This is the delta on y, the reduced argument

#log2(%);

# Slower reduction using three-part Cody and Waite, up to |k|=2^31

bitsCh_0:=23: # 22 or 23
Ch:= round(evalf(  C * 2^(bitsCh_0-expC-1))) / (2^(bitsCh_0-expC-1)):
# recompute bitsCh in case we are lucky
bitsCh:=1+log2(op(2,ieeedouble(Ch)[3])) :  # this means the log of the denominator

r := C-Ch:
Cmed := round(evalf(  r * 2^(2*bitsCh-expC-1))) / (2^(2*bitsCh-expC-1)):
bitsCmed:=1+log2(op(2,ieeedouble(Cmed)[3])) :

Cl:=nearest(C - Ch - Cmed):

kmax_cw3 := 2^min(53-bitsCh, 53-bitsCmed, 31):# Otherwise we have integer overflow

# The constants to move to the .h file
RR_CW3_CH  := Ch:
RR_CW3_CM  := Cmed:
RR_CW3_MCL := -Cl:
XMAX_CODY_WAITE_3 := nearest(kmax_cw3*C):

# The error in this case (we need absolute error)
delta_repr_C_cw3   := abs(C - Ch - Cmed - Cl):
delta_round_cw3    := kmax_cw3 * 1/2 * ulp(Cl) :
delta_cody_waite_3 := kmax_cw3 * delta_repr_C_cw3 + delta_round_cw3:
# This is the delta on y, the reduced argument

#log2(%);

# Third range reduction, using double-double arithmetic, for |k| up to 2^51-1

# This max int value can be produced by DOUBLE2LONGINT
kmax:=2^46-1:
XMAX_DDRR:=nearest(kmax*C):

#in this case we have C stored as 3 doubles
Ch   := nearest(C):
Cmed := nearest(C-Ch):
Cl   := nearest(C-Ch-Cmed):

RR_DD_MCH := -Ch:
RR_DD_MCM := -Cmed:
RR_DD_CL  := Cl:

delta_repr_C := abs(C - Ch - Cmed - Cl):

# and we have only exact Add12 and Mul12  operations. The only place
# with possible rounding errors is:
#       Add22 (pyh, pyl,    (x + kch_h) , (kcm_l - kd*RR_DD_CL),   th, tl) ;
# where (x + kch_h) is exact (Sterbenz) with up to kmax bits of cancellation
# and the error is simply the error in  (kcm_l - kd*RR_DD_CL)
# At the very worst :
delta_round :=
              kmax * 1/2 * ulp(Cl) # for   kd*RR_DD_CL
              + kmax*ulp(Cl)         # for the subtraction
              + 2^(-100) * Pi/512 :    # for the Add22
delta_RR_DD :=  kmax * delta_repr_C + delta_round:

#  the last case, Payne and Hanek reduction, gives a very small delta:
#  red arg is on 9*30 bits, then rounded to a double-double (106 bits)
# This should, of course, be optimized some day
delta_PayneHanek := 2^(-100):

# Finally the max delta on the reduced argument is
delta_ArgRed := max(delta_cody_waite_2, delta_cody_waite_3,
                    delta_RR_DD, delta_PayneHanek):

#print("delta_ArgRed to move to the .gappa file = ",  evalf(delta_ArgRed)):
#log2(delta_ArgRed);

# Now we use the above absolute error when k mod 256 <> 0
# otherwise we need to worry about relative accuracy of the result.

# First, what is the worst case for cancellation ?

emax := ieeedouble(XMAX_DDRR)[2] +1 :
# above emax, we will use Payne and Hanek so we do not worry

(wcn, wce, wceps) := WorstCaseForAdditiveRangeReduction(2,53,-8, emax, C):
wcx := wcn * 2^wce:
wck := round(wcx/C):
wcy := wcx - wck*C:

#log2(wcy);   # y > 2^(-67);

# In these cases we use the double-double range reduction, for |k|<kmax_cw3
# and the relative precision in the worst case is for wcy

delta_round := kmax_cw3 * 1/2 * ulp(Cl)      # for   kd*RR_DD_CL
              + kmax_cw3 * ulp(Cl) :         # for the subtraction

delta_RR_DD :=  kmax_cw3 * delta_repr_C + delta_round:

eps_ArgRed := (1+delta_RR_DD/wcy)*(1+2^(-100)) -1:

#log2(eps_ArgRed);

# In all the other cases we use Payne and Hanek, and eps_ArgRed is
# much smaller, so this is the max.

###########
# Polynomials for do_sine and do_cos, and for the case 2
degreeSin := 8:
degreeCos := 7:

maxepsk := (1+epsinvC)*(1+2^(-53))-1:

ymaxCase3  := evalf(Pi/512 + XMAX_DDRR*maxepsk):
#print("ymaxCase3 to move to the .gappa file = ",  ymaxCase3):

y2maxCase3 := ymaxCase3^2:
# These are the parameters that can be varied to fine-tune performance
#   (they should always be larger than Pi/512
xmaxCosCase2   := Pi/256:
xmaxSinCase2   := Pi/256:

x2maxSinCase2:= xmaxSinCase2^2:
x2maxCosCase2:= xmaxCosCase2^2:

# We had the difficulty here to find minimax polynomials which are good for case 2
# as well as for case 3. A simple solution was to set xmaxCosCase2 =ymaxCase3...
# However we found another answer: in the future we intend to use these polynomials for second
# step, too. Therefore, no minimax, only Taylor.

polySin:=  poly_exact(convert( series(sin(x), x=0, degreeSin+1), polynom)):
polyTs := expand(polySin/x-1):
polyTs2 := subs(x=sqrt(y), polyTs):

polyCos  := poly_exact (convert( series(cos(x), x=0, degreeCos+1), polynom)):
polyTc2 := subs(x=sqrt(y), polyCos - 1):

epsApproxSinCase2 := numapprox[infnorm]((x*polyTs+x -sin(x))/sin(x), x=0..xmaxSinCase2):
epsApproxSinCase3 := numapprox[infnorm]((x*polyTs +x -sin(x))/sin(x), x=0..ymaxCase3):

deltaApproxSinCase2 := numapprox[infnorm]((x*polyTs+x -sin(x)), x=0..xmaxSinCase2):
deltaApproxSinCase3 := numapprox[infnorm]((x*polyTs +x -sin(x)), x=0..ymaxCase3):

deltaApproxCosCase2:= numapprox[infnorm](polyCos -  cos(x), x=0..xmaxCosCase2):
deltaApproxCosCase3:= numapprox[infnorm](polyCos -  cos(x), x=0..ymaxCase3):

#print("deltaApproxSinCase3 to move to the .gappa file = ",  deltaApproxSinCase3):
#print("deltaApproxCosCase3 to move to the .gappa file = ",  deltaApproxCosCase3):

########################## Case 2 for sine  ###########################

# evaluate this polynomial in double. The error on x*x is at most half an ulp
errlist:=errlist_quickphase_horner(degree(polyTs2),0,0,2^(-53), 0):
(eps_rounding_Ts, delta_rounding_Ts, minTs, maxTs):=
	compute_horner_rounding_error(polyTs2,y,x2maxSinCase2, errlist, true):

eps_poly_TsCase2 := numapprox[infnorm]((x*polyTs)/(sin(x)-x) -1, x=0..xmaxSinCase2):
maxeps2 := (1+eps_poly_TsCase2)*(1+eps_rounding_Ts)*(1+2^(-53))-1:

maxepstotalSinCase2 := maxeps2 * numapprox[infnorm](1-x/sin(x), x=0..xmaxSinCase2):
rnconstantSinCase2 := evalf(compute_rn_constant(maxepstotalSinCase2)):

##################################### Case2 cos ###########################

# evaluate this polynomial in double. The error on x*x is at most half an ulp
errlist        := errlist_quickphase_horner(degree(polyTc2),0,0,2**(-53), 0):
(eps_rounding_Tc, delta_rounding_Tc, minTc, maxTc):=
              compute_horner_rounding_error(polyTc2,y,x2maxCosCase2, errlist, true):

# Then we have an Add12 which is exact. The result is greater then cos(xmaxCosCase2):
miny := cos(xmaxCosCase2):
maxepstotalCosCase2 :=  (delta_rounding_Tc + deltaApproxCosCase2) / miny :
#log2(%);
rnconstantCosCase2 := evalf(compute_rn_constant(maxepstotalCosCase2)):

######################## Case2 Tangent #########################
#
# Compute the Taylor series
degreeTanCase2 := 12:
xmaxTanCase2   := 2**(-4):
xminTanCase2   := 2**(-30):

Poly_P := convert(series(tan(sqrt(x))/(x^(3/2))-1/x, x=0, degreeTanCase2*4),polynom):
Poly_cheb := numapprox[chebpade](Poly_P, x=xminTanCase2..xmaxTanCase2^2, [degreeTanCase2/2-2,0]):
polyTanCase2 :=  poly_exact2(expand(x + x^3 * subs(x=x^2, Poly_cheb)), 4):

#polyTanCase2 :=  poly_exact2(convert(series(tan(x), x=0, degreeTanCase2), polynom), 4):

maxepsApproxTanCase2:=numapprox[infnorm](1 - polyTanCase2 / tan(x), x=xminTanCase2..xmaxTanCase2):

maxepsOverXTanCase2 :=numapprox[infnorm]((1 - polyTanCase2 / tan(x))/x, x=xminTanCase2..xmaxTanCase2):
# Now we pass these values to Gappa

filename:="TEMPTRIG/TanCase2.sed":
fd:=fopen(filename, WRITE, TEXT):
  t3h, t3l := hi_lo(coeff(polyTanCase2,x,3)):
  fprintf(fd, " s/_t3h/%1.40e/g\n", t3h):
  fprintf(fd, " s/_t3l/%1.40e/g\n", t3l):
  fprintf(fd, " s/_t5/%1.40e/g\n",  coeff(polyTanCase2,x,5)):
  fprintf(fd, " s/_t7/%1.40e/g\n",  coeff(polyTanCase2,x,7)):
  fprintf(fd, " s/_t9/%1.40e/g\n",  coeff(polyTanCase2,x,9)):
  fprintf(fd, " s/_t11/%1.40e/g\n", coeff(polyTanCase2,x,11)):
  fprintf(fd, " s/_xmax/%1.40e/g\n", xmaxTanCase2):
  fprintf(fd, " s/_maxEpsApproxOverX/%1.40e/g\n", maxepsOverXTanCase2*1.00001):
  fprintf(fd, " s/_maxEpsApprox/%1.40e/g\n", maxepsApproxTanCase2*1.00001):
fclose(fd):

printf("\n\n************ DONE TEMPTRIG/TanCase2.sed ************\n");
printf("Now you should use  \n    sed -f TEMPTRIG/TanCase2.sed trigoTanCase2.gappa | gappa  > /dev/null \n");

maxepstotalTanCase2:=4.59602e-19:  # Cut from Gappa output

log2(maxepstotalTanCase2): # almost 61 bits

###############################################################################
#   Computing errors for Case3 : now we have an error due to arg red

# First DoSinZero. The notations are those of the paper proof

# Approximation error already computed above as epsApproxSinCase3;

#  polynomial evaluation in double, with an error on  y*y of epsilonArgRed
errlist:=errlist_quickphase_horner(degree(polyTs2),0,0,eps_ArgRed, 0):
(epsRoundingTsSinZero, deltaRoundingTsSinZero, minTs, maxTs):=
	compute_horner_rounding_error(polyTs2,y,y2maxCase3, errlist, true):

# just as in the paper proof
maxepsSinZero1 := (1+epsApproxSinCase3)*(1+epsRoundingTsSinZero)*(1+2^(-53))*(1+2^(-53))  - 1:

# just as in the paper proof. For x>0 the absolute values are as given

epstotalSinZero :=  ( (x-sin(x))*maxepsSinZero1 + x*eps_ArgRed + 2^(-53)*x^3/3 )  / sin(x):
maxepstotalSinZero := numapprox[infnorm]( epstotalSinZero , x=2^(-30)..ymaxCase3):

printf("\nMax rel error for DoSinZero is %1.5e, if it's smaller than 2^(-66) (%1.5e) then the proof is OK\n\n", maxepstotalSinZero, 2^(-66)):

##############################SinCosCase3############################
  SinCosSize:= 128: # size f the table

# The Gappa files in TEMPTRIG
for i from 1 to SinCosSize/2 do
    filename:=cat("TEMPTRIG/SinACosA_",i,".sed"):
    fd:=fopen(filename, WRITE, TEXT):

    # The table values
    s:=hi_lo(sin(i*Pi/(2*SinCosSize))):
    c:=hi_lo(cos(i*Pi/(2*SinCosSize))):
    fprintf(fd, " s/_cah/%1.40e/g\n", c[1]):
    fprintf(fd, " s/_cal/%1.40e/g\n", c[2]):
    fprintf(fd, " s/_sah/%1.40e/g\n", s[1]):
    fprintf(fd, " s/_sal/%1.40e/g\n", s[2]):

    # The polynomial coefficients
    fprintf(fd, " s/_s3/%1.40e/g\n", coeff(polySin,x,3)):
    fprintf(fd, " s/_s5/%1.40e/g\n", coeff(polySin,x,5)):
    fprintf(fd, " s/_s7/%1.40e/g\n", coeff(polySin,x,7)):
    fprintf(fd, " s/_s9/%1.40e/g\n", coeff(polySin,x,9)):
    fprintf(fd, " s/_c2/%1.40e/g\n", coeff(polyCos,x,2)):
    fprintf(fd, " s/_c4/%1.40e/g\n", coeff(polyCos,x,4)):
    fprintf(fd, " s/_c6/%1.40e/g\n", coeff(polyCos,x,6)):
    fprintf(fd, " s/_c8/%1.40e/g\n", coeff(polyCos,x,8)):

    # The approximation errors
    fprintf(fd, " s/_ymaxCase3/%1.40e/g\n", ymaxCase3*1.00001):
    fprintf(fd, " s/_delta_ArgRed/%1.40e/g\n", delta_ArgRed*1.00001):
    fprintf(fd, " s/_delta_approx_Sin_Case3/%1.40e/g\n", deltaApproxSinCase3*1.00001):
    fprintf(fd, " s/_delta_approx_Cos_Case3/%1.40e/g\n", deltaApproxCosCase3*1.00001):

    fclose(fd):
od:

printf("************ DONE TEMPTRIG/*.sed ************\n"):

# A shell script to use them
filename:="../gappa/run-trigo-proof.sh":
fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, "#!/bin/sh\n"):
fprintf(fd, "# You probably need to edit the path to the gappa executable\n"):
fprintf(fd, "for file in ../maple/TEMPTRIG/SinACosA*.sed  \n"):
fprintf(fd, "do\n"):
fprintf(fd, "  echo $file:\n"):
fprintf(fd, "  sed -f  $file  trigoSinCosCase3.gappa | ~/gappa/src/gappa\n"):
fprintf(fd, "  echo\n"):
fprintf(fd, "done\n"):
fclose(fd):

printf("************ DONE trigo_test.sh ************\n"):
printf("To run the Gappa proof, you should go to the gappa directory and run\n"):
printf(" sh run-trigo-proof.sh 2> ../maple/TEMPTRIG/Gappa.out\n"):

printf("Then look at maple/TEMPTRIG/Gappa.out. It shouldn't contain 'No proof'.\n This means that everything is OK and the rounding constants in TEMPTRIG/trigo_fast.h are proven upper bounds.\n\n"):

#  This value has been validated by Gappa (using all the previous)
maxepstotalSinCosCase3:=3*2^(-66):
rnconstantSinCosCase3 := evalf(compute_rn_constant(maxepstotalSinCosCase3)):

# The error of sin, the error of cos, then the error of Div22
maxepstotalTanCase3:= 2.1*maxepstotalSinCosCase3:
rnconstantTanCase3 := evalf(compute_rn_constant(maxepstotalTanCase3)):

##############################################
## Compute constants for SCS arg red
oldDigits:=Digits:
Digits:=1000:
# for 2/Pi:
n:=round(2^(30*48)*evalf(2/Pi)):
digitlist:=[]:
for i from 1 to 48 do
  r:=n mod (2^30):
  n:=floor(n/(2^30)):
  hexstring:= convert(convert(r,hex),string):
  digitlist:=[hexstring, op(digitlist)]:
end:
digitlist:

# for 256/Pi:
n:=round(2^(30*47)*evalf(256/Pi)):
digitlist:=[]:
for i from 1 to 48 do
  r:=n mod (2^30):
  n:=floor(n/(2^30)):
  hexstring:= convert(convert(r,hex),string):
  digitlist:=[hexstring, op(digitlist)]:
end:
digitlist:
Digits:=oldDigits:

# an auxiliary output function:
# Outputs the high part of a double, and the double in comment.
# As all these high parts are used in code as
# if(absxhi < XMAX_COS_CASE2)
# we have to remove one LSB to the high part, or, divide var by
# (1+2^(-20))
# Now we have absxhi<highpart(var/(1+2^(-20))
# => absxhi*(1+2^(-20))<var
# => x<var
outputHighPart:=proc(cvarname, var)
  local varMinusLSB:
  Digits:=8:
  varMinusLSB:=var/(1+2^(-20)):
  ("#define " || cvarname || " 0x" || (ieeehexa(varMinusLSB)[1])
    ||  "        /* " || (convert(evalf(varMinusLSB),string)) ||  " */" )
end proc:

# Output:

filename:="TEMPTRIG/trigo_fast.h":
fd:=fopen(filename, WRITE, TEXT):

fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):

fprintf(fd, "\n/*File generated by maple/trigo.pl*/\n"):
fprintf(fd, "\n"):

fprintf(fd,  "%s\n",  outputHighPart("XMAX_RETURN_X_FOR_SIN", xmax_return_x_for_sin) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_SIN_CASE2        ", xmaxSinCase2) ):

fprintf(fd,  "%s\n",  outputHighPart("XMAX_RETURN_1_FOR_COS_RN", xmax_return_1_for_cos_RN) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_RETURN_1_FOR_COS_RDIR", xmax_return_1_for_cos_RDIR) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_COS_CASE2        ", xmaxCosCase2) ):

fprintf(fd,  "%s\n",  outputHighPart("XMAX_RETURN_X_FOR_TAN", xmax_return_x_for_tan) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_TAN_CASE2        ", xmaxTanCase2) ):

fprintf(fd, "\n"):
fprintf(fd, "#define ONE_ROUNDED_DOWN  %1.25e \n", one_rounded_down):
fprintf(fd, "\n"):

fprintf(fd, "#define EPS_SIN_CASE2     %1.25e \n", maxepstotalSinCase2):
fprintf(fd, "#define RN_CST_SIN_CASE2  %1.25f \n", rnconstantSinCase2):
fprintf(fd, "\n"):
fprintf(fd, "#define EPS_COS_CASE2     %1.25e \n", maxepstotalCosCase2):
fprintf(fd, "#define RN_CST_COS_CASE2  %1.25f \n", rnconstantCosCase2):
fprintf(fd, "\n"):
fprintf(fd, "#define EPS_SINCOS_CASE3     %1.25e \n", maxepstotalSinCosCase3):
fprintf(fd, "#define RN_CST_SINCOS_CASE3  %1.25f \n", rnconstantSinCosCase3):
fprintf(fd, "\n"):
fprintf(fd, "#define EPS_TAN_CASE2     %1.25e \n", maxepstotalTanCase2):
fprintf(fd, "#define EPS_TAN_CASE3     %1.25e \n", maxepstotalTanCase3):
fprintf(fd, "#define RN_CST_TAN_CASE3  %1.25f \n", rnconstantTanCase3):

fprintf(fd, "\n"):

fprintf(fd, "#define INV_PIO256        %1.25f \n", 1/C):
fprintf(fd, "\n"):

fprintf(fd,  "%s\n",  outputHighPart("XMAX_CODY_WAITE_2", XMAX_CODY_WAITE_2) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_CODY_WAITE_3", XMAX_CODY_WAITE_3) ):
fprintf(fd,  "%s\n",  outputHighPart("XMAX_DDRR        ", XMAX_DDRR) ):
#fprintf(fd,  "%s\n",  outputHighPart("", ) ):
fprintf(fd, "\n"):

fprintf(fd, "#define RR_CW2_CH  %1.25e\n", RR_CW2_CH):
fprintf(fd, "#define RR_CW2_MCL %1.25e\n", RR_CW2_MCL):
fprintf(fd, "\n"):

fprintf(fd, "#define RR_CW3_CH  %1.25e\n", RR_CW3_CH):
fprintf(fd, "#define RR_CW3_CM  %1.25e\n", RR_CW3_CM):
fprintf(fd, "#define RR_CW3_MCL %1.25e\n", RR_CW3_MCL):
fprintf(fd, "\n"):

fprintf(fd, "#define RR_DD_MCH  %1.25e\n", RR_DD_MCH):
fprintf(fd, "#define RR_DD_MCM  %1.25e\n", RR_DD_MCM):
fprintf(fd, "#define RR_DD_CL   %1.25e\n", RR_DD_CL):
fprintf(fd, "\n"):

fprintf(fd, "\n"):

  fprintf(fd, "\n\n"):

  # The 256/Pi SCS array
  fprintf(fd, "static const int digits_256_over_pi[] = \n{"):
  for i from 0 to 11 do
    for j from 1 to 4 do
      fprintf(fd, " 0x%s,  \t",digitlist[4*i+j]):
    end:
    fprintf(fd, "\n "):
  end:
  fprintf(fd, "};\n\n"):

  # The Pi/256 SCS constant
  fprintf(fd, "static const scs Pio256=\n"):
  WriteSCS(fd, evalf(C)):
  fprintf(fd, ";\n#define Pio256_ptr  (scs_ptr)(& Pio256)\n\n"):

fprintf(fd,"#ifdef WORDS_BIGENDIAN\n"):
for isbig from 1 to 0 by -1 do

  if(isbig=0) then
    fprintf(fd,"#else\n"):
  fi:

  # The sine polynomial

  fprintf(fd, "static db_number const s3 = "):
  printendian(fd, coeff(polySin,x,3), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const s5 = "):
  printendian(fd, coeff(polySin,x,5), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const s7 = "):
  printendian(fd, coeff(polySin,x,7), isbig):
  fprintf(fd, ";\n\n"):

  # the cos polynomial

  fprintf(fd, "static db_number const c2 = "):
  printendian(fd, coeff(polyCos,x,2), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const c4 = "):
  printendian(fd, coeff(polyCos,x,4), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const c6 = "):
  printendian(fd, coeff(polyCos,x,6), isbig):
  fprintf(fd, ";\n\n"):

  # the tan polynomial

  t3h, t3l := hi_lo(coeff(polyTanCase2,x,3)):
  fprintf(fd, "static db_number const t3h = "):
  printendian(fd, t3h, isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const t3l = "):
  printendian(fd, t3l, isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const t5 = "):
  printendian(fd, coeff(polyTanCase2,x,5), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const t7 = "):
  printendian(fd, coeff(polyTanCase2,x,7), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const t9 = "):
  printendian(fd, coeff(polyTanCase2,x,9), isbig):
  fprintf(fd, ";\n"):
  fprintf(fd, "static db_number const t11 = "):
  printendian(fd, coeff(polyTanCase2,x,11), isbig):
  fprintf(fd, ";\n\n"):

  # The sincos table

  fprintf(fd, "\n/*  sine and cos of kPi/256 in double-double */\n"):
  fprintf(fd, "static db_number const sincosTable[%d] =\n{\n",  4*(SinCosSize/2+1)):
  for i from 0 to SinCosSize/2 do
      s:=hi_lo(sin(i*Pi/(2*SinCosSize))):
      c:=hi_lo(cos(i*Pi/(2*SinCosSize))):
      printendian(fd,s[1],isbig):
      fprintf(fd," ,\n"):
      printendian(fd,s[2],isbig):
      fprintf(fd," ,\n"):
      printendian(fd,c[1],isbig):
      fprintf(fd," ,\n"):
      printendian(fd,c[2],isbig):
    if i<SinCosSize-1 then fprintf(fd," ,\n"): fi:
  od:
  fprintf(fd, "\n};\n\n"):

od:
fprintf(fd,"#endif /* WORDS_BIGENDIAN */\n\n\n"):

fclose(fd):

printf("\n\n************ DONE TEMPTRIG/trigo_fast.h ************\n Copy it to the crlibm source directory.\n\n");

#################################################
# Stuff for the accurate phase,
#
# Should very bad cases be found in the future, only the degrees below
# should  be increased, then this script rerun. Everything should
# compile OK with the newly generated trigo_accurate.h

printf("\n-------------------------------------------------------\n"):
printf("--------------------Accurate phase---------------------\n"):
printf("-------------------------------------------------------\n"):

xminSCS   := 0:
xmaxSCS   := Pi/4:

#-----------------Sine-----------------------
degreeSinSCS := 26:
Poly_P := convert(series(sin(sqrt(x))/(x^(3/2))-1/x, x=0, degreeSinSCS*4),polynom):
Poly_cheb := numapprox[chebpade](Poly_P, x=xminSCS..xmaxSCS^2, [degreeSinSCS/2-2,0]):
polySinSCS :=  poly_exact_SCS(expand(x + x^3 * subs(x=x^2, Poly_cheb))):

DEGREE_SIN_SCS := degree(polySinSCS):
maxepsApproxSinSCS:=numapprox[infnorm](1 - polySinSCS/sin(x), x=xminSCS..xmaxSCS):
printf("The sine polynomial for the second phase (degree %d) is accurate to %f bits on O..Pi/4\n", DEGREE_SIN_SCS, -log2(maxepsApproxSinSCS)):
# For the proof of the very bad cases, near zero
maxepsApproxSinSCSNearZero:=numapprox[infnorm](1 - polySinSCS/sin(x), x=0..2^(-17)):
printf("  ... and  accurate to %f bits on O..2^(-17)\n", -log2(maxepsApproxSinSCSNearZero)):

#-----------------Cos-----------------------
degreeCosSCS := 28:
Poly_P := convert(series((cos(sqrt(x)) -1)/x, x=0, degreeCosSCS*4),polynom):
Poly_cheb := numapprox[chebpade](Poly_P, x=xminSCS..xmaxSCS^2, [degreeCosSCS/2-2,0]):
polyCosSCS :=  poly_exact_SCS(expand(1 + x^2 * subs(x=x^2, Poly_cheb))):
DEGREE_COS_SCS := degree(polyCosSCS):

maxepsApproxCosSCS:=numapprox[infnorm](1 - polyCosSCS/cos(x), x=xminSCS..xmaxSCS):
printf("The cos polynomial for the second phase (degree %d) is accurate to %f bits on O..Pi/4\n", degree(polyCosSCS), -log2(maxepsApproxCosSCS)):
# For the proof of the very bad cases, near zero
maxepsApproxCosSCSNearZero:=numapprox[infnorm](1 - polyCosSCS/cos(x), x=0..2^(-18)):
printf("  ... and  accurate to %f bits on O..2^(-18)\n", -log2(maxepsApproxCosSCSNearZero)):

#-----------------Tan-----------------------
degreeTanSCS := 70:
Poly_P := convert(series(tan(sqrt(x))/(x^(3/2))-1/x, x=0, degreeTanSCS*4),polynom):
Poly_cheb := numapprox[chebpade](Poly_P, x=xminSCS..xmaxSCS^2, [degreeTanSCS/2-2,0]):
polyTanSCS :=  poly_exact_SCS(expand(x + x^3 * subs(x=x^2, Poly_cheb))):

DEGREE_TAN_SCS := degree(polyTanSCS):
maxepsApproxTanSCS:=numapprox[infnorm](1 - polyTanSCS/tan(x), x=xminSCS..xmaxSCS):
printf("The tan polynomial for the second phase (degree %d) is accurate to %f bits on O..Pi/4\n", DEGREE_TAN_SCS, -log2(maxepsApproxTanSCS)):
# For the proof of the very bad cases, near zero
maxepsApproxTanSCSNearZero:=numapprox[infnorm](1 - polyTanSCS/tan(x), x=0..2^(-17)):
printf("  ... and  accurate to %f bits on O..2^(-17)\n", -log2(maxepsApproxTanSCSNearZero)):

filename:="TEMPTRIG/trigo_accurate.h":
fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, "/*File generated by maple/trigo.mpl*/\n\n"):
fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n\n"):
fprintf(fd, "#define DEGREE_SIN_SCS  %d \n", DEGREE_SIN_SCS):
fprintf(fd, "#define DEGREE_COS_SCS  %d \n", DEGREE_COS_SCS):
fprintf(fd, "#define DEGREE_TAN_SCS  %d \n", DEGREE_TAN_SCS):
fprintf(fd, "\n"):
fprintf(fd, "#define  sin_scs_poly_ptr (scs_ptr)&sin_scs_poly \n"):
fprintf(fd, "#define  cos_scs_poly_ptr (scs_ptr)&cos_scs_poly \n"):
fprintf(fd, "#define  tan_scs_poly_ptr (scs_ptr)&tan_scs_poly \n"):
fprintf(fd, "\n"):
Write_SCS_poly(fd, sin_scs_poly, polySinSCS):
fprintf(fd, "\n\n"):
Write_SCS_poly(fd, cos_scs_poly, polyCosSCS):
fprintf(fd, "\n\n"):
Write_SCS_poly(fd, tan_scs_poly, polyTanSCS):
fprintf(fd, "\n"):
fclose(fd):

printf("\n\n************ DONE TEMPTRIG/trigo_accurate.h ************\n Copy it to the crlibm source directory.\n\n");

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               2020 Susi Lehtola

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_c_p86_params *params;

  assert(p->params != NULL);
  params = (gga_c_p86_params * )(p->params);
*)

$define lda_c_pz_params
$include "lda_c_pz.mpl"

(* Equation (4) *)
p86_DD := z  -> sqrt(opz_pow_n(z,5/3) + opz_pow_n(-z,5/3))/sqrt(2):

(* Equation (6) *)
p86_CC := rs ->
  + params_a_aa
  + (params_a_bb + params_a_malpha*rs + params_a_mbeta*rs^2)/(1 + params_a_mgamma*rs + params_a_mdelta*rs^2 + 1.0e4*params_a_mbeta*rs^3):
p86_CCinf := params_a_aa + params_a_bb:

(* Equation (9) *)
p86_x1   := (rs, xt) -> xt/sqrt(rs/RS_FACTOR):
p86_mPhi := (rs, xt) -> params_a_ftilde*(p86_CCinf/p86_CC(rs))*p86_x1(rs, xt):

(* Equation (8) *)
p86_H := (rs, z, xt) -> p86_x1(rs, xt)^2*exp(-p86_mPhi(rs, xt))*p86_CC(rs)/p86_DD(z):

f_p86 := (rs, z, xt, xs0, xs1) ->
  f_pz(rs, z) + p86_H(rs, z, xt):

f := (rs, z, xt, xs0, xs1) ->
  f_p86(rs, z, xt, xs0, xs1):


text:
(* type: mgga_exc *)
(* prefix:
  mgga_c_ltapw_params *params;

  assert(p->params != NULL);
  params = (mgga_c_ltapw_params * )(p->params);
*)

$define lda_c_pw_params
$include "lda_c_pw.mpl"

(* kinetic energy density to electron density *)
nt_tau := t -> (t/K_FACTOR_C)^(3*params_a_ltafrac/5):

(* effective density *)
n_eff_s := (rs, z, t) -> n_spin(rs, z) * nt_tau(t):
n_eff   := (rs, z, ts0, ts1) -> n_eff_s(rs, z, ts0) + n_eff_s(rs, -z, ts1):

(* recompute rs and zeta *)
eff_rs := (rs, z, ts0, ts1) -> r_ws(n_eff(rs, z, ts0, ts1)):
eff_z  := (rs, z, ts0, ts1) ->
  (n_eff_s(rs, z, ts0) - n_eff_s(rs, -z, ts1))/n_eff(rs, z, ts0, ts1):

f := (rs, z, xt, xs0, xs1, us0, us1, ts0, ts1) ->
  f_pw(eff_rs(rs, z, ts0, ts1), eff_z(rs, z, ts0, ts1)):


text:
`is_element/colour_profiles` := (C::set) -> (n::nonnegint) -> proc(c)
 type(c,list) and nops(c) = n and {op(c)} minus C = {};
end:

`is_equal/colour_profiles` := (C::set) -> (n::nonnegint) -> (c,b) -> evalb(c = b);

`list_elements/colour_profiles` := (C::set) -> proc(n::nonnegint)
 local L,i,x,c;
 L := [[]];
 for i from 1 to n do
  L := [seq(seq([op(c),x],x in C),c in L)];
 od:
 return L;
end:

`count_elements/colour_profiles` := (C::set) -> (n::nonnegint) -> nops(C)^n;

`random_element/colour_profiles` := (C::set) -> (n::nonnegint) -> proc()
 local i;
 [seq(random_element_of(C),i=1..n)];
end:

`o/colour_profiles` := proc() map(op,[args]); end:

`gamma/colour_profiles` := (C::set) -> (L::list(nonnegint)) -> proc(c,d)
 map(op,d);
end:

`circ/colour_profiles` := (C::set) -> (i,m,n) -> proc(c,b)
 [op(1..i-1,c),op(b),op(i+1..m,c)];
end:


text:
program nwhilest;  {sample27}
var i, j, k    : integer;
begin
i := 1;
while i <10 do begin
j := 1;
while j < 10 do begin
k := 1;
while k < 10 do begin
if (k div 2)*2=k then begin k := k+1 end
else begin k:=k+1end
end;
j := j+1;
end;
i:=i+1
end;
writeln                    ('All End')end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
var $g ptr
var $h i32
func $multiwayfunc ( var %p ptr) ptr {
 multiway (dread ptr %p) @labdft {
   (dread ptr $g): goto @lab1
   (add ptr (dread i32 $g, dread i32 $h)): goto @lab0
   (conststr ptr "world"): goto @lab9 }
@lab0
 return (dread ptr $g)
@labdft
 return (conststr ptr "hello")
@lab9
 return (add ptr (conststr ptr "foo", dread i32 $h))
@lab1
 return (addrof ptr $h) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
################################
# Derivation of BSSN equations
# Using GR-Tensor
# See bssn_spher.mw for interactive
# version of this script
################################
grtw();
grload(g,"metrics/spher_3.mpl");
grdisplay(g(dn,dn));
grcalc(g(up,up));
grdisplay(g(up,up));
grdef(`g_f{a b} := kdelta{a $x}*kdelta{b $x}*1 + kdelta{a $theta}*kdelta{b $theta}*x^2 + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2`); grcalc(g_f(dn,dn)); grdisplay(g_f(dn,dn));
grdef(`g_f_inv{^a ^b} := kdelta{^a $x}*kdelta{^b $x}*1 + kdelta{^a $theta}*kdelta{^b $theta}*1/x^2 + kdelta{^a $phi}*kdelta{^b $phi}*1/(x^2*(sin(theta))^2)`); grcalc(g_f_inv(up,up)); grdisplay(g_f_inv(up,up));
grdef(`ChL_f{d a b} := 1/2*( g_f{b d ,a} + g_f{a d,b} - g_f{a b,d})`);
grdef(`MCh_f{^c a b} := g_f_inv{^c ^d} * ChL_f{d a b}`):
grcalc(MCh_f(up,dn,dn)); grdisplay(MCh_f(up,dn,dn));
grdef(`ChL{d a b} := 1/2*( g{b d ,a} + g{a d,b} - g{a b,d})`):
grdef(`MCh{^c a b} := g{^c ^d} * ChL{d a b}`):
grcalc(MCh(up,dn,dn));
gralter(MCh(up,dn,dn),simplify,expand);
grdisplay(MCh(up,dn,dn));
grdef(`TG{^i} := g{^j ^k}*MCh{^i j k}`);
grcalc(TG(up));
grdisplay(TG(up));
grdef(`TG2{^i} := g{^j ^k}*(MCh{^i j k}-MCh_f{^i j k})`); grcalc(TG2(up)); gralter(TG2(up),simplify,expand); grdisplay(TG2(up));
grdef(`Lam{^a} := [Lamx(t,x),0,0]`);
grcalc(Lam(up)); grdisplay(Lam(up));
grdef(`GA{^i} := Lam{^i} + g{^j ^k}*MCh_f{^i j k}`);  grcalc(GA(up)); grdisplay(GA(up));
grdef(`SS{i j} := 2*g{^l ^m}*MCh{^k l (i}*ChL{j) k m} + g{^l ^m}*MCh{^k i m}*ChL{k l j}`);
grcalc(SS(dn,dn));
gralter(SS(dn,dn),simplify);
grdisplay(SS(dn,dn));
grdef(`RB{i j} := -1/2*g{^l ^m}*g{i j ,l ,m} + g{k (i}*GA{^k ,j)} + GA{^k}*ChL{(i j) k} + SS{i j}`);
grcalc(RB(dn,dn));
gralter(RB(dn,dn),simplify,expand);
grdisplay(RB(dn,dn));
grdef(`CF{} := phi(t,x)`);
grcalc(CF);
grdef(`DCF{i} := CF{;i}`);
grcalc(DCF(dn));
grdisplay(DCF(dn));
grdef(`Rphi{i j} := -2*DCF{j ;i} - 2*g{i j}*DCF{^k ;k} + 4*DCF{i}*DCF{j} - 4*g{i j}*DCF{^k}*DCF{k}`);
grcalc(Rphi(dn,dn));
gralter(Rphi(dn,dn),expand,simplify,expand);
grdisplay(Rphi(dn,dn));
grcalc(R(dn,dn)); grdisplay(R(dn,dn));
grdef(`RBtest{i j} := -1/2*g{^l ^m}*g{i j ,l ,m} + g{k (i}*TG{^k ,j)} + TG{^k}*ChL{(i j) k} + SS{i j}`);
grdef(`Rresid{i j} := RBtest{i j} - R{i j}`);
grcalc(Rresid(dn,dn));
gralter(Rresid(dn,dn),simplify);
grdisplay(Rresid(dn,dn));
grdef(`Sh{^a}:=[beta(t,x),0,0]`);
grcalc(Sh(up));
grdisplay(Sh(up));
grdef(`Q{i j} := Sh{^k}*g{i j ,k} + g{i k}*Sh{^k ,j} + g{k j}*Sh{^k ,i} -2/3*g{i j}*vee*divbeta(t,x)`);
grcalc(Q(dn,dn));
grdisplay(Q(dn,dn));
grdef(`A{a b} := kdelta{a $x}*kdelta{b $x}*Axx(t,x)/1 + kdelta{a $theta}*kdelta{b $theta}*x^2*Athth(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*Athth(t,x)`);
grcalc(A(dn,dn));grdisplay(A(dn,dn));
grdef(`P{i j} := Sh{^k}*A{i j,k} + A{i k}*Sh{^k ,j} + A{k j}*Sh{^k ,i} - 2/3*A{i j}*vee*divbeta(t,x)`):
grcalc(P(dn,dn));
grdisplay(P(dn,dn));
grdef(`Extc{} := K(t,x)`);
grcalc(Extc);
grdef(`Lapse{} := alpha(t,x)`);
grcalc(Lapse);
grdef(`S{i} :=[Sx(t,x),0,0]`);
grcalc(S(dn));

grdef(`DBS{} := divbeta(t,x)`);
grdef(`DGA{^i} := g{^l ^j}*Sh{^i ,l ,j} - 2*A{^i ^j}*Lapse{,j} + 2*Lapse*(MCh{^i j k}*A{^k ^j} - 2/3*g{^i ^j}*Extc{,j} - g{^i ^j}*S{j} + 6*A{^i ^j}*DCF{j}) + Sh{^j}*GA{^i ,j} - GA{^j}*Sh{^i ,j} +vee/3*(2*GA{^i}*divbeta(t,x)+g{^i ^l}*DBS{,l})`);
grcalc(DGA(up));

gralter(DGA(up),simplify,expand,trigsin,factor,expand);
grdisplay(DGA(up));
grdef(`met{a b} := exp(4*CF)*(kdelta{a $x}*kdelta{b $x}*A(t,x) + kdelta{a $theta}*kdelta{b $theta}*B(t,x)*x^2 + kdelta{a $phi}*kdelta{b $phi}*B(t,x)*x^2*sin(theta)^2)`);
grcalc(met(dn,dn));
grdisplay(met(dn,dn));

grdef(`met_s{a b} := kdelta{a $x}*kdelta{b $x}*metxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*metthetatheta(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*metthetatheta(t,x)`);
grcalc(met_s(dn,dn));grdisplay(met_s(dn,dn));
grdef(`invmet_s{^a ^b}:=kdelta{^a $x}*kdelta{^b $x}*1/metxx(t,x) + kdelta{^a $theta}*kdelta{^b $theta}*1/(x^2*metthetatheta(t,x)) + kdelta{^a $phi}*kdelta{^b $phi}*1/(x^2*(sin(theta))^2*metthetatheta(t,x))`);
grcalc(invmet_s(up,up));
grdisplay(invmet_s(up,up));
grdef(`rChL{d a b} := 1/2*( met_s{b d ,a} + met_s{a d,b} - met_s{a b,d})`);
grdef(`rMCh{^c a b} := invmet_s{^c ^d} * rChL{d a b}`);
grdef(`DDL{b c} :=   Lapse{,c ,b} - rMCh{^d b c}*Lapse{,d}`);
grcalc(DDL(dn,dn));
grdisplay(DDL(dn,dn));
grdef(`rR{u p} := rMCh{^v u p ,v} - rMCh{^v v p ,u} + rMCh{^a u p}*rMCh{^v a v} - rMCh{^a v p}*rMCh{^v a u}`);
grdef(`RHO{} := rho(t,x)`);
grcalc(RHO);
grdef(`TraceS{} := TS(t,x)`);
grcalc(TraceS);
grdef(`DDL_s{a b} := kdelta{a $x}*kdelta{b $x}*DDLxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*DDLthetatheta(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*DDLthetatheta(t,x)*(sin(theta))^2`);
grcalc(DDL_s(dn,dn)); grdisplay(DDL_s(dn,dn));
grdef(`TD{} := invmet_s{^i ^j} * DDL_s{i j}`);
grcalc(TD);
grdisplay(TD);
grdef(`DDLTF{i j} := DDL_s{i j} - 1/3*TD*met_s{i j}`);
grcalc(DDLTF(dn,dn));
grdisplay(DDLTF(dn,dn));

grdef(`DTK{} := -TD +Lapse*(A{i j}*A{^i ^j} +1/3*Extc^2) +1/2*Lapse*(RHO+TraceS) +Sh{^i}*Extc{,i}`);
grcalc(DTK);
gralter(DTK,simplify,expand);
grdisplay(DTK);
grdef(`JS{a b} := kdelta{a $x}*kdelta{b $x}*JSxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*JSthth(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*JSthth(t,x)`);
grcalc(JS(dn,dn));
grdisplay(JS(dn,dn));
grdef(`TJS := invmet_s{^i ^j}*JS{i j}`);
grcalc(TJS);
gralter(TJS,simplify,expand);
grdisplay(TJS);
grdef(`U{i j} := (Extc*A{i j} - 2*A{i l}*A{^l j})`);
grcalc(U(dn,dn));
gralter(U(dn,dn),simplify,expand);
grdisplay(U(dn,dn));
grdef(`dg{i j} := -2*Lapse*A{i j} + Q{i j}`);
grcalc(dg(dn,dn));
gralter(dg(dn,dn),simplify,factor);
grdisplay(dg(dn,dn));
grdef(`dgu{^i ^j} := -g{^i ^l}*dg{l m}*g{^m ^j}`); grcalc(dgu(up,up)); grdisplay(dgu(up,up));
grdef(`deltaDGA{^i} := -dgu{^j ^k}*MCh_f{^i j k}`); grcalc(deltaDGA(up)); grdisplay(deltaDGA(up));
grdef(`DLAM{^i} := DGA{^i} + deltaDGA{^i}`); grcalc(DLAM(up)); gralter(DLAM(up),simplify,expand); grdisplay(DLAM(up));
grdef(`DPH{} := -1/6*Lapse*Extc + Sh{^i}*CF{,i} + 1/6*divbeta(t,x)`);
grcalc(DPH); grdisplay(DPH);
grdef(`TRB{} := g{^i ^j}*RB{i j}`);
grcalc(TRB);
gralter(TRB,simplify);
grdisplay(TRB);
grdef(`PSI{} := psi(t,x)`);
grcalc(PSI);
grdef(`C1 {} := -1/8*TRB`); grcalc(C1); gralter(C1,simplify,expand,factor,expand); grdisplay(C1);
grdef(`C5{} := 1/8*A{i j}*A{^i ^j} -1/12*Extc^2 + 1/4*rho(t,x)`);
grcalc(C5); gralter(C5,simplify,expand); grdisplay(C5);
grdef(`HS3{} := g{^i ^j}*PSI{;j ;i} + C5psi(t,x) + C1s(t,x)*PSI`);
grdef(`HS{} := g{^i ^j}*PSI{;j ;i} + C5s(t,x)*PSI^5 + C1s(t,x)*PSI`); grdef(`HSR{} := g{^i ^j}*PSI{;i ;j} - 1/8*TRB*PSI + 1/8*A{i j}*A{^i ^j}*PSI^5 -1/12*PSI^5*Extc^2 + 1/4*rho(t,x)*PSI^5`); grcalc(HSR);gralter(HSR,simplify,expand,factor,expand);grdisplay(HSR);
grcalc(HS);
gralter(HS,simplify,expand,factor,expand);
grdisplay(HS);
grcalc(HS3); grdisplay(HS3);gralter(HS3,expand);
grdef(`RR{i j} := RB{i j} + Rphi{i j}`);
grcalc(RR(dn,dn));
gralter(RR(dn,dn),simplify,expand,factor,expand);
grdisplay(RR(dn,dn));
grdef(`RR_s{a b} := kdelta{a $x}*kdelta{b $x}*RRxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*RRthetatheta(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*RRthetatheta(t,x)`);
grcalc(RR_s(dn,dn)); grdisplay(RR_s(dn,dn));
grdef(`TR{} := invmet_s{^i ^j}*RR_s{i j}`);
grcalc(TR);
grdisplay(TR);

grdef(`RRTF{i j} := RR_s{i j} - 1/3*met_s{i j}*TR`);
grcalc(RRTF(dn,dn));
grdisplay(RRTF(dn,dn));
grdef(`P_s{a b} := kdelta{a $x}*kdelta{b $x}*Pxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*Pthetatheta(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*Pthetatheta(t,x)`);
grcalc(P_s(dn,dn)); grdisplay(P_s(dn,dn));
grdef(`U_s{a b} := kdelta{a $x}*kdelta{b $x}*Uxx(t,x) + kdelta{a $theta}*kdelta{b $theta}*x^2*Uthetatheta(t,x) + kdelta{a $phi}*kdelta{b $phi}*x^2*(sin(theta))^2*Uthetatheta(t,x)`);
grcalc(U_s(dn,dn));
grdisplay(U_s(dn,dn));
grdef(`DAA{i j} := em4phi(t,x)*(-DDLTF{i j} +  Lapse*(RRTF{i j} - ( JS{i j} - 1/3*met_s{i j}*TJS)    )) + Lapse*U_s{i j} + P_s{i j}`);
grcalc(DAA(dn,dn));
gralter(DAA(dn,dn),simplify,expand);
grdisplay(DAA(dn,dn));
grdef(`EMPH{} := exp(-4*CF)`);grcalc(EMPH);grdisplay(EMPH);
grdef(`MKTR1{^j ^i} := PSI^6*A{^j ^i}`);
grdef(`MK{^i} := MKTR1{^j ^i ;j} - 2/3*PSI^6*g{^i ^j}*Extc{;j} - PSI^6*g{^i ^j}*S{j}`); grcalc(MK(up)); gralter(MK(up),simplify,expand); grdisplay(MK(up));
grdef(`DB:= Sh{^i ;i}`); grcalc(DB); grdisplay(DB);gralter(DB,simplify,expand); grdisplay(DB);