EVAlgebra script
#*******************************************************************
#                           EVA 1.2  rev 7                         *      
#*******************************************************************
#
# contact at beyhgm@gmail.com
# EVA site beyh.free.fr
#
# installing EVA :
# download Eigenmath from eigenmath.sourceforge.net
# run Eigenmath
# press "Edit Script" and copy/paste this script
# press "Run Script"
#
# to do : optimize product inp
#         check domain validity for math functions  
#         define asin acos atan                      
#         isVersor, isBlade, versor factorization
#         enhance meet, join
#         paravectors specific functions
#         add more examples
#
# eigenmath control instructions :
# do( expression1, expression2,..., last_expression)   return last_expression
# test( predicate1, do(...), predicate2, do(...),..., do(...))
# for(k,1,n, do(...))
# if only one expression, do(...) not necessary.
#
# last revision : 15/02/2012
#
#(c) Copyright, Bernard Eyheramendy,
#  author of this script. Permission is given for free use,
#  including copying, and modifying, provided only that the
#  original program (EVA or EVAlgebra) is referenced.
#  This file is provided AS IS with NO WARRANTY OF ANY KIND.
#**********************************************************************
# global variables : Clmax=p+q, Cldim=2^Clmax, Clsig = signature, Clindex = basis vectors index
# e0, e1, ... basis vectors
# Gtab Otab Itab multiplication tables for   gp outp  inp
# j  oriented volume ( j=e1..k   with k=p+q )
trace=0
svtty=1
fixnum=4     # number of decimals (rnd(u) round results to fixnum else rnd(u,n))
maxTaylor=12
INF=-log(0)
# Internal functions  Base(u)  Index(u)  consGtab(u)  consOtab(u)  consItab(u) find(u)  sign(u)   prod(u)
# Gtable(u) Otable(u) Itable(u)
# local variables : names ending with 00
 
# controls *************************************************************
                                             
# test versor(u),  blade(u)                                  # versor = a1 a2 ... ak product of k vectors
versor(u)=do(R00=1,u00=rnd(gp(invol(u),inverse(u))),            # use : test(versor(u)) --> 1 if versor else 0
test(not(grade0(u00)),R00=0,                                 # grade(u' 1/u) = 0 ?
    u00==rnd(gp(inverse(u),invol(u)),preservGrd(u),R00=0),R00))# u' 1/u = 1/u u'   ?
preservGrd(u)=for(k,2,Clmax+1,                              
    test(grade1(gp(gp(invol(u),unit(Cldim)[k]),rev(u))),1,R00=0))    # grade(u' ei u~) = 1 ?
 
blade(u)=do(G00=0,for(k,1,Clmax+1,test(grade(u,k-1)==0,1,G00=G00+1)), # use : test(blade(u)) --> 1 if blade else 0
    test(G00=1,R00=1,R00=0),R00*test(versor(u)))                     # blade = versor with one grade only
 
grade0(u)=test(u==e0*u[1],1,0)  # use : test(grade0(u)) --> 1 if grade(u)=0
grade1(u)=do(u00=zero(Cldim),for(k,2,Clmax+1,u00[k]=u[k]),test(u==u00,1,0))
grade2(u)=do(u00=zero(Cldim),Bv=(0,0,4,7,11,16),for(k,Clmax+2,Bv[Clmax+1],u00[k]=u[k]),test(u==u00,1,0))
oddGrade(u) =test(u==1/2(u-invol(u)))                        #  return 1 if true        
evenGrade(u)=test(u==1/2(u+invol(u)))                        #       "
isVersor(u)=test(versor(u),1,             errmsg(msg6,u))    #
isInteger(u)=test(isinteger(u),1,         errmsg(msg5,u))    #  7
isNumber(u)=test(number(u),1,             errmsg(msg4,u))    #  23.9
isScalar(u)=test(u==grade(u,0), 1,u=0,1   errmsg(msg3,u))    # (3,0,0,....,0)
isVector(u)=test(u==grade(u,1), 1,u=0,1   errmsg(msg2,u))    # (0,2.......,0)
isMvector(u)=test(dim(u)==Cldim,1,u=0,1,  errmsg(msg1,u))    # (1,0,5,..-4,.)
isGradeover(n)=do(isNumber(n),
                 test(n>Clmax,           errmsg(msg7,u),
                      n<0,               errmsg(msg7,u)))   #  n < 0 or > Clmax
isGT1(u)=test(u[1] >= 1, 1,               errmsg(msg8,u))    #  u[1]  >=  1
isPos(u)=test(u[1] >= 0, 1,               errmsg(msg9,u))    #  u[1]  >=  0
isZero(u)=test(u[1] = 0,                  errmsg(msg10,u),1) #  u[1]   =  0
isRangeLT1(u)= test(abs(u[1])< 1, 1,      errmsg(msg11,u))   # |u[1]|  <  1
isRangeLE1(u)= test(abs(u[1])<=1, 1,      errmsg(msg12,u))   # |u[1]| <=  1
 
