diff --git a/N2A/src/gov/sandia/n2a/backend/c/JobC.java b/N2A/src/gov/sandia/n2a/backend/c/JobC.java index 5ed7f650..239b8019 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/JobC.java +++ b/N2A/src/gov/sandia/n2a/backend/c/JobC.java @@ -12,6 +12,7 @@ import gov.sandia.n2a.eqset.EquationEntry; import gov.sandia.n2a.eqset.EquationSet; import gov.sandia.n2a.eqset.EquationSet.Conversion; +import gov.sandia.n2a.eqset.EquationSet.ExponentContext; import gov.sandia.n2a.host.Host; import gov.sandia.n2a.host.Remote; import gov.sandia.n2a.host.Windows; @@ -78,8 +79,9 @@ public class JobC extends Thread { protected static Map> runtimeBuilt = new HashMap> (); // collection of Hosts for which runtime has already been checked/built during this session - public MNode job; - protected EquationSet digestedModel; + public MNode job; + protected EquationSet digestedModel; + protected ExponentContext exponentContext; public Host env; public Path localJobDir; @@ -503,6 +505,7 @@ public void rebuildRuntime () throws Exception runtimeBuilt.put (env, runtimes); } CompilerFactory factory = BackendC.getFactory (env); + supportsUnicodeIdentifiers = factory.supportsUnicodeIdentifiers (); String runtimeName = factory.prefixLibrary (shared) + runtimeName () + factory.suffixLibrary (shared); Path runtimeLib = runtimeDir.resolve (runtimeName); for (ProvideOperator pf : extensions) @@ -537,7 +540,6 @@ public void rebuildRuntime () throws Exception if (runtimes.contains (runtimeName)) return; // Compile runtime - supportsUnicodeIdentifiers = factory.supportsUnicodeIdentifiers (); List sources = new ArrayList (); // List of source names sources.add ("runtime"); sources.add ("holder"); @@ -902,7 +904,11 @@ public void digestModel () throws Exception digestedModel.determineTypes (); digestedModel.determineDuration (); digestedModel.assignParents (); - if (fixedPoint) digestedModel.determineExponents (); + if (fixedPoint) + { + exponentContext = new ExponentContext (digestedModel); + digestedModel.determineExponents (exponentContext); + } digestedModel.findConnectionMatrix (); analyzeEvents (digestedModel); analyze (digestedModel); @@ -1246,7 +1252,7 @@ public void generateCode (Path source) throws Exception SIMULATOR = "Simulator<" + T + ">::instance" + (tls ? "->" : "."); - result.append ("#include \"math.h\"\n"); // math.h must always come first, because it messes with mode in which is included. + result.append ("#include \"mymath.h\"\n"); // math.h must always come first, because it messes with mode in which is included. result.append ("#include \"runtime.h\"\n"); if (kokkos) { @@ -1963,7 +1969,7 @@ public void generateMainInitializers (RendererC context) for (Input i : mainInput) { result.append (" " + i.name + " = inputHelper<" + T + "> (\"" + i.operands[0].getString () + "\""); - if (fixedPoint) result.append (", " + i.exponent); + if (fixedPoint) result.append (", " + i.exponent + ", " + i.exponentRow); result.append (");\n"); boolean smooth = i.getKeywordFlag ("smooth"); @@ -2778,7 +2784,7 @@ else if (bed.index != null) result.append (" " + type (bed.n) + " " + mangle (bed.n) + ";\n"); } s.simplify ("$init", bed.globalInit); - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.globalInit); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.globalInit); for (Variable v : bed.globalInit) { multiconditional (v, context, " "); @@ -2799,7 +2805,7 @@ else if (bed.index != null) if (bed.n != null) // and not singleton, so trackN is true { result.append (" resize (" + resolve (bed.n.reference, context, bed.nInitOnly)); - if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent - Operator.MSB)); + if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent)); result.append (");\n"); } } @@ -2828,9 +2834,9 @@ else if (bed.index != null) { result.append (" " + resolve (v.reference, context, false) + " = preserve->" + mangle (v) + " + "); // For fixed-point: - // raw result = exponentDerivative+exponentTime-MSB - // shift = raw-exponentVariable = exponentDerivative+exponentTime-MSB-exponentVariable - int shift = v.derivative.exponent + bed.dt.exponent - Operator.MSB - v.exponent; + // raw result = exponentDerivative+exponentTime + // shift = raw-exponentVariable = exponentDerivative+exponentTime-exponentVariable + int shift = v.derivative.exponent + bed.dt.exponent - v.exponent; if (shift != 0 && fixedPoint) { result.append ("(int) ((int64_t) " + resolve (v.derivative.reference, context, false) + " * dt" + RendererC.printShift (shift) + ");\n"); @@ -2846,7 +2852,7 @@ else if (bed.index != null) for (Variable v : bed.globalIntegrated) { result.append (" " + resolve (v.reference, context, false) + " += "); - int shift = v.derivative.exponent + bed.dt.exponent - Operator.MSB - v.exponent; + int shift = v.derivative.exponent + bed.dt.exponent - v.exponent; if (shift != 0 && fixedPoint) { result.append ("(int) ((int64_t) " + resolve (v.derivative.reference, context, false) + " * dt" + RendererC.printShift (shift) + ");\n"); @@ -2875,7 +2881,7 @@ else if (bed.index != null) result.append (" " + type (v) + " " + mangle ("next_", v) + ";\n"); } s.simplify ("$live", bed.globalUpdate); - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.globalUpdate); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.globalUpdate); for (Variable v : bed.globalUpdate) { multiconditional (v, context, " "); @@ -2899,7 +2905,7 @@ else if (bed.index != null) { // $n may be assigned during the regular update cycle, so we need to monitor it. result.append (" if (" + mangle ("$n") + " != " + mangle ("next_", "$n") + ") " + SIMULATOR + "resize (this, max (0, " + mangle ("next_", "$n")); - if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent - Operator.MSB)); + if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent)); result.append ("));\n"); result.append (" else " + SIMULATOR + "resize (this, -1);\n"); } @@ -2966,7 +2972,7 @@ else if (bed.index != null) returnN = false; } result.append (" " + SIMULATOR + "resize (this, max (0, " + mangle ("$n")); - if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent - Operator.MSB)); + if (context.useExponent) result.append (RendererC.printShift (bed.n.exponent)); result.append ("));\n"); } } @@ -2978,7 +2984,7 @@ else if (bed.index != null) returnN = false; } result.append (" int floorN = max (0, "); - if (context.useExponent) result.append (mangle ("$n") + RendererC.printShift (bed.n.exponent - Operator.MSB)); + if (context.useExponent) result.append (mangle ("$n") + RendererC.printShift (bed.n.exponent)); else result.append ("(int) " + mangle ("$n")); result.append (");\n"); result.append (" if (n != floorN) " + SIMULATOR + "resize (this, floorN);\n"); @@ -3050,7 +3056,7 @@ else if (bed.index != null) result.append (" " + type (v) + " " + mangle ("next_", v) + ";\n"); } s.simplify ("$live", bed.globalDerivativeUpdate); // This is unlikely to make any difference. Just being thorough before call to multiconditional(). - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.globalDerivativeUpdate); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.globalDerivativeUpdate); for (Variable v : bed.globalDerivativeUpdate) { multiconditional (v, context, " "); @@ -3232,7 +3238,7 @@ else if (bed.index != null) result.append (" " + T + " p = c->getP ();\n"); result.append (" if (p <= 0) continue;\n"); result.append (" if (p < 1 && p < uniform<" + T + "> ()"); - if (fixedPoint) result.append (" >> 16"); + if (fixedPoint) result.append (RendererC.printShift (-1 - Operator.MSB - bed.p.exponent)); // shift = exponentUniform - p.exponent = (-1-MSB) - p.exponent result.append (") continue;\n"); result.append (" if (poll && pollSorted.count (c)) continue;\n"); result.append ("\n"); @@ -3830,7 +3836,7 @@ public void generateDefinitionsLocal (RendererC context) throws Exception result.append (" lastT = " + SIMULATOR + "currentEvent->t;\n"); } s.simplify ("$init", bed.localInit); - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.localInit); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.localInit); if (bed.localInit.contains (bed.dt)) { result.append (" EventStep<" + T + "> * event = getEvent ();\n"); @@ -3934,7 +3940,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant for (Variable v : bed.localIntegrated) { result.append (" " + resolve (v.reference, context, false) + " = preserve->" + mangle (v) + " + "); - int shift = v.derivative.exponent + bed.dt.exponent - Operator.MSB - v.exponent; + int shift = v.derivative.exponent + bed.dt.exponent - v.exponent; if (shift != 0 && fixedPoint) { if (v.type instanceof Matrix) @@ -3958,7 +3964,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant for (Variable v : bed.localIntegrated) { result.append (pad + " " + resolve (v.reference, context, false) + " += "); - int shift = v.derivative.exponent + bed.dt.exponent - Operator.MSB - v.exponent; + int shift = v.derivative.exponent + bed.dt.exponent - v.exponent; if (shift != 0 && fixedPoint) { if (v.type instanceof Matrix) @@ -4003,7 +4009,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant result.append (" " + type (v) + " " + mangle ("next_", v) + ";\n"); } s.simplify ("$live", bed.localUpdate); - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.localUpdate); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.localUpdate); for (Variable v : bed.localUpdate) { multiconditional (v, context, " "); @@ -4180,12 +4186,12 @@ else if (bed.setDt) // implies that bed.dt exists and is constant result.append (" if (pow (" + resolve (bed.p.reference, context, false) + ", " + resolve (bed.dt.reference, context, false)); if (context.useExponent) { - result.append (RendererC.printShift (bed.dt.exponent - 15)); // second operand must have exponent=15 + result.append (RendererC.printShift (bed.dt.exponent + Operator.MSB / 2)); // Second operand must have exponent=-MSB/2. shift = dt.exponent - -MSB/2 result.append (", " + bed.p.exponent); // exponentA result.append (", " + bed.p.exponent); // exponentResult } result.append (") < uniform<" + T + "> ()"); - if (context.useExponent) result.append (RendererC.printShift (-1 - bed.p.exponent)); // -1 is hard-coded from the Uniform function. + if (context.useExponent) result.append (RendererC.printShift (-1 - Operator.MSB - bed.p.exponent)); // shift = exponentUniform - p.exponent = (-1-MSB) - p.exponent result.append (")\n"); } } @@ -4198,7 +4204,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant for (Variable t : s.ordered) if (t.hasAttribute ("temporary") && bed.p.dependsOn (t) != null) list.add (t); list.add (bed.p); s.simplify ("$live", list, bed.p); - if (fixedPoint) EquationSet.determineExponentsSimplified (list); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, list); for (Variable v : list) { multiconditional (v, context, " "); @@ -4208,12 +4214,12 @@ else if (bed.setDt) // implies that bed.dt exists and is constant result.append (" if (" + mangle ("$p") + " <= 0 || " + mangle ("$p") + " < " + context.print (1, bed.p.exponent) + " && pow (" + mangle ("$p") + ", " + resolve (bed.dt.reference, context, false)); if (context.useExponent) { - result.append (RendererC.printShift (bed.dt.exponent - 15)); + result.append (RendererC.printShift (bed.dt.exponent + Operator.MSB / 2)); result.append (", " + bed.p.exponent); result.append (", " + bed.p.exponent); } result.append (") < uniform<" + T + "> ()"); - if (context.useExponent) result.append (RendererC.printShift (-1 - bed.p.exponent)); + if (context.useExponent) result.append (RendererC.printShift (-1 - Operator.MSB - bed.p.exponent)); result.append (")\n"); } result.append (" {\n"); @@ -4269,7 +4275,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant result.append (" " + type (v) + " " + mangle ("next_", v) + ";\n"); } s.simplify ("$live", bed.localDerivativeUpdate); - if (fixedPoint) EquationSet.determineExponentsSimplified (bed.localDerivativeUpdate); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, bed.localDerivativeUpdate); for (Variable v : bed.localDerivativeUpdate) { multiconditional (v, context, " "); @@ -4575,7 +4581,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant } list.add (project); s.simplify ("$connect", list, project); - if (fixedPoint) EquationSet.determineExponentsSimplified (list); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, list); for (Variable v : list) { multiconditional (v, context, " "); @@ -4712,7 +4718,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant } list.add (bed.p); s.simplify ("$connect", list, bed.p); - if (fixedPoint) EquationSet.determineExponentsSimplified (list); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, list); for (Variable v : list) { multiconditional (v, context, " "); @@ -4740,7 +4746,7 @@ else if (bed.setDt) // implies that bed.dt exists and is constant for (Variable t : s.ordered) if (t.hasAttribute ("temporary") && bed.xyz.dependsOn (t) != null) list.add (t); list.add (bed.xyz); s.simplify ("$live", list, bed.xyz); // evaluate in $live phase, because endpoints already exist when connection is evaluated. - if (fixedPoint) EquationSet.determineExponentsSimplified (list); + if (fixedPoint) EquationSet.determineExponentsSimplified (exponentContext, list); for (Variable v : list) { multiconditional (v, context, " "); @@ -4888,9 +4894,9 @@ else if (bed.setDt) // implies that bed.dt exists and is constant case Variable.DIVIDE: { // The current and buffered values of the variable have the same exponent. - // raw = exponentV + exponentV - MSB - // shift = raw - exponentV = exponentV - MSB - int shift = v.exponent - Operator.MSB; + // raw = exponentV + exponentV + // shift = raw - exponentV = exponentV + int shift = v.exponent; if (shift != 0 && fixedPoint) { result.append (" = (int64_t) " + current + " * " + buffered + RendererC.printShift (shift) + ";\n"); @@ -5339,9 +5345,9 @@ public void renderEquation (RendererC context, EquationEntry e) result.append (" += "); break; case Variable.MULTIPLY: - // raw exponent = exponentV + exponentExpression - MSB - // shift = raw - exponentV = expnentExpression - MSB - shift = e.expression.exponentNext - Operator.MSB; + // raw exponent = exponentV + exponentExpression + // shift = raw - exponentV = expnentExpression + shift = e.expression.exponentNext; if (shift != 0 && fixedPoint) { if (shift < 0) result.append (" = (int64_t) " + LHS + " * "); @@ -5353,9 +5359,9 @@ public void renderEquation (RendererC context, EquationEntry e) } break; case Variable.DIVIDE: - // raw = exponentV - exponentExpression + MSB - // shift = raw - exponentV = MSB - exponentExpression - shift = Operator.MSB - e.expression.exponentNext; + // raw = exponentV - exponentExpression + // shift = raw - exponentV = -exponentExpression + shift = -e.expression.exponentNext; if (shift != 0 && fixedPoint) { if (shift > 0) result.append (" = ((int64_t) " + LHS + RendererC.printShift (shift) + ") / "); @@ -5565,7 +5571,7 @@ public boolean visit (Operator op) bed.defined.add (v); } context.result.append (pad + "InputHolder<" + T + "> * " + i.name + " = inputHelper<" + T + "> (" + i.fileName); - if (fixedPoint) context.result.append (", " + i.exponent); + if (fixedPoint) context.result.append (", " + i.exponent + ", " + i.exponentRow); context.result.append (");\n"); boolean smooth = i.getKeywordFlag ("smooth"); @@ -5756,7 +5762,7 @@ public boolean visit (Operator op) context.result.append (pad + name + " = "); value.render (context); context.result.append (";\n"); - if (fixedPoint) context.result.append (pad + name + " /= powf (2, FP_MSB - " + value.exponent + ");\n"); + if (fixedPoint) context.result.append (pad + name + " /= powf (2, - " + value.exponent + ");\n"); continue; // The model matrix will get coerced to float in the call to draw. @@ -5799,7 +5805,7 @@ else if (material) context.result.append (" = "); if (fixedPoint && convertToFloat) context.result.append ("("); value.render (context); - if (fixedPoint && convertToFloat) context.result.append (") / powf (2, FP_MSB - " + value.exponent + ")"); + if (fixedPoint && convertToFloat) context.result.append (") / powf (2, - " + value.exponent + ")"); context.result.append (";\n"); } } diff --git a/N2A/src/gov/sandia/n2a/backend/c/RendererC.java b/N2A/src/gov/sandia/n2a/backend/c/RendererC.java index 8a362de9..0cade8f2 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/RendererC.java +++ b/N2A/src/gov/sandia/n2a/backend/c/RendererC.java @@ -164,8 +164,7 @@ public boolean render (Operator op) } else if (useExponent) { - int shiftX = Operator.MSB - y.exponent; - int x = shiftX >= 0 ? 0x1 << shiftX : 0; + int x = y.exponent >= 0 ? 0x1 << y.exponent : 0; result.append (", " + x); } } @@ -249,14 +248,14 @@ else if (useExponent) if (op instanceof Draw3D) { boolean needModel = ((Draw3D) d).needModelMatrix (); - int exponentP = Operator.MSB; + int exponentP = 0; if (useExponent) { // position and model matrix share same exponent Operator model = d.getKeyword ("model"); if (model != null) exponentP = model.exponentNext; else if (d.operands.length > 1) exponentP = d.operands[1].exponentNext; - // else model has been initialized to identity, presumably integer 1. Operator.MSB set above is suitable default. + // else model has been initialized to identity, presumably integer 1. Zero set above is suitable default. } result.append (d.name + "->" + op + " (" + job.SIMULATOR + "currentEvent->t, material"); diff --git a/N2A/src/gov/sandia/n2a/backend/c/RendererCfp.java b/N2A/src/gov/sandia/n2a/backend/c/RendererCfp.java index 839c8877..0688e575 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/RendererCfp.java +++ b/N2A/src/gov/sandia/n2a/backend/c/RendererCfp.java @@ -192,9 +192,8 @@ public boolean render (Operator op) OperatorBinary b = (OperatorBinary) op; // Explanation of shift -- The exponent of the result will be the sum of the exponents - // of the two operands. That new exponent will be associated with bit position 2*MSB. - // We want the exponent at bit position MSB. - int exponentRaw = b.operand0.exponentNext + b.operand1.exponentNext - Operator.MSB; // Exponent at MSB position after a direct integer multiply. + // of the two operands. + int exponentRaw = b.operand0.exponentNext + b.operand1.exponentNext; // Exponent at LSB position after a direct integer multiply. int shift = exponentRaw - b.exponentNext; if (shift == 0) @@ -266,8 +265,8 @@ public boolean render (Operator op) // Explanation of shift -- In a division, the quotient is effectively down-shifted by // the number of bits in the denominator, and its exponent is the difference between // the exponents of the numerator and denominator. In A/B, exponentA-exponentB will - // be the exponent at bit position 0. We want the exponent at bit position MSB. - int exponentRaw = d.operand0.exponentNext - d.operand1.exponentNext + Operator.MSB; // Exponent in MSB from a direct integer division. + // be the exponent at bit position 0. + int exponentRaw = d.operand0.exponentNext - d.operand1.exponentNext; // Exponent in LSB from a direct integer division. int shift = exponentRaw - d.exponentNext; if (shift == 0) @@ -395,13 +394,13 @@ public boolean render (Operator op) { Ceil c = (Ceil) op; Operator a = c.operands[0]; - if (a.exponent >= Operator.MSB) // LSB is above decimal, so ceil() operation is impossible. This test must be on a.exponent rather than a.exponentNext because ... + if (a.exponent >= 0) // LSB is above decimal, so ceil() operation is impossible. This test must be on a.exponent rather than a.exponentNext because ... { // a.exponentNext has been set to f.exponentNext by Round.determineExponentNextStatic(), // so no shift is necessary. a.render (this); } - else // a.exponentNext is in [0, MSB), as enforced by Round.determineExponentNextStatic() + else // a.exponentNext is in [-MSB, -1], as enforced by Round.determineExponentNextStatic() { // Create a mask for bits below the decimal point. // When this mask is added to the number, it will add 1 to the first bit position @@ -409,8 +408,7 @@ public boolean render (Operator op) // any residual bits below the decimal. // This works for both positive and negative numbers. int shift = a.exponentNext - c.exponentNext; - int decimalPlaces = Operator.MSB - a.exponentNext; - int wholeMask = 0xFFFFFFFF << decimalPlaces; + int wholeMask = 0xFFFFFFFF << -a.exponentNext; int decimalMask = ~wholeMask; if (shift != 0) result.append ("("); @@ -456,8 +454,7 @@ public boolean render (Operator op) if (op.parent == null || op.parent instanceof OperatorLogicalInput) return super.render (op); Event e = (Event) op; - int exponentRaw = Operator.MSB - e.eventType.valueIndex; - int shift = exponentRaw - e.exponentNext; + int shift = -e.eventType.valueIndex - e.exponentNext; if (shift != 0) result.append ("("); result.append ("(flags & (" + bed.localFlagType + ") 0x1 << " + e.eventType.valueIndex + ")"); @@ -473,7 +470,7 @@ public boolean render (Operator op) // See Ceil above for similar code. Floor f = (Floor) op; Operator a = f.operands[0]; - if (a.exponent >= Operator.MSB) + if (a.exponent >= 0) { a.render (this); } @@ -481,8 +478,7 @@ public boolean render (Operator op) { // Mask off bits below the decimal point. This works for both positive and negative numbers. int shift = a.exponentNext - f.exponentNext; - int decimalPlaces = Operator.MSB - a.exponentNext; - int wholeMask = 0xFFFFFFFF << decimalPlaces; + int wholeMask = 0xFFFFFFFF << -a.exponentNext; if (shift != 0) result.append ("("); result.append ("("); @@ -528,14 +524,14 @@ public boolean render (Operator op) // See Ceil above for similar code. Round r = (Round) op; Operator a = r.operands[0]; - if (a.exponent >= Operator.MSB) + if (a.exponent >= 0) { a.render (this); } - else // a.exponentNext in [0, MSB) + else // a.exponentNext in [-MSB, -1] { int shift = a.exponentNext - r.exponentNext; - int decimalPlaces = Operator.MSB - a.exponentNext; + int decimalPlaces = -a.exponentNext; int mask = 0xFFFFFFFF << decimalPlaces; int half = 0x1 << decimalPlaces - 1; @@ -555,7 +551,7 @@ public boolean render (Operator op) { Signum s = (Signum) op; Operator a = s.operands[0]; - int one = 0x1 << Operator.MSB - s.exponentNext; + int one = 0x1 << -s.exponentNext; boolean needParens = a.precedence () >= new LT ().precedence (); if (needParens) result.append ("("); a.render (this); @@ -609,14 +605,14 @@ public boolean provides64bit (Operator op) if (op instanceof Multiply) { OperatorBinary b = (OperatorBinary) op; - int exponentRaw = b.operand0.exponentNext + b.operand1.exponentNext - Operator.MSB; + int exponentRaw = b.operand0.exponentNext + b.operand1.exponentNext; int shift = exponentRaw - b.exponentNext; return shift < 0; // only use 64-bit when we downshift the raw result } if (op instanceof Divide) { OperatorBinary b = (OperatorBinary) op; - int exponentRaw = b.operand0.exponentNext - b.operand1.exponentNext + Operator.MSB; + int exponentRaw = b.operand0.exponentNext - b.operand1.exponentNext; int shift = exponentRaw - b.exponentNext; return shift > 0; // only use 64-bit when we upshift the raw result } @@ -665,7 +661,7 @@ public String print (double d, int exponent) bits |= 0x10000000000000l; // set implied msb of mantissa (bit 52) to 1 bits &= 0x1FFFFFFFFFFFFFl; // clear sign and exponent bits if (negate) bits = -bits; - bits >>= 52 - Operator.MSB + exponent - e; + bits >>= 52 + exponent - e; return Integer.toString ((int) bits); } } diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/fixedpoint.cc b/N2A/src/gov/sandia/n2a/backend/c/runtime/fixedpoint.cc index 90e2f9b3..24e4871f 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/fixedpoint.cc +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/fixedpoint.cc @@ -70,75 +70,95 @@ identity (const MatrixStrided & A, int one) int norm (const MatrixStrided & A, int n, int exponentA, int exponentResult) { - const int exponentN = 15; + const int exponentN = -FP_MSB2; int * a = A.base (); int * end = a + A.rows () * A.columns (); // Assumes MatrixFixed, so that elements are dense in memory. - int result = 0; - if (n == INFINITY) + if (n == 0 || n == INFINITY) { - while (a < end) result = std::max (std::abs (*a++), result); - int shift = exponentA - exponentResult; - if (shift > 0) return result << shift; - if (shift < 0) return result >> -shift; - return result; - } - if (n == 0) - { - while (a < end) if (*a++) result++; - int shift = FP_MSB - exponentResult; + int result = 0; + int shift = -exponentResult; + if (n) // INFINITY + { + while (a < end) result = std::max (std::abs (*a++), result); + shift += exponentA; + } + else // 0 + { + while (a < end) if (*a++) result++; + } if (shift > 0) return result << shift; if (shift < 0) return result >> -shift; return result; } - if (n == 0x1 << exponentN) + + uint64_t sum = 0; + if (n == 0x1 << -exponentN) { - while (a < end) result += std::abs (*a++); + while (a < end) sum += std::abs (*a++); int shift = exponentA - exponentResult; - if (shift > 0) return result << shift; - if (shift < 0) return result >> -shift; - return result; + if (shift > 0) sum <<= shift; + if (shift < 0) sum >>= -shift; + if (sum > INFINITY) return INFINITY; + return sum; } // Fully general form - // "result" will hold the sum, and exponentA will hold exponentSum when done. - int root; // exponent=15 - if (n == 0x2 << exponentN) + // exponentA will hold exponentSum when done. + int root; // exponent=-MSB/2 + if (n == 0x2 << -exponentN) { - root = 0x4000; // 0.5 - exponentA = exponentA * 2 - FP_MSB; // raw result of squaring elements of A - register uint64_t sum = 0; + root = 0x1 << FP_MSB2 - 1; // 0.5 + exponentA = exponentA * 2; // raw result of squaring elements of A while (a < end) { int t = *a++; sum += (int64_t) t * t; } - while (sum > INFINITY) - { - sum >>= 1; - exponentA++; - } - result = sum; // truncate to 32 bits } - else + else // fractional power or power > 2 { // for root: - // raw division = exponentOne-exponentN+MSB = MSB-MSB/2+MSB - // want exponentN, so shift = raw-exponentN = (MSB-MSB/2+MSB)-MSB/2 = MSB + // raw division = exponentOne - exponentN = 0 - -MSB/2 = MSB/2 + // want exponentN, so shift = raw - exponentN = MSB/2 - -MSB/2 = MSB root = (0x1 << FP_MSB) / n; + // Estimate center bit position + // Ideally this would be the median of MSB positions, but that + // requires sorting. Instead, we compute the average. + int count = 0; + int center = 0; + int * t = a; + while (t < end) + { + int temp = std::abs (*t++); + if (! temp) continue; + count++; + while (temp) + { + temp >>= 1; + center++; + } + } + if (count) center /= count; + else center = FP_MSB2; // Though, in this case it doesn't matter because matrix is all zeroes. + // for exponentSum: - // assume center of A = MSB/2 - // center power of A = centerA = exponentA - MSB/2 - // center power of one term = centerTerm = centerA*n - // want center of term at MSB/2, so exponentSum = centerTerm+MSB/2 = (exponentA-MSB/2)*n+MSB/2 = exponentA*n+(1-n)*MSB/2 - int exponentSum = ((exponentA - FP_MSB2) * n >> exponentN) + FP_MSB2; + // centerA = center power of A = exponentA + center + // centerTerm = center power of one term = centerA*n + // want center of term at MSB/2, so exponentSum = centerTerm - MSB/2 = (exponentA+center)*n - MSB/2 + int exponentSum = ((exponentA + center) * n >> -exponentN) - FP_MSB2; - while (a < end) result += pow (std::abs (*a++), n, exponentA, exponentSum); + while (a < end) sum += pow (std::abs (*a++), n, exponentA, exponentSum); exponentA = exponentSum; } - return pow (result, root, exponentA, exponentResult); + while (sum > INFINITY) + { + sum >>= 1; + exponentA++; + } + return pow ((int) sum, root, exponentA, exponentResult); } Matrix @@ -146,15 +166,16 @@ normalize (const MatrixStrided & A, int exponentA) { // Calculate 2-norm of A // Allow for magnitude of "scale" to be larger than the magnitude of individual elements. - int count = norm (A, 0, exponentA, FP_MSB); // Number of nonzero elements + int count = norm (A, 0, exponentA, 0); // Number of nonzero elements int bits = 0; while (count >>= 1) bits++; int exponentScale = exponentA + bits; - int scale = norm (A, 0x2 << 15, exponentA, exponentScale); // 2-norm + int scale = norm (A, 0x2 << FP_MSB2, exponentA, exponentScale); // 2-norm // Divide A - // Goal is for result to be at exponent=0 - // See comments on Divide in RendererCfp + // raw = exponentA - exponentScale + // goal = -MSB = everything in [0,1] + // shift = raw - goal = exponentA - exponentScale - -MSB int shift = exponentA - exponentScale + FP_MSB; return divide (A, scale, shift); } @@ -466,17 +487,17 @@ glFrustum (int left, int right, int bottom, int top, int near, int far, int expo Matrix result (4, 4); clear (result); - // After division, raw = exponent - exponent + MSB = MSB - // Goal is to shift back to original exponent. shift = raw - exponent = MSB - exponent - int shift = FP_MSB - exponent; + // After division, raw = exponent - exponent = 0 + // Goal is to shift back to original exponent. shift = raw - exponent = -exponent + int shift = -exponent; result(0,0) = ( (int64_t) 2 * near << shift) / (right - left); result(1,1) = ( (int64_t) 2 * near << shift) / (top - bottom); result(0,2) = ( (int64_t) right + left << shift) / (right - left); result(1,2) = ( (int64_t) top + bottom << shift) / (top - bottom); result(2,2) = (-((int64_t) far + near) << shift) / (far - near); - // shift = MSB - exponent + // shift = exponentOne - exponent = 0 - exponent result(3,2) = -1 << shift; - // raw = (exponent + exponent - MSB) - exponent + MSB = exponent + // raw = (exponent + exponent) - exponent = exponent // shift = raw - exponent = 0 result(2,3) = (int64_t) -2 * far * near / (far - near); @@ -489,19 +510,19 @@ glOrtho (int left, int right, int bottom, int top, int near, int far, int expone Matrix result (4, 4); clear (result); - // raw = MSB - exponent + MSB = 2*MSB - exponent - // shift = raw - exponent = 2*MSB - 2*exponent - int shift = 2 * (FP_MSB - exponent); + // raw = exponentOne - exponent = 0 - exponent = -exponent + // shift = raw - exponent = -2*exponent + int shift = -2 * exponent; result(0,0) = ((int64_t) 2 << shift) / (right - left); result(1,1) = ((int64_t) 2 << shift) / (top - bottom); result(2,2) = ((int64_t) -2 << shift) / (far - near); - // raw = exponent - exponent + MSB = MSB - // shift = raw - exponent = MSB - exponent - shift = FP_MSB - exponent; + // raw = exponent - exponent = 0 + // shift = raw - exponent = -exponent + shift = -exponent; result(0,3) = (-((int64_t) right + left ) << shift) / (right - left); result(1,3) = (-((int64_t) top + bottom) << shift) / (top - bottom); result(2,3) = (-((int64_t) far + near ) << shift) / (far - near); - // shift = MSB - exponent + // shift = exponentOne - exponent = -exponent result(3,3) = 1 << shift; return result; @@ -512,13 +533,13 @@ glLookAt (const MatrixFixed & eye, const MatrixFixed & center, { // Create an orthonormal frame Matrix f = center - eye; - f = normalize (f, exponent); // f exponent=0 - Matrix u = normalize (up, exponent); // u exponent=0 - Matrix s = cross (f, u, FP_MSB); // s exponent=0; but s is not necessarily unit length - s = normalize (s, 0); + f = normalize (f, exponent); // f exponent=-MSB + Matrix u = normalize (up, exponent); // u exponent=-MSB + Matrix s = cross (f, u, FP_MSB); // s exponent=-MSB; but s is not necessarily unit length + s = normalize (s, -FP_MSB); u = cross (s, f, FP_MSB); - Matrix R (4, 4); // R exponent=0 + Matrix R (4, 4); // R exponent=-MSB clear (R); R(0,0) = s[0]; R(0,1) = s[1]; @@ -532,12 +553,12 @@ glLookAt (const MatrixFixed & eye, const MatrixFixed & center, R(3,3) = 1 << FP_MSB; Matrix Tr (4, 4); // Tr has the exponent passed to this function - identity (Tr, 1 << FP_MSB - exponent); + identity (Tr, 1 << -exponent); Tr(0,3) = -eye[0]; Tr(1,3) = -eye[1]; Tr(2,3) = -eye[2]; - // raw = 0 + exponent - MSB + // raw = -MSB + exponent // goal = exponent // shift = raw - goal = -MSB return multiply (R, Tr, FP_MSB); @@ -546,29 +567,29 @@ glLookAt (const MatrixFixed & eye, const MatrixFixed & center, Matrix glPerspective (int fovy, int aspect, int near, int far, int exponent) { - // raw = (exponent + 1 - MSB) - MSB + MSB = exponent + 1 - MSB - // goal = 1, same as M_PI - // shift = raw - goal = exponent - MSB - int shift = exponent - FP_MSB; + // raw = (exponent + 1-MSB) - 0 = exponent+1-MSB + // goal = 1-MSB, same as M_PI + // shift = raw - goal = exponent+1-MSB - (1-MSB) = exponent + int shift = exponent; fovy = ::shift ((int64_t) fovy * M_PI / 180, shift); - // raw = MSB - 3 + MSB = 2*MSB - 3 + // raw = exponentOne - exponentTan = 0 - 3-MSB = MSB-3 // goal = exponent - // shift = raw - goal = 2*MSB - 3 - exponent - shift = 2 * FP_MSB - 3 - exponent; - int f = ((int64_t) 1 << shift) / tan (fovy / 2, 1, 3); // tan() goes to infinity, but 8 (2^3) should be sufficient for almost all cases. + // shift = raw - goal = MSB-3 - exponent + shift = FP_MSB - 3 - exponent; + int f = ((int64_t) 1 << shift) / tan (fovy / 2, 1-FP_MSB, 3-FP_MSB); // tan() goes to infinity, but 8 (2^3) should be sufficient for almost all cases. Matrix result (4, 4); clear (result); - // raw = exponent - exponent + MSB = MSB + // raw = exponent - exponent = 0 // goal = exponent - // shift = raw - goal = MSB - exponent - shift = FP_MSB - exponent; + // shift = raw - goal = -exponent + shift = -exponent; result(0,0) = ((int64_t) f << shift) / aspect; result(1,1) = f; result(2,2) = ((int64_t) far + near << shift) / (near - far); result(3,2) = -1 << FP_MSB - exponent; - // raw = (exponent + exponent - MSB) - exponent + MSB = exponent + // raw = (exponent + exponent) - exponent = exponent result(2,3) = (int64_t) 2 * far * near / (near - far); return result; @@ -583,42 +604,43 @@ glRotate (int angle, const MatrixFixed & axis, int exponent) Matrix glRotate (int angle, int x, int y, int z, int exponent) { - // raw = (exponent + 1 - MSB) - MSB + MSB = exponent + 1 - MSB - // goal = 1, same as M_PI - // shift = raw - goal = exponent - MSB - int shift = exponent - FP_MSB; + // raw = (exponent + 1-MSB) - 0 = exponent + 1 - MSB + // goal = 1-MSB, same as M_PI + // shift = raw - goal = exponent + int shift = exponent; angle = ::shift ((int64_t) angle * M_PI / 180, shift); - // c, s and c1 all have exponent 1 + // c, s and c1 all have exponent 1-MSB int c = cos (angle, 1); int s = sin (angle, 1); int c1 = (1 << FP_MSB - 1) - c; // normalize([x y z]) - // raw = exponent + exponent - MSB + // raw = exponent + exponent // result = exponent + 2 bits of headroom for additions - int l = sqrt ((int64_t) x * x + (int64_t) y * y + (int64_t) z * z, 2 * exponent - FP_MSB, exponent + 2); - // raw = exponent - (exponent + 2) + MSB = MSB - 2 - // goal = 0 + int l = sqrt ((int64_t) x * x + (int64_t) y * y + (int64_t) z * z, 2 * exponent, exponent + 2); + // raw = exponent - (exponent + 2) = -2 + // goal = -MSB + // shift = raw - goal = -2 - -MSB = MSB-2 shift = FP_MSB - 2; x = ((int64_t) x << shift) / l; y = ((int64_t) y << shift) / l; z = ((int64_t) z << shift) / l; - // exponentResult=0 + // exponentResult=-MSB Matrix result (4, 4); clear (result); - // raw = (0 + 0 - MSB) + 1 - MSB = -2*MSB + 1 - // goal = 1 to match c - // shift = -2 * MSB, applied in two stages + // raw = -MSB + -MSB + 1-MSB = -3*MSB + 1 + // goal = 1-MSB to match c + // shift = -2*MSB, applied in two stages // Then we need one bit upshift to match exponentResult result(0,0) = (((int64_t) x * x >> FP_MSB) * c1 >> FP_MSB) + c << 1; result(1,1) = (((int64_t) y * y >> FP_MSB) * c1 >> FP_MSB) + c << 1; result(2,2) = (((int64_t) z * z >> FP_MSB) * c1 >> FP_MSB) + c << 1; result(3,3) = 1 << FP_MSB; // For second term: - // raw = 0 + 1 - MSB - // goal = 1 + // raw = -MSB + 1-MSB + // goal = 1-MSB result(1,0) = (((int64_t) y * x >> FP_MSB) * c1 >> FP_MSB) + ((int64_t) z * s >> FP_MSB) << 1; result(2,0) = (((int64_t) x * z >> FP_MSB) * c1 >> FP_MSB) - ((int64_t) y * s >> FP_MSB) << 1; result(0,1) = (((int64_t) x * y >> FP_MSB) * c1 >> FP_MSB) - ((int64_t) z * s >> FP_MSB) << 1; @@ -643,7 +665,7 @@ glScale (int sx, int sy, int sz, int exponent) result(0,0) = sx; result(1,1) = sy; result(2,2) = sz; - result(3,3) = 1 << FP_MSB - exponent; + result(3,3) = 1 << -exponent; return result; } @@ -657,22 +679,22 @@ Matrix glTranslate (int x, int y, int z, int exponent) { Matrix result (4, 4); - identity (result, 1 << FP_MSB - exponent); + identity (result, 1 << -exponent); result(0,3) = x; result(1,3) = y; result(2,3) = z; return result; } -// exponentResult = 1, to accommodate [-pi, pi] +// exponentResult = 1-MSB, to accommodate [-pi, pi] // exponentY == exponentX, but it doesn't matter what the exponent is; only ratio matters int atan2 (int y, int x) { - // Using CORDIC algorithm. See https://www.mathworks.com/help/fixedpoint/ug/calculate-fixed-point-arctangent.html + // Using the CORDIC algorithm. See https://www.mathworks.com/help/fixedpoint/ug/calculate-fixed-point-arctangent.html // Look-up table for values of atan(2^-i), i=0,1,2,... - // Converted to fixed-point with exponent=1 (same as result of this function). + // Converted to fixed-point with exponent=1-MSB (same as result of this function). // Limited to 12 terms, as a compromise between accuracy and time+space cost. // The Mathworks article discusses these tradeoffs. static const int lut[] = {421657428, 248918914, 131521918, 66762579, 33510843, 16771757, 8387925, 4194218, 2097141, 1048574, 524287, 262143}; //, 131071, 65535, 32767, 16383, 8191, 4095, 2047, 1023, 511, 255, 127, 63, 31, 15, 7, 4, 2, 1}; @@ -717,6 +739,9 @@ atan2 (int y, int x) if (x >> 4 >= y) // Use small-angle formula. { + // raw = 0 + // goal = 1-MSB + // shift = 0 - (1-MSB) = MSB-1 result += (int) (((int64_t) y << FP_MSB - 1) / x); } else // Use CORDIC @@ -758,10 +783,9 @@ int ceil (int a, int exponentA, int exponentResult) { int result; - if (exponentA >= 0 && exponentA < FP_MSB) + if (exponentA >= -FP_MSB && exponentA < 0) { - int decimalPlaces = FP_MSB - exponentA; - int wholeMask = 0xFFFFFFFF << decimalPlaces; + int wholeMask = 0xFFFFFFFF << -exponentA; int decimalMask = ~wholeMask; result = a + decimalMask & wholeMask; } @@ -779,32 +803,37 @@ ceil (int a, int exponentA, int exponentResult) int cos (int a, int exponentA) { - // We want to add PI/2 to a. M_PI exponent=1. To induce down-shift, claim it is 0. - // Thus, shift is exactly exponentA. - if (exponentA >= 0) return sin (a + (M_PI >> exponentA), exponentA); - // If exponentA is negative, then a is too small to use as-is. - if (exponentA < -FP_MSB) return 0x20000000; // one, with exponent=1 - return sin ((a >> -exponentA) + M_PI, 0); + // We want to add PI/2 to a. M_PI exponent=1-MSB. To induce down-shift, claim it is -MSB. + // current = -MSB (formerly 1-MSB) + // goal = exponentA + // shift = current - goal = exponentA + MSB + if (exponentA >= -FP_MSB) return sin (a + (M_PI >> exponentA + FP_MSB), exponentA); + + // If a is too small, then result is basically 1. + if (exponentA < -2 * FP_MSB) return 0x20000000; // one, with exponent=1-MSB + + // Shift a to match M_PI. Strategy for PI/2 remains the same. + return sin ((a >> -exponentA - FP_MSB) + M_PI, -FP_MSB); } int exp (int a, int exponentResult) { - const int exponentA = 7; // Hard-coded value established in gov.sandia.n2a.language.function.Exp class. + const int exponentA = 7 - FP_MSB; // Hard-coded value established in gov.sandia.n2a.language.function.Exp class. if (a == 0) { - int shift = FP_MSB - exponentResult; + int shift = -exponentResult; if (shift < 0) return 0; return 1 << shift; } - const int one = 1 << FP_MSB - exponentA; - if (a == one) + const int one = 1 << -exponentA; + if (a == one) // special case for returning e, the natural logarithm base. { - int shift = 1 - exponentResult; // M_E exponent=1 + int shift = 1 - FP_MSB - exponentResult; // M_E exponent=1-MSB + if (shift == 0) return M_E; if (shift < 0) return M_E >> -shift; - if (shift > 0) return INFINITY; // Up-shifting M_E is nonsense, since it already uses all the bits. - return M_E; + return INFINITY; // Up-shifting M_E is nonsense, since it already uses all the bits. } // Algorithm: @@ -820,13 +849,13 @@ exp (int a, int exponentResult) int exponentWork = exponentA; // Shift for inner loop: - // i has exponent=MSB - // a has exponent=7 per above comment + // i has exponent=0 + // a has exponent=7-MSB per above comment // term and result have exponentWork - // raw multiply = exponentA+exponentWork at bit 60 - // raw divide = (exponentA+exponentWork)-MSB at bit 30 - // We want exponentWork at bit 30, so shift = raw-exponentWork = (exponentA+exponentWork-MSB)-exponentWork = exponentA-MSB - const int shift = FP_MSB - exponentA; // preemtively flip the sign, since this is always used in the positive form + // raw = (exponentA + exponentWork) - 0 + // goal = exponentWork + // shift = raw - goal = exponentA + const int shift = -exponentA; // preemtively flip the sign, since this is always used in the positive form const int round = 1 << shift - 1; const int maximum = 1 << FP_MSB; @@ -848,13 +877,15 @@ exp (int a, int exponentResult) if (negate) { - // Let 1 have exponent=0 at bit 60 (2*MSB) - // raw result of inversion = 0-exponentWork at bit 30 + // Let 1 have exponent=-2*MSB + // raw result of inversion = -2*MSB - exponentWork + // goal = exponentResult + // shift = raw - goal = -2*MSB - exponentWork - exponentResult uint64_t temp = ((uint64_t) 1 << 2 * FP_MSB) / result; - int shift = -exponentWork - exponentResult; + int shift = -2 * FP_MSB - exponentWork - exponentResult; if (shift < 0) { - if (shift < -60) return 0; // Prevent weird effects from modulo arithmetic on size of shift. + if (shift < -2 * FP_MSB) return 0; // Prevent weird effects from modulo arithmetic on size of shift. return temp >> -shift; } if (shift > 0) @@ -877,8 +908,9 @@ exp (int a, int exponentResult) if (shift > 0) { if (shift > FP_MSB) return INFINITY; - // Don't bother trapping overflow with 32-bit math. Our fixed-point analysis should keep numbers in range in any case. - return result <<= shift; + uint64_t temp = (uint64_t) result << shift; + if (temp > INFINITY) return INFINITY; + return temp; } return result; } @@ -888,10 +920,9 @@ int floor (int a, int exponentA, int exponentResult) { int result; - if (exponentA >= 0 && exponentA < FP_MSB) + if (exponentA >= -FP_MSB && exponentA < 0) { - int decimalPlaces = FP_MSB - exponentA; - int wholeMask = 0xFFFFFFFF << decimalPlaces; + int wholeMask = 0xFFFFFFFF << -exponentA; result = a & wholeMask; } else @@ -909,11 +940,13 @@ int log (int a, int exponentA, int exponentResult) { // We use the simple identity log_e(a) = log_2(a) / log_2(e) - // exponentRaw = exponentResult - 0 + MSB + // exponentRaw = exponentResult - -MSB // shift = exponentRaw - exponentResult = MSB return ((int64_t) log2 (a, exponentA, exponentResult) << FP_MSB) / M_LOG2E; } +// This implementation tries to keep everything within a single 32-bit word, +// so it throws away some precision. int log2 (int a, int exponentA, int exponentResult) { @@ -922,7 +955,7 @@ log2 (int a, int exponentA, int exponentResult) // If a<1, then the result is -log2(1/a) bool negate = false; - if (exponentA < 0 || (exponentA < FP_MSB && a < 1 << FP_MSB - exponentA)) + if (exponentA < -FP_MSB || (exponentA < 0 && a < 1 << -exponentA)) { negate = true; @@ -936,45 +969,48 @@ log2 (int a, int exponentA, int exponentResult) } // compute 1/a - // Let the numerator 1 have center power of 0, with center at MSB/2, for exponent=MSB/2 - // Center of a is presumably at MSB/2 - // Center of inverse should be at MSB/2 - // Center power of a is exponentA-MSB/2 - // Center power of inverse is 0-(exponentA-MSB/2) = MSB/2-exponentA - // Exponent of inverse = center power + MSB/2 = MSB-exponentA - // Exponent of unshifted division = exponentRaw = MSB/2-exponentA+MSB = 3*MSB/2-exponentA - // Needed shift for exponent of inverse = exponentRaw - (MSB-exponentA) = MSB/2 - // If shifting 1 up to necessary position: (1 << MSB/2) << MSB/2 = 1 << MSB + // Let the numerator 1 have exponent=-MSB, to maximize significant bits in result. + // raw division exponent = -MSB - exponentA a = (1 << FP_MSB) / a; - exponentA = FP_MSB - exponentA; + exponentA = -FP_MSB - exponentA; } // At this point a >= 1 // Using the identity log(ab)=log(a)+log(b), we put a into normal form: // operand = a*2^exponentA // log2(operand) = log2(a)+log2(2^exponentA) = log2(a)+exponentA - int exponentWork = 15; - int one = 1 << FP_MSB - exponentWork; + // Our remaining work will be only on the mantissa of a, not its exponent. + // Goal is to shift the mantissa be in [1,2). Alternately, we define "1" + // with an exponent. + int exponentOne = -FP_MSB2; + int one = 1 << -exponentOne; while (a < one) { one >>= 1; - exponentWork++; + exponentOne++; } - int result = exponentA - exponentWork; // Represented as a pure integer, with exponent=MSB - int two = 2 * one; + int result = exponentA - exponentOne; // Represented as a pure integer, with exponent=0. Later this will have to be shifted to exponentResult. + int two = 2 * one; // Same as upshift by 1. while (a >= two) // This could also be done with a bit mask that checks for any bits in the twos position or higher. { result++; a = (a >> 1) + (a & 1); // divide-by-2 with rounding } + // "result" is now the exponent of the bit in a at power 0. + // That gives us the whole part of log_2(operand). We now need to compute the fractional part. + // TODO: Guard against large shifts. - int shift = FP_MSB - exponentResult; + int shift = -exponentResult; if (a > one) // Otherwise a==one, in which case the following algorithm will do nothing. { - while (shift > 0) + while (shift > 0) // The requested result has some fractional bits, so fill them. { - a = multiplyRound (a, a, exponentWork - FP_MSB); // exponentRaw - exponentWork = (2*exponentWork-MSB) - exponentWork + // exponentOne is the working exponent of mantissa a + // raw exponent of a^2 = 2*exponentOne + // goal = exponentOne + // shift = raw - goal = exponentOne + a = multiplyRound (a, a, exponentOne); result <<= 1; shift--; if (a >= two) @@ -983,7 +1019,7 @@ log2 (int a, int exponentA, int exponentResult) a = (a >> 1) + (a & 1); } } - a = multiplyRound (a, a, exponentWork - FP_MSB); + a = multiplyRound (a, a, exponentOne); if (a >= two) result++; } @@ -1056,15 +1092,15 @@ modFloor (int a, int b, int exponentA, int exponentB) int pow (int a, int b, int exponentA, int exponentResult) { - // exponentB = 15 + // exponentB = -MSB/2 // Use the identity: a^b = e^(b*ln(a)) // Most of the complexity of this function is in trapping special cases. // For details, see man page for floating-point pow(). // We don't have signed zero, so ignore all distinctions based on that. bool negate = false; - int blna = 1; // exponent=7, as required by exp(); Nonzero indicates that blna needs to be calculated. - int shift = FP_MSB - exponentA; + int blna = 1; // exponent=7-MSB, as required by exp(); Nonzero indicates that blna needs to be calculated. + int shift = -exponentA; int one; if (shift < 0) one = 0; else one = 1 << shift; @@ -1120,9 +1156,10 @@ pow (int a, int b, int exponentA, int exponentResult) if (blna) { - // raw multiply = exponentB+7-MSB at bit 30 - // shift = (exponentB+7-MSB)-7 = exponentB-MSB = -15 - int64_t temp = (int64_t) b * log (a, exponentA, 7) >> 15; + // raw multiply = exponentB + exponentBLNA + // goal = exponentBLNA + // shift = raw - goal = exponentB = -MSB/2 + int64_t temp = (int64_t) b * log (a, exponentA, 7-FP_MSB) >> FP_MSB2; if (temp > INFINITY) return INFINITY; if (temp < -INFINITY) return 0; blna = temp; @@ -1137,9 +1174,9 @@ int round (int a, int exponentA, int exponentResult) { int result; - if (exponentA >= 0 && exponentA < FP_MSB) + if (exponentA >= -FP_MSB && exponentA < 0) { - int decimalPlaces = FP_MSB - exponentA; + int decimalPlaces = -exponentA; int mask = 0xFFFFFFFF << decimalPlaces; int half = 0x1 << decimalPlaces - 1; result = a + half & mask; @@ -1159,7 +1196,7 @@ int sgn (int a, int exponentResult) { if (a == 0) return 0; - int result = 0x1 << FP_MSB - exponentResult; // This breaks for exponentResult outside [0, MSB], but the calling code is already meaningless in that case. + int result = 0x1 << -exponentResult; // This breaks for exponentResult outside [-MSB, 0], but the calling code is already meaningless in that case. if (a < 0) return -result; return result; } @@ -1170,8 +1207,8 @@ sqrt (int a, int exponentA, int exponentResult) if (a < 0) return NAN; // Simple approach: apply the identity a^0.5=e^(ln(a^0.5))=e^(0.5*ln(a)) - //int l = log (a, exponentA, MSB/2) >> 1; - //return exp (l, MSB/2, exponentResult); + //int l = log (a, exponentA, -MSB/2) >> 1; + //return exp (l, -MSB/2, exponentResult); // More efficient approach: Use digit-by-digit method described in // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots @@ -1186,13 +1223,12 @@ sqrt (int a, int exponentA, int exponentResult) // sqrt(2m) >> -(n-1)/2 uint32_t m = a; // "m" for mantissa - int exponent0 = exponentA - FP_MSB; // exponent at bit position 0 - if (exponent0 % 2) // Odd, so leave remainder inside sqrt() + if (exponentA % 2) // Odd, so leave remainder inside sqrt() { m <<= 1; // equivalent to "2m" in the comments above - exponent0--; + exponentA--; } - int exponentRaw = exponent0 / 2 + FP_MSB; // exponent of raw result at MSB + int exponentRaw = exponentA / 2; // exponent of raw result uint32_t bit; if (m & 0xFFFE0000) bit = 1 << 30; @@ -1243,13 +1279,12 @@ sqrt (int64_t a, int exponentA, int exponentResult) if (a < 0) return NAN; uint64_t m = a; // "m" for mantissa - int exponent0 = exponentA - FP_MSB; // exponent at bit position 0 - if (exponent0 % 2) // Odd, so leave remainder inside sqrt() + if (exponentA % 2) // Odd, so leave remainder inside sqrt() { m <<= 1; // equivalent to "2m" in the comments above - exponent0--; + exponentA--; } - int exponentRaw = exponent0 / 2 + FP_MSB; // exponent of raw result at MSB + int exponentRaw = exponentA / 2; uint64_t bit; if (m & 0xFFFFFFFF80000000) bit = (int64_t) 1 << 60; @@ -1292,12 +1327,14 @@ sqrt (int64_t a, int exponentA, int exponentResult) int sin (int a, int exponentA) { + const int exponentResult = 1 - FP_MSB; + // Limit a to [0,pi/2) // To create 2PI, we lie about the exponent of M_PI, increasing it by 1. - a = modFloor (a, M_PI, exponentA, 2); // exponent = min (exponentA, 2) - int shift = exponentA - 2; - if (shift < 0) a >>= -shift; - const int PIat2 = M_PI >> 1; // M_PI with exponent=2 rather than 1 + a = modFloor (a, M_PI, exponentA, 2-FP_MSB); // exponent = min (exponentA, 2-MSB) + int shift = exponentA + FP_MSB - 2; // shift = exponentA - (2-MSB); If negative, then result above has exponentA. + if (shift < 0) a >>= -shift; // Shift it to match exponent of 2*pi. + const int PIat2 = M_PI >> 1; // M_PI with exponent=2-MSB rather than 1-MSB bool negate = false; if (a > PIat2) { @@ -1305,27 +1342,27 @@ sin (int a, int exponentA) negate = true; } if (a > PIat2 >> 1) a = PIat2 - a; - a <<= 1; // Now exponent=1, which matches our promised exponentResult. + a <<= 1; // Now exponent=1-MSB, which matches our promised exponentResult. // Use power-series to compute sine, similar to exp() // sine(a) = sum_0^inf (-1)^n * x^(2n-1) / (2n+1)! = x - x^3/3! + x^5/5! - x^7/7! ... int term = a; int result = a; // zeroth term - int n1 = 0; // exponent=MSB + int n1 = 0; // exponent=0 int n2 = 1; for (int n = 1; n < 7; n++) { n1 = n2 + 1; n2 = n1 + 1; // raw exponent of operations below, in evaluation order: - // 2*exponentResult at bit 60 - // 2*exponentResult-MSB at bit 30 - // shift = (2*exponenetResult-MSB)-exponentResult = exponentResult-MSB = -29 - // 2*exponentResult at bit 60 - // 2*exponentResult-MSB at bit 30 + // multiply: 2*exponentResult + // divide: 2*exponentResult + // shift = 2*exponenetResult - exponentResult = exponentResult = 1-MSB + // multiply: 2*exponentResult + // divide : 2*exponentResult // same shift again - term = ((int64_t) -term * a / n1 >> 29) * a / n2 >> 29; + term = ((int64_t) -term * a / n1 >> -exponentResult) * a / n2 >> -exponentResult; if (term == 0) break; result += term; } @@ -1340,32 +1377,37 @@ tan (int a, int exponentA, int exponentResult) // See http://mathworld.wolfram.com/MaclaurinSeries.html // However, to save space we simply compute sin()/cos(). // raw division exponent=0 at bit 0 - return ((int64_t) sin (a, exponentA) << exponentResult) / cos (a, exponentA); // Don't do any saturation checks. We are not really interested in infinity. + return ((int64_t) sin (a, exponentA) << -exponentResult) / cos (a, exponentA); // Don't do any saturation checks. We are not really interested in infinity. } int tanh (int a, int exponentA) { // result = (exp(2a) - 1) / (exp(2a) + 1) - // exponentResult = 0 + // exponentResult = -MSB // tanh() is symmetric around 0, so only deal with one sign bool negate = a < 0; if (negate) a = -a; if (a == 0) return 0; // This also traps NAN - // Determine exponent desired from exp(2a). The result from exp(2a) is never smaller than 1, so exponent>=0 + // Determine "exponentOne" desired from exp(2a). + // The result from exp(2a) is never smaller than 1, so exponentOne>=-MSB + // 1 must be representable, so exponentOne<=0. // We want enough bits to contain the msb of the output. This is - // exponent = log2(exp(2a)) = ln(exp(2a)) / ln(2) = 2a / (log2(2) / log2(e)) = 2a * log2(e) - // "exponent" should be expressed as a simple integer - // raw = exponentA + 1 + exponentLOG2E - MSB = exponentA + 1 - MSB - // shift = raw - MSB = exponentA + 1 - 2 * MSB - int exponent = 0; - if (exponentA >= -1) + // exponentMSB = log2(exp(2a)) = ln(exp(2a)) / ln(2) = 2a / (log2(2) / log2(e)) = 2a * log2(e) + // Claiming that exponentA is one higher has the effect of multiplying by 2. + // exponents of the multiplication itself are: + // raw = (exponentA+1) + exponentLOG2E = exponentA + 1 + -MSB + // goal = 0 (integer count of bits) + // shift = raw - goal = exponentA + 1 - MSB + int exponentOne = -FP_MSB; + if (exponentA >= -1 - FP_MSB) // Otherwise 2a is less than 1, so no need for the following calculation. { - exponent = multiplyCeil (a, M_LOG2E, exponentA + 1 - 2 * FP_MSB); - // If exponent gets too large, the result will always be 1 or -1 - if (exponent > FP_MSB) return negate ? -0x40000000 : 0x40000000; + exponentOne = multiplyCeil (a, M_LOG2E, exponentA + 1 - FP_MSB); // exponentOne now refers to MSB + // If exponentOne gets too large, the result will always be 1 or -1 + if (exponentOne > FP_MSB || exponentOne == 0) return negate ? -0x40000000 : 0x40000000; + exponentOne -= FP_MSB; // exponentOne now refers to LSB } // Find true magnitude of a @@ -1377,21 +1419,31 @@ tanh (int a, int exponentA) // Require at least 16 bits for exp() after downshifting. // Otherwise answer is not as accurate as simple linear. - if (exponentA < 22 - FP_MSB) // Includes offset of 6 for the downshift. + // See the downshift after this if-block. This is equivalent to ensuring that + // total amount of shift (6-MSB-exponentA) < 16. + // Solving for exponentA: exponentA > -10 - MSB + // Negating logic ... + if (exponentA <= -10 - FP_MSB) { // Return linear answer, if possible. - if (exponentA < -FP_MSB) return 0; // Can't return result with correct magnitude, so only option is zero. - int result = a >> -exponentA; + // goal = exponentResult = -MSB + // shift = exponentA - goal = exponentA + MSB + if (exponentA < -2 * FP_MSB) return 0; // Can't return result with correct magnitude, so only option is zero. + int result = a >> -exponentA - FP_MSB; if (negate) return -result; return result; } // Set correct magnitude for exp(). - // exp(a) expects exponentA=7, but we want exp(2a), so shift to exponentA=6 and lie about it - a >>= 6 - exponentA; + // exp(a) expects exponentA=7-MSB, but we want exp(2a), so shift to exponentA=6 and lie about it + // goal = 6-MSB + // shift = exponentA - goal = exponentA - 6 + MSB + // This implies that exponentA <= 6-MSB. + if (exponentA > 6 - FP_MSB) return negate ? -0x40000000 : 0x40000000; + a >>= 6 - FP_MSB - exponentA; // call exp() and complete calculation - int result = exp (a, exponent); - int one = 1 << FP_MSB - exponent; + int result = exp (a, exponentOne); + int one = 1 << -exponentOne; result = ((int64_t) result - one << FP_MSB) / ((int64_t) result + one); if (negate) return -result; diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.cc b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.cc index 8f7b369e..86afdefa 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.cc +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.cc @@ -120,6 +120,7 @@ template class MatrixStrided; template class Matrix; template class MatrixFixed; +template SHARED void clear ( MatrixStrided & A, const float scalar); template SHARED float norm (const MatrixStrided & A, float n); template SHARED Matrix normalize (const MatrixAbstract & A); template SHARED Matrix operator * (const MatrixStrided & A, const MatrixAbstract & B); @@ -150,16 +151,16 @@ template class InputHolder; template class OutputHolder; template SHARED IteratorNonzero * getIterator (MatrixAbstract * A); #ifdef n2a_FP -template SHARED MatrixInput * matrixHelper (const String & fileName, int exponent, MatrixInput * oldHandle); -template SHARED InputHolder * inputHelper (const String & fileName, int exponent, InputHolder * oldHandle); +template SHARED MatrixInput * matrixHelper (const String & fileName, int exponent, MatrixInput * oldHandle); +template SHARED InputHolder * inputHelper (const String & fileName, int exponent, int exponentRow, InputHolder * oldHandle); #else -template SHARED MatrixInput * matrixHelper (const String & fileName, MatrixInput * oldHandle); -template SHARED InputHolder * inputHelper (const String & fileName, InputHolder * oldHandle); +template SHARED MatrixInput * matrixHelper (const String & fileName, MatrixInput * oldHandle); +template SHARED InputHolder * inputHelper (const String & fileName, InputHolder * oldHandle); #endif -template SHARED Mfile * MfileHelper (const String & fileName, Mfile * oldHandle); -template SHARED OutputHolder * outputHelper (const String & fileName, OutputHolder * oldHandle); -template SHARED ImageInput * imageInputHelper (const String & fileName, ImageInput * oldHandle); -template SHARED ImageOutput * imageOutputHelper(const String & fileName, ImageOutput * oldHandle); +template SHARED Mfile * MfileHelper (const String & fileName, Mfile * oldHandle); +template SHARED OutputHolder * outputHelper (const String & fileName, OutputHolder * oldHandle); +template SHARED ImageInput * imageInputHelper (const String & fileName, ImageInput * oldHandle); +template SHARED ImageOutput * imageOutputHelper(const String & fileName, ImageOutput * oldHandle); #ifdef HAVE_JNI @@ -398,6 +399,14 @@ putUnique (std::vector & vertices, float x, float y, float z) void icosphere (std::vector & vertices, std::vector & indices) { + // This function is always defined as float. + // When compiling the fixed-point runtime, it is necessary to restore the + // floating-point definition of pi. +# ifdef n2a_FP +# undef M_PI +# define M_PI 3.14159265359f +# endif + float angleH = 2 * M_PI / 5; // 72 degrees float angleV = atan (0.5); // elevation = 26.565 degree diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.h b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.h index c2b784b2..76c7cbd1 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.h +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.h @@ -379,7 +379,8 @@ class SHARED InputHolder : public Holder bool delimiterSet; T epsilon; ///< for time values # ifdef n2a_FP - int exponent; ///< of value returned by get() + int exponent; ///< of value returned by get() + int exponentRow; ///< of row value passed into get() # endif InputHolder (const String & fileName); @@ -391,9 +392,9 @@ class SHARED InputHolder : public Holder Matrix get (T row); }; #ifdef n2a_FP -template SHARED InputHolder * inputHelper (const String & fileName, int exponent, InputHolder * oldHandle = 0); +template SHARED InputHolder * inputHelper (const String & fileName, int exponent, int exponentRow, InputHolder * oldHandle = 0); #else -template SHARED InputHolder * inputHelper (const String & fileName, InputHolder * oldHandle = 0); +template SHARED InputHolder * inputHelper (const String & fileName, InputHolder * oldHandle = 0); #endif template diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.tcc b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.tcc index 162732bb..1f8c8c22 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.tcc +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/holder.tcc @@ -21,6 +21,7 @@ the U.S. Government retains certain rights in this software. #include #ifdef _MSC_VER # define stat _stat +# define timegm _mkgmtime #else # include #endif @@ -274,7 +275,7 @@ convert (String input, int exponent) bits |= 0x10000000000000l; // set implied msb of mantissa (bit 52) to 1 bits &= 0x1FFFFFFFFFFFFFl; // clear sign and exponent bits if (negate) bits = -bits; - return bits >> 52 - FP_MSB + exponent - e; + return bits >> 52 + exponent - e; } #endif @@ -594,7 +595,7 @@ ImageInput::get (String channelName, T now, bool step) double nextPTS = atof (video->get ("nextPTS").c_str ()); if (framePeriod) nextPTS /= framePeriod; # ifdef n2a_FP - t = (T) (nextPTS * pow (2.0, FP_MSB - Event::exponent)); + t = (T) (nextPTS * pow (2.0, -Event::exponent)); # else t = nextPTS; # endif @@ -605,7 +606,7 @@ ImageInput::get (String channelName, T now, bool step) if (pattern.size ()) { # ifdef n2a_FP - index = (int) (now / pow (2.0f, FP_MSB - Event::exponent)); + index = (int) (now / pow (2.0f, -Event::exponent)); # else index = (int) now; # endif @@ -662,9 +663,9 @@ ImageInput::get (String channelName, T now, bool step) # ifdef n2a_FP // Convert buffer to int. - float conversion = pow (2.0f, FP_MSB - exponent); + float conversion = pow (2.0f, -exponent); int count = image.width * image.height * 3; - Pointer q (count * sizeof (T)); + n2a::Pointer q (count * sizeof (T)); float * from = (float *) p; T * to = (T *) q; T * end = to + count; @@ -870,7 +871,7 @@ template<> void setColor (float target[], const Matrix & color, bool withAlpha) { - float conversion = powf (2, FP_MSB); // exponent=0 + float conversion = powf (2, FP_MSB); // exponent=MSB target[0] = color[0] / conversion; target[1] = color[1] / conversion; target[2] = color[2] / conversion; @@ -896,11 +897,11 @@ ImageOutput::setClearColor (const Matrix & color) if (count == 1) count = color.columns (); bool hasAlpha = count > 3; - uint32_t r = std::min (1.0f, std::max (0.0f, color[0])) * 255; - uint32_t g = std::min (1.0f, std::max (0.0f, color[1])) * 255; - uint32_t b = std::min (1.0f, std::max (0.0f, color[2])) * 255; + uint32_t r = std::min (1.0f, std::max (0.0f, (float) color[0])) * 255; + uint32_t g = std::min (1.0f, std::max (0.0f, (float) color[1])) * 255; + uint32_t b = std::min (1.0f, std::max (0.0f, (float) color[2])) * 255; uint32_t a = 0xFF; - if (hasAlpha) a = std::min (1.0f, std::max (0.0f, color[3])) * 255; + if (hasAlpha) a = std::min (1.0f, std::max (0.0f, (float) color[3])) * 255; clearColor = r << 24 | g << 16 | b << 8 | a; # ifdef HAVE_GL @@ -971,7 +972,7 @@ ImageOutput::drawDisc (T now, bool raw, const MatrixFixed & center, next (now); # ifdef n2a_FP - double conversion = pow (2.0, FP_MSB - exponent); + double conversion = pow (2.0, -exponent); MatrixFixed center = centerFP; center /= conversion; double radius = radiusFP / conversion; @@ -1007,7 +1008,7 @@ ImageOutput::drawSquare (T now, bool raw, const MatrixFixed & center, next (now); # ifdef n2a_FP - double conversion = pow (2.0, FP_MSB - exponent); + double conversion = pow (2.0, -exponent); n2a::Point center = centerFP; center /= conversion; double w = wFP / conversion; @@ -1045,7 +1046,7 @@ ImageOutput::drawSegment (T now, bool raw, const MatrixFixed & p1, c next (now); # ifdef n2a_FP - double conversion = pow (2.0, FP_MSB - exponent); + double conversion = pow (2.0, -exponent); n2a::Point p1 = p1FP; n2a::Point p2 = p2FP; p1 /= conversion; @@ -1146,7 +1147,7 @@ ImageOutput::writeImage () if (timeScale) { # ifdef n2a_FP - canvas.timestamp = timeScale * t / pow (2.0, FP_MSB - Event::exponent); + canvas.timestamp = timeScale * t / pow (2.0, -Event::exponent); # else canvas.timestamp = timeScale * t; # endif @@ -1424,7 +1425,7 @@ T #ifdef n2a_FP ImageOutput::drawCube (T now, const Material & material, const Matrix & model, int exponentP) { - const_cast &> (model) /= powf (2, FP_MSB - exponentP); + const_cast &> (model) /= powf (2, -exponentP); #else ImageOutput::drawCube (T now, const Material & material, const Matrix & model) { @@ -1515,8 +1516,8 @@ T #ifdef n2a_FP ImageOutput::drawCylinder (T now, const Material & material, const MatrixFixed & p1, int exponentP, float r1, int exponentR, const MatrixFixed & p2, float r2, int cap1, int cap2, int steps, int stepsCap) { - float scaleP = powf (2, FP_MSB - exponentP); - float scaleR = powf (2, FP_MSB - exponentR); + float scaleP = powf (2, -exponentP); + float scaleR = powf (2, -exponentR); const_cast &> (p1) /= scaleP; const_cast &> (p2) /= scaleP; r1 /= scaleR; @@ -1573,9 +1574,17 @@ ImageOutput::drawCylinder (T now, const Material & material, const MatrixFixe // This also works when r1 is shorter than r2. float slopeZ = (r2 - r1) / length; +# ifdef n2a_FP +# pragma push_macro("M_PI") +# undef M_PI +# define M_PI 3.14159265359f +# endif float angleStep = M_PI * 2 / steps; float angleStep2 = M_PI / steps; // half of a regular step float angleStepCap = M_PI / 2 / (stepsCap + 1); +# ifdef n2a_FP +# pragma pop_macro("M_PI") +# endif int rowsBegin = 0; // Index of first vertex in first ring. int rowsEnd = 0; // Index of first vertex after last ring. @@ -1789,7 +1798,7 @@ T #ifdef n2a_FP ImageOutput::drawPlane (T now, const Material & material, const Matrix & model, int exponentP) { - const_cast &> (model) /= powf (2, FP_MSB - exponentP); + const_cast &> (model) /= powf (2, -exponentP); #else ImageOutput::drawPlane (T now, const Material & material, const Matrix & model) { @@ -1838,7 +1847,7 @@ T #ifdef n2a_FP ImageOutput::drawSphere (T now, const Material & material, const Matrix & model, int exponentP, int steps) { - const_cast &> (model) /= powf (2, FP_MSB - exponentP); + const_cast &> (model) /= powf (2, -exponentP); #else ImageOutput::drawSphere (T now, const Material & material, const Matrix & model, int steps) { @@ -2235,7 +2244,7 @@ InputHolder::getRow (T row) bool valid = false; if (index == timeColumn) { - int year = 1970; // will be adjusted below for mktime() + int year = 1970; // will be adjusted below for timegm() int month = 1; // ditto int day = 1; int hour = 0; @@ -2296,12 +2305,12 @@ InputHolder::getRow (T row) struct tm date; date.tm_isdst = 0; // time is strictly UTC, with no DST - // ignoring tm_wday and tm_yday, as mktime() doesn't do anything with them + // ignoring tm_wday and tm_yday, as timegm() doesn't do anything with them - // Hack to adjust for mktime() that can't handle dates before posix epoch (1970/1/1). + // Hack to adjust for timegm() that can't handle dates before posix epoch (1970/1/1). // This simple hack only works for years after ~1900. // Solution comes from https://bugs.php.net/bug.php?id=17123 - // Alternate solution would be to implement a simple mktime() right here. + // Alternate solution would be to implement a simple timegm() right here. // Since we don't care about DST or timezones, all it has to do is handle Gregorion leap-years. time_t offset = 0; if (year <= 70) // Yes, that includes 1970 itself. @@ -2314,7 +2323,7 @@ InputHolder::getRow (T row) date.tm_hour = 0; date.tm_min = 0; date.tm_sec = 0; - offset = mktime (&date); + offset = timegm (&date); } date.tm_year = year; @@ -2324,10 +2333,10 @@ InputHolder::getRow (T row) date.tm_min = minute; date.tm_sec = second; - nextValues[index] = mktime (&date) - offset; // Unix time; an integer, so exponent=MSB + nextValues[index] = timegm (&date) - offset; // Unix time; an integer, so exponent=0 # ifdef n2a_FP // Need to put value in expected exponent. - int shift = FP_MSB - (time ? Event::exponent : exponent); + int shift = -(time ? exponentRow : exponent); if (shift >= 0) nextValues[index] <<= shift; else nextValues[index] >>= -shift; # endif @@ -2337,7 +2346,7 @@ InputHolder::getRow (T row) if (! valid) // Not a date, so general case ... { # ifdef n2a_FP - nextValues[index] = convert (field, time && index == timeColumn ? Event::exponent : exponent); + nextValues[index] = convert (field, time && index == timeColumn ? exponentRow : exponent); # else nextValues[index] = (T) atof (field.c_str ()); # endif @@ -2516,9 +2525,9 @@ InputHolder::get (T row) template InputHolder * #ifdef n2a_FP -inputHelper (const String & fileName, int exponent, InputHolder * oldHandle) +inputHelper (const String & fileName, int exponent, int exponentRow, InputHolder * oldHandle) #else -inputHelper (const String & fileName, InputHolder * oldHandle) +inputHelper (const String & fileName, InputHolder * oldHandle) #endif { InputHolder * handle = (InputHolder *) SIMULATOR getHolder (fileName, oldHandle); @@ -2527,7 +2536,8 @@ inputHelper (const String & fileName, InputHolder * oldHandle) handle = new InputHolder (fileName); SIMULATOR holders.push_back (handle); # ifdef n2a_FP - handle->exponent = exponent; + handle->exponent = exponent; + handle->exponentRow = exponentRow; # endif } return handle; @@ -2602,7 +2612,7 @@ OutputHolder::trace (T now) { columnMap["$t"] = 0; # ifdef n2a_FP - columnValues.push_back ((float) t / pow (2.0f, FP_MSB - Event::exponent)); + columnValues.push_back ((float) t / pow (2.0f, -Event::exponent)); # else columnValues.push_back (t); # endif @@ -2611,7 +2621,7 @@ OutputHolder::trace (T now) else { # ifdef n2a_FP - columnValues[0] = (float) t / pow (2.0f, FP_MSB - Event::exponent); + columnValues[0] = (float) t / pow (2.0f, -Event::exponent); # else columnValues[0] = t; # endif @@ -2670,7 +2680,7 @@ OutputHolder::trace (T now, const String & column, T value, c if (valueFP == NAN) value = std::numeric_limits::quiet_NaN (); else if (valueFP == INFINITY) value = std::numeric_limits::infinity (); else if (valueFP == -INFINITY) value = -std::numeric_limits::infinity (); - else value = (float) valueFP / pow (2.0f, FP_MSB - exponent); + else value = (float) valueFP / pow (2.0f, -exponent); # endif std::unordered_map::iterator result = columnMap.find (column); @@ -2757,27 +2767,24 @@ OutputHolder::trace (T now, const String & column, const Matrix & A, const template T #ifdef n2a_FP -OutputHolder::trace (T now, T columnFP, T valueFP, int exponent, const char * mode) +OutputHolder::trace (T now, T index, T valueFP, int exponent, const char * mode) #else -OutputHolder::trace (T now, T column, T value, const char * mode) +OutputHolder::trace (T now, T column, T value, const char * mode) #endif { trace (now); # ifdef n2a_FP - float column = (float) columnFP / pow (2.0f, FP_MSB - 15); // column has fixed exponent of 15 float value; if (valueFP == NAN) value = std::numeric_limits::quiet_NaN (); else if (valueFP == INFINITY) value = std::numeric_limits::infinity (); else if (valueFP == -INFINITY) value = -std::numeric_limits::infinity (); - else value = (float) valueFP / pow (2.0f, FP_MSB - exponent); + else value = (float) valueFP / pow (2.0f, -exponent); +# else + int index = (int) column; // truncates floating-point # endif - String columnName; - int index; // Only used for "raw" mode. - if (raw) columnName = index = (int) round (column); - else columnName = column; - + String columnName = index; std::unordered_map::iterator result = columnMap.find (columnName); if (result == columnMap.end ()) { diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/matrix.h b/N2A/src/gov/sandia/n2a/backend/c/runtime/matrix.h index ef56c477..d1b0a582 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/matrix.h +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/matrix.h @@ -55,7 +55,7 @@ class SHARED MatrixAbstract virtual int columns () const = 0; }; -template class Matrix; +template class SHARED Matrix; template SHARED void clear (const MatrixAbstract & A, const T scalar = (T) 0); ///< Set all elements to given value. const on A is a lie to allow transient matrix regions. template SHARED void identity (const MatrixAbstract & A); ///< Set diagonal to 1 and all other elements to 0. Does not have to be a square matrix. const on A is a lie to allow transient matrix regions. @@ -180,8 +180,8 @@ template SHARED Matrix operator - (const T scalar, const // Fixed-point operations on MatrixStrided // These are not templates. Their implementations are in fixedpoint.cc SHARED void identity (const MatrixStrided & A, int one); // "one" is passed explicitly, already scaled with the right exponent -SHARED int norm (const MatrixStrided & A, int n, int exponentA, int exponentResult); // exponentN=15 -SHARED Matrix normalize (const MatrixStrided & A, int exponentA); // result has exponent=0 +SHARED int norm (const MatrixStrided & A, int n, int exponentA, int exponentResult); // exponentN=-MSB/2 +SHARED Matrix normalize (const MatrixStrided & A, int exponentA); // result has exponent=-MSB SHARED Matrix cross (const MatrixStrided & A, const MatrixStrided & B, int shift); SHARED Matrix visit (const MatrixStrided & A, int (*function) (int, int), int exponent1); diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/mymath.h b/N2A/src/gov/sandia/n2a/backend/c/runtime/mymath.h index 37ad55c5..37996a32 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/mymath.h +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/mymath.h @@ -80,9 +80,9 @@ namespace n2a #define FP_MSB 30 #define FP_MSB2 15 -#define M_LOG2E 1549082004 // log_2(e) = 1.4426950408889634074; exponent=0 -#define M_E 1459366444 // exponent=1 -#define M_PI 1686629713 // exponent=1 +#define M_LOG2E 1549082004 // log_2(e) = 1.4426950408889634074; exponent=-MSB +#define M_E 1459366444 // exponent=1-MSB +#define M_PI 1686629713 // exponent=1-MSB #define NAN 0x80000000 #define INFINITY 0x7FFFFFFF @@ -110,22 +110,22 @@ namespace std } } -SHARED int atan2 (int y, int x); // returns angle in [-pi,pi], exponentResult=1; exponentA must match exponentB, but it may be arbitrary. +SHARED int atan2 (int y, int x); // returns angle in [-pi,pi], exponentResult=1-MSB; exponentA must match exponentB, but it may be arbitrary. SHARED int ceil (int a, int exponentA, int exponentResult); // Only needed for matrices. For scalars, an immediate implementation is emitted. -SHARED int cos (int a, int exponentA); // exponentResult=1 -SHARED int exp (int a, int exponentResult); // exponentA=7 +SHARED int cos (int a, int exponentA); // exponentResult=1-MSB +SHARED int exp (int a, int exponentResult); // exponentA=7-MSB SHARED int floor (int a, int exponentA, int exponentResult); // Only needed for matrices. SHARED int log (int a, int exponentA, int exponentResult); -SHARED int log2 (int a, int exponentA, int exponentResult); // Use as a subroutine. Not directly available to user. +SHARED int log2 (int a, int exponentA, int exponentResult); // Used as a subroutine. Not directly available to user. SHARED int modFloor (int a, int b, int exponentA, int exponentB); // exponentResult is promised to be min(exponentA,exponentB) -SHARED int pow (int a, int b, int exponentA, int exponentResult); // exponentB=15 +SHARED int pow (int a, int b, int exponentA, int exponentResult); // exponentB=-MSB/2 SHARED int round (int a, int exponentA, int exponentResult); // Only needed for matrices. SHARED int sgn (int a, int exponentResult); // Only needed for matrices. -SHARED int sin (int a, int exponentA); // exponentResult=1 +SHARED int sin (int a, int exponentA); // exponentResult=1-MSB SHARED int sqrt (int a, int exponentA, int exponentResult); -SHARED int sqrt (int64_t a, int exponentA, int exponentResult); // 64-bit version of sqrt(). exponentA still refers to power in bit position FP_MSB, not any bit in the upper half of a. +SHARED int sqrt (int64_t a, int exponentA, int exponentResult); // 64-bit version of sqrt(). SHARED int tan (int a, int exponentA, int exponentResult); -SHARED int tanh (int a, int exponentA); // exponentResult=0 +SHARED int tanh (int a, int exponentA); // exponentResult=-MSB #endif // n2a_FP diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.cc b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.cc index ed408e31..d7fe0491 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.cc +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.cc @@ -28,10 +28,18 @@ template SHARED n2a_T pulse (n2a_T t, n2a_T width, n2a_T period, n2a_T rise, n2a template SHARED n2a_T unitmap (const MatrixAbstract & A, n2a_T row, n2a_T column); -#ifdef n2a_FP +#define float 1 // This is an evil hack to get around lack of string comparison in c preprocessor. +#define double 2 +#define int 3 +#if n2a_T != float +#undef float template SHARED Matrix glOrtho (float left, float right, float bottom, float top, float near, float far); -// Non-templated integer gl functions are defined in fixedpoint.cc -#else +#endif +#undef float // Undo the hack +#undef double +#undef int + +#ifndef n2a_FP template SHARED Matrix glFrustum (n2a_T left, n2a_T right, n2a_T bottom, n2a_T top, n2a_T near, n2a_T far); template SHARED Matrix glOrtho (n2a_T left, n2a_T right, n2a_T bottom, n2a_T top, n2a_T near, n2a_T far); template SHARED Matrix glLookAt (const MatrixFixed & eye, const MatrixFixed & center, const MatrixFixed & up); @@ -43,6 +51,7 @@ template SHARED Matrix glScale (n2a_T sx, n2a_T sy, n2a_T sz); template SHARED Matrix glTranslate (const MatrixFixed & position); template SHARED Matrix glTranslate (n2a_T x, n2a_T y, n2a_T z); #endif +// else n2a_FP is defined -- Non-templated integer gl functions are defined in fixedpoint.cc template SHARED void removeMonitor (std::vector *> & partList, Part * part); diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.h b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.h index 188b80fb..5dbcf5d1 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.h +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.h @@ -68,7 +68,7 @@ template MatrixFixed uniform (const MatrixFixed (); # ifdef n2a_FP - return multiply (sigma, temp, 31); // See uniform(int) in runtime.tcc + return multiply (sigma, temp, 1 + FP_MSB); // See uniform(int) in runtime.tcc # else return sigma * temp; # endif @@ -92,7 +92,7 @@ template MatrixFixed gaussian (const MatrixFixed (); # ifdef n2a_FP - return multiply (sigma, temp, 28); // See gaussian(int) in runtime.tcc + return multiply (sigma, temp, FP_MSB - 2); // See gaussian(int) in runtime.tcc # else return sigma * temp; # endif @@ -109,11 +109,11 @@ template MatrixFixed sphere (const MatrixFixed & s // directions on a sphere of dimension R, then scale that vector so that it fills // the volume evenly between radius 0 and 1. This requires less density near the center. # ifdef n2a_FP - int n = norm (result, 2 << 15, 2, 4); // exponentA is result of gaussian(). exponentResult should be roughly 2+log2(R) (exponent of gaussian, plus log of the number of additions). Currently hardcoded for up to 4 rows. - int R1 = (1 << 30) / R >> 15; // result of division is exponent 0. R1 must be exponent 15. - int scale = pow (uniform (), R1, -1, -1); - result = divide (result, n, FP_MSB-2); // result of division is 2-4+MSB = MSB-2. Result is always in [0,1], so exponentResult should be 0. - result = multiply (result, scale, FP_MSB+1); // result of multiply is 0+(-1)-MSB=-MSB-1. Result is always in [0,1]. + int n = norm (result, 2 << FP_MSB2, 2-FP_MSB, 4-FP_MSB); // exponentA is result of gaussian(). exponentResult should be roughly 2+log2(R)-MSB (exponent of gaussian, plus log of the number of additions). Currently hardcoded for up to 4 rows. + int R1 = (1 << FP_MSB2) / R; // R1 exponent is -MSB/2, as required by pow() + int scale = pow (uniform (), R1, -1-FP_MSB, -1-FP_MSB); + result = divide (result, n, FP_MSB-2); // result of division is (2-MSB) - (4-MSB) = -2. Result is always in [0,1], so exponentResult should be -MSB. + result = multiply (result, scale, FP_MSB+1); // result of multiply is -MSB + -1-MSB = -2*MSB-1. Result is always in [0,1]. return multiplyElementwise (sigma, result, FP_MSB); // Result should have same exponent as sigma # else T scale = pow (uniform (), (T) 1 / R) / norm (result, (T) 2); @@ -131,9 +131,9 @@ template MatrixFixed sphere (const MatrixFixed (); # ifdef n2a_FP - int n = norm (temp, 2 << 15, 2, 4); - int C1 = (1 << 30) / C >> 15; - int scale = pow (uniform (), C1, -1, -1); + int n = norm (temp, 2 << FP_MSB2, 2-FP_MSB, 4-FP_MSB); + int C1 = (1 << FP_MSB2) / C; + int scale = pow (uniform (), C1, -1-FP_MSB, -1-FP_MSB); temp = divide (temp, n, FP_MSB-2); temp = multiply (temp, scale, FP_MSB+1); return multiply (sigma, temp, FP_MSB); @@ -289,7 +289,7 @@ class SHARED Part : public Simulatable // Accessors for $variables virtual T getLive (); ///< @return 1 if we are in normal simulation. 0 if we have died. Default is 1. - virtual T getP (); ///< Default is 1 (always create) + virtual T getP (); ///< Default is 1 (always create). exponent=-MSB virtual void getXYZ (MatrixFixed & xyz); ///< Default is [0;0;0]. // Events diff --git a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.tcc b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.tcc index 77e81bda..5f01a494 100644 --- a/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.tcc +++ b/N2A/src/gov/sandia/n2a/backend/c/runtime/runtime.tcc @@ -346,7 +346,7 @@ template<> int uniform () { - // exponent=-1; We promise the semi-open interval [0,1), so must never actaully reach 1. + // exponent=-1-MSB; We promise the semi-open interval [0,1), so must never actaully reach 1. #if RAND_MAX == 0x7FFFFFFF return rand (); #elif RAND_MAX == 0x7FFF @@ -360,7 +360,10 @@ template<> int uniform (int sigma) { - return (int64_t) uniform () * sigma >> 31; // shift = -1 - MSB + // raw = (-1-MSB) + exponentSigma + // goal = exponentSigma + // shift = raw - goal = -1-MSB + return (int64_t) uniform () * sigma >> 1 + FP_MSB; } template<> @@ -374,7 +377,7 @@ uniform (int lo, int hi, int step) // Box-Muller method (polar variant) for Gaussian random numbers. // Although this method can return very large values, we limit it to strictly -// less than 8 std (3 bits above the decimal point). Result exponent=2. +// less than 8 std (3 bits above the decimal point). Result exponent=2-MSB. template<> int gaussian () @@ -389,28 +392,29 @@ gaussian () } else { - const int half = 0x40000000; // 0.5, with exponent=-1 - const int big = 0xFFFF; // Too large for log(). exponent=14 - const int small = 0x8; // Too small for the division that creates multiplier. exponent=14 + const int half = 0x40000000; // 0.5, with exponent=-1-MSB + const int big = 0xFFFF; // Too large for log(). exponent=-16 + const int small = 0x8; // Too small for the division that creates multiplier. exponent=-16 int v1, v2, s; do { - v1 = uniform () - half; // u-0.5; Then implicitly double by treating exponent as 0 rather than -1. + v1 = uniform () - half; // u-0.5; Then implicitly double by treating exponent as -MSB rather than -1-MSB. v2 = uniform () - half; - // Squaring v puts exponent=0 at bit MSB*2 - // Down-shift puts exponent=14 at bit MSB. + // raw after multiply = -2*MSB + // goal = -16 + // shift = raw - goal = -2*MSB + 16 // We could keep more bits, but this approach is better conditioned. - s = (int64_t) v1 * v1 + (int64_t) v2 * v2 >> FP_MSB + 14; + s = (int64_t) v1 * v1 + (int64_t) v2 * v2 >> 2 * FP_MSB - 16; } while (s >= big || s <= small); - // log (s, 14, 14) / s -- Raw result of division has exponent=MSB - // Median absolute value of result is near 1 (ln(0.5)/0.5~=-1.4), so we want center power of 0, for exponent=15. - // Ideal shift is 15(=MSB-15), to put exponent=15 at MSB. - // We also multiply by 2, so claim exponent=16. - int multiplier = sqrt (((int64_t) log (s, 14, 14) << FP_MSB - 15) / -s, 16, 14); // multiplier has exponent=14; v1 and v2 have exponent=0 - nextGaussian = (int64_t) v2 * multiplier >> FP_MSB - 12; // product has exponent=14 at bit MSB*2; shift so exponent=2 at MSB + // log (s, -16, -16) / s -- Raw result of division has exponent= -16 - -16 = 0 + // Median absolute value of result is near 1 (ln(0.5)/0.5~=-1.4), so we want center power of 0. + // Want center bit in the middle of the word, position 15. Resulting exponent = -15. + // We also multiply by 2, so claim exponent=-14. + int multiplier = sqrt (((int64_t) log (s, -16, -16) << 15) / -s, -14, -16); // multiplier has exponent=-16; v1 and v2 have exponent=-MSB + nextGaussian = (int64_t) v2 * multiplier >> 18; // product has exponent=-16-MSB; shift so exponent=2-MSB haveNextGaussian = true; - return (int64_t) v1 * multiplier >> FP_MSB - 12; + return (int64_t) v1 * multiplier >> 18; } } @@ -418,7 +422,10 @@ template<> int gaussian (int sigma) { - return (int64_t) gaussian () * sigma >> 28; // ones bit of gaussian draw is at position MSB - 2 + // raw = 2-MSB + X + // goal = X + // shift = raw - goal = 2-MSB + return (int64_t) gaussian () * sigma >> FP_MSB - 2; } template<> @@ -426,7 +433,7 @@ MatrixFixed grid (int i, int nx, int ny, int nz) { MatrixFixed result = gridRaw (i, nx, ny, nz); - result[0] = (((int64_t) result[0] << 1) + 1 << FP_MSB) / nx; // exponentResult = -1 + result[0] = (((int64_t) result[0] << 1) + 1 << FP_MSB) / nx; // exponentResult = -1-MSB result[1] = (((int64_t) result[1] << 1) + 1 << FP_MSB) / ny; result[2] = (((int64_t) result[2] << 1) + 1 << FP_MSB) / nz; return result; @@ -448,16 +455,16 @@ pulse (int t, int width, int period, int rise, int fall) template<> int -unitmap (const MatrixAbstract & A, int row, int column) // row and column have exponent=0 +unitmap (const MatrixAbstract & A, int row, int column) // row and column have exponent=-MSB { // Just assume handle is good. - int rows = A.rows (); // exponent=MSB + int rows = A.rows (); // exponent=0 int columns = A.columns (); int lastRow = rows - 1; int lastColumn = columns - 1; - int64_t scaledRow = row * rows - (0x1 << FP_MSB - 1); // raw exponent = 0+MSB-MSB = 0 + int64_t scaledRow = row * rows - (0x1 << FP_MSB - 1); // raw exponent = 0 + -MSB = -MSB int64_t scaledColumn = column * columns - (0x1 << FP_MSB - 1); - int r = scaledRow >> FP_MSB; // to turn raw result into integer, shift = 0-MSB = -MSB + int r = scaledRow >> FP_MSB; // to turn raw result into integer, shift = -MSB - 0 = -MSB int c = scaledColumn >> FP_MSB; if (r < 0) { @@ -465,7 +472,7 @@ unitmap (const MatrixAbstract & A, int row, int column) // row and column else if (c >= lastColumn) return A(0,lastColumn); else { - int b = scaledColumn & 0x3FFFFFFF; // fractional part, with exponent = 0 (same as raw exponent) + int b = scaledColumn & 0x3FFFFFFF; // fractional part, with exponent = -MSB (same as raw exponent) int b1 = (1 << FP_MSB) - b; return (int64_t) b1 * A(0,c) + (int64_t) b * A(0,c+1) >> FP_MSB; } @@ -698,7 +705,11 @@ template T Part::getP () { +# ifdef n2a_FP + return 1 << FP_MSB2; +# else return 1; +# endif } template @@ -988,8 +999,10 @@ ConnectPopulation::reset (bool newOnly) if (newOnly) count = std::max (0, size - firstborn); else count = size; # ifdef n2a_FP - // raw multiply = -1+MSB-MSB = -1 - // shift = -1 - MSB + // raw multiply = -1-MSB + 0 = -1-MSB + // goal = 0 + // shift = -1-MSB - 0 = -1-MSB + // Could instead call multiplyRound() if (count > 1) i = (int) (((int64_t) uniform () * (count - 1) >> FP_MSB) + 1 >> 1); // Add 1 just below the decimal point to cause rounding. Total shift is -(MSB+1) # else if (count > 1) i = (int) round (uniform () * (count - 1)); @@ -1275,13 +1288,13 @@ Population::connect () outer->setProbe (c); while (outer->next ()) { - T create = c->getP (); + T p = c->getP (); // Yes, we need all 3 conditions. If create is 0 or 1, we do not do a random draw, since it would have no effect. - if (create <= 0) continue; + if (p <= 0) continue; # ifdef n2a_FP - if (create < 1 && create < uniform () >> 16) continue; + if (p < 1 && p < uniform () >> 1 + FP_MSB2) continue; // p exponent is -MSB/2; uniform() exponent is -1-MSB. shift = (-1-MSB) - -MSB/2 = -1 - MSB/2 # else - if (create < 1 && create < uniform ()) continue; + if (p < 1 && p < uniform ()) continue; # endif c->enterSimulation (); @@ -1512,7 +1525,7 @@ void Simulator::init (WrapperBase * wrapper) { # ifdef n2a_FP - EventStep * event = new EventStep (0, (1 << FP_MSB) / 10000); // Works for exponentTime=0. For any other case, it is necessary for top-level part to call setPeriod(). + EventStep * event = new EventStep (0, ((int64_t) 1 << -Event::exponent) / 1000); // 1ms step sizes, rather than 0.1ms like below. We are only promised (weakly) 10 bits below the decimal. # else EventStep * event = new EventStep ((T) 0, (T) 1e-4); # endif @@ -1778,7 +1791,7 @@ RungeKutta::run (Event & event) event.visit ([](Visitor * visitor) { visitor->part->finalizeDerivative (); - visitor->part->multiplyAddToStack (2 << FP_MSB - 1); // exponent=1, just enough to hold the values used by RK4 + visitor->part->multiplyAddToStack (2 << FP_MSB - 1); // exponent=1-MSB, just enough to hold the values used by RK4 }); } es.dt = dt; // restore original values diff --git a/N2A/src/gov/sandia/n2a/backend/vensim/Spreadsheet.java b/N2A/src/gov/sandia/n2a/backend/vensim/Spreadsheet.java index 46a5a261..b6fa8cfa 100644 --- a/N2A/src/gov/sandia/n2a/backend/vensim/Spreadsheet.java +++ b/N2A/src/gov/sandia/n2a/backend/vensim/Spreadsheet.java @@ -80,13 +80,13 @@ public void determineExponent (ExponentContext context) if (getKeyword ("info") == null) // normal mode. This includes "prefix" mode, but in that case we return a string, so don't care about exponent. { int centerNew = MSB / 2; - int exponentNew = getExponentHint (0) + MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } else // info mode { if (getType () instanceof Text) return; // If we return a string, leave exponent as unknown. - updateExponent (context, MSB, 0); // Return an integer + updateExponent (context, 0, 0); // Return an integer } } diff --git a/N2A/src/gov/sandia/n2a/eqset/EquationSet.java b/N2A/src/gov/sandia/n2a/eqset/EquationSet.java index 8e8aff4b..6f9b1bb6 100644 --- a/N2A/src/gov/sandia/n2a/eqset/EquationSet.java +++ b/N2A/src/gov/sandia/n2a/eqset/EquationSet.java @@ -3272,9 +3272,8 @@ public void assignParents () for (EquationSet s : parts) s.assignParents (); } - public void determineExponents () + public void determineExponents (ExponentContext context) { - ExponentContext context = new ExponentContext (this); determineExponentsInit (context); int limit = context.depth * 2 + 2; // One cycle for each variable to get initial exponent, and another cycle for each var to influence itself, plus a couple more for good measure. @@ -3303,7 +3302,7 @@ public void determineExponents () // List results of analysis to error stream PrintStream ps = Backend.err.get (); - ps.println ("Results of fixed-point analysis. Column 1 is expected median absolute value as a power of 10. Column 2 is binary power of MSB."); + ps.println ("Results of fixed-point analysis. Column 1 is expected median absolute value as a power of 10. Column 2 is binary power of LSB."); if (dumpMedians (ps, context)) throw new AbortRun (); } @@ -3328,15 +3327,10 @@ public ExponentContext (EquationSet root) { double duration = root.metadata.getOrDefault (0.0, "duration"); haveDuration = duration != 0; - if (haveDuration) exponentTime = (int) Math.floor (Math.log (duration) / Math.log (2)); + if (haveDuration) exponentTime = (int) Math.floor (Math.log (duration) / Math.log (2)) - Operator.MSB; else exponentTime = Operator.UNKNOWN; } - public ExponentContext (int exponentTime) - { - this.exponentTime = exponentTime; - } - public void updateTime () { if (haveDuration) return; // Don't override exponent based on duration. It is the most accurate. @@ -3348,12 +3342,12 @@ public void updateTime () { if (e.expression != null && e.expression.exponent != Operator.UNKNOWN) { - min = Math.min (min, e.expression.centerPower ()); + min = Math.min (min, e.expression.centerPower ()); // Mainly to accommodate constants, whose center encompasses all their significant bits. } } } - if (min == Integer.MAX_VALUE) exponentTime = 0; // If no value of $t' has been set yet, estimate duration as 1s, and $t has exponent=0. - else exponentTime = min + 20; // +20 allows one million minimally-sized timesteps, each with 10 bit resolution + if (min == Integer.MAX_VALUE) exponentTime = -Operator.MSB; // If no value of $t' has been set yet, estimate duration as 1s, and $t has exponent=-MSB. + else exponentTime = min - 10; // 10 bits to represent value of smallest timestep. Assuming MSB==30, this allows one million minimally-sized timesteps. } public boolean updateInputs () @@ -3361,27 +3355,40 @@ public boolean updateInputs () boolean changed = false; for (ArrayList list : inputs.values ()) { - int count = 0; - int pow = 0; + int count = 0; + int pow = 0; + int countRow = 0; + int powRow = 0; for (Input i : list) { + if (i.exponent != Operator.UNKNOWN) + { + count++; + pow += i.exponent; + } + if (i.operands.length <= 1) continue; Operator operand1 = i.operands[1]; if (operand1.exponent != Operator.UNKNOWN) { - count++; - pow += operand1.exponent; + countRow++; + powRow += operand1.exponent; } } - if (count > 0) + if (count > 0) pow /= count; + else pow = Operator.UNKNOWN; + if (countRow > 0) powRow /= countRow; + else powRow = Operator.UNKNOWN; + for (Input i : list) { - pow /= count; - for (Input i : list) + if (i.exponent != pow) { - if (i.exponentTime != pow) - { - i.exponentTime = pow; - changed = true; - } + i.exponent = pow; + changed = true; + } + if (i.exponentRow != powRow) + { + i.exponentRow = powRow; + changed = true; } } } @@ -3472,22 +3479,8 @@ public void determineExponentNext () /** Update the newly-created operators in a subset of variables that has been simplified for a specific phase. **/ - public static void determineExponentsSimplified (List list) + public static void determineExponentsSimplified (ExponentContext context, List list) { - // Determine exponentTime. - // This value is only used to update $t and $t', and it is available from either $t or $t', - // if one of them exists in the list. - int exponentTime = Operator.UNKNOWN; - for (Variable v : list) - { - if (v.name.equals ("$t")) - { - exponentTime = v.exponent; - break; - } - } - ExponentContext context = new ExponentContext (exponentTime); // Could be UNKNOWN, but if so it won't hurt anything, because it is only needed for $t. - // Update any newly-created operators. for (Variable v : list) { @@ -3577,7 +3570,7 @@ else if (op instanceof Power) // No need to warn about large v.exponent, because the user will be able to view the list. // Convert center power to an approximate decimal value. - int centerPower = v.exponent - Operator.MSB + v.center; + int centerPower = v.exponent + v.center; String base10 = Integer.toString ((int) Math.floor (centerPower / b2d)); String exponent = Integer.toString (v.exponent); if (v.exponent == Operator.UNKNOWN) diff --git a/N2A/src/gov/sandia/n2a/eqset/Variable.java b/N2A/src/gov/sandia/n2a/eqset/Variable.java index ca505730..4c96f16d 100644 --- a/N2A/src/gov/sandia/n2a/eqset/Variable.java +++ b/N2A/src/gov/sandia/n2a/eqset/Variable.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -910,14 +910,14 @@ public void determineExponent (ExponentContext context) } else if (name.equals ("$index") || name.equals ("$type")) // Variables stored as a pure integer. { - centerNew = 0; - exponentNew = Operator.MSB; + centerNew = 0; // Unless we have an estimate of population size, this is the best we can do. + exponentNew = 0; } else if (name.equals ("$connect") || name.equals ("$init") || name.equals ("$live")) { - // Booleans are stored as regular floats. - centerNew = Operator.MSB / 2; - exponentNew = Operator.MSB - centerNew; + // Booleans are stored as simple integers where possible. + centerNew = 0; + exponentNew = 0; } else { @@ -975,10 +975,10 @@ else if (name.equals ("$connect") || name.equals ("$init") || name.equals (" } else if (v.assignment == DIVIDE) { - int centerPower = v.exponent - Operator.MSB + v.center; + int centerPower = v.exponent + v.center; centerPower *= -1; // Negate to account for division 1/v cent += Operator.MSB / 2; // Fix quotient center at MSB/2 - pow += centerPower + Operator.MSB / 2; + pow += centerPower - Operator.MSB / 2; multiplicative = true; } else @@ -1021,7 +1021,7 @@ else if (derivative != null && derivative.exponent != Operator.UNKNOWN) else { int b = bound.exponent; // For expressions, which have variable magnitude. - if (bound instanceof Constant) b = bound.centerPower (); // Constants have exactly one magnitude (and center is always shifted to hold it), so use directly. + if (bound instanceof Constant) b = bound.centerPower () - Operator.MSB; // Constants have exactly one magnitude (and center is always shifted to hold it), so use directly. if (b > exponentNew) // Ensure we can accommodate max magnitude. { centerNew -= b - exponentNew; @@ -1049,23 +1049,51 @@ else if (derivative != null && derivative.exponent != Operator.UNKNOWN) } exponentNew = context.exponentTime; } + if (name.equals ("$p") && order == 0) + { + // The exponent of $p is hard-coded in the C runtime: getP() and connect(). + // If we know that both functions are overridden for this part, then we can let exponent be + // determined by equations writing $p. However, it is not worth the effort to check that. + // The hard-coded value (assuming 32-bit words) provides both sufficient precision and head-room. + centerNew = Operator.MSB / 2; + exponentNew = -Operator.MSB; + } } // User-specified exponent overrides any calculated value. - String magnitude = ""; - if (metadata != null) magnitude = metadata.get ("median"); - if (! magnitude.isEmpty ()) + if (metadata != null) { - if (exponent != Operator.UNKNOWN) return; // Already processed the hint. - try + String medianHint = metadata.get ("median"); + String centerHint = metadata.get ("center"); + boolean haveMedian = ! medianHint.isBlank (); + boolean haveCenter = ! centerHint.isBlank (); + if (haveMedian || haveCenter) { - double value = Operator.parse (magnitude).getDouble (); - int centerPower = 0; - if (value > 0) centerPower = (int) Math.floor (Math.log (value) / Math.log (2)); // log2 (value) - centerNew = Operator.MSB / 2; - exponentNew = centerPower + Operator.MSB - centerNew; + if (exponent != Operator.UNKNOWN) return; // Already processed the hint. + + if (haveCenter) + { + try + { + centerNew = Integer.valueOf (centerHint); + if (! haveMedian) exponentNew = -centerNew; // assume center power 0 + } + catch (NumberFormatException e) {} + } + + if (haveMedian) + { + try + { + double value = Operator.parse (medianHint).getDouble (); + int centerPower = 0; + if (value > 0) centerPower = (int) Math.floor (Math.log (value) / Math.log (2)); // log2 (value) + if (! haveCenter) centerNew = Operator.MSB / 2; + exponentNew = centerPower - centerNew; + } + catch (Exception e) {} + } } - catch (Exception e) {} } if (exponentNew != exponent || centerNew != center) changed = true; diff --git a/N2A/src/gov/sandia/n2a/language/AccessElement.java b/N2A/src/gov/sandia/n2a/language/AccessElement.java index fee88db8..e6e30729 100644 --- a/N2A/src/gov/sandia/n2a/language/AccessElement.java +++ b/N2A/src/gov/sandia/n2a/language/AccessElement.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -139,7 +139,7 @@ public void determineExponentNext () for (int i = 1; i < operands.length; i++) { Operator op = operands[i]; - op.exponentNext = MSB; // forces pure integer + op.exponentNext = 0; // forces pure integer op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/BuildMatrix.java b/N2A/src/gov/sandia/n2a/language/BuildMatrix.java index 06bcb94d..2dc49dee 100644 --- a/N2A/src/gov/sandia/n2a/language/BuildMatrix.java +++ b/N2A/src/gov/sandia/n2a/language/BuildMatrix.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -143,7 +143,7 @@ public Operator simplify (Variable from, boolean evalOnly) Matrix A = new MatrixDense (rows, cols); // potential constant to replace us boolean isConstant = true; // any element that is not constant will change this to false int cent = 0; // for fixed-point analysis - int pow = 0; + int pow = 0; // power of the bit at "cent" int count = 0; for (int c = 0; c < cols; c++) { @@ -202,12 +202,12 @@ public Operator simplify (Variable from, boolean evalOnly) cent /= count; pow /= count; result.center = cent; - result.exponent = pow + MSB - cent; + result.exponent = pow - cent; } else { result.center = MSB / 2; - result.exponent = MSB - result.center; + result.exponent = -result.center; } return result; } @@ -224,7 +224,7 @@ public void determineExponent (ExponentContext context) for (Operator e : c) { e.determineExponent (context); - if (! (e instanceof Constant) || e.getDouble () != 0) // avoid counting zeros + if (! (e instanceof Constant) || e.getDouble () != 0) // Avoid counting zeros. Operators that are not Constant are assumed to be nonzero. { cent += e.center; pow += e.exponent; @@ -240,7 +240,7 @@ public void determineExponent (ExponentContext context) else { cent = MSB / 2; - pow = MSB - cent; + pow = -cent; } updateExponent (context, pow, cent); } diff --git a/N2A/src/gov/sandia/n2a/language/Comparison.java b/N2A/src/gov/sandia/n2a/language/Comparison.java index 4fb973d7..ca5291c1 100644 --- a/N2A/src/gov/sandia/n2a/language/Comparison.java +++ b/N2A/src/gov/sandia/n2a/language/Comparison.java @@ -1,5 +1,5 @@ /* -Copyright 2017-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2017-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -47,7 +47,7 @@ public void determineExponent (ExponentContext context) operand0.determineExponent (context); operand1.determineExponent (context); if (operand0.exponent != UNKNOWN || operand1.exponent != UNKNOWN) alignExponent (context); - updateExponent (context, MSB, 0); // Output is 1 or 0 + updateExponent (context, 0, 0); // Output is 1 or 0, pure integer. // Any time a variable is compared to a value, it is a clue about the expected range of the variable. // The following two blocks apply this heuristic. diff --git a/N2A/src/gov/sandia/n2a/language/Constant.java b/N2A/src/gov/sandia/n2a/language/Constant.java index cc19370a..4dcf530e 100644 --- a/N2A/src/gov/sandia/n2a/language/Constant.java +++ b/N2A/src/gov/sandia/n2a/language/Constant.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -88,7 +88,7 @@ public void determineExponent (ExponentContext context) int e = 0; if (v != 0 && Double.isFinite (v)) e = Math.getExponent (v); - int exponentNew = e + MSB - centerNew; + int exponentNew = e - centerNew; updateExponent (context, exponentNew, centerNew); } // Matrix constants are built by BuildMatrix with their exponent and center values set correctly. @@ -104,11 +104,11 @@ public void determineExponent (ExponentContext context, int exponentOther) if (shift == 0) return; if (shift < 0) // down-shift { - // The mantissa of a float is 24 bits (1 implicit + 23 explicit). - // If this were aligned at MSB, we would have an extra (MSB+1)-24=MSB-23 zero bits beyond any zeros in the mantissa. - // Since the mantissa is actually aligned with center, we must subtract MSB-center bits from that count. - // available zero bits = zeros(mantissa)+MSB-23-(MSB-center) = zeros(mantissa)-23+center - int z = trailingZeros () - 23 + center; + // To avoid throwing away precision, we only shift until we run out of trailing zeros. + // How many trailing zeros do we have? + // Assuming 24 bits (1 implicit + 23 explicit) of a float, we have (center+1)-24 zeros + // that are just padding, plus any trailing zeros in the value itself. + int z = trailingZeros () + center - 23; z = Math.max (z, 0); // Don't allow negative z. This could happen if we are truncating some bits (center is less than 23, but there are no trailing zeros). shift = Math.max (shift, -z); } diff --git a/N2A/src/gov/sandia/n2a/language/Operator.java b/N2A/src/gov/sandia/n2a/language/Operator.java index a41e0c22..e247b519 100644 --- a/N2A/src/gov/sandia/n2a/language/Operator.java +++ b/N2A/src/gov/sandia/n2a/language/Operator.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -116,24 +116,23 @@ public class Operator implements Cloneable **/ public static int MSB = 30; /** - The power of bit that occupies the MSB position, before any shift to prepare value for use by the next operator. + The power of bit that occupies the LSB position, before any shift to prepare value for use by the next operator. In the fixed-point analysis implemented by Operator, all bits are fractional just like IEEE floats, - and we keep track of the power of the most significant bit, just like the IEEE float exponent. - In the case of a simple integer, exponent is equal to MSB. - We expect that the fixed-point implementation will do saturation checks, so we don't accommodate - the largest possible output of each operation, only the median. - Ideally, only the lower half of the available bits will be occupied. + but unlike IEEE float we keep track of the power of the least significant bit. That simplifies some calculations. + This value is equivalent to the amount of shift needed to make the bit in position zero have power zero. + In the case of a simple integer, exponent is 0. **/ public int exponent = UNKNOWN; /** Zero-based index of median magnitude. + We expect that the fixed-point implementation will do saturation checks, so we don't accommodate + the largest possible output of each operation, only the median. For numbers with a range of magnitudes, this will generally be MSB/2 (equivalent to Q16.15 format). It is expected that about half the time, all nonzero bits are at or below this position in the word. - Simple integers, on the other hand, will set this to 0, and all their nonzero bits are at or above this position. **/ public int center = MSB / 2; /** - The power of bit that occupies the MSB position, as required by the operator that contains us. + The power of bit that occupies the LSB position, as required by the operator that contains us. Used to determine shifts. **/ public int exponentNext = UNKNOWN; @@ -282,7 +281,7 @@ public void updateExponent (ExponentContext context, int exponentNew, int center public int centerPower () { - return exponent - MSB + center; + return exponent + center; } /** diff --git a/N2A/src/gov/sandia/n2a/language/Split.java b/N2A/src/gov/sandia/n2a/language/Split.java index 34764649..60166a81 100644 --- a/N2A/src/gov/sandia/n2a/language/Split.java +++ b/N2A/src/gov/sandia/n2a/language/Split.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -38,7 +38,7 @@ public void getOperandsFrom (SimpleNode node) throws ParseException public void determineExponent (ExponentContext context) { - updateExponent (context, MSB, 0); // integer + updateExponent (context, 0, 0); // integer } public void determineUnit (boolean fatal) throws Exception diff --git a/N2A/src/gov/sandia/n2a/language/function/Atan.java b/N2A/src/gov/sandia/n2a/language/function/Atan.java index 24b5a2fb..7c7f66d4 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Atan.java +++ b/N2A/src/gov/sandia/n2a/language/function/Atan.java @@ -1,5 +1,5 @@ /* -Copyright 2019-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2019-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -38,7 +38,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { for (int i = 0; i < operands.length; i++) operands[i].determineExponent (context); - updateExponent (context, 1, MSB - 1); // in [-pi,pi] + updateExponent (context, 1 - MSB, MSB - 1); // in [-pi,pi] } public void determineExponentNext () @@ -55,8 +55,8 @@ public void determineExponentNext () { // atan(y) = atan2(y,1), so treat x as 1 // If y is so small that all its bits are lost, then angle can be treated as zero. - // If y is so large that all the bits of x are lost, then angle can be treated as pi. - next = Math.max (0, y.exponent); + // If y is so large that all the bits of x are lost, then angle can be treated as pi/2 (1.57...). + next = Math.max (-MSB, y.exponent); } } else diff --git a/N2A/src/gov/sandia/n2a/language/function/Columns.java b/N2A/src/gov/sandia/n2a/language/function/Columns.java index a9e29007..ced6e36f 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Columns.java +++ b/N2A/src/gov/sandia/n2a/language/function/Columns.java @@ -1,5 +1,5 @@ /* -Copyright 2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2020-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -34,7 +34,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { - updateExponent (context, MSB, 0); // small integer + updateExponent (context, 0, 0); // small integer } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/Cosine.java b/N2A/src/gov/sandia/n2a/language/function/Cosine.java index 2e3a60c8..330ad29a 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Cosine.java +++ b/N2A/src/gov/sandia/n2a/language/function/Cosine.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -38,7 +38,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { operands[0].determineExponent (context); - updateExponent (context, 1, MSB - 2); // Largest absolute value is 1, but we allow one extra bit above the decimal for the convenience of the C implementation. + updateExponent (context, 1 - MSB, MSB - 2); // Largest absolute value is 1, but we allow one extra bit above the decimal for the convenience of the C implementation. } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/Draw.java b/N2A/src/gov/sandia/n2a/language/function/Draw.java index 2a79210a..d7473f05 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Draw.java +++ b/N2A/src/gov/sandia/n2a/language/function/Draw.java @@ -95,7 +95,7 @@ public void determineExponent (ExponentContext context) { op.determineExponent (context); } - updateExponent (context, MSB, 0); // Our output is always an integer (1). + updateExponent (context, 0, 0); // Our output is always an integer (1). if (keywords == null) return; for (Operator op : keywords.values ()) @@ -139,11 +139,11 @@ public void determineExponentNext () case "emission": case "specular": // exponent depends on whether color is an integer or a matrix - if (op.getType () instanceof Matrix) op.exponentNext = 0; - else op.exponentNext = MSB; + if (op.getType () instanceof Matrix) op.exponentNext = -MSB; + else op.exponentNext = 0; break; default: - op.exponentNext = MSB; // Suitable for integer or boolean values. + op.exponentNext = 0; // Suitable for integer or boolean values. } op.determineExponentNext (); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Draw2D.java b/N2A/src/gov/sandia/n2a/language/function/Draw2D.java index a11f4975..837d90a9 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Draw2D.java +++ b/N2A/src/gov/sandia/n2a/language/function/Draw2D.java @@ -17,7 +17,7 @@ public void determineExponentNext () // Last arg is color, which is always a raw integer. int last = operands.length - 1; Operator c = operands[last]; - c.exponentNext = MSB; + c.exponentNext = 0; c.determineExponentNext (); // All pixel-valued operands must agree on exponent. diff --git a/N2A/src/gov/sandia/n2a/language/function/Equal.java b/N2A/src/gov/sandia/n2a/language/function/Equal.java index fdb18332..15cf51df 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Equal.java +++ b/N2A/src/gov/sandia/n2a/language/function/Equal.java @@ -1,5 +1,5 @@ /* -Copyright 2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2023-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -71,7 +71,7 @@ public void determineExponent (ExponentContext context) if (op0 instanceof Constant) ((Constant) op0).determineExponent (context, op1.exponent); else if (op1 instanceof Constant) ((Constant) op1).determineExponent (context, op0.exponent); } - updateExponent (context, MSB, 0); // Output is 1 or 0 + updateExponent (context, 0, 0); // Output is 1 or 0 // Any time a variable is compared to a value, it is a clue about the expected range of the variable. // The following two blocks apply this heuristic. diff --git a/N2A/src/gov/sandia/n2a/language/function/Event.java b/N2A/src/gov/sandia/n2a/language/function/Event.java index 2f792235..ba133be1 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Event.java +++ b/N2A/src/gov/sandia/n2a/language/function/Event.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -47,7 +47,7 @@ public void determineExponent (ExponentContext context) // Our output is boolean, so use standard values for that. int cent = MSB / 2; - int pow = MSB - cent; + int pow = -cent; updateExponent (context, pow, cent); // Need to record the exponent for time, since we avoid passing context to determineExponentNext(). diff --git a/N2A/src/gov/sandia/n2a/language/function/Exp.java b/N2A/src/gov/sandia/n2a/language/function/Exp.java index ce65255d..c363ece0 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Exp.java +++ b/N2A/src/gov/sandia/n2a/language/function/Exp.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -49,8 +49,7 @@ public void determineExponent (ExponentContext context) // Currently, op is always signed. The best we can do is assume output is centered around 1. int centerNew = MSB / 2; - int exponentNew = getExponentHint (0); - exponentNew += MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } @@ -61,7 +60,7 @@ public void determineExponentNext () // Our goal with fixed-point is to roughly match the performance of single-precision float. // The largest number representable in float is 2^127, so allocating 8 bits above the decimal // should be more than enough. - op.exponentNext = 7; + op.exponentNext = 7 - MSB; op.determineExponentNext (); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Gaussian.java b/N2A/src/gov/sandia/n2a/language/function/Gaussian.java index 48111c45..3e3a6cd8 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Gaussian.java +++ b/N2A/src/gov/sandia/n2a/language/function/Gaussian.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -68,7 +68,7 @@ public void determineExponent (ExponentContext context) // exceed 6.66 standard deviations when using 32-bit uniform numbers. Thus, 7 std is safe. // log2(7)~=2.81, so magnitude of msb is 2 // Since about 68% of all results are less than 1 sigma, center can point to bit holding 2^-1. - updateExponent (context, 2, MSB - 3); + updateExponent (context, 2 - MSB, MSB - 3); } else { diff --git a/N2A/src/gov/sandia/n2a/language/function/Grid.java b/N2A/src/gov/sandia/n2a/language/function/Grid.java index 35de620f..6d9932ce 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Grid.java +++ b/N2A/src/gov/sandia/n2a/language/function/Grid.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -43,8 +43,8 @@ public void determineExponent (ExponentContext context) } boolean raw = getKeywordFlag ("raw"); - if (raw) updateExponent (context, MSB, 0); // integer - else updateExponent (context, -1, MSB - 1); // Since output never quite reaches 1, all bits can be fractional. + if (raw) updateExponent (context, 0, 0); // integer + else updateExponent (context, -1 - MSB, MSB - 1); // Since output never quite reaches 1, all bits can be fractional. } public void determineExponentNext () @@ -52,7 +52,7 @@ public void determineExponentNext () for (int i = 0; i < operands.length; i++) { Operator op = operands[i]; - op.exponentNext = MSB; // grid() requires integers + op.exponentNext = 0; // grid() requires integers op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/HyperbolicTangent.java b/N2A/src/gov/sandia/n2a/language/function/HyperbolicTangent.java index 87513e6a..ebdb2420 100644 --- a/N2A/src/gov/sandia/n2a/language/function/HyperbolicTangent.java +++ b/N2A/src/gov/sandia/n2a/language/function/HyperbolicTangent.java @@ -1,5 +1,5 @@ /* -Copyright 2020-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2020-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -37,7 +37,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { operands[0].determineExponent (context); - updateExponent (context, 0, MSB - 1); // result is always in [-1,1]. TODO: select center more carefully, based on center of operand + updateExponent (context, -MSB, MSB - 1); // result is always in [-1,1]. TODO: select center more carefully, based on center of operand } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/Input.java b/N2A/src/gov/sandia/n2a/language/function/Input.java index c71ab610..2d4ff06e 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Input.java +++ b/N2A/src/gov/sandia/n2a/language/function/Input.java @@ -36,7 +36,7 @@ public class Input extends Function { public boolean timeWarning; - public int exponentTime = UNKNOWN; // For C backend with integer math. The exponent used to convert time values to integer, both from the input file and from the caller. + public int exponentRow = UNKNOWN; // For C backend with integer math. The exponent used to convert row or time values into a floating-point number. public String name; // For C backend, the name of the InputHolder object. public String fileName; // For C backend, the name of the string variable holding the file name, if any. @@ -71,7 +71,7 @@ public void determineExponent (ExponentContext context) for (Operator op : operands) op.determineExponent (context); int centerNew = MSB / 2; - int exponentNew = getExponentHint (0) + MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } @@ -87,14 +87,14 @@ public void determineExponentNext () if (operands.length > 1) { Operator op = operands[1]; - if (usesTime ()) op.exponentNext = exponentTime; - else op.exponentNext = MSB; // We expect an integer. + if (usesTime ()) op.exponentNext = exponentRow; + else op.exponentNext = 0; // We expect an integer. op.determineExponentNext (); } if (operands.length > 2) { Operator op = operands[2]; - op.exponentNext = MSB; // We expect an integer. + op.exponentNext = 0; // We expect an integer. op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/Log.java b/N2A/src/gov/sandia/n2a/language/function/Log.java index 3f4ab79f..e74ab271 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Log.java +++ b/N2A/src/gov/sandia/n2a/language/function/Log.java @@ -1,5 +1,5 @@ /* -Copyright 2017-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2017-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -51,7 +51,7 @@ public void determineExponent (ExponentContext context) int p = 0; if (o != 0) p = (int) (Math.signum (o) * Math.round (Math.log (Math.abs (o)) / Math.log (2))); int centerNew = MSB / 2; - int exponentNew = p + MSB - centerNew; + int exponentNew = p - centerNew; updateExponent (context, exponentNew, centerNew); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Mcount.java b/N2A/src/gov/sandia/n2a/language/function/Mcount.java index 14e373bc..3c3dd8ac 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Mcount.java +++ b/N2A/src/gov/sandia/n2a/language/function/Mcount.java @@ -34,7 +34,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { - updateExponent (context, MSB, 0); // small integer + updateExponent (context, 0, 0); // small integer } public void determineExponentNext () @@ -44,7 +44,7 @@ public void determineExponentNext () for (int i = 1; i < operands.length; i++) { Operator op = operands[i]; - op.exponentNext = MSB; + op.exponentNext = 0; op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/Mfile.java b/N2A/src/gov/sandia/n2a/language/function/Mfile.java index d33f4f68..9e10d9df 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Mfile.java +++ b/N2A/src/gov/sandia/n2a/language/function/Mfile.java @@ -52,7 +52,7 @@ public void determineExponent (ExponentContext context) // Set exponent based on hint. Suitable for numerics outputs. Does no harm for string outputs. int centerNew = MSB / 2; - int exponentNew = getExponentHint (0) + MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } @@ -65,7 +65,7 @@ public void determineExponentNext () for (int i = 1; i < operands.length; i++) { Operator op = operands[i]; - op.exponentNext = MSB; + op.exponentNext = 0; op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/Norm.java b/N2A/src/gov/sandia/n2a/language/function/Norm.java index 4f5a35d4..49b1e829 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Norm.java +++ b/N2A/src/gov/sandia/n2a/language/function/Norm.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -85,7 +85,7 @@ public Type get (VariableReference r) throws EvaluationException { // Result is an integer centerNew = 0; - exponentNew = MSB; + exponentNew = 0; } else if (Double.isInfinite (n)) { @@ -105,7 +105,7 @@ public void determineExponentNext () op0.determineExponentNext (); Operator op1 = operands[1]; // n - op1.exponentNext = Operator.MSB / 2; // required by fixedpoint in C backend + op1.exponentNext = -MSB / 2; // required by fixedpoint in C backend op1.determineExponentNext (); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Output.java b/N2A/src/gov/sandia/n2a/language/function/Output.java index 08ce1df3..1ca34024 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Output.java +++ b/N2A/src/gov/sandia/n2a/language/function/Output.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -100,7 +100,7 @@ public void determineExponentNext () if (operands.length > 2) { Operator op2 = operands[2]; - op2.exponentNext = 15; // To avoid passing extra exponent parameter, we use a general-purpose value. + op2.exponentNext = 0; // If index, it will be truncated to an integer. op2.determineExponentNext (); } @@ -386,7 +386,12 @@ else if (cols == 1) public String getColumnName (Instance context) { // Explicit column name - if (hasColumnName) return operands[2].eval (context).toString (); + if (hasColumnName) + { + Type name = operands[2].eval (context); + if (name instanceof Scalar) return String.valueOf ((long) ((Scalar) name).value); + return name.toString (); + } // Auto-generate column name if (context instanceof InstanceTemporaries) context = ((InstanceTemporaries) context).wrapped; diff --git a/N2A/src/gov/sandia/n2a/language/function/Pulse.java b/N2A/src/gov/sandia/n2a/language/function/Pulse.java index 6fcee9ce..5715ca50 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Pulse.java +++ b/N2A/src/gov/sandia/n2a/language/function/Pulse.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -40,7 +40,7 @@ public void determineExponent (ExponentContext context) { operands[i].determineExponent (context); } - updateExponent (context, 0, MSB - 2); // Output reaches exactly 1. + updateExponent (context, -MSB, MSB - 2); // Output reaches exactly 1. } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/ReadImage.java b/N2A/src/gov/sandia/n2a/language/function/ReadImage.java index a447973a..e0d53277 100644 --- a/N2A/src/gov/sandia/n2a/language/function/ReadImage.java +++ b/N2A/src/gov/sandia/n2a/language/function/ReadImage.java @@ -1,5 +1,5 @@ /* -Copyright 2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2022-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -67,7 +67,7 @@ public boolean canBeInitOnly () public void determineExponent (ExponentContext context) { // Color channels are generally in [0,1] - updateExponent (context, 0, MSB - 1); + updateExponent (context, -MSB, MSB - 1); // operands[0] is the file string // operands[1] is the color string diff --git a/N2A/src/gov/sandia/n2a/language/function/ReadMatrix.java b/N2A/src/gov/sandia/n2a/language/function/ReadMatrix.java index c210207f..b3e41659 100644 --- a/N2A/src/gov/sandia/n2a/language/function/ReadMatrix.java +++ b/N2A/src/gov/sandia/n2a/language/function/ReadMatrix.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -57,7 +57,7 @@ public boolean isMatrixInput () public void determineExponent (ExponentContext context) { int centerNew = MSB / 2; - int exponentNew = getExponentHint (0) + MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Round.java b/N2A/src/gov/sandia/n2a/language/function/Round.java index 7e5c58a7..b3b189a0 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Round.java +++ b/N2A/src/gov/sandia/n2a/language/function/Round.java @@ -1,5 +1,5 @@ /* -Copyright 2016-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2016-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -51,10 +51,10 @@ public static void determineExponentStatic (Function f, ExponentContext context) int centerPower = Math.max (0, op.centerPower ()); // because we always output an integer int pow = op.exponent; - int cent = MSB - (pow - centerPower); // We trust that op has its center positioned within the range [0,MSB], so (pow - centerPower) <= MSB. - if (pow < 0) + int cent = centerPower - pow; // We trust that op has its center positioned within the range [0,MSB], so (centerPower - pow) >= 0. + if (pow < -MSB) { - pow = 0; + pow = -MSB; cent = MSB; } f.updateExponent (context, pow, cent); @@ -71,9 +71,9 @@ Shared by round(), ceil() and floor() public static void determineExponentNextStatic (Function f, int exponentNext) { Operator op = f.operands[0]; - if (op.exponent < 0) op.exponentNext = 0; // Must have at least one bit above the decimal point in order to round. - else if (op.exponent < MSB) op.exponentNext = op.exponent; // Decimal point is visible, so we can process this. - else op.exponentNext = exponentNext; // Otherwise, just pass through. + if (op.exponent < -MSB) op.exponentNext = -MSB; // Must have at least one bit above the decimal point in order to round. + else if (op.exponent < 0 ) op.exponentNext = op.exponent; // Decimal point is visible, so we can process this. + else op.exponentNext = exponentNext; // Otherwise, just pass through. op.determineExponentNext (); f.exponent = op.exponentNext; } diff --git a/N2A/src/gov/sandia/n2a/language/function/Rows.java b/N2A/src/gov/sandia/n2a/language/function/Rows.java index 125ed5b5..176e9f19 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Rows.java +++ b/N2A/src/gov/sandia/n2a/language/function/Rows.java @@ -1,5 +1,5 @@ /* -Copyright 2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2020-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -34,7 +34,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { - updateExponent (context, MSB, 0); // small integer + updateExponent (context, 0, 0); // small integer } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/Signum.java b/N2A/src/gov/sandia/n2a/language/function/Signum.java index 7f990c73..bb9290f6 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Signum.java +++ b/N2A/src/gov/sandia/n2a/language/function/Signum.java @@ -1,5 +1,5 @@ /* -Copyright 2017-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2017-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -40,7 +40,7 @@ public void determineExponent (ExponentContext context) Operator op = operands[0]; op.determineExponent (context); int centerNew = MSB / 2; - int exponentNew = MSB - centerNew; + int exponentNew = -centerNew; updateExponent (context, exponentNew, centerNew); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Sine.java b/N2A/src/gov/sandia/n2a/language/function/Sine.java index 0c6f9937..c1ebd46a 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Sine.java +++ b/N2A/src/gov/sandia/n2a/language/function/Sine.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -38,7 +38,7 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { operands[0].determineExponent (context); - updateExponent (context, 1, MSB - 2); // Largest absolute value is 1, but we allow one extra bit above the decimal for the convenience of the C implementation. + updateExponent (context, 1 - MSB, MSB - 2); // Largest absolute value is 1, but we allow one extra bit above the decimal for the convenience of the C implementation. } public void determineExponentNext () diff --git a/N2A/src/gov/sandia/n2a/language/function/Sphere.java b/N2A/src/gov/sandia/n2a/language/function/Sphere.java index 41113fa7..cee40445 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Sphere.java +++ b/N2A/src/gov/sandia/n2a/language/function/Sphere.java @@ -1,5 +1,5 @@ /* -Copyright 2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2022-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -54,7 +54,7 @@ public Operator simplify (Variable from, boolean evalOnly) Operator sigma = new Constant (new MatrixDense (3, 1, 1)); sigma.parent = this; sigma.center = MSB / 2; // Constant matrix should always carry fixed-point information. - sigma.exponent = MSB - sigma.center; + sigma.exponent = -sigma.center; operands[0] = sigma; return this; } @@ -134,7 +134,7 @@ public void determineExponent (ExponentContext context) { if (operands.length == 0) { - updateExponent (context, -1, MSB); // Current implementation never quite reaches 1, only [0,1). + updateExponent (context, -1 - MSB, MSB); // Current implementation never quite reaches 1, only [0,1). } else { diff --git a/N2A/src/gov/sandia/n2a/language/function/SquareRoot.java b/N2A/src/gov/sandia/n2a/language/function/SquareRoot.java index a6c151cd..32646935 100644 --- a/N2A/src/gov/sandia/n2a/language/function/SquareRoot.java +++ b/N2A/src/gov/sandia/n2a/language/function/SquareRoot.java @@ -1,5 +1,5 @@ /* -Copyright 2018-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2018-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -46,7 +46,7 @@ public void determineExponent (ExponentContext context) if (op.exponent == UNKNOWN) return; int pow = op.centerPower () / 2; int centerNew = MSB / 2; - int exponentNew = pow + MSB - centerNew; + int exponentNew = pow - centerNew; updateExponent (context, exponentNew, centerNew); } diff --git a/N2A/src/gov/sandia/n2a/language/function/SumSquares.java b/N2A/src/gov/sandia/n2a/language/function/SumSquares.java index a6ebc493..acb9dd53 100644 --- a/N2A/src/gov/sandia/n2a/language/function/SumSquares.java +++ b/N2A/src/gov/sandia/n2a/language/function/SumSquares.java @@ -1,5 +1,5 @@ /* -Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2021-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -57,7 +57,7 @@ public Type get (VariableReference r) throws EvaluationException { int shift = (int) Math.floor (Math.log (Asize) / Math.log (2)); // log2(Asize) int centerNew = MSB / 2; - int exponentNew = op0.centerPower () * 2 + shift + MSB - centerNew; + int exponentNew = op0.centerPower () * 2 + shift - centerNew; updateExponent (context, exponentNew, centerNew); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/Tangent.java b/N2A/src/gov/sandia/n2a/language/function/Tangent.java index b19aadc5..b38dbacf 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Tangent.java +++ b/N2A/src/gov/sandia/n2a/language/function/Tangent.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -42,7 +42,7 @@ public void determineExponent (ExponentContext context) op.determineExponent (context); int centerNew = MSB / 2; - int exponentNew = getExponentHint (0) + MSB - centerNew; + int exponentNew = getExponentHint (0) - centerNew; updateExponent (context, exponentNew, centerNew); } diff --git a/N2A/src/gov/sandia/n2a/language/function/Uniform.java b/N2A/src/gov/sandia/n2a/language/function/Uniform.java index 3b428135..32fa9ab1 100644 --- a/N2A/src/gov/sandia/n2a/language/function/Uniform.java +++ b/N2A/src/gov/sandia/n2a/language/function/Uniform.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -63,7 +63,7 @@ public void determineExponent (ExponentContext context) { if (operands.length == 0) { - updateExponent (context, -1, MSB); // Current implementation never quite reaches 1, only [0,1). + updateExponent (context, -1 - MSB, MSB); // Current implementation never quite reaches 1, only [0,1). return; } diff --git a/N2A/src/gov/sandia/n2a/language/function/UnitMap.java b/N2A/src/gov/sandia/n2a/language/function/UnitMap.java index 28e5b407..74210bc2 100644 --- a/N2A/src/gov/sandia/n2a/language/function/UnitMap.java +++ b/N2A/src/gov/sandia/n2a/language/function/UnitMap.java @@ -1,5 +1,5 @@ /* -Copyright 2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2021-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -64,7 +64,7 @@ public void determineExponentNext () for (int i = 1; i < operands.length; i++) { Operator op = operands[i]; - op.exponentNext = 0; // a number in [0,1] with some provision for going slightly out of bounds + op.exponentNext = -MSB; // a number in [0,1] with some provision for going slightly out of bounds op.determineExponentNext (); } } diff --git a/N2A/src/gov/sandia/n2a/language/function/glFrustum.java b/N2A/src/gov/sandia/n2a/language/function/glFrustum.java index 736e39cf..1a65504c 100644 --- a/N2A/src/gov/sandia/n2a/language/function/glFrustum.java +++ b/N2A/src/gov/sandia/n2a/language/function/glFrustum.java @@ -63,8 +63,8 @@ public static void determineExponent (Function f, ExponentContext context) // Impose bounds on exponent. This will potentially make the resulting matrix lose significance. // There's nothing we can do about that. Fixed-point has limited capacity to represent a projection matrix. int shift = 0; - if (pow < 0 ) shift = pow; // Must be able to represent 1 for rotation matrix. - else if (pow > MSB - 8) shift = pow - MSB + 8; // Must have at least 8 bits below decimal to have a crude rotation matrix. + if (pow < -MSB) shift = pow + MSB; // Must be able to represent 1 for rotation matrix. + else if (pow > -8 ) shift = pow + 8; // Must have at least 8 bits below decimal to have a crude rotation matrix. cent += shift; pow -= shift; diff --git a/N2A/src/gov/sandia/n2a/language/function/glRotate.java b/N2A/src/gov/sandia/n2a/language/function/glRotate.java index 1f114ba1..86a10bcb 100644 --- a/N2A/src/gov/sandia/n2a/language/function/glRotate.java +++ b/N2A/src/gov/sandia/n2a/language/function/glRotate.java @@ -37,13 +37,13 @@ public Operator createInstance () public void determineExponent (ExponentContext context) { for (Operator op : operands) op.determineExponent (context); - updateExponent (context, 0, MSB - 1); // A pure rotation matrix only needs to represent 1. + updateExponent (context, -MSB, MSB - 1); // A pure rotation matrix only needs to represent 1. // No relevant keyword arguments. } public void determineExponentNext () { - exponentNext = exponent; // exponent is assigned 0 above + exponentNext = exponent; // exponent is assigned -MSB above // Determine average exponent of operands, then force them to agree. int avg = 0; diff --git a/N2A/src/gov/sandia/n2a/language/operator/AND.java b/N2A/src/gov/sandia/n2a/language/operator/AND.java index 48b88e6c..88c3703e 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/AND.java +++ b/N2A/src/gov/sandia/n2a/language/operator/AND.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -76,7 +76,7 @@ public void determineExponent (ExponentContext context) { operand0.determineExponent (context); operand1.determineExponent (context); - updateExponent (context, MSB, 0); + updateExponent (context, 0, 0); } public void determineUnit (boolean fatal) throws Exception diff --git a/N2A/src/gov/sandia/n2a/language/operator/Add.java b/N2A/src/gov/sandia/n2a/language/operator/Add.java index 032221a1..d14f357c 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Add.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Add.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -87,6 +87,18 @@ public void determineExponent (ExponentContext context) int c0 = operand0.center - (pow - operand0.exponent); int c1 = operand1.center - (pow - operand1.exponent); int cent = Math.max (c0, c1); + int min = Math.min (c0, c1); + if (cent >= MSB) + { + // The most optimistic thing we could hope for is that value only exceeds center by 1 bit. + pow += cent - (MSB - 1); + cent = MSB - 1; + } + else if (min < 0) + { + pow += min; // decreases pow + cent -= min; // increases cent + } updateExponent (context, pow, cent); } else if (operand0.exponent != UNKNOWN) diff --git a/N2A/src/gov/sandia/n2a/language/operator/Divide.java b/N2A/src/gov/sandia/n2a/language/operator/Divide.java index 1da13c1c..806863af 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Divide.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Divide.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -67,7 +67,7 @@ public void determineExponent (ExponentContext context) { int cent = MSB / 2; int pow = operand0.centerPower () - operand1.centerPower (); - pow += MSB - cent; + pow -= cent; updateExponent (context, pow, cent); } // else don't to propagate a bad guess. Instead, hope that better information arrives in next cycle. diff --git a/N2A/src/gov/sandia/n2a/language/operator/Modulo.java b/N2A/src/gov/sandia/n2a/language/operator/Modulo.java index 5f05c60e..9bc43fd1 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Modulo.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Modulo.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -47,7 +47,7 @@ public void determineExponent (ExponentContext context) // The most precise answer is smaller than the magnitude of either operand. int pow = Math.min (operand0.exponent, operand1.exponent); int centerPower = Math.min (operand0.centerPower (), operand1.centerPower ()); - int cent = MSB - (pow - centerPower); + int cent = centerPower - pow; updateExponent (context, pow, cent); } else if (operand0.exponent != UNKNOWN) diff --git a/N2A/src/gov/sandia/n2a/language/operator/Multiply.java b/N2A/src/gov/sandia/n2a/language/operator/Multiply.java index 8de62101..1d1d5245 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Multiply.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Multiply.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -91,7 +91,7 @@ public void determineExponent (ExponentContext context) { int cent = MSB / 2; int pow = operand0.centerPower () + operand1.centerPower (); - pow += MSB - cent; + pow -= cent; updateExponent (context, pow, cent); } // else don't to propagate a bad guess. Instead, hope that better information arrives in next cycle. diff --git a/N2A/src/gov/sandia/n2a/language/operator/MultiplyElementwise.java b/N2A/src/gov/sandia/n2a/language/operator/MultiplyElementwise.java index 1890cfa9..1694a481 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/MultiplyElementwise.java +++ b/N2A/src/gov/sandia/n2a/language/operator/MultiplyElementwise.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -66,7 +66,7 @@ public void determineExponent (ExponentContext context) { int cent = MSB / 2; int pow = operand0.centerPower () + operand1.centerPower (); - pow += MSB - cent; + pow -= cent; updateExponent (context, pow, cent); } // else don't to propagate a bad guess. Instead, hope that better information arrives in next cycle. diff --git a/N2A/src/gov/sandia/n2a/language/operator/NOT.java b/N2A/src/gov/sandia/n2a/language/operator/NOT.java index 32c741b1..f6d95332 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/NOT.java +++ b/N2A/src/gov/sandia/n2a/language/operator/NOT.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -53,10 +53,12 @@ public Operator simplify (Variable from, boolean evalOnly) Constant c = (Constant) result; if (c.value instanceof Matrix) { - // See determineExponent() below + // Divide 1/operand + // center power = 0 - (operand center power) + // shift so center is at MSB/2 if (operand.exponent == UNKNOWN) return result; result.center = MSB / 2; - result.exponent = 0 - operand.centerPower () + MSB - result.center; + result.exponent = 0 - operand.centerPower () - result.center; } } return result; @@ -74,12 +76,12 @@ public void determineExponent (ExponentContext context) if (operand.exponent == UNKNOWN) return; int cent = MSB / 2; int pow = 0 - operand.centerPower (); // See Divide class. We're treating this as 1/A, where 1 has center power 0. - pow += MSB - cent; + pow -= cent; updateExponent (context, pow, cent); } else // Logical not { - updateExponent (context, MSB, 0); + updateExponent (context, 0, 0); } } diff --git a/N2A/src/gov/sandia/n2a/language/operator/OR.java b/N2A/src/gov/sandia/n2a/language/operator/OR.java index 7f791868..952514fc 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/OR.java +++ b/N2A/src/gov/sandia/n2a/language/operator/OR.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -76,7 +76,7 @@ public void determineExponent (ExponentContext context) { operand0.determineExponent (context); operand1.determineExponent (context); - updateExponent (context, MSB, 0); + updateExponent (context, 0, 0); } public void determineUnit (boolean fatal) throws Exception diff --git a/N2A/src/gov/sandia/n2a/language/operator/Power.java b/N2A/src/gov/sandia/n2a/language/operator/Power.java index 61ed0682..e80c6408 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Power.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Power.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2022 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -151,7 +151,7 @@ public void determineExponent (ExponentContext context) } if (exponentNew != UNKNOWN) { - exponentNew += MSB - centerNew; + exponentNew -= centerNew; updateExponent (context, exponentNew, centerNew); } } @@ -159,7 +159,7 @@ public void determineExponent (ExponentContext context) public void determineExponentNext () { operand0.exponentNext = operand0.exponent; - operand1.exponentNext = MSB / 2; // Exponentiation is very sensitive, so no benefit in allowing arbitrary size of input. + operand1.exponentNext = -MSB / 2; // Exponentiation is very sensitive, so no benefit in allowing arbitrary size of input. operand0.determineExponentNext (); operand1.determineExponentNext (); } diff --git a/N2A/src/gov/sandia/n2a/language/operator/Subtract.java b/N2A/src/gov/sandia/n2a/language/operator/Subtract.java index 1e922d0d..29b17c26 100644 --- a/N2A/src/gov/sandia/n2a/language/operator/Subtract.java +++ b/N2A/src/gov/sandia/n2a/language/operator/Subtract.java @@ -1,5 +1,5 @@ /* -Copyright 2013-2021 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +Copyright 2013-2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. */ @@ -106,6 +106,18 @@ public void determineExponent (ExponentContext context) int c0 = operand0.center - (pow - operand0.exponent); int c1 = operand1.center - (pow - operand1.exponent); int cent = Math.max (c0, c1); + int min = Math.min (c0, c1); + if (cent >= MSB) + { + // The most optimistic thing we could hope for is that value only exceeds center by 1 bit. + pow += cent - (MSB - 1); + cent = MSB - 1; + } + else if (min < 0) + { + pow += min; // decreases pow + cent -= min; // increases cent + } updateExponent (context, pow, cent); } else if (operand0.exponent != UNKNOWN)