-- RealVecImpl.mesa, an implementation for dealing with vectors of real numbers
-- Last Modified On July 1, 1982 4:47 pm By Paul Rovner
DIRECTORY
PriorityQueue USING[Ref, Predict, Insert, Empty, Remove],
Real USING[FixI, RoundI],
RealVec USING[Handle, Object],
SafeStorage USING[NewZone];
RealVecImpl: MONITOR -- protects sortee
IMPORTS PriorityQueue, Real, SafeStorage
EXPORTS RealVec =
{ OPEN RealVec;
-- PUBLIC SIGNALS
LengthFault: PUBLIC ERROR[vec: Handle, expectedLength: NAT] = CODE;
-- module variables
vecZone: ZONE = SafeStorage.NewZone[];
vecQZone: ZONE = SafeStorage.NewZone[quantized];
-- PUBLIC PROCs
-- nElements <= h.length.
-- Result has averaged values as per some reasonable algorithm
Smooth: PUBLIC PROC[h: Handle, nElements: NAT] RETURNS[ans: Handle] =
{ IF nElements = h.length THEN RETURN[h];
--a, b, c... are in the old space, i, j, k ... in the new
ans ← vecZone.NEW[Object[nElements]];
FOR j: NAT IN [0..nElements)
DO
x: REAL = j;
startT: REAL ← x / nElements;
endT: REAL ← (x + 1.0) / nElements;
a: NAT ← Real.FixI[startT * h.length];
b: NAT ← Real.RoundI[endT * h.length]; -- may be too small by 1
dtOld: REAL ← 1.0/h.length;
dtNew: REAL ← endT - startT;
IF b*dtOld < endT THEN b ← b + 1;
IF a = b THEN ERROR;
ans.elements[j] ← 0.0;
-- now add up the contributions of each old bucket that the new bucket
-- overlaps
FOR z: NAT IN [a..b) -- left edge index (from 0) identifies the bucket
-- value
DO --add in the value in the zth bucket times the fraction of the time
-- covered by the zth bucket in the new bucket
ans.elements[j] ← ans.elements[j]
+ h.elements[z] -- the value in the zth bucket
* (MIN[dtOld, MIN[endT, (z + 1)*dtOld]
- MAX[startT, z*dtOld]])/dtNew
ENDLOOP;
ENDLOOP};
-- convenient constructors
All: PUBLIC PROC[length: NAT, value: REAL ← 0.0] RETURNS[ans: Handle] =
{ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length) DO ans.elements[i] ← value ENDLOOP};
-- h[0] = 0, h[n] = n
IndexVec: PUBLIC PROC[length: NAT] RETURNS[ans: Handle] =
{ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length) DO ans.elements[i] ← i ENDLOOP};
-- convenient composition and extraction operations
SubVec: PUBLIC PROC[h: Handle, firstIndex, lastIndex: NAT]
RETURNS[ans: Handle] =
{length: NAT ← (IF h = NIL THEN 0 ELSE h.length);
IF lastIndex < firstIndex OR h = NIL OR lastIndex >= h.length
THEN ERROR LengthFault[h, length];
length ← lastIndex - firstIndex + 1;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h.elements[firstIndex + i] ENDLOOP};
ConCat: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF LOOPHOLE[LOOPHOLE[h1.length, CARDINAL]
+ LOOPHOLE[h2.length, CARDINAL], INTEGER] < 0
THEN ERROR LengthFault[NIL, 0]; -- overflow check
length ← h1.length + h2.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..h1.length)
DO ans.elements[i] ← h1.elements[i] ENDLOOP;
FOR i: NAT IN [0..h2.length)
DO ans.elements[h1.length + i] ← h2.elements[i] ENDLOOP};
-- h.elements[FixI[perms.elements[i]]] is put into ans.elements[i]
-- e.g., Permute[h, SortOrder[h]] produces a sorted version of h
Permute: PUBLIC PROC[h: Handle, perms: Handle] RETURNS[ans: Handle] =
{ length: NAT ← perms.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO x: NAT = Real.FixI[perms.elements[i]];
IF x >= h.length THEN ERROR LengthFault[h, x+1];
ans.elements[i] ← h.elements[x];
ENDLOOP;
};
sortee: Handle;
sortPred: SAFE PROC[x, y: REF ANY, data: REF ANY]
RETURNS[BOOLEAN] = TRUSTED
{i: NAT = LOOPHOLE--NARROW--[x, REF NAT]↑;
j: NAT = LOOPHOLE--NARROW--[y, REF NAT]↑;
RETURN[sortee.elements[i] < sortee.elements[j]]};
-- perms[i] is the index in h of the ith smallest element
-- i.e. (h.elements[perms[i]] <= h.elements[perms[i+1]])
SortOrder: PUBLIC ENTRY PROC[h: Handle] RETURNS[perms: Handle] =
{ENABLE UNWIND => NULL;
length: NAT ← h.length;
q: PriorityQueue.Ref;
sortee ← h;
q ← PriorityQueue.Predict[pred: sortPred, size: h.length];
perms ← vecQZone.NEW[Object[length]];
FOR i: NAT IN [0..h.length)
DO PriorityQueue.Insert[q, vecQZone.NEW[NAT ← i]] ENDLOOP;
FOR i: NAT IN [0..h.length)
DO perms.elements[i] ← NARROW[PriorityQueue.Remove[q], REF NAT]↑
ENDLOOP;
IF NOT PriorityQueue.Empty[q] THEN ERROR;
};
-- element-by-element operations
Add: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];
length ← h1.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h1.elements[i] + h2.elements[i] ENDLOOP};
Subtract: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];
length ← h1.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h1.elements[i] - h2.elements[i] ENDLOOP};
Multiply: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];
length ← h1.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h1.elements[i] * h2.elements[i] ENDLOOP};
Divide: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];
length ← h1.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h1.elements[i] / h2.elements[i] ENDLOOP};
DotProduct: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: REAL] =
{ length: NAT;
IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];
IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];
length ← h1.length;
ans ← 0.0;
FOR i: NAT IN [0..length)
DO ans ← ans + h1.elements[i] * h2.elements[i] ENDLOOP};
-- scalar element-by-element operations
ScalarAdd: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =
{ length: NAT;
IF h = NIL THEN RETURN[NIL];
length ← h.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h.elements[i] + s ENDLOOP};
ScalarSubtract: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =
{ length: NAT;
IF h = NIL THEN RETURN[NIL];
length ← h.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h.elements[i] - s ENDLOOP};
ScalarMultiply: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =
{ length: NAT;
IF h = NIL THEN RETURN[NIL];
length ← h.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h.elements[i] * s ENDLOOP};
ScalarDivide: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =
{ length: NAT;
IF h = NIL THEN RETURN[NIL];
length ← h.length;
ans ← vecZone.NEW[Object[length]];
FOR i: NAT IN [0..length)
DO ans.elements[i] ← h.elements[i] / s ENDLOOP};
Equal: PUBLIC PROC[h1, h2: Handle] RETURNS[BOOLEAN] =
{ IF h1 = NIL OR h2 = NIL THEN RETURN[h1 = h2];
FOR i: NAT IN [0..h1.length)
DO IF NOT almost[h1.elements[i], h2.elements[i]]
THEN RETURN[FALSE]
ENDLOOP;
RETURN[TRUE]};
almost: PROC[p, q: REAL] RETURNS[BOOLEAN] = INLINE
{RETURN[MAX[p, q] = 0.0 OR ABS[p - q]/MAX[p, q] < 0.00001]};
}.