(no commit message)
[libreriscv.git] / simple_v_extension / vector_ops.mdwn
index edcd88ce74d1d3d18f19402c682ed9f202af8ccf..5294cbdb5012e5dfde9f2c84ae89f871553beb50 100644 (file)
@@ -1,16 +1,75 @@
+[[!tag oldstandards]]
+
+**OBSOLETE**, see [[openpower/sv/3d_vector_ops]]
+
 # Vector Operations Extension to SV
 
+This extension defines vector operations that would otherwise take several cycles to complete in software. With 3D priorities being to compute as many pixels per clock as possible, the normal RISC rules (reduce opcode count and make heavy use of macro op fusion) do not necessarily apply. 
+
 This extension is usually dependent on SV SUBVL being implemented. When SUBVL is set to define the length of a subvector the operations in this extension interpret the elements as a single vector.
 
 Normally in SV all operations are scalar and independent, and the operations on them may inherently be independently parallelised, with the result being a vector of length exactly equal to the input vectors.
 
 In this extension, the subvector itself is typically the unit, although some operations will work on scalars or standard vectors as well, or the result is a scalar that is dependent on all elements within the vector arguments.
 
+However given that some of the parameters are vectors (with and without SUBVL set), and some are scalars (where SUBVL will not apply), some clear rules need to be defined as to how the operations work.
+
 Examples which can require SUBVL include cross product and may in future involve complex numbers.
 
+## CORDIC
+
+6 opcode options (fmt3):
+
+* CORDIC.lin.rot vd, vs, beta
+* CORDIC.cir.rot vd, vs, beta
+* CORDIC.hyp.rot vd, vs, beta
+* CORDIC.lin.vec vd, vs, beta
+* CORDIC.cir.vec vd, vs, beta
+* CORDIC.hyp.vec vd, vs, beta
+
+
+| Instr | result | src1 | src2 | SUBVL | VL | Notes |
+| ------------------ | ------ | ---- | ---- | ----- | -- | ------ |
+| CORDIC.x.t vd, vs1, rs2 | vec2 | vec2 | scal | 2 | any | src2 ignores SUBVL |
+
+SUBVL must be set to 2 and applies to vd and vs. SUBVL is *ignored* on beta.  vd and vs must be marked as vectors.
+
+VL may be applied.  beta as a scalar is ok (applies across all vectors vd and vs). Predication is also ok (single predication) sourced from vd. Use of swizzle is also ok.
+
+Non vector args vd, vs are reserved encodings.
+
+CORDIC is an extremely general-purpose algorithm useful for a huge number
+of diverse purposes.  In its full form it does however require quite a
+few parameters, one of which is a vector, making it awkward to include in
+a standard "scalar" ISA.  Additionally the coordinates can be set to circular,
+linear or hyperbolic, producing three different modes, and the algorithm
+may also be run in either "vector" mode or "rotation" mode.  See [[discussion]]
+
+CORDIC can also be used for performing DCT.  See
+<https://arxiv.org/abs/1606.02424>
+
+CORDIC has several RADIX-4 papers for efficient pipelining.  Each stage requires its own ROM tables which can get costly.  Two combinatorial blocks may be chained together to double the RADIX and halve the pipeline depth, at the cost of doubling the latency.
+
+Also, to get good accuracy, particularly at the limits of CORDIC input range, requires double the bitwidth of the output in internal computations. This similar to how MUL requires double the bitwidth to compute.
+
+Links:
+
+* <http://www.myhdl.org/docs/examples/sinecomp/>
+* <https://www.atlantis-press.com/proceedings/jcis2006/232>
+
 ## Vector cross product
 
-Result is the cross product of x and y, i.e., the resulting components are, in order:
+* VCROSS vd, vs1, vs1
+
+Result is the cross product of x and y.
+
+SUBVL must be set to 3, and all regs must be vectors. VL nonzero produces multiple results in vd.
+
+| Instr | result | src1 | src2 | SUBVL | VL |
+| ------------------ | ------ | ---- | ---- | ----- | -- |
+| VCROSS vd, vs1, vs2 | vec3 | vec3 | vec3 | 3 | any |
+
+The resulting components are, in order:
 
     x[1] * y[2] - y[1] * x[2]
     x[2] * y[0] - y[2] * x[0]
@@ -18,23 +77,101 @@ Result is the cross product of x and y, i.e., the resulting components are, in o
 
 All the operands must be vectors of 3 components of a floating-point type.
 
