@@ -34,16 +34,35 @@ IEEEEvaluateWithAbsoluteError::usage =
3434"intervals (assumed to not carry any error), or a list of two intervals " <>
3535"for a value and its absolute error."
3636IEEEEvaluateWithAbsoluteError ::argnum =
37- "IEEEEvaluate called with `1` arguments; 1 argument is expected." ;
37+ "IEEEEvaluateWithAbsoluteError called with `1` arguments; 1 argument is expected." ;
3838IEEEEvaluateWithAbsoluteError ::badass =
3939"IEEEEvaluateWithAbsoluteError does not support associativity, expressions must be " <>
4040"parenthesized."
4141
4242
43+ IEEEEvaluateWithRelativeError ;
44+ IEEEEvaluateWithRelativeError ::usage =
45+ "IEEEEvaluateWithRelativeeError[\!\(\* StyleBox[\" x\" ,\n FontSlant->\" Italic\" ]\) ] " <>
46+ "evaluates \!\(\* StyleBox[\" x\" ,\n FontSlant->\" Italic\" ]\) using the rules of " <>
47+ "IEEE arithmetic, and propagates the relative error bounds at each stage of " <>
48+ "the computation. Returns a list of two intervals: the first one is the " <>
49+ "interval of the result evaluated with proper IEEE rounding, the second one is " <>
50+ "the range of the relative error on the result. The argument is an expression " <>
51+ "which can include numbers (assumed to be exact after correct rounding), " <>
52+ "intervals (assumed to not carry any error), or a list of two intervals " <>
53+ "for a value and its relative error." ;
54+ IEEEEvaluateWithRelativeError ::argnum =
55+ "IEEEEvaluateWithRelativeError called with `1` arguments; 1 argument is expected." ;
56+ IEEEEvaluateWithRelativeError ::badass =
57+ "IEEEEvaluateWithRelativeError does not support associativity, expressions must be " <>
58+ "parenthesized."
59+
60+
4361UseFMA ;
4462UseFMA ::usage =
45- "UseFMA is an option for IEEEEvaluate and IEEEEvaluateWithAbsoluteError that " <>
46- "specifies whether to use FMA for expressions of the form \!\(\*
63+ "UseFMA is an option for IEEEEvaluate, IEEEEvaluateWithAbsoluteError, and " <>
64+ "IEEEEvaluateWithRelativeError that specifies whether to use FMA for " <>
65+ "expressions of the form \!\(\*
4766StyleBox[\" a\" ,\n FontSlant->\" Italic\" ]\) * \!\(\*
4867StyleBox[\" b\" ,\n FontSlant->\" Italic\" ]\) + \!\(\*
4968StyleBox[\" c\" ,\n FontSlant->\" Italic\" ]\) ." ;
@@ -91,37 +110,37 @@ errorBelow1:=If[
91110
92111
93112(* ::Text:: *)
94- (*Returns the error bound for the largest element (in absolute value) of its argument. The returned value is positive, regardless of the sign of the argument. The argument is an interval or an unbound variable . Note that if the largest element is a power of two, the error bound is the one below that power of two.*)
113+ (*Returns the error bound for the largest element (in absolute value) of its argument. The returned value is positive, regardless of the sign of the argument. The argument is an interval. Note that if the largest element is a power of two, the error bound is the one below that power of two.*)
95114
96115
97116(* Would want to write Max[Log2[...]] below but that doesn't work somehow. *)
98117absoluteErrorBound [x_ ]:= Block [{exponent = Log2 [Max [Abs [x ]]]},errorBelow1 2 ^ Ceiling [exponent ]];
99118
100119
101120(* Special case because for an interval x*x^2 is not x^3. *)
102- applyOp [cube ,{va_ ,\[Delta ]a_ }]:= Block [
121+ applyOpWithAbsoluteError [cube ,{va_ ,\[Delta ]a_ }]:= Block [
103122 {va2 ,\[Delta ]a2 ,h ,r ,\[Delta ]r },
104- {va2 ,\[Delta ]a2 }= applyOp [# ^ 2 & ,{va ,\[Delta ]a }];
123+ {va2 ,\[Delta ]a2 }= applyOpWithAbsoluteError [# ^ 2 & ,{va ,\[Delta ]a }];
105124 r = va ^ 3 ; (* Wrong! *)
106125 r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
107126 h = absoluteErrorBound [r ];
108127 \[Delta ]r = ReplaceAll [Expand [(v2 + \[Delta ]2 )(v + \[Delta ])- (v2 v )],{v -> va ,\[Delta ]-> \[Delta ]a ,v2 -> va2 ,\[Delta ]2 -> \[Delta ]a2 }]+ Interval [{- h ,h }];
109128 {r ,\[Delta ]r }];
110- applyOp [op_ ,{va_ ,\[Delta ]a_ }]:= Block [
129+ applyOpWithAbsoluteError [op_ ,{va_ ,\[Delta ]a_ }]:= Block [
111130 {h ,r ,\[Delta ]r },
112131 r = op [va ];
113132 r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
114133 h = absoluteErrorBound [r ];
115134 \[Delta ]r = ReplaceAll [Expand [op [v + \[Delta ]]- op [v ]],{v -> va ,\[Delta ]-> \[Delta ]a }]+ Interval [{- h ,h }];
116135 {r ,\[Delta ]r }];
117- applyOp [op_ ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ }]:= Block [
136+ applyOpWithAbsoluteError [op_ ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ }]:= Block [
118137 {h ,r ,\[Delta ]r },
119138 r = op [va ,vb ];
120139 r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
121140 h = absoluteErrorBound [r ];
122141 \[Delta ]r = ReplaceAll [Expand [op [v1 + \[Delta ]1 ,v2 + \[Delta ]2 ]- op [v1 ,v2 ]],{v1 -> va ,\[Delta ]1 -> \[Delta ]a ,v2 -> vb ,\[Delta ]2 -> \[Delta ]b }]+ Interval [{- h ,h }];
123142 {r ,\[Delta ]r }];
124- applyOp [op_ ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ },{vc_ ,\[Delta ]c_ }]:= Block [
143+ applyOpWithAbsoluteError [op_ ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ },{vc_ ,\[Delta ]c_ }]:= Block [
125144 {h ,r ,\[Delta ]r },
126145 r = op [va ,vb ,vc ];
127146 r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
@@ -141,27 +160,140 @@ Block[
141160SetAttributes [evae ,HoldAll ];
142161evae [a_ * b_ + c_ ]:= If [
143162usefma ,
144- applyOp [#1 #2 + #3 & ,evae [a ],evae [b ],evae [c ]],
145- applyOp [Plus ,applyOp [Times ,evae [a ],evae [b ]],evae [c ]]];
146- evae [a_ + b_ ]:= applyOp [Plus ,evae [a ],evae [b ]];
163+ applyOpWithAbsoluteError [#1 #2 + #3 & ,evae [a ],evae [b ],evae [c ]],
164+ applyOpWithAbsoluteError [Plus ,applyOpWithAbsoluteError [Times ,evae [a ],evae [b ]],evae [c ]]];
165+ evae [a_ + b_ ]:= applyOpWithAbsoluteError [Plus ,evae [a ],evae [b ]];
147166evae [a_ + b__ ]:= (Message [IEEEEvaluateWithAbsoluteError ::badass ]; $Failed );
148- evae [a_ * b_ ]:= applyOp [Times ,evae [a ],evae [b ]];
167+ evae [a_ * b_ ]:= applyOpWithAbsoluteError [Times ,evae [a ],evae [b ]];
149168evae [a_ * b__ ]:= (Message [IEEEEvaluateWithAbsoluteError ::badass ]; $Failed );
150- evae [a_ / b_ ]:= applyOp [Divide ,evae [a ],evae [b ]];
169+ evae [a_ / b_ ]:= applyOpWithAbsoluteError [Divide ,evae [a ],evae [b ]];
151170(* Negation is exact. *)
152171evae [- a_ ]:= - evae [a ];
153172(* Squaring an interval is not the same as multiplying two identical intervals. *)
154- evae [a_ ^ 2 ]:= applyOp [# ^ 2 & ,evae [a ]];
155- evae [a_ ^ 3 ]:= applyOp [cube ,evae [a ]];
156- evae [a_ ^ 4 ]:= applyOp [# ^ 2 & ,applyOp [# ^ 2 & ,evae [a ]]];
173+ evae [a_ ^ 2 ]:= applyOpWithAbsoluteError [# ^ 2 & ,evae [a ]];
174+ evae [a_ ^ 3 ]:= applyOpWithAbsoluteError [cube ,evae [a ]];
175+ evae [a_ ^ 4 ]:= applyOpWithAbsoluteError [# ^ 2 & ,applyOpWithAbsoluteError [# ^ 2 & ,evae [a ]]];
157176evae [a_ ? NumberQ ]:= Block [{cra = CorrectlyRound [a ]},evae [Interval [{cra ,cra }]]];
158177evae [{v_ Interval ,\[Delta ]_ Interval }]:= {v ,\[Delta ]};
159178evae [a_ Interval ]:= {a ,Interval [{0 ,0 }]};
160179evae [a_ ? ValueQ ]:= evae [Evaluate [a ]];
161180evae [a_ ]:= a ;
162181evae [x ]];
163182IEEEEvaluateWithAbsoluteError [_ , args__ ]:=
164- (Message [IEEEEvaluate ::argnum , Length [{args }] + 1 ]; $Failed );
183+ (Message [IEEEEvaluateWithAbsoluteError ::argnum , Length [{args }] + 1 ]; $Failed );
184+
185+
186+ (* ::Text:: *)
187+ (*The relative error bound on an IEEE computation.*)
188+
189+
190+ relativeErrorBound := If [
191+ RoundingMode []== NearestTiesToEven ,
192+ FromRepresentation [Representation [1 ]+ 1 / 2 ]- 1 ,
193+ FromRepresentation [Representation [1 ]+ 1 ]- 1 ];
194+
195+
196+ (* ::Text:: *)
197+ (*Returns the relative error bound for a computation. The relative error is "small", i.e., close to 2^-53, not to (1+2^-53), and is an interval. The arguments are pairs of {value interval, relative error interval}.*)
198+
199+
200+ applyOpWithRelativeError [square ,{va_ ,\[Delta ]a_ }]:= Block [
201+ {h ,r ,\[Delta ]r },
202+ r = va ^ 2 ;
203+ r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
204+ h = relativeErrorBound ;
205+ \[Delta ]r = (1 + \[Delta ]a )^ 2 Interval [{1 - h ,1 + h }]- 1 ;
206+ {r ,\[Delta ]r }];
207+ applyOpWithRelativeError [cube ,{va_ ,\[Delta ]a_ }]:= Block [
208+ {h ,r ,\[Delta ]r },
209+ r = va ^ 3 ;
210+ h = relativeErrorBound ;
211+ \[Delta ]r = (1 + \[Delta ]a )^ 3 Interval [{1 - h ,1 + h }] Interval [{1 - h ,1 + h }]- 1 ;
212+ {r ,\[Delta ]r }];
213+ applyOpWithRelativeError [fourth ,{va_ ,\[Delta ]a_ }]:= Block [
214+ {h ,r ,\[Delta ]r },
215+ r = va ^ 4 ;
216+ h = relativeErrorBound ;
217+ \[Delta ]r = (1 + \[Delta ]a )^ 4 Interval [{1 - h ,1 + h }]Interval [{1 - h ,1 + h }]- 1 ;
218+ {r ,\[Delta ]r }];
219+ applyOpWithRelativeError [Plus ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ }]:= Block [
220+ {corners ,err ,h ,\[Theta ]2 ,\[Theta ]\[Prime ]2 ,r ,\[Delta ]r },
221+ r = va + vb ;
222+ r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
223+ If [
224+ Min [r ]< 0 && Max [r ]> 0 ,
225+ \[Delta ]r = Interval [{- \[Infinity ],+ \[Infinity ]}],
226+ h = relativeErrorBound ;
227+ \[Theta ]2 = (1 + \[Delta ]a )Interval [{1 - h ,1 + h }]- 1 ;
228+ \[Theta ]\[Prime ]2 = (1 + \[Delta ]b )Interval [{1 - h ,1 + h }]- 1 ;
229+ (* At the origin the function err has an indeterminate value
230+ bounded by \[Delta]wa and \[Delta]wb *)
231+ err [0 ,\[Delta ]wa_ ,0 ,\[Delta ]wb_ ]:= Interval [Min [{\[Delta ]wa ,\[Delta ]wb }],Max [{\[Delta ]wa ,\[Delta ]wb }]];
232+ err [wa_ ,\[Delta ]wa_ ,wb_ ,\[Delta ]wb_ ]:= (wa \[Delta ]wa + wb \[Delta ]wb )/ (wa + wb );
233+ (* The function err is monotonic, and therefore its extrema
234+ are reached at the corners of its domain. *)
235+ corners = Outer [
236+ err ,
237+ {Min [va ],Max [va ]},
238+ {Min [\[Theta ]2 ],Max [\[Theta ]2 ]},
239+ {Min [vb ],Max [vb ]},
240+ {Min [\[Theta ]\[Prime ]2 ],Max [\[Theta ]\[Prime ]2 ]}];
241+ \[Delta ]r = Interval [{Min [corners ],Max [corners ]}]];
242+ {r ,\[Delta ]r }];
243+ applyOpWithRelativeError [Times ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ }]:= Block [
244+ {h ,r ,\[Delta ]r },
245+ r = va vb ;
246+ r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
247+ h = relativeErrorBound ;
248+ \[Delta ]r = (1 + \[Delta ]a )(1 + \[Delta ]b )Interval [{1 - h ,1 + h }]- 1 ;
249+ {r ,\[Delta ]r }];
250+ applyOpWithRelativeError [Divide ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ }]:= Block [
251+ {h ,r ,\[Delta ]r },
252+ r = va / vb ;
253+ r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
254+ h = relativeErrorBound ;
255+ \[Delta ]r = Interval [{1 - h ,1 + h }](1 + \[Delta ]a )/ (1 + \[Delta ]b )- 1 ;
256+ {r ,\[Delta ]r }];
257+ applyOpWithRelativeError [fma ,{va_ ,\[Delta ]a_ },{vb_ ,\[Delta ]b_ },{vc_ ,\[Delta ]c_ }]:= Block [
258+ {h ,r ,\[Delta ]r },
259+ r = va vb + vc ;
260+ r = Interval [{CorrectlyRound [Min [r ]],CorrectlyRound [Max [r ]]}];
261+ h = relativeErrorBound ;
262+ \[Theta ]3 = (1 + \[Delta ]a )(1 + \[Delta ]b )Interval [{1 - h ,1 + h }]- 1 ;
263+ \[Theta ]2 = (1 + \[Delta ]c )Interval [{1 - h ,1 + h }]- 1 ;
264+ \[Delta ]r = Interval [{Min [\[Theta ]3 ,\[Theta ]2 ],Max [\[Theta ]3 ,\[Theta ]2 ]}];
265+ {r ,\[Delta ]r }];
266+
267+
268+ SetAttributes [IEEEEvaluateWithRelativeError ,HoldAll ];
269+ Options [IEEEEvaluateWithRelativeError ]= {UseFMA -> True };
270+ IEEEEvaluateWithRelativeError [x_ ,OptionsPattern []]:=
271+ Block [
272+ {Plus ,Times ,evre ,usefma = OptionValue [UseFMA ]},
273+ SetAttributes [evre ,HoldAll ];
274+ evre [a_ * b_ + c_ ]:= If [
275+ usefma ,
276+ applyOpWithRelativeError [fma ,evre [a ],evre [b ],evre [c ]],
277+ applyOpWithRelativeError [Plus ,applyOpWithRelativeError [Times ,evre [a ],evre [b ]],evre [c ]]];
278+ evre [a_ + b_ ]:= applyOpWithRelativeError [Plus ,evre [a ],evre [b ]];
279+ evre [a_ + b__ ]:= (Message [IEEEEvaluateWithRelativeError ::badass ]; $Failed );
280+ evre [a_ * b_ ]:= applyOpWithRelativeError [Times ,evre [a ],evre [b ]];
281+ evre [a_ * b__ ]:= (Message [IEEEEvaluateWithRelativeError ::badass ]; $Failed );
282+ evre [a_ / b_ ]:= applyOpWithRelativeError [Divide ,evre [a ],evre [b ]];
283+ (* Negation is exact. *)
284+ evre [- a_ ]:= - evre [a ];
285+ (* Squaring an interval is not the same as multiplying two identical intervals. *)
286+ evre [a_ ^ 2 ]:= applyOpWithRelativeError [square ,evre [a ]];
287+ evre [a_ ^ 3 ]:= applyOpWithRelativeError [cube ,evre [a ]];
288+ evre [a_ ^ 4 ]:= applyOpWithRelativeError [fourth ,evre [a ]];
289+ evre [a_ ? NumberQ ]:= Block [{cra = CorrectlyRound [a ]},evre [Interval [{cra ,cra }]]];
290+ evre [{v_ Interval ,\[Delta ]_ Interval }]:= {v ,\[Delta ]};
291+ evre [a_ Interval ]:= {a ,Interval [{0 ,0 }]};
292+ evre [a_ ? ValueQ ]:= evre [Evaluate [a ]];
293+ evre [a_ ]:= a ;
294+ evre [x ]];
295+ IEEEEvaluateWithRelativeError [_ , args__ ]:=
296+ (Message [IEEEEvaluateWithRelativeError ::argnum , Length [{args }] + 1 ]; $Failed );
165297
166298
167299End []
0 commit comments