How do I define a new numeric constant in Mathematica?

Question

What is the best way to define a numerical constant in Mathematica?

For example, say I want g to be the approximate acceleration due to gravity on the surface of the Earth. I give it a numerical value (in m/s^2), tell Mathematica it's numeric, positive and a constant using

Unprotect[g];
ClearAll[g]
N[g] = 9.81;
NumericQ[g] ^= True;
Positive[g] ^= True;
SetAttributes[g, Constant];
Protect[g];

Then I can use it as a symbol in symbolic calculations that will automatically evaluate to 9.81 when numerical results are called for. For example 1.0 g evaluates to 9.81.

This does not seem as well tied into Mathematica as built in numerical constants. For example Pi > 0 will evaluate to True, but g > 0 will not. (I could add g > 0 to the global $Assumptions but even then I need a call to Simplify for it to take effect.) Also, Positive[g] returns True, but Positive[g^2] does not evaluate - compare this with the equivalent statements using Pi.

So my question is, what else should I do to define a numerical constant? What other attributes/properties can be set? Is there an easier way to go about this? Etc...

Answer 1

I'd recommend using a zero-argument "function". That way it can be given both the NumericFunction attribute and a numeric evaluation rule. that latter is important for predicates such as Positive.

SetAttributes[gravUnit, NumericFunction]
N[gravUnit[], prec_: $MachinePrecision] := N[981/100, prec]

In[121]:= NumericQ[gravitUnit[]]
Out[121]= True

In[122]:= Positive[gravUnit[]^2 - 30]
Out[122]= True

Daniel Lichtblau

Answer 2

May be I am naive, but to my mind your definitions are a good start. Things like g > 0->True can be added via UpValues. For Positive[g^2] to return True, you probably have to overload Positive, because of the depth-1 limitation for UpValues. Generally, I think the exact set of auto-evaluated expressions involving a constant is a moving target, even for built-in constants. In other words, those extra built-in rules seem to be determined from convenience and frequent uses, on a case-by-case basis, rather than from the first principles. I would just add new rules as you go, whenever you feel that you need them. You probably can not expect your constants to be as well integrated in the system as built-ins, but I think you can get pretty close. You will probably have to overload a number of built-in functions on these symbols, but again, which ones those will be, will depend on what you need from your symbol.

EDIT

I was hesitating to include this, since the code below is a hack, but it may be useful in some circumstances. Here is the code:

Clear[evalFunction];
evalFunction[fun_Symbol, HoldComplete[sym_Symbol]] := False;

Clear[defineAutoNValue];
defineAutoNValue[s_Symbol] :=
  Module[{inSUpValue},
    s /: expr : f_[left___, s, right___] :=
      Block[{inSUpValue = True},
        With[{stack = Stack[_]},
          If[
            expr === Unevaluated[expr] &&
               (evalFunction[f, HoldComplete[s]] ||
                  MemberQ[
                    stack,
                    HoldForm[(op_Symbol /; evalFunction[op, HoldComplete[s]])
                       [___, x_ /; ! FreeQ[Unevaluated[x], HoldPattern@expr], ___]],
                    Infinity
                  ]
               ),
            f[left, N[s], right],
            (* else *)
            expr
      ]]] /; ! TrueQ[inSUpValue]];

ClearAll[substituteNumeric];
SetAttributes[substituteNumeric, HoldFirst];
substituteNumeric[code_, rules : {(_Symbol :> {__Symbol}) ..}] :=
  Internal`InheritedBlock[{evalFunction},
     MapThread[
       Map[Function[f, evalFunction[f, HoldComplete[#]] = True], #2] &,
       Transpose[List @@@ rules]
     ];
     code]

With this, you may enable a symbol to auto-substitute its numerical value in places where we indicate some some functions surrounding those function calls may benefit from it. Here is an example:

ClearAll[g, f];
SetAttributes[g, Constant];
N[g] = 9.81;
NumericQ[g] ^= True;
defineAutoNValue[g];
f[g] := "Do something with g";

Here we will try to compute some expressions involving g, first normally:

In[391]:= {f[g],g^2,g^2>0, 2 g, Positive[2 g+1],Positive[2g-a],g^2+a^2,g^2+a^2>0,g<0,g^2+a^2<0}
Out[391]= {Do something with g,g^2,g^2>0,2 g,Positive[1+2 g],
   Positive[-a+2 g],a^2+g^2,a^2+g^2>0,g<0,a^2+g^2<0}

And now inside our wrapper (the second argument gives a list of rules, to indicate for which symbols which functions, when wrapped around the code containing those symbols, should lead to those symbols being replaced with their numerical values):

In[392]:= 
substituteNumeric[{f[g],g^2,g^2>0, 2 g, Positive[2 g+1],Positive[2g-a],g^2+a^2,g^2+a^2>0,
      g<0,g^2+a^2<0},
     {g:>{Positive,Negative,Greater}}]

Out[392]= {Do something with g,g^2,True,2 g,True,Positive[19.62\[VeryThinSpace]-a],
a^2+g^2,96.2361\[VeryThinSpace]+a^2>0,g<0,a^2+g^2<0}

Since the above is a hack, I can not guarantee anything about it. It may be useful in some cases, but that must be decided on a case-by-case basis.