+Pseudocode:
+
+    vec3 a, b; // elements in order a.x, a.y, a.z
+    // compute a cross b:
+    vec3 t1 = a.yzx; // produce vector [a.y, a.z, a.x]
+    vec3 t2 = b.zxy;
+    vec3 t3 = a.zxy;
+    vec3 t4 = b.yzx;
+    vec3 p = t3 * t4;
+    vec3 cross = t1 * t2 - p;
+
+Assembler:
+
+    fswizzlei,2130 F4, F1
+    fswizzlei,1320 F5, F1
+    fswizzlei,2130 F6, F2
+    fswizzlei,1320 F7, F2
+    fmul F8, F5, F6
+    fmulsub F3, F4, F7, F8
+
 ## Vector dot product
 
-Computes the dot product of two vectors. Internal accuracy must be greater than the
-input vectors and the result.
+* VDOT rd, vs1, vs2
+
+Computes the dot product of two vectors. Internal accuracy must be
+greater than the input vectors and the result.
+
+There are two possible argument options:
+
+* SUBVL=2,3,4 vs1 and vs2 set as vectors,  multiple results are generated. When VL is set, only the first (unpredicated) SUBVector is used to create a result, if rd is scalar (standard behaviour for single predication). Otherwise, if rd is a vector, multiple scalar results are calculated (i.e. SUBVL is always ignored for rd). Swizzling may be applied.
+* When rd=scalar, SUBVL=1 and vs1=vec, vs2=vec, one scalar result is generated from the entire src vectors.  Predication is allowed on the src vectors.
+
+
+| Instr | result | src1 | src2 | SUBVL | VL |
+| ------------------ | ------ | ---- | ---- | ----- | -- |
+| VDOT rd, vs1, vs2 | scal | vec  | vec | 2-4 | any |
+| VDOT rd, vs1, vs2 | scal | vec  | vec | 1 | any |
+
+Pseudocode in python:
+
+    from operator import mul
+    sum(map(mul, A, B))
+
+Pseudocode in c:
+
+    double dot_product(float v[], float u[], int n)
+    {
+        double result = 0.0;
+        for (int i = 0; i < n; i++)
+            result += v[i] * u[i];
+        return result;
+    }
+
+## Vector Normalisation (not included)
+
+Vector normalisation may be performed through dot product, recip square root and multiplication:
+
+    fdot F3, F1, F1 # vector dot with self
+    rcpsqrta F3, F3
+    fscale,0 F2, F3, F1
+
+Or it may be performed through VLEN (Vector length) and division.
 
 ## Vector length
 
+* rd=scalar, vs1=vec (SUBVL=1)
+* rd=scalar, vs1=vec (SUBVL=2,3,4) only 1 (predication rules apply)
+* rd=vec, SUBVL ignored; vs1=vec, SUBVL=2,3,4
+* rd=vec, SUBVL ignored; vs1=vec, SUBVL=1: reserved encoding.
+
+* VLEN rd, vs1
+
 The scalar length of a vector:
 
     sqrt(x[0]^2 + x[1]^2 + ...).
 
+One option is for this to be a macro op fusion sequence, with inverse-sqrt also being a second macro op sequence suitable for normalisation.
+
 ## Vector distance
 
-The scalar distance between two vectors. Subtracts one vector from the other and returns length
+* VDIST rd, vs1, vs2
+
+The scalar distance between two vectors. Subtracts one vector from the
+other and returns length:
+
+    length(v0 - v1)
 
 ## Vector LERP
 
+* VLERP vd, vs1, rs2 # SUBVL=2: vs1.v0 vs1.v1
+
+| Instr | result | src1 | src2 | SUBVL | VL |
+| ------------------ | ------ | ---- | ---- | ----- | -- |
+| VLERP vd, vs1, rs2 | vec2 | vec2 | scal | 2 | any |
+
 Known as **fmix** in GLSL.
 
 <https://en.m.wikipedia.org/wiki/Linear_interpolation>
@@ -56,6 +193,11 @@ Pseudocode:
 
 ## Vector SLERP
 
+* VSLERP vd, vs1, vs2, rs3
+
+Not recommended as it is not commonly used and has several trigonometric
+functions, although CORDIC in vector rotate circular mode is designed for this purpose. Also a costly 4 arg operation.
+
 <https://en.m.wikipedia.org/wiki/Slerp>
 
 Pseudocode:
@@ -100,6 +242,116 @@ Pseudocode:
         return (s0 * v0) + (s1 * v1);
     }
 
