1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
% Macros for dealing with very large and very small numbers in `Mlog form'.
% A number x in Mlog form represents mexp(x) if x is an even multiple of
% epsilon and -mexp(x) if x is an odd multiple of epsilon, where epsilon=1/65536
% is the basic unit for MetaPost's fixpoint numbers. Such numbers can represent
% values large as 3.8877e+55 or as small as 1.604e-28. (Anything less than that
% is treated as zero.)
% Mlog_str <string> convert a string like "6.02e23" to Mlog form
% Mexp_str <M_numeric> convert from Mlog form to a string like "6.02e23"
% Meform(q) find (x,y) such that q is the Mlog form of x*10^y
% Mlog <numeric> convert a number to Mlog form
% Mexp <M_numeric> convert from Mlog form into ordinary numeric form
% Mabs <M_numeric> absolute value
% Mlog_Str <string or numeric> convert to Mlog form; unknown results in unknown
% <M_numeric> Madd <M_numeric> add
% <M_numeric> Msub <M_numeric> subtract
% <M_numeric> Mmul <M_numeric> multiply
% <M_numeric> Mdiv <M_numeric> divide
% <M_numeric> Mleq <M_numeric> the boolean operator <=
% Mzero constant that represents zero
% Mten constant that represents 10.0
% All other externally visible names start with `M' and end with `_'.
warningcheck := 0; % Need warningcheck:=0 anytime these macros are in use
if unknown Mzero: % If we have not already seen this file
newinternal Mzero;
Mzero := -16384; % Anything at least this small is treated as zero
input string.mp
% Ideal value of mlog(10) is 589.4617838 or 38630967.46/65536. To get an even
% numerator, we round to 38630968/65536 or 589.4617920.
newinternal Mten;
Mten := 589.46179;
% Convert x*10^y to Mlog form, assuming x is already in Mlog form.
primarydef x M_e_ y = (x+Mten*y) enddef;
def Mgobble_ text t = enddef;
pair M_iv_; M_iv_=(0,4);
% String s is a number in scientific notatation with no leading spaces.
% Convert it to Mlog form.
vardef Mlog_str primary s =
save k, t, e, r;
string t;
if substring(0,1) of s="-":
r=-epsilon; t=substring (1,infinity) of s;
else:
r=0; t=s;
fi
let e = Mgobble_;
if begingroup scantokens substring M_iv_ of t endgroup >= 1000:
k = cspan(t,isdigit);
t := substring M_iv_ of t & "." & substring(4,k) of t &
substring (k if substring(k,k+1) of t=".": +1 fi, infinity) of t;
r := r + (k-4)*Mten;
fi
let e = M_e_;
r + Mlog scantokens t
enddef;
% Convert x from Mlog form into a pair whose xpart gives a mantissa and whose
% ypart gives a power of ten. The constant 1768.38985 is slightly more than
% 3Mten as is required to ensure that 3Mten-epsilon<=q+e*Mten<4Mten-epsilon.
% This forces the xpart of the result to be between 1000 and 10000.
vardef Meform(expr q) =
if q<=Mzero: (0,0)
else:
save e; e=floor((q-1768.38985)/Mten);
(Mexp(q-e*Mten), e)
fi
enddef;
% Perform the inverse of Mlog_str, converting from Mlog form to a string.
vardef Mexp_str primary x =
save p;
pair p; p=Meform(x);
decimal xpart p
if ypart p<>0: & "e" & decimal ypart p fi
enddef;
vardef Mabs primary x = x*.5*2 enddef;
% Convert a number to Mlog form
vardef Mlog primary x =
if x>0: Mabs mlog x
elseif x<0: epsilon + Mabs mlog(-x)
else: Mzero
fi
enddef;
% Convert a number from Mlog form
vardef Mexp primary x =
if x=Mabs x: mexp x
else: -mexp x
fi
enddef;
primarydef a Mmul b =
if (a<=Mzero) or (b<=Mzero): Mzero
else: (a+b)
fi
enddef;
primarydef a Mdiv b =
if a<=Mzero: Mzero
else: (a-b)
fi
enddef;
% 256ln(1579)=123556596.0003/65536=1885.3240356445
secondarydef a Madd b =
if a>=b: (Mlog(1579 + Mexp(b Mmul (1885.32404-a))) + a-1885.32404)
else: b Madd a
fi
enddef;
secondarydef a Msub b = a Madd (b-epsilon) enddef;
tertiarydef a Mleq b =
(if a=Mabs a:
if b=Mabs b: a<=b else: false fi
elseif b=Mabs b: true else: b<=a fi
)
enddef;
vardef Mlog_Str primary x =
if unknown x: whatever
elseif string x: Mlog_str x
else: Mlog x
fi
enddef;
fi % end the if..fi that encloses this file
|
