#example functions used later on
f(x)=x^3-x

#tangent line at point p of function f(x)
t(x,p)=f'(p)(x-p)+f(p)
#so f(x) # t(x,1) # t(x,-2) will show the tangent line of f(x) at 1, which intersects at -2p so -2

#perpendicular to the tangent at given point p
pt(x,p)=(p-x)/f'(p)+f(p)

#vertical line from top down reflecting at point p of f(x), to graph original line run y0+p
r(x,p)=fm=f'(p):(x-p)(fm^2-1)/(2fm)+f(p)
#to only show the reflection
#r(x,p)=pw({x(f'(p)^2-1)/(2f'(p))+f(p)-p*(f'(p)^2-1)/(2f'(p)),(sgn(f'(p))==-1&&x>=p)||(sgn(f'(p))==1&&x<=p)})

#r(x,p) extended to a line mx+b reflecting off of f(x) at point p
r(x,p,m)=(x-p)tan(2atan(f'(p))-atan(m))+f(p)
#in reference to the vars for the function the original line would be mx+f(p)-m*p
rl(x,p,m)=m(x-p)+f(p)

#newtons method to find 0, off of f(x)
n(x)=x-f(x)/f'(x)

#trapazoidal method of approximenting area under the curve as n goes towards infinity
ar(a,b,n)=sum(k,(b-a)/n(f(a+k(b-a)/n)+f(a+(k+1)(b-a)/n))/2,0,n-1)

#basic method of approximenting arc length as n goes towards infinity
al(a,b,n)=sum(k,sqrt(((b-a)/n)^2+(f(a+(k+1)(b-a)/n)-f(a+k(b-a)/n))^2),0,n-1)
#gl(a,b,n) shows the lines al(a,b,n) uses to calculate the arc length, so vrange=6;f(x);gl(0,2,4) shows 4 lines that approximate the arc length of f(x) from 0 to 2
gl(a,b,n)=mat(k,{a+k(b-a)/n,f(a+k(b-a)/n)},0,n)

#nth derivitive, precision should be around n*128, so this fails at the 5th derivitive at default precision
#dn(x,n)=sum(k,(-1)^k C(n,k) f(x+(n-k)2^-100),0,n)/2^(-100n)
#dn(x,n)=pw({slope(p,dn(p,n-1),x),n>=1},{f(x),1})

#f(x) rotated by θ
fr(x,θ)=φ=atan(x,f(x))+θ:{cos(φ),sin(φ)}sqrt(x^2+f(x)^2)

#the following is a function to estimate the path of a moving positive charge under the effects of any number of stationary positive or negitive charges, loss of accuracy when close to a stationary charge
#c is the starting {x,y,vx,vy} of the moving charge
#t is the time intervals that will be computed
#a is the cordonates and charge of the stationary charges, organized in a matrix like {{x1,y1,c1},{x2,y2,c2}...{xN,yN,cN}}
#example graph function with the moving charge stationary at the origin with 1/10th second time intervals and a positive chage at {0,1} and negitive charge at {2,0}
#part(iter(g,ef(g,1/10,{{0,1,1},{2,0,-1}}),{0,0,0,0},100,1),-1,0..1))#{{0,1}}#{{2,0}}
ef(c,t,a)=x=part(c,0):y=part(c,1):vx=part(c,2):vy=part(c,3):r=vec(i,(part(a,i,0)-x)^2+(part(a,i,1)-y)^2,0,len(a)-1):θ=mat(i,{part(a,i,0)-x,part(a,i,1)-y}/sqrt(part(r,i)),0,len(a)-1):xp=sum(i,-part(a,i,2)part(θ,i,0)/part(r,i),0,len(a)-1):yp=sum(i,-part(a,i,2)part(θ,i,1)/part(r,i),0,len(a)-1):{x+t*vx+t^2/2*xp,y+t*vy+t^2/2*yp,vx+t*xp,vy+t*yp}

#refraction of light using snells law, x=1 to x=2 is considered a different medium with a refraction index of 'b', 'm' is starting slope of line, 'a' is refraction index outside of x=1 to x=2
#graphing s(x,2,1.5,1.75) will show a refraction example with an index of refraction of 1.5 going to 1.75 with a initial slope of 2
s(x,m,a,b)=d=(am/sqrt(m^2+1)/b)/(sqrt(1-(am/sqrt(m^2+1)/b)^2)):r=(bd/sqrt(d^2+1)/a)/(sqrt(1-(bd/sqrt(d^2+1)/a)^2)):pw({rx+d-r,x>2},{dx+m-d,2>x>1},mx)

#turns a complex number into coordinates on a riemann sphere
rs(z)=θ=atan(re(z),im(z)):c=cos(θ):s=sin(θ):n=1-2/(|z|+1):r=sqrt((1-n^2)/(c^2+s^2)):{c r,s r,n}

#nth deritivive of x^p at x=b, negitive p values corrospond to the anti derititive. cant graph negitive integers p likely do to logarithms from integration, however where logarithms for integration do not apply the whole numbers will work, ie g(-1,-2,3)=-1/3 but g(-2,-2,3) will have a logarithm so it fails
g(n,p,b)=P(p,n) b^(p-n)
#below works for -p values but does not integrate properly for -p values
#g(n,p,b)=lim(a,re(P(p+ai,n)),0) b^(p-n)
#another continuation that can graph negitive p values, however it outputs complex numbers at non whole n values
#g(n,p,b)=(-1)^n b^(p-n) ph(-p,n)
#for p=-1 fixes negitive integer n values
#g(n,p,b)=lim(a,(-1)^(n+a) b^(p-n-a) ph(-p,n+a) + b^(p-n) / (a(p-n)!),0)


#functions used for tangent planes
z(x,y)=x^3-x+y^3-y
zx(x,y)=slope(p,z(p,y),x)
zy(x,y)=slope(p,z(x,p),y)
#tangent plane at point px, py
zt(x,y,px,py)=(zx(px,py))*x+(zy(px,py))*y+(z(px,py)-zx(px,py)*px-zy(px,py)*py)
#so z(x,y) # zt(x,y,2,3) will have the tangent plane of z(x,y) at (2,3)

#z(x,y) rotated by θ and φ
zr(x,y,θ,φ)=α=atan(z(x,y),sqrt(x^2+y^2))+φ:β=atan(x,y)+θ:{sin(α)cos(β),sin(α)sin(β),cos(α)}sqrt(x^2+y^2+z(x,y)^2)


#gets how many digits match, ch(22/7,pi)=3 as 3.142.. equals 3.141.. upto the third decimal
ch(a,b)=floor(log(10,max{a,b}))-floor(log(10,|a-b|))


#i,j,k unit vectors if you like those
#i={1,0,0}
#j={0,1,0}
#k={0,0,1}

#a cube given by 3d points
cube={{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,0},{0,0,1},{0,1,1},{0,1,0},{0,1,1},{1,1,1},{1,1,0},{1,1,1},{1,0,1},{1,0,0},{1,0,1},{0,0,1}}
#the following is graphing this cube and it rotated pi/4 rad in yaw a roll directions
#ticks=3;xr=-1,2;yr=-1,2;zr=1.5;cube#cube*rotate(pi/4,0,pi/4)

#minimum distance between b^x and log(b,x) over b where b>=root(e,e)
minexp(x)=sqrt(2)log(x,eln(x))


#minimum distance between left and right side of x^-(2n+1) over n(n is natural)
minpow(x)=2sqrt(root(x,2x+2)^2+root(x,2x+2)^(-2*x))


#minimum distance between sec(a*x) and its inverse, asec(x)/a
minsec(x)=sqrt(2)(sec(asin((sqrt(x^2+4)-x)/2))-(asin((sqrt(x^2+4)-x)/2)/x))


#following is minimum distance between the center curve and curve right of tan(a*x)
#distance between (x,tan(ax)) and (pi/a-x,-tan(ax))
disttan(x,a)=sqrt((2x-pi/a)^2+4tan(ax)^2)

#derivitive of above set to 0 simplified
mintan(x,a)=x+a*sin(ax)/cos(ax)^3-pi/(2a)

#derivitive of above to find 0 by newtons method
mintanp(x,a)=a^2((1+2sin(ax)^2)/cos(ax)^4)+1

#modified newtons method to allow graphing along 'a'
n(x,a)=x-mintan(x,a)/mintanp(x,a)

#disttan(n(n(n(n(1/x,x),x),x),x),x) for example will graph it fairly nicely, we start with 1/x to scale properly and not get trapped in a asymptote
#the following graphically shows that it finds the minimum for just tan(c*x) by inputing the following into the calculator as the circles only intersect the other curve once with 1 line going between the 2 points of interest and the other should intersect the 2 circles at the same points while being perpendicular to the other line, you may need to increase the amount of 'n' to increase precision
# c=2 ; a=n(n(n(n(n(1/(2c),c),c),c),c),c) ; b=disttan(a,c) ; tan(c*x) # piecewise({+-sqrt(b^2-(x-a)^2)+tan(a*c),x<=(a+b)&&x>=(a-b)}) # piecewise({+-sqrt(b^2-(x-(pi/c-a))^2)-tan(c*a),x<=((pi/c-a)+b)&&x>=((pi/c-a)-b)}) # (-2tan(a*c)/(pi/c-2a))x+tan(a*c)-(-2tan(a*c)/(pi/c-2a))a # x*c*sec(a*c)^2-sec(a*c)^2*pi/2


#{year,month,day,hour,min,second} derived from epoch time, with o as hour offset
#date(o)=t=epoch:{1970+floor(t->year),floor((fract(t->year)year->day)/30.4375)+1,floor(fract((fract(t->year)year->day)/30.4375)day30.4375->day),floor(fract(fract(t->year)year->day)day->hour)+o,floor(fract(fract(fract(t->year)year->day)day->hour)hour->min),floor(fract(fract(fract(fract(t->year)year->day)day->hour)hour->min)min->second)}


#gets the distribution data for the sum of k amounts of n sided die starting at 0, dice(6,2,2) will give the frequency that 2 6 sided die will give a sum of 4 or (p+k) which is 3 so 3/6^2 will give the chance of rolling a sum of 4 from 2 6 sided die, than you can use 'vec(b,dice(6,2,b),0,(6-1)2)' to get the entire distribution data for 2 6 sided die
#dice(n,k,p)=piecewise({1,p==0},{0,(0>p)||(p>(k(n-1)))},{sum(a,dice(n,k-1,p-a),0,n-1),1})


#hyperoperation function
H(a,n,b)=piecewise({b+1,n==0},{a,n==1&&b==0},{0,n==2&&b==0},{1,n>=3&&b==0},{H(a,n-1,H(a,n,b-1)),1})


#graph of the chance for the sum of 2 die to be less then value of a card picked
# a=dice{6,6} ; b=1..13 ; vec(c,part(a,c)*part(b,c),0,10) # vec(c,part(a,c)*(13-part(b,c)),0,10)


#vector input example, m{1,2,3,4}=1^2+2^2+3^2+4^2
m(v)=sum(a,part(v,a)^2,0,len(v)-1)


#x^3+px+q=0, depressed cubic :c
dc(p,q)=cbrt(-q/2+sqrt(q^2/4+p^3/27))+cbrt(-q/2-sqrt(q^2/4+p^3/27))


#complex expansion of fibonacci sequence
fibonacci(x)=(φ^x-(-1/φ)^x)/sqrt(5)