+However this algorithm does not involve transcendentals except in
+the computation of the tables: <https://en.wikipedia.org/wiki/CORDIC#Rotation_mode>
+
+    function v = cordic(beta,n)
+        % This function computes v = [cos(beta), sin(beta)] (beta in radians)
+        % using n iterations. Increasing n will increase the precision.
+
+        if beta < -pi/2 || beta > pi/2
+            if beta < 0
+                v = cordic(beta + pi, n);
+            else
+                v = cordic(beta - pi, n);
+            end
+            v = -v; % flip the sign for second or third quadrant
+            return
+        end
+
+        % Initialization of tables of constants used by CORDIC
+        % need a table of arctangents of negative powers of two, in radians:
+        % angles = atan(2.^-(0:27));
+        angles =  [  ...
+            0.78539816339745   0.46364760900081   
+            0.24497866312686   0.12435499454676 ...
+            0.06241880999596   0.03123983343027   
+            0.01562372862048   0.00781234106010 ...
+            0.00390623013197   0.00195312251648   
+            0.00097656218956   0.00048828121119 ...
+            0.00024414062015   0.00012207031189   
+            0.00006103515617   0.00003051757812 ...
+            0.00001525878906   0.00000762939453   
+            0.00000381469727   0.00000190734863 ...
+            0.00000095367432   0.00000047683716   
+            0.00000023841858   0.00000011920929 ...
+            0.00000005960464   0.00000002980232   
+            0.00000001490116   0.00000000745058 ];
+        % and a table of products of reciprocal lengths of vectors [1, 2^-2j]:
+        % Kvalues = cumprod(1./abs(1 + 1j*2.^(-(0:23))))
+        Kvalues = [ ...
+            0.70710678118655   0.63245553203368   
+            0.61357199107790   0.60883391251775 ...
+            0.60764825625617   0.60735177014130   
+            0.60727764409353   0.60725911229889 ...
+            0.60725447933256   0.60725332108988   
+            0.60725303152913   0.60725295913894 ...
+            0.60725294104140   0.60725293651701   
+            0.60725293538591   0.60725293510314 ...
+            0.60725293503245   0.60725293501477   
+            0.60725293501035   0.60725293500925 ...
+            0.60725293500897   0.60725293500890   
+            0.60725293500889   0.60725293500888 ];
+        Kn = Kvalues(min(n, length(Kvalues)));
+
+        % Initialize loop variables:
+        v = [1;0]; % start with 2-vector cosine and sine of zero
+        poweroftwo = 1;
+        angle = angles(1);
+
+        % Iterations
+        for j = 0:n-1;
+            if beta < 0
+                sigma = -1;
+            else
+                sigma = 1;
+            end
+            factor = sigma * poweroftwo;
+            % Note the matrix multiplication can be done using scaling by 
+            % powers of two and addition subtraction
+            R = [1, -factor; factor, 1];
+            v = R * v; % 2-by-2 matrix multiply
+            beta = beta - sigma * angle; % update the remaining angle
+            poweroftwo = poweroftwo / 2;
+            % update the angle from table, or eventually by just dividing by two
+            if j+2 > length(angles)
+                angle = angle / 2;
+            else
+                angle = angles(j+2);
+            end
+        end
+
+        % Adjust length of output vector to be [cos(beta), sin(beta)]:
+        v = v * Kn;
+        return
+
+    endfunction
+
+2x2 matrix multiply can be done with shifts and adds:
+
+    x = v[0] - sigma * (v[1] * 2^(-j));
+    y = sigma * (v[0] * 2^(-j)) + v[1];
+    v = [x; y];
+
+The technique is outlined in a paper as being applicable to 3D:
+<https://www.atlantis-press.com/proceedings/jcis2006/232>
+
 # Expensive 3-operand OP32 operations
 
 3-operand operations are extremely expensive in terms of OP32 encoding space.  A potential idea is to embed 3 RVC register formats across two out of three 5-bit fields rs1/rs2/rd
+
+Another is to overwrite one of the src registers.
+
+# Opcode Table
+
+TODO
+
+# Links
+
+* <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-September/002736.html>
+* <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-September/002733.html>
+* <http://bugs.libre-riscv.org/show_bug.cgi?id=142>
+
+Research Papers
+
+* <https://www.researchgate.net/publication/2938554_PLX_FP_An_Efficient_Floating-Point_Instruction_Set_for_3D_Graphics>