# error messages
msg1="not a multivector : "
msg2="not a vector      : "
msg3="not a scalar      : "
msg4="not a number      : "
msg5="not integer       : "
msg6="not a versor      : "
msg7=" grade from 0 to p+q  : "
# error if        acos1(|u| > 1)  asin1(|u| > 1)  acosh1(u < 1)  atanh1(|u| > 1)  log1(u < 0)  sqrt1(u < 0)
msg8=" u[1] must be >= 1 : "
msg9=" u[1] must be > 0    : "
msg10=" u[1] must be different from 0 : "
msg11=" |u[1]| must be < 1 : "
msg12=" |u[1]| must be <= 1 : "
errmsg(msg,object)=do(print(msg,object),stop())               # print msg and stop
 
# end controls  *****************************************************
"   EVA version :  1.27 "
" exemples :"
" Cl(3)   --> define space Cl3 : basis = e0,e1,e2,e3,e12,e23,e13,e123"
" e0 is a scalar "
" e1, e2, e3 are vectors, e12, e23, e13 are bivectors"
" e123 is a pseudoscalar "
""
" define vector a=(3,2,-1), b=(3,0,-5)     :   a=3e1+2e2-e3   b=3e1-5e3"
" geometric product  a b                   :   gp(a,b)"
" inner product      a.b                   :   inp(a,b)"
" outer product      a^b                   :   outp(a,b)"
" define multivector B=(3,1,-5,0,1,1,0,0)  :   B =       3e0+e1-5e2+e12+2e13"
" grade projection : scalar                :   grade(B,0) = e0"
"                    vector                :   grade(B,1) = e1-5e2"
"                    bivector              :   grade(B,2) = e12+2e13"
"                    pseudoscalar          :   grade(B,3) = 0"
" involutions  :     reversal              :   rev(B)   = 3e0+e1-5e2-e12-2e13"
"                    grade involution      :   invol(B) = 3e0-e1+5e2+e12+2e13"
"                    clifford conjugation  :   cj(B)    = 3e0-e1+5e2+e12+2e13"
" inverse 1/B                              :   inverse(B)"
" dual B                                   :   dual(B)"
" magnitude  |B|                           :   magnitude(B)"
" normalize  B/|B|                         :   normalize(B)"
" math functions :"
" exp1, log1, sqrt1, pow1, sin1, cos1, tan1, sinh1, cosh1, tanh1"
" asin1, acos1, atan1, asinh1, acosh1, atanh1"
""
""
"**********************************************************************"
" ENTER Cl(p) or Cl(p,q) :  p+q number basis vectors up to 5 "
" p positive squares, q negative squares"
" e.g.  Cl(3) for a 3D space  or  Cl(3,1) for a 4D Minkowski space"
cl(p,q1)=
do(tty=1,test(number(q1),q=q1,q=0),Clmax=p+q,Cldim=2^Clmax,test(Clmax=1,Clsig=(0,0),Clsig=zero(Clmax)),
check(p>=0),check(q>=0),for(i,1,Clmax,Clsig[i]=1),
test(q>0,for(i,p+1,Clmax,Clsig[i]=-1)),
Base(Clmax),Index(Clmax),consGtab(Cldim),consOtab(Cldim),consItab(Cldim),
test(Clmax=1,Clsig=Clsig[1]),               #  (0,0) --> (0) sinon zero(1) = 0 et non (0)
print("signature : "),Clsig)
 
Cl(p,q)=cl(p,q)
tty=svtty
signature(u)=print("  signature : ", Clsig)
 

# tables for cl p+q = 1 to 5  **************************************
 
Base(u)=do(Jtab=("e1","e12","e123","e1234","e12345"),         # Call : Base(Clmax)
test(u=1,do(B=unit(2),e0=B[1],e1=B[2],j=e1),
    u=2,do(B=unit(4),e0=B[1],e1=B[2],e2=B[3],e12=B[4],j=e12),
    u=3,do(B=unit(8),e0=B[1],e1=B[2],e2=B[3],e3=B[4],e12=B[5],e13=B[6],e23=B[7],e123=B[8],j=e123),
    u=4,do(B=unit(16),e0=B[1],e1=B[2],e2=B[3],e3=B[4],e4=B[5],e12=B[6],e13=B[7],e14=B[8],e23=B[9],e24=B[10],e34=B[11],e123=B[12],e124=B[13],e134=B[14],e234=B[15],e1234=B[16],j=e1234),
    u=5,do(B=unit(32),e0=B[1],e1=B[2],e2=B[3],e3=B[4],e4=B[5],e5=B[6],e12=B[7],e13=B[8],e14=B[9],e15=B[10],e23=B[11],
           e24=B[12],e25=B[13],e34=B[14],e35=B[15],e45=B[16],e123=B[17],e124=B[18],e125=B[19],e134=B[20],e135=B[21],
           e145=B[22],e234=B[23],e235=B[24],e245=B[25],e345=B[26],e1234=B[27],e1235=B[28],e1245=B[29],e1345=B[30],
           e2345=B[31],e12345=B[32],j=e12345)),B=quote(B),print("oriented volume : "," j ="Jtab[Clmax]))
 
