From 191f348318e85711f9c8c39eae593310f0a375bc Mon Sep 17 00:00:00 2001 From: Fred Rothganger Date: Mon, 13 May 2024 11:29:07 -0600 Subject: [PATCH] C: fixed-point analysis now tracks exponent of LSB rather than MSB This is in preparation for working with different word sizes in same model. Would also allow the whole runtime to be compiled for a different word size. --- N2A/src/gov/sandia/n2a/backend/c/JobC.java | 92 ++-- .../gov/sandia/n2a/backend/c/RendererC.java | 7 +- .../gov/sandia/n2a/backend/c/RendererCfp.java | 38 +- .../n2a/backend/c/runtime/fixedpoint.cc | 446 ++++++++++-------- .../sandia/n2a/backend/c/runtime/holder.cc | 25 +- .../gov/sandia/n2a/backend/c/runtime/holder.h | 7 +- .../sandia/n2a/backend/c/runtime/holder.tcc | 91 ++-- .../gov/sandia/n2a/backend/c/runtime/matrix.h | 6 +- .../gov/sandia/n2a/backend/c/runtime/mymath.h | 22 +- .../sandia/n2a/backend/c/runtime/runtime.cc | 15 +- .../sandia/n2a/backend/c/runtime/runtime.h | 22 +- .../sandia/n2a/backend/c/runtime/runtime.tcc | 77 +-- .../n2a/backend/vensim/Spreadsheet.java | 4 +- N2A/src/gov/sandia/n2a/eqset/EquationSet.java | 73 ++- N2A/src/gov/sandia/n2a/eqset/Variable.java | 68 ++- .../sandia/n2a/language/AccessElement.java | 4 +- .../gov/sandia/n2a/language/BuildMatrix.java | 12 +- .../gov/sandia/n2a/language/Comparison.java | 4 +- N2A/src/gov/sandia/n2a/language/Constant.java | 14 +- N2A/src/gov/sandia/n2a/language/Operator.java | 19 +- N2A/src/gov/sandia/n2a/language/Split.java | 4 +- .../sandia/n2a/language/function/Atan.java | 8 +- .../sandia/n2a/language/function/Columns.java | 4 +- .../sandia/n2a/language/function/Cosine.java | 4 +- .../sandia/n2a/language/function/Draw.java | 8 +- .../sandia/n2a/language/function/Draw2D.java | 2 +- .../sandia/n2a/language/function/Equal.java | 4 +- .../sandia/n2a/language/function/Event.java | 4 +- .../gov/sandia/n2a/language/function/Exp.java | 7 +- .../n2a/language/function/Gaussian.java | 4 +- .../sandia/n2a/language/function/Grid.java | 8 +- .../language/function/HyperbolicTangent.java | 4 +- .../sandia/n2a/language/function/Input.java | 10 +- .../gov/sandia/n2a/language/function/Log.java | 4 +- .../sandia/n2a/language/function/Mcount.java | 4 +- .../sandia/n2a/language/function/Mfile.java | 4 +- .../sandia/n2a/language/function/Norm.java | 6 +- .../sandia/n2a/language/function/Output.java | 11 +- .../sandia/n2a/language/function/Pulse.java | 4 +- .../n2a/language/function/ReadImage.java | 4 +- .../n2a/language/function/ReadMatrix.java | 4 +- .../sandia/n2a/language/function/Round.java | 14 +- .../sandia/n2a/language/function/Rows.java | 4 +- .../sandia/n2a/language/function/Signum.java | 4 +- .../sandia/n2a/language/function/Sine.java | 4 +- .../sandia/n2a/language/function/Sphere.java | 6 +- .../n2a/language/function/SquareRoot.java | 4 +- .../n2a/language/function/SumSquares.java | 4 +- .../sandia/n2a/language/function/Tangent.java | 4 +- .../sandia/n2a/language/function/Uniform.java | 4 +- .../sandia/n2a/language/function/UnitMap.java | 4 +- .../n2a/language/function/glFrustum.java | 4 +- .../n2a/language/function/glRotate.java | 4 +- .../gov/sandia/n2a/language/operator/AND.java | 4 +- .../gov/sandia/n2a/language/operator/Add.java | 14 +- .../sandia/n2a/language/operator/Divide.java | 4 +- .../sandia/n2a/language/operator/Modulo.java | 4 +- .../n2a/language/operator/Multiply.java | 4 +- .../operator/MultiplyElementwise.java | 4 +- .../gov/sandia/n2a/language/operator/NOT.java | 12 +- .../gov/sandia/n2a/language/operator/OR.java | 4 +- .../sandia/n2a/language/operator/Power.java | 6 +- .../n2a/language/operator/Subtract.java | 14 +- 63 files changed, 715 insertions(+), 573 deletions(-) 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)