Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 21 additions & 4 deletions lagoon/desk/lib/lagoon.hoon
Original file line number Diff line number Diff line change
Expand Up @@ -903,10 +903,20 @@
?: ?=(%fixp kind.meta.a)
(scalar-to-ray meta.a (fixp-fdp ;;([a=@ b=@] tail.meta.a) (ravel a) (ravel b)))
(cumsum (mul a b))
:: +dotc: Hermitian (conjugated) dot product, Sum conj(a_i)*b_i. This is
:: +dotc: [$ray $ray] -> $ray
::
:: Returns the Hermitian (conjugated) dot product Sum conj(a_i)*b_i. This is
:: the inner product whose norm <a,a> = Sum |a_i|^2 is real and nonnegative
:: (what Saloon's eig/QR want). For real kinds conj is identity, so dotc
:: coincides with dot.
:: (what Saloon's eig/QR want). On real kinds conj is identity, so dotc
:: coincides with +dot.
:: Examples
:: > =a (en-ray:la [[~[1 2] 6 %cplx ~] ~[~[0x4000.0000.3f80.0000 0x4080.0000.4040.0000]]]) :: [1+2i 3+4i]
:: > =b (en-ray:la [[~[1 2] 6 %cplx ~] ~[~[0x40c0.0000.40a0.0000 0x4100.0000.40e0.0000]]]) :: [5+6i 7+8i]
:: > (get-item:la (dotc:la a b) ~[0 0])
:: 0xc100.0000.428c.0000 :: Sum conj(x)*y = 70-8i
:: > (get-item:la (dot:la a b) ~[0 0])
:: 0x4288.0000.c190.0000 :: Sum x*y = -18+68i
:: Source
++ dotc
~/ %dotc
|= [a=ray b=ray]
Expand Down Expand Up @@ -1011,7 +1021,14 @@
^- ray
(el-wise-op a (trans-scalar bloq.meta.a kind.meta.a %abs))
::
:: +conj: elementwise complex conjugate (identity on real kinds).
:: +conj: $ray -> $ray
::
:: Returns the elementwise complex conjugate of a ray. Identity on real
:: kinds, which is what makes +dotc reduce to +dot on reals.
:: Examples
:: > (get-item:la (conj:la (fill:la [~[1 1] 6 %cplx ~] 0x4000.0000.3f80.0000)) ~[0 0])
:: 0xc000.0000.3f80.0000 :: conj(1+2i)=1-2i
:: Source
++ conj
~/ %conj
|= a=ray
Expand Down
251 changes: 242 additions & 9 deletions libmath/desk/lib/complex.hoon
Original file line number Diff line number Diff line change
Expand Up @@ -22,32 +22,120 @@
::
|%
+$ rounding-mode ?(%n %u %d %z)
:: +cs: complex-single (@cs), two @rs (32-bit) components.
:: +cs: complex-single (@cs), two @rs (32-bit) components.
::
:: Door keyed on the rounding mode. In the examples below 1+2i packs to
:: 0x4000.0000.3f80.0000 and 3+4i to 0x4080.0000.4040.0000 (real in the low
:: 32 bits).
::
++ cs
|_ rnd=rounding-mode
:: component accessors / constructor (real low, imag high)
::
:: Components
::
:: +re: @cs -> @rs
::
:: Returns the real component (low 32 bits) of a packed @cs value.
:: Examples
:: > (~(re cs:complex %n) 0x4000.0000.3f80.0000) :: re(1+2i)
:: .1
:: Source
++ re |=(p=@ `@rs`(end [0 32] p))
:: +im: @cs -> @rs
::
:: Returns the imaginary component (high 32 bits) of a packed @cs value.
:: Examples
:: > (~(im cs:complex %n) 0x4000.0000.3f80.0000) :: im(1+2i)
:: .2
:: Source
++ im |=(p=@ `@rs`(rsh [0 32] p))
:: +pak: [@rs @rs] -> @cs
::
:: Packs real and imaginary @rs components into one @cs atom (real low).
:: Examples
:: > (~(pak cs:complex %n) .1 .2)
:: 0x4000.0000.3f80.0000
:: Source
++ pak |=([r=@rs i=@rs] ^-(@ (con `@`r (lsh [0 32] `@`i))))
:: component helpers
:: component helpers (internal)
++ fneg |=(x=@rs (~(sub rs rnd) .0 x))
++ fabs |=(x=@rs ?:((~(lth rs rnd) x .0) (~(sub rs rnd) .0 x) x))
:: constants
::
:: Constants
::
:: +zero: @cs
::
:: The additive identity 0+0i.
:: Examples
:: > ~(zero cs:complex %n)
:: 0x0
:: Source
++ zero ^-(@ 0)
:: +one: @cs
::
:: The multiplicative identity 1+0i.
:: Examples
:: > ~(one cs:complex %n)
:: 0x3f80.0000
:: Source
++ one (pak .1 .0)
:: arithmetic
::
:: Arithmetic
::
:: +add: [@cs @cs] -> @cs
::
:: Returns the sum of two complex values (component-wise).
:: Examples
:: > (~(add cs:complex %n) 0x4000.0000.3f80.0000 0x4080.0000.4040.0000)
:: 0x40c0.0000.4080.0000
:: Source
++ add |=([p=@ q=@] (pak (~(add rs rnd) (re p) (re q)) (~(add rs rnd) (im p) (im q))))
:: +sub: [@cs @cs] -> @cs
::
:: Returns the difference of two complex values (component-wise).
:: Examples
:: > (~(sub cs:complex %n) 0x4000.0000.3f80.0000 0x4080.0000.4040.0000)
:: 0xc000.0000.c000.0000
:: Source
++ sub |=([p=@ q=@] (pak (~(sub rs rnd) (re p) (re q)) (~(sub rs rnd) (im p) (im q))))
:: +neg: @cs -> @cs
::
:: Returns the additive inverse -z.
:: Examples
:: > (~(neg cs:complex %n) 0x4000.0000.3f80.0000) :: -(1+2i)
:: 0xc000.0000.bf80.0000
:: Source
++ neg |=(p=@ (pak (fneg (re p)) (fneg (im p))))
:: +conj: @cs -> @cs
::
:: Returns the complex conjugate (negate the imaginary component).
:: Examples
:: > (~(conj cs:complex %n) 0x4000.0000.3f80.0000) :: conj(1+2i)=1-2i
:: 0xc000.0000.3f80.0000
:: Source
++ conj |=(p=@ (pak (re p) (fneg (im p))))
:: +mul: [@cs @cs] -> @cs
::
:: Returns the complex product (ar+ai i)(br+bi i), rounded per component.
:: Examples
:: > (~(mul cs:complex %n) 0x4000.0000.3f80.0000 0x4080.0000.4040.0000)
:: 0x4120.0000.c0a0.0000 :: (1+2i)(3+4i)=-5+10i
:: Source
++ mul
|= [p=@ q=@]
=/ ar (re p) =/ ai (im p)
=/ br (re q) =/ bi (im q)
%+ pak
(~(sub rs rnd) (~(mul rs rnd) ar br) (~(mul rs rnd) ai bi))
(~(add rs rnd) (~(mul rs rnd) ar bi) (~(mul rs rnd) ai br))
:: +div: Smith's algorithm (scale by the larger denominator component).
:: +div: [@cs @cs] -> @cs
::
:: Returns the complex quotient via Smith's algorithm (scale by the larger
:: denominator component for numerical stability).
:: Examples
:: > (~(div cs:complex %n) 0x4000.0000.4000.0000 0x3f80.0000.3f80.0000)
:: 0x4000.0000 :: (2+2i)/(1+1i)=2
:: Source
++ div
|= [p=@ q=@]
=/ ar (re p) =/ ai (im p)
Expand All @@ -63,7 +151,14 @@
%+ pak
(~(div rs rnd) (~(add rs rnd) (~(mul rs rnd) ar r) ai) dn)
(~(div rs rnd) (~(sub rs rnd) (~(mul rs rnd) ai r) ar) dn)
:: +abs: modulus |z| as a real-valued complex [|z|, 0], via hypot.
:: +abs: @cs -> @cs
::
:: Returns the modulus |z| as a real-valued complex [|z|, 0], computed via
:: hypot (scale by the larger component to avoid overflow).
:: Examples
:: > (~(abs cs:complex %n) 0x4080.0000.4040.0000) :: |3+4i|=5
:: 0x40a0.0000
:: Source
++ abs
|= p=@
=/ xr (fabs (re p)) =/ xi (fabs (im p))
Expand All @@ -73,31 +168,142 @@
(pak (~(mul rs rnd) xr (~(sqt rs rnd) (~(add rs rnd) .1 (~(mul rs rnd) t t)))) .0)
=/ t (~(div rs rnd) xr xi)
(pak (~(mul rs rnd) xi (~(sqt rs rnd) (~(add rs rnd) .1 (~(mul rs rnd) t t)))) .0)
:: +equ/+neq: bit equality of both components (matches the %i754 convention).
::
:: Comparison (no total order on complex -- only equality)
::
:: +equ: [@cs @cs] -> ?
::
:: Loobean: bit-equality of both components (matches the %i754 convention).
:: Examples
:: > (~(equ cs:complex %n) 0x4000.0000.3f80.0000 0x4000.0000.3f80.0000)
:: %.y
:: > (~(equ cs:complex %n) 0x4000.0000.3f80.0000 0x4080.0000.4040.0000)
:: %.n
:: Source
++ equ |=([p=@ q=@] &(=((re p) (re q)) =((im p) (im q))))
:: +neq: [@cs @cs] -> ?
::
:: Loobean negation of +equ.
:: Examples
:: > (~(neq cs:complex %n) 0x4000.0000.3f80.0000 0x4080.0000.4040.0000)
:: %.y
:: Source
++ neq |=([p=@ q=@] =(| (equ p q)))
--
:: +cd: complex-double (@cd), two @rd (64-bit) components.
:: +cd: complex-double (@cd), two @rd (64-bit) components.
::
:: The @rd/@cd analogue of +cs over double-precision components. Here 1+2i
:: packs to 0x4000.0000.0000.0000.3ff0.0000.0000.0000 and 3+4i to
:: 0x4010.0000.0000.0000.4008.0000.0000.0000 (real in the low 64 bits).
::
++ cd
|_ rnd=rounding-mode
::
:: Components
::
:: +re: @cd -> @rd
::
:: Returns the real component (low 64 bits) of a packed @cd value.
:: Examples
:: > (~(re cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000) :: re(1+2i)
:: .~1
:: Source
++ re |=(p=@ `@rd`(end [0 64] p))
:: +im: @cd -> @rd
::
:: Returns the imaginary component (high 64 bits) of a packed @cd value.
:: Examples
:: > (~(im cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000) :: im(1+2i)
:: .~2
:: Source
++ im |=(p=@ `@rd`(rsh [0 64] p))
:: +pak: [@rd @rd] -> @cd
::
:: Packs real and imaginary @rd components into one @cd atom (real low).
:: Examples
:: > (~(pak cd:complex %n) .~1 .~2)
:: 0x4000.0000.0000.0000.3ff0.0000.0000.0000
:: Source
++ pak |=([r=@rd i=@rd] ^-(@ (con `@`r (lsh [0 64] `@`i))))
:: component helpers (internal)
++ fneg |=(x=@rd (~(sub rd rnd) .~0 x))
++ fabs |=(x=@rd ?:((~(lth rd rnd) x .~0) (~(sub rd rnd) .~0 x) x))
::
:: Constants
::
:: +zero: @cd
::
:: The additive identity 0+0i.
:: Examples
:: > ~(zero cd:complex %n)
:: 0x0
:: Source
++ zero ^-(@ 0)
:: +one: @cd
::
:: The multiplicative identity 1+0i.
:: Examples
:: > ~(one cd:complex %n)
:: 0x3ff0.0000.0000.0000
:: Source
++ one (pak .~1 .~0)
::
:: Arithmetic
::
:: +add: [@cd @cd] -> @cd
::
:: Returns the sum of two complex values (component-wise).
:: Examples
:: > (~(add cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4010.0000.0000.0000.4008.0000.0000.0000)
:: 0x4018.0000.0000.0000.4010.0000.0000.0000
:: Source
++ add |=([p=@ q=@] (pak (~(add rd rnd) (re p) (re q)) (~(add rd rnd) (im p) (im q))))
:: +sub: [@cd @cd] -> @cd
::
:: Returns the difference of two complex values (component-wise).
:: Examples
:: > (~(sub cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4010.0000.0000.0000.4008.0000.0000.0000)
:: 0xc000.0000.0000.0000.c000.0000.0000.0000
:: Source
++ sub |=([p=@ q=@] (pak (~(sub rd rnd) (re p) (re q)) (~(sub rd rnd) (im p) (im q))))
:: +neg: @cd -> @cd
::
:: Returns the additive inverse -z.
:: Examples
:: > (~(neg cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000) :: -(1+2i)
:: 0xc000.0000.0000.0000.bff0.0000.0000.0000
:: Source
++ neg |=(p=@ (pak (fneg (re p)) (fneg (im p))))
:: +conj: @cd -> @cd
::
:: Returns the complex conjugate (negate the imaginary component).
:: Examples
:: > (~(conj cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000) :: conj(1+2i)=1-2i
:: 0xc000.0000.0000.0000.3ff0.0000.0000.0000
:: Source
++ conj |=(p=@ (pak (re p) (fneg (im p))))
:: +mul: [@cd @cd] -> @cd
::
:: Returns the complex product (ar+ai i)(br+bi i), rounded per component.
:: Examples
:: > (~(mul cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4010.0000.0000.0000.4008.0000.0000.0000)
:: 0x4024.0000.0000.0000.c014.0000.0000.0000 :: (1+2i)(3+4i)=-5+10i
:: Source
++ mul
|= [p=@ q=@]
=/ ar (re p) =/ ai (im p)
=/ br (re q) =/ bi (im q)
%+ pak
(~(sub rd rnd) (~(mul rd rnd) ar br) (~(mul rd rnd) ai bi))
(~(add rd rnd) (~(mul rd rnd) ar bi) (~(mul rd rnd) ai br))
:: +div: [@cd @cd] -> @cd
::
:: Returns the complex quotient via Smith's algorithm (scale by the larger
:: denominator component for numerical stability).
:: Examples
:: > (~(div cd:complex %n) 0x4000.0000.0000.0000.4000.0000.0000.0000 0x3ff0.0000.0000.0000.3ff0.0000.0000.0000)
:: 0x4000.0000.0000.0000 :: (2+2i)/(1+1i)=2
:: Source
++ div
|= [p=@ q=@]
=/ ar (re p) =/ ai (im p)
Expand All @@ -113,6 +319,14 @@
%+ pak
(~(div rd rnd) (~(add rd rnd) (~(mul rd rnd) ar r) ai) dn)
(~(div rd rnd) (~(sub rd rnd) (~(mul rd rnd) ai r) ar) dn)
:: +abs: @cd -> @cd
::
:: Returns the modulus |z| as a real-valued complex [|z|, 0], computed via
:: hypot (scale by the larger component to avoid overflow).
:: Examples
:: > (~(abs cd:complex %n) 0x4010.0000.0000.0000.4008.0000.0000.0000) :: |3+4i|=5
:: 0x4014.0000.0000.0000
:: Source
++ abs
|= p=@
=/ xr (fabs (re p)) =/ xi (fabs (im p))
Expand All @@ -122,7 +336,26 @@
(pak (~(mul rd rnd) xr (~(sqt rd rnd) (~(add rd rnd) .~1 (~(mul rd rnd) t t)))) .~0)
=/ t (~(div rd rnd) xr xi)
(pak (~(mul rd rnd) xi (~(sqt rd rnd) (~(add rd rnd) .~1 (~(mul rd rnd) t t)))) .~0)
::
:: Comparison (no total order on complex -- only equality)
::
:: +equ: [@cd @cd] -> ?
::
:: Loobean: bit-equality of both components (matches the %i754 convention).
:: Examples
:: > (~(equ cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4000.0000.0000.0000.3ff0.0000.0000.0000)
:: %.y
:: > (~(equ cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4010.0000.0000.0000.4008.0000.0000.0000)
:: %.n
:: Source
++ equ |=([p=@ q=@] &(=((re p) (re q)) =((im p) (im q))))
:: +neq: [@cd @cd] -> ?
::
:: Loobean negation of +equ.
:: Examples
:: > (~(neq cd:complex %n) 0x4000.0000.0000.0000.3ff0.0000.0000.0000 0x4010.0000.0000.0000.4008.0000.0000.0000)
:: %.y
:: Source
++ neq |=([p=@ q=@] =(| (equ p q)))
--
--
Loading