Index(u)=do(                                                   # call : Index(Clmax)
test(u=1,Clindex=((0,0),(1,0)),
    u=2,Clindex=((0,0),(1,0),(0,1),(1,1)),
    u=3,Clindex=((0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1)),
    u=4,Clindex=((0,0,0,0),(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0),(1,0,1,0),(1,0,0,1),(0,1,1,0),(0,1,0,1),(0,0,1,1),(1,1,1,0),(1,1,0,1),(1,0,1,1),(0,1,1,1),(1,1,1,1)),
    u=5,Clindex=((0,0,0,0,0),(1,0,0,0,0),(0,1,0,0,0),(0,0,1,0,0),(0,0,0,1,0),(0,0,0,0,1),(1,1,0,0,0),(1,0,1,0,0),(1,0,0,1,0),(1,0,0,0,1),(0,1,1,0,0),(0,1,0,1,0),(0,1,0,0,1),(0,0,1,1,0),(0,0,1,0,1),(0,0,0,1,1),
    (1,1,1,0,0),(1,1,0,1,0),(1,1,0,0,1),(1,0,1,1,0),(1,0,1,0,1),(1,0,0,1,1),(0,1,1,1,0),(0,1,1,0,1),
    (0,1,0,1,1),(0,0,1,1,1),(1,1,1,1,0),(1,1,1,0,1),(1,1,0,1,1),(1,0,1,1,1),(0,1,1,1,1),(1,1,1,1,1))))
find(u)=do(for(i,1,Cldim,test(u=Clindex[i],rang=i)),rang)
 
ds(u)=disp(u)
disp(v)=do(isMvector(v),u=rnd(v),
test(Clmax=1,do(print(" ",u[1]*"e0",u[2]*"e1")),
    Clmax=2,do(print(" ",u[1]*"e0",u[2]*"e1"+u[3]*"e2",u[4]*"e12")),
    Clmax=3,do(print(" ",u[1]*"e0",u[2]*"e1"+u[3]*"e2"+u[4]*"e3",u[5]*"e12"+u[6]*"e13"+u[7]*"e23",u[8]*"e123")),
    Clmax=4,do(print(" ",u[1]*"e0",u[2]*"e1"+u[3]*"e2"+u[4]*"e3"+u[5]*"e4",u[6]*"e12"+u[7]*"e13"+u[8]*"e14"+u[9]*"e23"+u[10]*"e24"+u[11]*"e34",u[12]*"e123"+u[13]*"e124"+u[14]*"e134"+u[15]*"e234",u[16]*"e1234")),
    Clmax=5,do(print(" ",u[1]*"e0",
           u[2]*"e1"+u[3]*"e2"+u[4]*"e3"+u[5]*"e4"+u[6]*"e5",u[7]*"e12"+u[8]*"e13"+u[9]*"e14"+u[10]*"e15"+u[11]*"e23"+u[12]*"e24"+u[13]*"e25"+u[14]*"e34"+u[15]*"e35"+u[16]*"e45",u[17]*"e123"+u[18]*"e124"+u[19]*"e125"+u[20]*"e134"+u[21]*"e135"+u[22]*"e145"+u[23]*"e234"+u[24]*"e235"+u[25]*"e245"+u[26]*"e345",+u[27]*"e1234"+u[28]*"e1235"+u[29]*"e1245"+u[30]*"e1345"+u[31]*"e2345" ,u[32]*"e12345")))," ")
 
#    end tables for Cl p+q = 1 to 5  *****************************************
 
# multiplication table Gtab =   e12..k x e12..k = p+q
consGtab(u)=                                                   # call : ConsGtab(Cldim)
do(Gtab=zero(u,u),for(i,1,u,for(k,1,u,Gtab[i,k]=prod(Clindex[i],Clindex[k]))),corrGtab(u))
corrGtab(u)=                                                   # diagonal must be +e0
for(i,1,u,test(Gtab[i,i] == -1,for(k,1,u,Gtab[k,i]=-Gtab[k,i])))
# evaluate sign change in products with base vectors permutations  e12 e13 --> -e23
sign(u,v)=do(s=0,for(i,2,Clmax,test(u[i] == 1, for(k,1,i-1,s=s+v[k]))),(-1)^(s))
 
# base elements product : e13 e124  --> u=(1,0,1,0)  v=(1,1,0,1) in Cl4
prod(u,v)=  do(sL0=1,Luv=u+v,test(Clmax=1,L0=(0,0) ,L0=zero(Clmax)),  
  for(i,1,Clmax,
        test(Luv[i] < 2,L0[i]=Luv[i],Luv[i] = 2,do(sL0=sL0*Clsig[i],L0[i]=0))),sL0=sL0*sign(u,v),sL0*find(L0))
 
range(u,a,b)=test(a<b,and(u>=a,u<=b),and(u>=b,u<=a))
consOtab(u)=do(Otab=unit(u),F00=(2,3,4,11,26),test(Clmax<4,consOtab3(u),Clmax=4,consOtab4(u),Clmax=5,consOtab5(u)))
 
