1# Math Library
2
3Two questions:
4* What does the math library look like, from a user perspective?
5* Which math operations should be primitive, and which are implemented in Curv?
6  Note that functions written in Curv are still compiled into optimized GL code.
7
8## Math Library Implementation
9* Some math operations are coded in Curv. That's best if there's no efficiency
10  reason to do otherwise, because the code is shorter.
11* Some math operations are built-in functions, for efficiency reasons.
12  * If the function is in the C library and is not something you'd want to
13    write in Curv, like trig and exponential functions.
14  * If the function is a GLSL built-in function. The goal is to generate
15    efficient GLSL code, and this is the most conservative option.
16* Some Curv primitive math operations are written as idioms that compile into
17  special internal operators. Eg, `(a*b)+c` compiles to FMA(a,b,c),
18  e^x compiles to EXP(x).
19
20A fluent and parsimonious math API:
21Instead of providing a large number of 'special case' numeric functions
22like expm1(x) == e^x-1, you instead write 'natural' code and the compiler
23will detect common patterns and rewrite them into more efficient/accurate
24code. For example:
25* e^x -> exp(x)
26* e^x - 1 -> expm1(x)
27* 1/sqrt(x) -> invertsqrt(x) // GLSL
28* log 10 x -> log10(x)
29* log 2 x -> log2(x)
30* ln(1+x) -> log1p(x)
31
32## Math Library API
33
34logarithm
35  Currently, just `log x` -- natural logarithm. But GLSL also has log2(x).
36  Add log base b. It needs a different name than natural log.
37    `ln x` for natural log, since `log x` is ambiguous (means natural log in
38      math, common log in chemistry).
39    `log b x` -- same arg order as $log_b x$ in math notation.
40      Use `x >> log b` for pipelining, map(log b) for combinator programming.
41      Use partial evaluation to convert `log 2` and `log 10` to internal
42      log2 and log10 primitives.
43      C, Javascript, etc `log2(x),log10(x)` => log 2 x, log 10 x
44  Notes:
45    Haskell: log x, logBase b x
46    Mathematica: log[x], log[b,x]
47    J:  ^.x,  b^.x
48    Python: log(x), log(x,b) -- `b` is second arg because it is optional
49    GLSL: log(x) is natural log, no 2 arg version. log2(x).
50    C, Julia, Javascript: log(x), log2(x), log10(x), logp1(x).
51
52Fractional part of a number.
53Mathworld names the function `frac(n)`, gives two competing definitions:
54* sawtooth wave: mod(n,1) or n - floor(n). Result is non-negative.
55  Eg: HLSL frac, GLSL fract, Mathematica SawtoothWave, MuPad frac,
56* rem(n,1) or n - trunc(n). Result has the same sign and digits after the
57  decimal point as the argument. Eg, C Julia modf, Mathematica FractionalPart
58  Maple frac(n).
59* Notationally, I think mod[n,1] and rem[n,1] are fine.
60
61Modulus: Currently `mod[a,m]`, m is the modulus.
62  pipeline: a `mod` m
63  tacit: (`mod` m) using a proposed extension. Eg frac = `mod` 1
64  Consistent with math notation 'a (mod m)' and 'mod(a,m)', seems good enough.
65
66mag(v) = sum(v^2);
67* Magnitude of a vector. More like high-school math than 'norm', which is also
68  a more abstract and general concept.
69
70phase[x,y] = atan2(y,x).
71* Convert a 2D vector to an angle. It's a standard engineering name for the
72  concept, as per MathWorld and Wikipedia. `arg(p)` is a pure math alternative,
73  but the word "argument" has a different meaning in programming.
74  angle(z) is "phase angle" of a complex number in MatLab.
75  Python cmath has "phase".
76
77cis(theta) = [cos theta, sin theta].
78* Convert an angle to a unit vector.
79  The `cis` function is available in many high performance math libraries,
80  where it is faster than calling cos and sin separately.
81  Python cmath has "cis".
82
83Complex numbers.
84[re,im] is the representation (no abstract type).
85* mag(z) -- complex magnitude -- like abs(z) in MatLab
86* phase(z) -- complex phase
87* conj[re,im] = [re,-im] -- complex conjugate.
88  Reflects a 2-vector across the X axis.
89  Is it useful in Curv? Dunno.
90  * cmul(z, conj(z)) == a^2 + b^2, which is real and non-negative // z==[a,b]
91    and is mag(z)^2.
92  * Used when computing complex division by hand to cancel out the i's.
93  * sqrt(z) * conj(z) == mag(z)
94  * To convert a real function to a positive value, you can square it.
95    The analog for complex functions is to multiply by the conjugate.
96* from polar coordinates: `r*cis(theta)`
97* z+w
98* z-w
99* cmul([a,b],[c,d]) = [a*c - b*d, b*c + a*d];
100  This adds the phases and multiplies the magnitudes.
101  eg, cmul(vec,cis(theta)) rotates a 2D vector by angle theta.
102* cdiv(z,w)
103  Subtract the phases and divide the magnitudes.
104* there are also exponentials, trig
105
106Normalize a vector (convert to a unit vector).
107Possible names: normalize(v) -- GLSL, unit(v) -- a popular name,
108    normalize(v,p) -- Julia,
109Implementation should avoid division by zero:
110  unit v =
111    let(magv = mag v)
112    if (magv > 0) // or EPS
113      v / magv
114    else
115      v;
116  unit v = (
117    magv = mag v;
118    if (magv > 0) // or EPS
119      v / magv
120    else
121      v
122  );
123
124Shadertoy.com classifies the following GLSL ES features as being useful.
125
126void bool int float vec2 vec3 vec4 bvec2 bvec3 bvec4 ivec2 ivec3 ivec4 mat2 mat3 mat4 sampler2D
127
128sampler2D represents a 2D texture
129
130type radians (type degrees)
131type degrees (type radians)
132type sin (type angle)
133type cos (type angle)
134type tan (type angle)
135type asin (type x)
136type acos (type x)
137type atan (type y, type x)
138type atan (type y_over_x)
139
140type pow (type x, type y)   x^y
141type exp (type x)           e^x
142type log (type x)
143type exp2 (type x)          2^x
144type log2 (type x)
145type sqrt (type x)
146type inversesqrt (type x)   1/sqrt(x)
147
148type abs (type x)
149type sign (type x)
150type floor (type x)
151type ceil (type x)
152type fract (type x)
153type mod (type x, float y)
154type mod (type x, type y)
155type min (type x, type y)
156type min (type x, float y)
157type max (type x, type y)
158type max (type x, float y)
159type clamp (type x, type minV, type maxV)
160type clamp (type x, float minV, float maxV)
161    x if x>minV && x<maxV, minV if x<=minV, maxV if x>=maxV
162type mix (type x, type y, type a)
163type mix (type x, type y, float a)
164    linear blend of x and y. x*(1-a)+y*a
165type step (type edge, type x)
166type step (float edge, type x)
167    0 if x<edge, else 1
168    bit(x<edge) or bit(edge>x)
169type smoothstep (type a, type b, type x)
170type smoothstep (float a, float b, type x)
171    0 if x<a, 1 if x>b, else interpolate between 0 and 1 using Hermite polynom.
172
173float length (type x)
174    euclidean norm
175float distance (type p0, type p1)
176    norm(p0-p1)
177float dot (type x, type y)
178vec3 cross (vec3 x, vec3 y)
179type normalize (type x)
180    x/norm(x)
181type faceforward (type N, type I, type Nref)
182    a vector that points in the same direction as Nref
183type reflect (type I, type N)
184type refract (type I, type N,float eta)
185float determinant(mat* m)
186type matrixCompMult (type x, type y)
187type inverse (type inverse)
188
189float intBitsToFloat (int value)
190uint uintBitsToFloat (uint value)
191int floatBitsToInt (float value)
192uint floatBitsToUint (float value)
193
194bvec lessThan(vec x, vec y)
195bvec lessThan(ivec x, ivec y)
196bvec lessThanEqual(vec x, vec y)
197bvec lessThanEqual(ivec x, ivec y)
198bvec greaterThan(vec x, vec y)
199bvec greaterThan(ivec x, ivec y)
200bvec greaterThanEqual(vec x, vec y)
201bvec greaterThanEqual(ivec x, ivec y)
202bvec equal(vec x, vec y)
203bvec equal(ivec x, ivec y)
204bvec equal(bvec x, bvec y)
205bvec notEqual(vec x, vec y)
206bvec notEqual(ivec x, ivec y)
207bvec notEqual(bvec x, bvec y)
208* shorter names?
209* <' <=' >' >=' ==' !='
210* less, greater, not_less, not_greater, equal, not_equal
211
212bool any(bvec x)
213bool all(bvec x)
214bvec not(bvec x)
215
216vec4 texture(sampler* sampler, vec2 coord [, float bias])
217vec4 textureLod(sampler* sampler, vec2 coord, float lod)
218vec4 textureGrad(sampler* sampler, vec2 P, vec2 dPdx, vec2 dPdy)
219vec4 textureProj(sampler* sampler, vec3 coord [, float bias])
220vec4 textureProjLod(sampler* sampler, vec3 coord, float lod)
221vec4 textureProjGrad(sampler* sampler, vec3 P, vec2 dPdx, vec2 dPdy)
222vec* textureSize(sampler* sampler, int lod)
223
224type dFdx( type x ), dFdy( type x )
225type fwidth( type p )
226