Numeric error propagation: Difference between revisions

=={{header|Factor}}==

This version defines a new type, <code>imprecise</code>, and also defines a custom syntax similar to the syntax for complex numbers. It uses multi-methods to handle various combinations of scalar and compound values. The <code>multi-methods</code> vocabulary is described as experimental, but aside from a (probably intentional) clunky interface which requires the programmer to disambiguate generic word definitions, and the fact that multi-methods don't show up in the help browser properly, I had no issues with them. In some stress tests, they don't appear to suffer from speed issues.

{{works with|Factor|0.99 2019-10-06}}

<lang factor>USING: accessors arrays fry kernel locals math math.functions

multi-methods parser prettyprint prettyprint.custom sequences ;

RENAME: GENERIC: multi-methods => MM-GENERIC:

FROM: syntax => M: ;

IN: imprecise

TUPLE: imprecise

{ value float read-only }

{ sigma float read-only } ;

C: <imprecise> imprecise

: >imprecise< ( imprecise -- value sigma )

[ value>> ] [ sigma>> ] bi ;

! Define a custom syntax for imprecise numbers.

<< SYNTAX: I{ \ } [ first2 <imprecise> ] parse-literal ; >>

M: imprecise pprint-delims drop \ I{ \ } ;

M: imprecise >pprint-sequence >imprecise< 2array ;

M: imprecise pprint* pprint-object ;

<PRIVATE

! Error functions

: f+-i ( float imprecise quot -- imprecise )

[ >imprecise< ] dip dip <imprecise> ; inline

: i+-i ( imprecise1 imprecise2 quot -- imprecise )

'[ [ value>> ] bi@ @ ]

[ [ sigma>> sq ] bi@ + sqrt <imprecise> ] 2bi ; inline

: f*/i ( float imprecise quot -- imprecise )

[ >imprecise< overd ] dip [ * abs ] 2bi* <imprecise> ;

inline

:: i*/i ( a b quot -- imprecise )

a b [ >imprecise< ] bi@ :> ( vala siga valb sigb )

vala valb quot call :> val

val sq siga sq * vala sq /f sigb sq + valb sq /f sqrt :> sig

val sig <imprecise> ; inline

PRIVATE>

MM-GENERIC: ~+ ( obj1 obj2 -- imprecise ) foldable flushable

METHOD: ~+ { float imprecise } [ + ] f+-i ;

METHOD: ~+ { imprecise float } swap ~+ ;

METHOD: ~+ { imprecise imprecise } [ + ] i+-i ;

MM-GENERIC: ~- ( obj1 obj2 -- imprecise ) foldable flushable

METHOD: ~- { float imprecise } [ - ] f+-i ;

METHOD: ~- { imprecise float } swap [ swap - ] f+-i ;

METHOD: ~- { imprecise imprecise } [ - ] i+-i ;

MM-GENERIC: ~* ( obj1 obj2 -- imprecise ) foldable flushable

METHOD: ~* { float imprecise } [ * ] f*/i ;

METHOD: ~* { imprecise float } swap ~* ;

METHOD: ~* { imprecise imprecise } [ * ] i*/i ;

MM-GENERIC: ~/ ( obj1 obj2 -- imprecise ) foldable flushable

METHOD: ~/ { float imprecise } [ /f ] f*/i ;

METHOD: ~/ { imprecise float } swap [ swap /f ] f*/i ;

METHOD: ~/ { imprecise imprecise } [ /f ] i*/i ;

:: ~^ ( a x -- imprecise )

a >imprecise< :> ( vala siga )

vala x ^ >rect drop :> val

val x * siga vala /f * abs :> sig

val sig <imprecise> ; foldable flushable

<PRIVATE

: imprecise-demo ( -- )

I{ 100 1.1 } I{ 200 2.2 } ~- 2. ~^

I{ 50 1.2 } I{ 100 2.3 } ~- 2. ~^ ~+ 0.5 ~^ . ;

PRIVATE>

MAIN: imprecise-demo</lang>

{{out}}

<pre>

I{ 111.8033988749895 2.487167063146342 }

</pre>

The following version is more idiomatic for Factor, where the convention for simple numeric types is to define regular words that work on sequences. For examples of this, check out words like <code>v+n</code> in the <code>math.vectors</code> vocabulary or <code>q+</code> in <code>math.quaternions</code>. Since these words perform no dispatch, all three forms are defined: <code>v-n</code> for subtracting a scalar from a vector, <code>n-v</code> for subtracting a vector from a scalar, and <code>v-</code> for subtracting two vectors. This convention has been used in the following example.

{{works with|Factor|0.99 2019-10-06}}

<lang factor>USING: arrays kernel locals math math.functions math.vectors

prettyprint sequences sequences.extras ;

IN: uncertain

<PRIVATE

: ubi* ( x y p q -- u )

[ [ first2 ] bi@ swapd ] 2dip 2bi* 2array ; inline

: err+ ( x y -- z ) [ sq ] bi@ + sqrt ;

:: (u*) ( u1 u2 quot -- u )

u1 u2 [ first2 ] bi@ :> ( v1 s1 v2 s2 )

v1 v2 quot call :> val

s1 v1 /f sq s2 v2 /f sq + val sq * sqrt :> sig

val sig 2array ; inline

PRIVATE>

: u+n ( u n -- u ) 0 2array v+ ;

: n+u ( n u -- u ) swap u+n ;

: u-n ( u n -- u ) 0 2array v- ;

: n-u ( n u -- u ) [ 0 2array ] dip v- ;

: u+ ( u u -- u ) [ + ] [ err+ ] ubi* ;

: u- ( u u -- u ) [ - ] [ err+ ] ubi* ;

: u*n ( u n -- u ) dup 2array [ * ] [ * abs ] ubi* ;

: n*u ( n u -- u ) swap u*n ;

: u/n ( u n -- u ) dup 2array [ /f ] [ * abs ] ubi* ;

: n/u ( n u -- u ) [ dup 2array ] dip [ /f ] [ * abs ] ubi* ;

: u* ( u u -- u ) [ * ] (u*) ;

: u/ ( u u -- u ) [ /f ] (u*) ;

:: u^n ( u n -- u )

u first2 :> ( v1 s1 )

v1 n ^ >rect drop :> val

val n * s1 v1 /f * abs :> sig

val sig 2array ;

<PRIVATE

: uncertain-demo ( -- )

{ 100 1.1 } { 200 2.2 } u- 2.0 u^n

{ 50 1.2 } { 100 2.3 } u- 2.0 u^n u+ 0.5 u^n . ;

PRIVATE>

MAIN: uncertain-demo</lang>

{{out}}

<pre>

{ 111.8033988749895 2.487167063146342 }

</pre>