consOtab3(u)=do(for(i,1,u,for(k,1,i,Otab[i,k]=Gtab[i,k])),consOtab31(u))
consOtab31(u)=do(for(i,2,u-1,for(k,2,i,test(abs(Otab[i,k]) > F00[Clmax],Otab[i,k]=0)))) # eij
consOtab4(u)=do(for(i,1,u,for(k,1,i,Otab[i,k]=Gtab[i,k])),consOtab41(u),consOtab42(u),consOtab43(u))
consOtab41(u)=do(for(i,2,u-1,for(k,2,i,test(abs(Otab[i,k]) > F00[Clmax],Otab[i,k]=0))))
consOtab42(u)=do(for(i,2,11,for(k,2,i,test(range(abs(Otab[i,k]),6,11),Otab[i,k]=0))))   # eij
consOtab43(u)=do(for(i,12,16,for(k,12,i,test(range(abs(Otab[i,k]),6,11),Otab[i,k]=0)))) # eij
 
consOtab5(u)=do(for(i,1,u,for(k,1,i,Otab[i,k]=Gtab[i,k])),consOtab51(u),consOtab52(u),consOtab53(u),consOtab54(u),consOtab55(u))
consOtab51(u)=do(for(i,2,u-1,for(k,2,i,test(abs(Otab[i,k]) > F00[Clmax],Otab[i,k]=0))))
consOtab52(u)=do(for(i,2,16,for(k,2,i,test(range(abs(Otab[i,k]),7,16),Otab[i,k]=0))))     # eij
consOtab53(u)=do(for(i,17,26,for(k,17,i,test(range(abs(Otab[i,k]),7,16),Otab[i,k]=0))))
consOtab54(u)=do(for(i,2,26,for(k,2,i,test(range(abs(Otab[i,k]),17,26),Otab[i,k]=0))))    # eijk
consOtab55(u)=do(for(i,27,32,for(k,17,i,test(range(abs(Otab[i,k]),17,26),Otab[i,k]=0))))
 
# end tables *****************************************************************
Gtable(u)=  # multiplication tables Gtab  Otab  Itab  filtered by u, some multivector.
do(d=dim(Gtab),tab=zero(d,d),
  for(i,1,d,for(k,1,d,do(indx=Gtab[i,k],aindx=abs(indx),tab[i,k]=(u[aindx])*sgn(indx)))),tab)
Otable(u)=do(d=dim(Otab),tab=zero(d,d),
  for(i,1,d,for(k,1,i,do(indx=Otab[i,k],aindx=abs(indx),test(aindx>0 ,tab[i,k]=(u[aindx])*sgn(indx))))),tab)
         
gp(u,v)=   # main products gp inverse outp inp : geometric, inverse, outer, inner product
do(isMvector(u),isMvector(v),rnd(dot(Gtable(float(u)),float(v)),6))
inverse(u)=test(number(u),1/u,rnd(do(isMvector(u),inv(transpose(Gtable(float(u))))[1]),6))
outp(u,v)=do(isMvector(u),isMvector(v),rnd(dot(Otable(float(u)),float(v)),6))
 
# involutions :  **************************************************************
rev(u)=sum(k,0,Clmax,grade(u,k)*(-1)^(k*(k-1)/2))  # reversion
cj(u)=sum(k,0,Clmax,grade(u,k)*(-1)^(k*(k+1)/2))   # Clifford conjugation
invol(u)=rev(cj(u))                                # grade involution
 
# grade operator : <u>n  ******************************************************
grdtab=((1,2,0,0,0,0),(1,3,4,0,0,0),(1,4,7,8,0,0),(1,5,11,15,16,0),(1,6,16,26,31,32))
grade(u,n)=do(isMvector(u),isInteger(n),isGradeover(n),T00=grdtab[Clmax],grd00=zero(Cldim),test(n=0,grd00=u[1]*e0,
for(k,T00[n]+1,T00[n+1],grd00[k]=u[k])),grd00)
 
# magnitude(u)=|u|  normalize(u)=u/|u|  dual1(u)=u/j  with j the oriented volume = e12..k   with k = p+q #
magnitude(u)=test(number(u),abs(u),sqrt(gp(u,rev(u))[1]))
dual1(u)=do(j=unit(Cldim)[Cldim],gp(u,inverse(j)))
dual(u)=dual1(u)
normalize(u)=test(magnitude(u)=0,u,u/magnitude(u))
 
# Right, Left contractions, doth, & inner product (slower) ********************
Re(u)=do(isMvector(u),u[1])
Pu(u)=do(isMvector(u),u[Cldim])
lc(u,v)=do(j=unit(Cldim)[Cldim],gp(outp(u,gp(v,j)),inverse(j)))
rc(u,v)=do(j=unit(Cldim)[Cldim],gp(inverse(j),outp(gp(j,u),v)))
doth(u,v)=inp(u,v)-Re(u)*e0-u*Re(v)-Re(u)*v+Re(u)*Re(v)*e0
Q(u)=test(u=0,0,gp(u,rev(u))[1])
Bl(u,v)=1/2(Q(u+v)-Q(u)-Q(v))
inp(u,v)=rc(u,v)+lc(u,v)-Bl(u,rev(v))*e0
 
real1(x)=1/2(x+rev(x))       # polar form for versors **************************
imag1(x)=1/2(x-rev(x))
arg1(x)=test(not(magnitude(real1(x))= 0),arctan(magnitude(imag1(x))/magnitude(real1(x))),
            magnitude(imag1(x))>0,pi/2,
            magnitude(imag1(x))<0,-pi/2)
