I’m a mechanical engineer and, although I studied for software engineering degrees right after I graduated, I did my fair share of numerical calculations during my five years of studies in engineering. Some of it in Turbo Pascal, but most of it in the recommended language, which was Fortran. It was very spartan but, one thing it did well, was get out of the way between engineers and their calculations. Around 2007, 17 full years away from purely numerical computing, I turned around to see what had changed in the domain of calculating with physical quantities. Believe me, any new version of Fortran that has come out was a huge improvement over Fortran 77, but I found that not much had changed in the calculations themselves. In particular, regarding calculations with physical units and proper dimensional analysis. The kind that could avoid disasters like these.
Before I tell you what I concocted in 2007 and then in 2014, let me stay a little longer on Fortran, or rather what was designed to be a modern-day Fortran: Fortress.
Fortress
Fortress, the brainchild of Guy Steele, was an experimental project in Sun to answer the HPCS challenge (cf. “Fwd: HPCS recap”). After it was not selected in subsequent rounds of the challenge, the project persisted for a total of ten years, before being axed by Oracle. People explain things like that because of the “pragmatic” approach of Oracle towards Java. Which is a polite way to say that Oracle makes Java walk the streets for money, and probably gives her a good beatin’ whenever she comes back with low earnings.
Fortress had nailed the dimension thing, offering dimensions (decribed as “type-like constructs that are separate from other types”) and units (described as “special constructs that modify types and values […] that are instances of dimensions”). Dimensionless and dimensioned numbers can use the whole gamut of language constructs. In addition to standard SI dimensions and units, many more are available in the standard library and more can be added.
Here is what some dimension-aware calculations would look in Fortress (Fortress is a language written to be read back as close to mathematics as possible, so this is how Fortress code would look like).
x:ℝ64 Length = 1.3m
t:ℝ64 Time = 5 s
v:ℝ64 Velocity = x/t
w:ℝ64 Velocity in nm/s = 17nm/s
Standard SI prefixes like (nano, micro etc) and their abbreviations (n, μ etc) can be used. In order to avoid misunderstandings, these prefixes cannot be used in front of units of dimension Information (bit and byte).
Both dimensions and units can be multiplied, divided, and raised to rational powers to produce new dimensions and units. Multiplication can be expressed using juxtaposition or · (yes, Fortress actually gives the natural meaning of multiplication to juxtaposition!) Division can be expressed using / or per. The syntactic operators square and cubic (or the special postfix syntactic operators squared and cubed) may be applied to a dimension or unit in order
to raise it to the second power, third power, respectively. The syntactic operator inverse may be applied to a dimension or unit to divide it into 1. Ex. grams per cubic centimeter, meter per second squared, inverse ohms.
All in all, a pretty solid treatment. Unfortunately, none of this was actually implemented… However, there is one existing language which does support dimensioned calculations in a very real way.
F#
F# is the only mainstream language, as far as I know, that contains specific constructs for working with units (Curl supposedly does also, but I’m not familiar with it and it’s a long stretch to call it mainstream).
It is not entirely pleasing as, there is no concept of dimension, only the concept of unit of measure, but all the machinery needed for safe dimensional calculations is there. As is common in the .Net and Java world, the units of measure are defined using an attribute (annotation in Java-speak). SI predefined “measures” can be found in the PowerPack.
[<Measure>] type kg
[<Measure>] type s
[<Measure>] type m
[<Measure>] type cm = m / 100
[<Measure>] type ml = cm^3
let gravityOnEarth = 9.81<m/s^2>
let heightOfFall = 15<m>
let speedOfImpact = sqrt(2.0 * gravityOnEarth * heightOfFall)
let someForce = 10<kg m s^-2>
As you see in the example, units can be defined by their own (in which case they can be thought as defining a dimension) or as derived from some unit expression. Such expressions can be suffixed to numeric literals between angle brackets (without any intervening space). As a nod to Fortress, who knows?, juxtaposition stands for multiplication in the numerator of unit expressions. A dimensionless quantity is denoted by the unit expression <1>. One can even have generic units! And, you know, .Net has real multidimensional array types, and value types and also arrays of value types. The whole makes a compelling proposal for a scientist or an engineer.
Dimenso
Back in 2007, none of that was there. So I set myself a goal of experimenting in that domain, leading to the open-source project 
. I focused on just SI quantities and then, just the MKS part together with current (no temperature, no amount of substance and no luminous intensity).
Dimenso: C++
The prototypical implementation was in C++. Believe it or not, the template language in C++ is pretty advanced, even with today’s standards. One could really make a case that C++ template metaprogramming is akin to using dependent types. Because templates in C++, unlike generics in the likes of Java or .Net, can actually be parameterized by integers, in addition to types. Here’s how the definition of the class that embodies dimensioned values would start: template<int L, int M, int T, int C> class DimensionedValue.
Defining what it means to multiply two dimensioned values together is this simple operator (yes! Java-lovers, most other languages do have operator overloads!).
template<int L2, int M2, int T2, int C2>
DimensionedValue<L+L2,M+M2,T+T2,C+C2> operator*(DimensionedValue<L2,M2,T2,C2> other) {
return DimensionedValue<L+L2,M+M2,T+T2,C+C2>(magnitude()*other.magnitude());
}
Given all the machinery available in C++, one can really make these objects look and behave like plain numbers. The challenge, then, was to go as close to that in C# and Java, which lacked support for dependent types.
Dimenso: C# and Java
The way I followed was: code generation. I generated all dimension combinations within specific bounds (but I also provided the Python script if one would need a more extended set). These combinations were named after the corresponding exponents (for example, type m_2k0s1a0 represents dimension L-2T) but actual names from physics were used, when available (but not in a right way… if you try to use Dimenso, stick to the ugly names).
This is how an arithmetic operation on “length” would look like in Java, if I hadn’t messed up the name thing.
public velocity divide(time v2) {
velocity ret = new velocity(this.value/v2.getValue());
return ret;
}
In C#, which does have operator overloads and value types, the whole is a little more palatable. But still, generating hundreds of classes feels… kind of wrong. Let’s see what the Java of today has to say on the matter.
JSR-275, JSR-363 and JScience
To make a long story short, JSR-275 “Measures and Units” was finally rejected and, its successor, JSR-363 “Units of Measurement API” has not really formed yet but, abandoning at that point would make it really short. That’s why I’m going to look at support for units in JScience, which claims support for such.
Let’s look at some sample code, which will highlight the point I was trying to make with Dimenso.
// Compiler Warnings
Unit<Acceleration> mps2 = Unit.valueOf("m/s²"); // Compiler Warning.
Unit<Power> WATT = JOULE.divide(SECOND); // Compiler Warning.
// Adding run-time check of dimension consistency removes the warning.
Unit<Acceleration> mps2 = Unit.valueOf("m/s²")
.asType(Acceleration.class);
Unit<Power> WATT = JOULE.divide(SECOND).asType(Power.class);
// Use wildcard parameterized type (no warning and no check).
Unit<?> kelvinPerSec = KELVIN.divide(SECOND);
JScience has a lot of capabilities for units, including arbitrary systems of dimensions, localizable unit parsing and output and compound units, like those for common time. But there is no way to have the compiler figure out what units two quantities multiplied or divided among themselves have. This is a shortcoming of Java, not of JScience. And, unless you’ve been ruined by too many years of Java, you can tell that no scientist or engineer would put up with so much unwarranted verbosity to do his job.
Do not despair, as Java is not the only choice you have on the JVM.
Scala
Let’s see what I’ve been able to do with some typelevel hacking in Scala. First of all, I’d like to point you to a better solution, Metascala Units, which did come up in my Google search, and a tab in my browser was open, waiting patiently for me to have a look at them, which I did, but after the fact. No matter, because what I really wanted was some practical experience with typelevel programming, with dimension-aware calculations as a good case to work on. The way that Metascala encodes integers in the typesystem is a lot more concise and a lot more elegant than my own, which is based on Shapeless’ encoding of naturals.
Here is some safe code using Metascala Units:
val dist : Length = m(2.3)
val time : Time = s(1.7)
val x = dist * time
val speed : Speed = dist / time
I hope you can appreciate the fact that this looks like using built-in features. It is as concise as it can get and offers full compile-time checking.
And here’s my own stab at this. It’s not part of Dimenso yet, because it’s not fully tested, I believe that if I build stronger typelevel kung-fu, I’ll be able to handle output of units at compile-time, as well, and there is probably some more magic to add with implicit value classes.
case class DimVal[M1 <: ZInt, K1 <: ZInt, S1 <: ZInt, A1 <: ZInt](val value : Double) {
type T = DimVal[M1,K1,S1,A1]
def +(b: DimVal[M1,K1,S1,A1]) = new DimVal[M1,K1,S1,A1](value + b.value)
def -(b: DimVal[M1,K1,S1,A1]) = new DimVal[M1,K1,S1,A1](value - b.value)
def *[M2 <: ZInt, M <: ZInt, K2 <: ZInt, K <: ZInt, S2 <: ZInt, S <: ZInt, A2 <: ZInt, A <: ZInt]
(b: DimVal[M2,K2,S2,A2])
(implicit ev: ZSumAux[M1,M2,M], ev2: ZSumAux[K1,K2,K], ev3: ZSumAux[S1,S2,S], ev4: ZSumAux[A1,A2,A]) : DimVal[M,K,S,A] =
new DimVal[M,K,S,A](value * b.value)
def /[M2 <: ZInt, M <: ZInt, K2 <: ZInt, K <: ZInt, S2 <: ZInt, S <: ZInt, A2 <: ZInt, A <: ZInt]
(b: DimVal[M2,K2,S2,A2])
(implicit ev: ZDiffAux[M1,M2,M], ev2: ZDiffAux[K1,K2,K], ev3: ZDiffAux[S1,S2,S], ev4: ZDiffAux[A1,A2,A]) : DimVal[M,K,S,A] =
new DimVal[M,K,S,A](value / b.value)
def asString () (implicit i1: ToZInt[M1], i2: ToZInt[K1], i3: ToZInt[S1], i4: ToZInt[A1]) = {
val x1 = i1.apply()
val x2 = i2.apply()
val x3 = i3.apply()
val x4 = i4.apply()
val u =
(x1,x2,x3,x4) match {
case (0,0,1,1) => " C"
case (1,1,-2,0) => " N"
// ... more special cases to handle derived units elided ...
case _ => {
val sb = new StringBuilder
// ... ugly, imperative code elided ...
if(sb.length > 0) " "+sb else ""
}
}
value.toString + u
}
}
The dimension calculations are done at the type level, using types ZSumAux and ZDiffAux, and type ToZInt handles conversion to language-level integers. Like I said, not all induction cases are tested, but here’s what they look like now.
import shapeless._
import shapeless.Nat._0
trait ZInt
case class ZPos[T <: Nat]() extends ZInt
case class ZNeg[T <: Nat]() extends ZInt
trait ZSumAux[A <: ZInt, B <: ZInt, C <: ZInt]
trait ZDiffAux[A <: ZInt, B <: ZInt, C <: ZInt]
object ZSumAux {
implicit def sum3[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[A, B, C]) = new ZSumAux[ZPos[A], ZPos[B], ZPos[C]] {}
implicit def sum4[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[A, B, C]) = new ZSumAux[ZNeg[A], ZNeg[B], ZNeg[C]] {}
implicit def sum5[A <: Nat] = new ZSumAux[ZPos[A], ZNeg[A], ZPos[_0]] {}
implicit def sum6[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[A, B, C], ev2: LT[B,A]) = new ZSumAux[ZPos[A], ZNeg[B], ZPos[C]] {}
implicit def sum7[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[B, A, C], ev2: LT[A,B]) = new ZSumAux[ZPos[A], ZNeg[B], ZNeg[C]] {}
implicit def sum8[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[A, B, C], ev2: LT[B,A]) = new ZSumAux[ZNeg[A], ZPos[B], ZNeg[C]] {}
implicit def sum9[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[B, A, C], ev2: LT[A,B]) = new ZSumAux[ZNeg[A], ZPos[B], ZPos[C]] {}
}
object ZDiffAux {
implicit def diff3[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[A, B, C], ev2: LT[B,A]) = new ZDiffAux[ZPos[A], ZPos[B], ZPos[C]] {}
implicit def diff3b[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[B, A, C], ev2: LT[A,B]) = new ZDiffAux[ZPos[A], ZPos[B], ZNeg[C]] {}
implicit def diff4[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[A, B, C], ev2: LT[B,A]) = new ZDiffAux[ZNeg[A], ZNeg[B], ZNeg[C]] {}
implicit def diff4b[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : DiffAux[B, A, C], ev2: LT[A,B]) = new ZDiffAux[ZNeg[A], ZNeg[B], ZPos[C]] {}
implicit def diff5[A <: Nat] = new ZDiffAux[ZPos[A], ZPos[A], ZPos[_0]] {}
implicit def diff5b[A <: Nat] = new ZDiffAux[ZNeg[A], ZNeg[A], ZPos[_0]] {}
implicit def diff6[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[A, B, C], ev2: LT[B,A]) = new ZDiffAux[ZPos[A], ZNeg[B], ZPos[C]] {}
implicit def diff7[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[B, A, C], ev2: LT[A,B]) = new ZDiffAux[ZPos[A], ZNeg[B], ZNeg[C]] {}
implicit def diff8[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[A, B, C], ev2: LT[B,A]) = new ZDiffAux[ZNeg[A], ZPos[B], ZNeg[C]] {}
implicit def diff9[A <: Nat, B <: Nat, C <: Nat]
(implicit ev : SumAux[B, A, C], ev2: LT[A,B]) = new ZDiffAux[ZNeg[A], ZPos[B], ZPos[C]] {}
}
trait ToZInt[Z <: ZInt] {
def apply() : Int
}
object ToZInt {
implicit def toPosInt[N <: Nat] (implicit i: ToInt[N]) = new ToZInt[ZPos[N]] {
def apply() = i.apply()
}
implicit def toNegInt[N <: Nat] (implicit i: ToInt[N]) = new ToZInt[ZNeg[N]] {
def apply() = - i.apply()
}
}
Some useful definitions are in object DimValOps.
object DimValOps {
type Z0 = ZPos[_0]
type Z1 = ZPos[Succ[_0]]
type Z2 = ZPos[Succ[Succ[_0]]]
type Z3 = ZPos[Succ[Succ[Succ[_0]]]]
type Z4 = ZPos[Succ[Succ[Succ[Succ[_0]]]]]
type Z_1 = ZNeg[Succ[_0]]
type Z_2 = ZNeg[Succ[Succ[_0]]]
type Z_3 = ZNeg[Succ[Succ[Succ[_0]]]]
type length = DimVal[Z1, Z0, Z0, Z0]
type mass = DimVal[Z0, Z1, Z0, Z0]
type time = DimVal[Z0, Z0, Z1, Z0]
type current = DimVal[Z0, Z0, Z0, Z1]
type charge = DimVal[Z0, Z0, Z1, Z1]
type area = DimVal[Z2, Z0, Z1, Z0]
type volume = DimVal[Z3, Z0, Z1, Z0]
type velocity = DimVal[Z1, Z0, Z_1, Z0]
type acceleration = DimVal[Z1, Z0, Z_2, Z0]
type momentum = DimVal[Z1, Z1, Z_1, Z0]
type angular_momentum = DimVal[Z2, Z1, Z_1, Z0]
type frequency = DimVal[Z0, Z0, Z_1, Z0]
type force = DimVal[Z1, Z1, Z_2, Z0]
type pressure = DimVal[Z_1, Z1, Z_2, Z0]
type viscosity = DimVal[Z_1, Z1, Z_1, Z0]
type energy = DimVal[Z2, Z1, Z_2, Z0]
type power = DimVal[Z2, Z1, Z_3, Z0]
type potential = DimVal[Z2, Z1, Z_3, Z_1]
type capacitance = DimVal[Z_2, Z_1, Z4, Z2]
type resistance = DimVal[Z2, Z1, Z_3, Z_2]
type conductance = DimVal[Z_2, Z_1, Z3, Z2]
type magnetic_flux = DimVal[Z2, Z1, Z_2, Z_1]
type magnetic_flux_density = DimVal[Z0, Z1, Z_2, Z_1]
type inductance = DimVal[Z2, Z1, Z_2, Z_2]
val m = new length(1) //meter
val kg = new mass(1) //kilogram
val s = new time(1) //second
val A = new current(1) //Ampere
val C = new charge(1) //Coulomb
val N = new force(1) //Newton
val Hz = new frequency(1) //Hertz
val Pa = new pressure(1) //Pascal
val ohm = new resistance(1)
val V = new potential(1) //Volt
val S = new conductance(1) //Siemens
val W = new power(1) //Watt
val J = new energy(1) //Joule
val F = new capacitance(1) //Farad
val H = new inductance(1) //Henry
val Wb = new magnetic_flux(1) //Weber
val T = new magnetic_flux_density(1) //Tesla
implicit def numToDim(d : Double) : DimVal[Z0,Z0,Z0,Z0] = new DimVal[Z0,Z0,Z0,Z0](d)
}
On the one hand, it looks frightening. On the other hand, it’s just a few lines of code to add a capability that is missing since forever. Here is some putting it through its paces.
val c0 = new charge(4)
println("c0="+c0.asString())
val v = 3 * m / s
println("v="+v.asString)
val mom = 10 * kg * v
println("mom="+mom.asString())
val e = mom * v
println("e="+e.asString())
val p = e / (3 * s)
println("p="+p.asString())
Which outputs:
c0=4.0 C
v=3.0 m·s-1
mom=30.0 m·kg·s-1
e=90.0 J
p=30.0 W
Look what happens when you try to compile an assignment with mismatched dimensions.
val p : pressure = e / (3 * s)
Error:(37, 26) type mismatch;
found : ShapelessExperiment.DimVal[ShapelessExperiment.ZPos[shapeless.Succ[shapeless.Succ[shapeless.Nat._0]]],ShapelessExperiment.ZPos[shapeless.Succ[shapeless.Nat._0]],ShapelessExperiment.ZNeg[shapeless.Succ[shapeless.Succ[shapeless.Succ[shapeless.Nat._0]]]],ShapelessExperiment.ZPos[shapeless.Nat._0]]
required: ShapelessExperiment.DimValOps.pressure
(which expands to) ShapelessExperiment.DimVal[ShapelessExperiment.ZNeg[shapeless.Succ[shapeless.Nat._0]],ShapelessExperiment.ZPos[shapeless.Succ[shapeless.Nat._0]],ShapelessExperiment.ZNeg[shapeless.Succ[shapeless.Succ[shapeless.Nat._0]]],ShapelessExperiment.ZPos[shapeless.Nat._0]]
val p : pressure = e / (3 * s)
^
Ok, the message is not entirely simple but, if you understand the type signatures, it is clear how to read it. In conclusion, even though Scala’s compile-time metaprogramming has been said to imitate, rather than implement, dependent types, you can see it is able to quack like a duck and walk like a duck. And it is perfectly capable of stepping in when Java proves too little, even though it cannot overcome fundamental limitations of the JVM.
dsouflis 