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