polar1(x)=do(isVersor(x),magnitude(x)*exp1(normalize(imag1(x))*arg1(x)))
rect1(x)=do(isVersor(x),test(arg1(x)=0,x,magnitude(x)*(normalize(real1(x))*cos(arg1(x))+normalize(imag1(x))*sin(arg1(x)))))
max(u)=                      # blade factorisation *****************************
do(isMvector(u),Umax=magnitude(u[1]),for(k,1,Cldim,test(magnitude(u[k])>Umax,do(Umax=magnitude(u[k]),Kmax=k))),Kmax)
factor1(u)=test(blade(u),do(Kmax=max(u),Bs=u/u[Kmax],b=zero(Clmax,Cldim),f=unit(Cldim),F1=inverse(f[Kmax]*u[Kmax]),
 for(k,1,Clmax,test( Clindex[Kmax,k]==0,1,b[k]=inp(inp(f[k+1],F1),Bs ))),b),print("not a blade")stop)
#   math functions *************************************************************
rnd(u,n)=do(test(number(n),Fixnum=n,Fixnum=fixnum),rnd00(u))
rnd00(u)= do(isMvector(u),zr00=zero(Cldim),for(k,1,Cldim,test(number(u[k]),zr00[k]=floor(u[k]*10^Fixnum+0.5)/10^Fixnum,zr00[k]=u[k])),float(zr00))  
   
 
exp0(u)=test(number(u),exp(u),                           # Padé coef (4,4)
do(u2=gp(u,u),u3=gp(u2,u),u4=gp(u2,u2),1.0*gp(e0+u/2+3u2/28+u3/84+u4/1680, inverse(e0-u/2+3u2/28-u3/84+u4/1680 ))))
exp1(u)= exp0(u-e0*Re(u))*float(exp(Re(u)))
qq(m,n,x)=sum(k,0,n,((m+n-k)!n!)/((m+n)!(n-k)!k! )(-x)^k)# compute the Padé coef. for exp = pp/qq  (ex: m=n=4 for exp1)
pp(m,n,x)=sum(k,0,m,((m+n-k)!m!)/((m+n)!(m-k)!k! )x^k)
exp11(v)=test(number(v),exp(v),do(u=float(v),un=e0,exp00=e0,for(k,1,maxTaylor,un=1.0*(gp(un,u)),exp00=exp00+un/k!),exp00))
 

power1(u,n)=test(number(u),u^n,do(isInteger(n),check(n>=0),pow00=e0,for(k,1,n,pow00=float((gp(pow00,u)))),pow00))
pow1(u,n)=power1(u,n)
 
sqrt0(u)=do(n=3,y(n) = test(n==0, u, (y(n-1) + inverse(z(n-1)))/2),  
z(n) = test(n==0, e0,(z(n-1) + inverse(y(n-1)))/2),1.0*y)
sqrt1(v)=do(u=normalize(v),rnd(sqrt0(u)*sqrt(magnitude(v)),6))        # square root  u^(1/2)
 
sqrtn(u,n)=test(n=1,sqrt1(u),sqrt1(sqrtn(u,n-1)))            # sq( sq( ...sq(u)  n times
sqrt10(u)=test(number(u),sqrt(u),exp1(log1(u)/2))            # alternate method
sqrt11(v)=do(u=normalize(v),rnd(sqrt10(u)*sqrt(magnitude(v)),6))
 
log0(v)=do(isPos(v),u=e0-v,u2=gp(u,u),u3=gp(u2,u),1.0*gp(-60u+60*u2-11*u3,inverse(60e0-90u+36*u2-3*u3)))
log1(v)=do(u=normalize(v),log0(u)+e0*log(magnitude(v)))  
### tester log(-v) # normalize
log01(u)=do(isPos(u),u1=gp(u-e0,inverse(u+e0)),u2=gp(u1,u1),un=u1,
                    test(u2==0,print("failed to calculate log",stop),1),                                
                    log00=u1,for(k,1,maxTaylor,un=float(gp(un,u2)), log00=log00+un/(2k+1)),
                    2*log00)+e0*log(magnitude(u))
log11(v)=do(u=normalize(v),log01(u)+e0*log(magnitude(v)))
#*********************************************************************************
# circular and hyperbolic functions
 
sin0(u)=test(number(u),sin(u),do(un=u,sin00=u,u2=gp(u,u),
           for(k,1,maxTaylor,un=1.0*(gp(un,u2)),
           sin00=sin00+(-1)^k*un/(2k+1)!),float(sin00)))
cos0(u)=test(number(u),cos(u),do(un=e0,cos00=e0,u2=gp(u,u),
            for(k,1,maxTaylor,un=gp(un,u2),
            cos00=cos00+(-1)^k*un/(2k)!),float(cos00)))
sin1(u)=float(sin(Re(u)))*cos0(u-e0*Re(u))+float(cos(Re(u)))*sin0(u-e0*Re(u))
cos1(u)=float(cos(Re(u)))*cos0(u-e0*Re(u))-float(sin(Re(u)))*sin0(u-e0*Re(u))
tan1(u)=test(number(u),tan(u),gp(sin1(u),inverse(cos1(u))))
asin1(u)=2*atan1(gp(u,inverse(e0+sqrt1(e0-gp(u,u)))))
acos1(u)=e0*float(pi/2)-asin1(a)
atan1(u)=do(j00=normalize(imag1(u)), test(number(u),arctan(u),0.5*gp(j00,log1(e0-gp(j00,u))-log1(e0+gp(j00,u)))))
 
sinh1(u)=test(number(u),sinh(u),0.5*(float(exp1(u))-1.0*exp1(-u)))
cosh1(u)=test(number(u),cosh(u),0.5*(exp1(u)+exp1(-u)))                          
tanh1(u)=test(number(u),tanh(u),gp(1.0*sinh1(u),inverse(1.0*cosh1(u))))
asinh1(u)=test(number(u),arcsinh(u),log1(u+sqrt1(power1(u,2)+e0)))
acosh1(u)=do(isGT1(u),test(number(u),arccosh(u),log1(u+sqrt1(power1(u,2)-e0))))
atanh1(u)=do(isRangeLT1(u),test(number(u),arctanh(u),0.5(log1(e0+u)-log1(e0-u))))
 
#********************************************************************************
#      outer algebra (product = outer product)  
                   
power2(u,n)=do(pow00=e0,for(k,1,n,pow00=outp(pow00,u)),pow00)
expi(u)=test(u=0,e0,do(un=e0,expi00=e0,for(k,1,maxTaylor,un=outp(un,u),expi00=expi00+un/k!),expi00))
exp2(u)=test(u=0,e0,do(u0=grade(u,0),expi(u-u0)*exp(u0[1])))
invi(u)=do(invi00=e0,for(k,1,maxTaylor,invi00=invi00+power2(u,k)*((-1)^k)),invi00)
 
inv2(u)=do(isZero(u[1]),float(invi((u-u[1]*e0)/u[1])/u[1]))
sqrti(u)=do( u2=outp(u,u),u3=outp(u2,u), u4=outp(u2,u2),u5=outp(u4,u),u6=outp(u4,u2),u7=outp(u6,u),u8=outp(u4,u4),u9=outp(u8,u),e0+u/2-u2/8+u3/16-u4*(5/128)+u5*(7/256)+u6*(21/1024)+u7*(33/2048)+u8*(429/32768)+u9*(715/65536))
sqrt2(u)=do(isPos(u[1]),float(sqrti((u-u[1]*e0)/u[1])*sqrt(u[1])))
sqrt22(u)=do(isPos(u[1]),exp2(0.5*log2(u)))
logi(u)=do(un=u,logi00=u,for(k,1,maxTaylor,un=outp(un,u),logi00=logi00+(-1)^k*un/(k+1)),1.0*logi00)
log2(u)=do(isPos(u[1]),test(u[1]=0,-INF*e0,logi((u-u[1]*e0)/u[1])+e0*log(u[1])))
sini(u)=test(u=0,0,do(un=u,sini00=u,u2=outp(u,u),for(k,1,maxTaylor,un=outp(un,u2),sini00=sini00+(-1)^k*un/(2k+1)!),sini00))
cosi(u)=test(u=0,e0,do(un=e0,cosi00=e0,u2=outp(u,u),for(k,1,maxTaylor,un=outp(un,u2),cosi00=cosi00+(-1)^k*un/(2k)!),cosi00))
sin2(u)=test(u=0,0,do(sin(u[1])*cosi(u-u[1]*e0)+cos(u[1])*sini(u-u[1]*e0)))
cos2(u)=test(u=0,e0,do(cos(u[1])*cosi(u-u[1]*e0)-sin(u[1])*sini(u-u[1]*e0)))
tan2(u)=test(u=0,0,do(outp(sin2(u),inv2(cos2(u)))))
asini(u)=test(u=0,0,do(u2=power2(u,2),un=u,asini00=u,
for(k,1,maxTaylor,un=outp(un,u2),asini00=asini00+((2k)!/(2^(2k)*((k!)^2)))*un/(2k+1)),1.0*asini00))
asin2(u)=test(u=0,0,do(isRangeLE1(u),arcsin(u[1])*e0+asini(u*sqrt(1-u[1]^2)-u[1]*sqrt2(e0-power2(u,2)))))
acos2(u)=e0*pi/2-asin2(u)
atani(u)=test(u=0,0,do(u2=power2(u,2),un=u,atani00=u,for(k,1,maxTaylor,un=outp(u2,un),atani00=atani00+((-1)^k)*un/(2k+1)),1.0*atani00))
atan2(u)=test(u=0,0,do(arctan(u[1])*e0+atani(outp((u-u[1]*e0),inv2(e0+u[1]*u)))))
sinh2(u)=do(0.5(exp2(u)-exp2(-u)))
cosh2(u)=do(0.5(exp2(u)+exp2(-u)))
tanh2(u)=outp(sinh2(u),inv2(cosh2(u)))
 
asinh2(u)=log2(u+sqrt2(power2(u,2)+e0))
acosh2(u)=log2(u+sqrt2(power2(u,2)-e0))
atanh2(u)=1/2*log2(outp(e0+u,inv2(e0-u)))
 
# differential operators     ****************************************************
 
nabla(u,x)=do(isMvector(u),isMvector(x),nabla00=zero(Cldim),ek=unit(Cldim),
  for(k,1,Cldim,test(not(x[k]==0),do(nabla00=nabla00+gp(ek[k],ek[k])[1]*gp(d(u,x[k]),ek[k])))),nabla00)
laplacian(u,x)=nabla(nabla(u,x),x)
 
# directional derivative of u with respect to x in direction of a
diff(u,a,x)=do(isMvector(u),isMvector(x),nabla00=zero(Cldim),ek=unit(Cldim),C(k)=gp(ek[k],ek[k])[1],
  for(k,1,Cldim,test(not(x[k]==0),do(nabla00=nabla00+ gp(a,ek[k])[1]*C[k]*d(u,x[k])))),nabla00)
 

# projection rotation reflection rotation   **************************************
project(u,v)=gp(inp(u,v),inverse(v))
reject(u,v)=u-project(u,v)
reflect(u,v)=do(test(evenGrade(v),gp(gp(-invol(v),u),inverse(v)),oddGrade(v),gp(gp(v,u),inverse(v))))
join(u,v)=do(isMvector(u),isMvector(v),H=outp(cunit(u),cunit(v)),cunit(H))
meet(u,v)=do(isMvector(u),isMvector(v),gp(outp(gp(u,-j),gp(v,-j)),j))
RR(R,u)=               # RR = R u R~  
do(test(versor(R),1,print("1st arg not a rotor",stop)),gp(gp(R,u),rev(R)))
R(B,theta)=            # R = exp(-B T/2)    B rotation plane   theta rotation angle in degree
do(test(grade2(B),1,print("1st arg not a bivector",stop)),test(number(theta),T=float(theta*pi/180),T=magnitude(B)),B1=normalize(B),
                  (exp1(-B1*T/2)))
M(R1)=do(test(Clmax=3,1,print("works only for Cl(3)",stop)),V00=zero(3,8),M00=zero(3,3),  # matrix form of rotor R1
V00[1]=RR(R1,e1),V00[2]=RR(R1,e2),V00[3]=RR(R1,e3),for(i,1,3,for(j,1,3,M00[i,j]=V00[i,j+1] )),M00)
 
# ortogonal  &  reciprocal frames  (3D)   *****************************************
orthogonal(u,v,w)=do(isVector(u),isVector(v),isVector(w),
test(outp(outp(u,v),w)=0,print("not a frame",stop)),
u0=u,v0=gp(outp(v,u0),inverse(u0)),w0=gp(outp(outp(w,u0),v0),inverse(outp(u0,v0))),
((u0[2],u0[3],u0[4]),(v0[2],v0[3],v0[4]),(w0[2],w0[3],w0[4])))
 
reciprocal(u,v,w)=do(isVector(u),isVector(v),isVector(w),
#  any set of linearly independant vectors [u,v,w] form a basis (frame)
#  reciprocal basis [r1,r2,r3]  with  r1.u=1 r2.v=1 r3.w=1
#  for orthonormal frame [e1,e2,e3] reciprocal is [e1,e2,e3]
test(outp(outp(u,v),w)=0,print("not a frame",stop)),
u123=outp(outp(u,v),w),
ru=gp(outp(v,w),inverse(u123)),rv=gp(outp(w,u),inverse(u123)),rw=gp(outp(u,v),inverse(u123)),
((ru[2],ru[3],ru[4]),(rv[2],rv[3],rv[4]),(rw[2],rw[3],rw[4])))
#***********************************************************************************
# examples :  add here tutorial scripts  and  delete after use to save script buffer space
#***********************************************************************************
"*******************************************************************************"
GAintro()=do( tty=1,print(" pls delete this tutorial when unnecessary           ",
"*******************************************************************************",
"",
" Clifford algebra Cl(V) is generated from an inner product vector space by ",
" a product of any vector with itself equal to their inner product :","",
"                                aa = a.a",
" (inner product defines the measurement of length and angle)","",
" in order to learn how to compute with Clifford numbers, we begin with a low",
" dimensional examples :","",
" a Euclidean plane R2 is spanned by    e1, e2 <- R2      with",
" e1.e1=1, e2.e2=1, e1.e2=0"," {e1, e2} is an orthonormal vector basis of R2",
""," under Clifford's associative geometric product we set",
"       e1² = e1e1 = e1.e1 = 1,               e2² = e2e2 = e2.e2 = 1 ",
" and  (e1 + e2)(e1 + e2) = e11 + e22 + e1e2 + e2e1 = 2 + e1e2 + e2e1",
"                         = (e1 + e2).(e1 + e2)     = 2",
" ==>   e1e2 + e2e1 = 0   <==>   e1e2 = -e2e1",
""," i.e. the geometric product of orthogonal vectors forms a new entity, called ",
" unit bi-vector :     e12 = e1e2, and is anti-symmetric. ","",
" for orthogonal vectors the geometric product equals exterior product, symbol ^",
"                              e12 = e1e2 = e1^e2 = -e2^e1 = -e2e1 = -e21","",
" using associativity, we can compute the products  e1e12 =  e1e1e2 =  e2",
"                                                   e2e12 = -e2e2e1 = -e1",
" which represent a mathematically positive (anti-clockwise) 90° rotation ","",
" the opposite order gives                          e12e1 = -e21e1 = -e2",
"                                                   e12e2 =  e1",
" which represents a mathematically negative (clockwise) 90° rotation ","",
" the bivector e12 acts like a rotation operator and we observe the general",
" anti-commutation property",
"                    ae12 = -e12a,   for any a = a1e1+a2e2 <- R2, a1,a2 <- R",
""," the square of the unit bivector is -1     e12² = e1e2e12 = e1(-e1) = -1",
" just like the imaginary unit i of complex numbers","",
" the general geometric product of two vectors a, b <- R2",
" ab = (a1e1 + a2e2)(b1e1 + b2e2) = a1b1 + a2b2 + (a1b2 - a2b1)e12",
"                                 = 1/2(ab + ba) + 1/2(ab - ba) = a.b + a^b",
"            ab = a.b + a^b",
" has therefore a scalar symmetric inner product part",
"            1/2(ab + ba) = a.b =  a1b1 + a2b2     = |a||b|cos(theta)",
" and a bi-vector skew-symmetric outer product part",
"            1/2(ab - ba) = a^b = (a1b2 - a2b1)e12 = |a||b|e12sin(theta)","",
"     parallel   vectors (theta(a,b) = 0°)  commute,         ab = a.b =  ba",
" and orthogonal vectors (theta(a,b) = 90°) anti-commute,    ab = a^b = -ba",
"", " the outer product part a^b represents the oriented area of",
" the parallelogram spanned by the vectors a and b in the plane of R2,",
" with oriented magnitude    det(a,b) = |a||b|sin(theta) = -(a^b)e12 ","",
" with the Euler formula we can rewrite the geometric product as",
" ab = |a||b|(cos(theta) + e12 sin(theta)) = |a||b|exp(e12 theta)",
" again because e12² = -1","",
 
" EVA example :",
" first, we have to define the space basis (Cl(2) in this example) :",
" {e0,e1,e2,e12} e0 for scalar, e1 and e2 for vector, e12 for bivector",
" a Cl(2) Clifford number representation is 3e0+2e1-e2+5e12 = (3, 2, -1, 5)",
" Cl(2)",Cl(2), a=a1*e1+a2*e2, b=b1*e1+b2*e2,
" we define two vectors a and b :  a=a1e1+a2e2, b=b1e1+b2e2",
" EVA functions for inner, outer and geometric products are :",
"  inp(a,b), outp(a,b), gp(a,b) ",tty=0,
" the geometric product a a is gp(a,a)", gp(a,a),
" the geometric product a b is gp(a,b)", gp(a,b),tty=1,
" result for Cl(2) is a 4D (= 2²) Clifford number (se0, 1ve1, 1ve2, 2ve12) with",
" s (or 0v)is scalar component,      grade0 or <ab>0 ",
" 1v 1v  are vector components,      grade1 or <ab>1 ",
" 2v     is bivector component,      grade2 or <ab>2 ","",
" to save space, we may output these results as line vectors :", gp(a,a), gp(a,b),
""," gp(a,a) gives a scalar a1^2+a2^2 as expected, other components are 0",
" gp(a,b) gives scalar component a1b1+a2b2 and a bivector component a1b2-a2b1",
" scalar component is the inner product   inp(a,b) :  ", inp(a,b),
" bivector component is the outer product out(a,b) :  ", outp(a,b),
" and gp(a,b) = inp(a,b) + outp(a,b) (for vectors) ","",
" we may project components in Cl(2) space with the EVA function grade(u,n) :",
" grade(gp(a,b),0) = ", grade(gp(a,b),0)," and grade(gp(a,b),2) = ", grade(gp(a,b),2),
""," the geometric product of vectors is invertible for all vectors",
" with non-zero square a² not= 0",
"         inverse(a)     =   a/a²",                    
"         a inverse(a)   =  a a/a²         = 1",
"         inverse(a) a   = (a/a²)a = a²/a² = 1",
" the inverse vector a/a² is a rescaled version (reflected at the unit circle)",
" of the vector a. This invertibility leads to enormous simplifications and ease",
" of practical computations.","",
" EVA function for inverse is inverse(x).  e.g. with  a=2e1+e2, b=3e1-4e2 :",
a=2e1+e2, b=3e1-4e2,
" inverse(a) =", inverse(a)," inverse(b) =", inverse(b),
" gp(inverse(a),a) =", gp(inverse(a),a), " gp(inverse(b),b) =", gp(inverse(b),b),
"","      projection :   project(a,b)    reject(a,b)  with a=2e1+e2, b=3e1-4e2 :",
"          projection of vector a onto b ",
" prab = (a.b/|b|)(b/|b|) = (a.b)(b/|b|² = (a.b)/b    ","",
" prab = project(a,b)  ", project(a,b),
"          rejection of vector a from b (perpendicular part) :",
" rjab = a - prab = a - (a.b)/b = ab/b - (a.b)/b = 1/b(ab - a.b) = (a^b)/b","",
" rjab = reject(a,b)  ", reject(a,b),""," prab + rjab = ",project(a,b)+reject(a,b),"",
" we can now use prab, rjab to compute the reflection of a = prab + rjab at the line",
" (hyperplane) with normal vector b, which means to reverse prab -> -prab",
" reflect(a,b) = ",reflect(a,b), " reflect(reflect(a,b),b) = ",reflect(reflect(a,b),b)))