Sorry that it took so long.
I hope that I did the changes to your liking.
The last 3 commits of
https://github.com/fricas/fricas/compare/master...hemmecke:fricas:wip/improve-stream-division
are the changes that I added on top of what I had in July.
All 4 commits in one patch is attached.
Ralf
On 7/31/24 23:15, Waldek Hebisch wrote:
On Wed, Jul 31, 2024 at 05:55:30PM +0200, Ralf Hemmecke wrote:
On 7/31/24 16:49, Waldek Hebisch wrote:>> ((*$sTO)(rs1, s3) - rs2)$sTO
gives the same "wrong" result.
It should be
((*$sTO)(s3, rs1) - rs2)$sTO
and this works without the patch.
I know, but why do you claim that it must be the way you say in a
non-commutative ring? I do not think that it is only my opinion that this is
ambiguous.
Well, people use various weird conventions. But mainstream practice
to is chose conventions in a way that make computations simpler.
Normal rule is s*(x/y) = (s*x)/y. Another rule is (x/y)*y = x.
Both rules would be broken if you define division differently.
Or to put is differenly, if you want simple rules do not move things
to the other side on the whim.
I have nothing against introducing right-/left- division, but / should be
reserved for the case when x * y^(-1) = y^(-1) * x otherwise I would
leave
it undefined.
Well, division solves equation x = d*y, so we have d = x * y^(-1).
How can you say so in a non-commutative ring? I guess some people claim that
it solves x=y*d.
This definition is necessary to have sane interaction with multiplication.
To do it on the other side is a different operation. If you look
at definitions in the algebra, FreeDivisionAlgebra, Group and
XPolynomialRing are quite explicit and define x/y as x * y^(-1).
In most other cases order does not matter as corresponding
multiplication is commutative. Possibly the only unclear definition
is the one in StreamTaylorSeriesOperations.
OK, you convinced me with our definitions in FriCAS. I still see danger that
our formatters transform x/y to \frac{x}{y} also for a non-commutative
domain (I haven't actually checked) and that would be certainly wrong.
If they do, we would need to fix this. But it is likely that formatters
do not produce 'x/y'.
And yes, the docstring of / in StreamTaylorSeriesOperations must definitely
reflect compatibility with FreeDivisionAlgebra and all other
(non-commutative) places.
Yes.
In other words we should have
if A has CommutativeRing then
"exquo" : (ST A,ST A) -> Union(ST A,"failed")
"/" : (ST A,ST A) -> ST A
recip : ST A -> UN
Maybe we can a bit weaker for recip, but I would still require that
x * recip(x) = recip(x) * x. Otherwise recip should fail.
We have:
recip : % -> Union(%,"failed")
++ recip(a) returns an element, which is both a left and a right
++ inverse of \spad{a},
++ or \spad{"failed"} if such an element doesn't exist or cannot
++ be determined (see unitsKnown).
So yes, x * recip(x) = recip(x) * x = e where e is the unit element
for multiplication.
Yes, but then we should also update the docstring in
StreamTaylorSeriesOperations
recip : ST A -> UN
++ recip(a) returns the power series reciprocal of \spad{a}, or
++ "failed" if not possible.
and make this explicit.
I do not like that users (even though pretty
intuitive) have to guess that the word "reciprocal" actually refers to the
docstring of recip: % -> Union(%,"failed") (somewhere else) as you cited
above, i.e. that is a left+right inverse.
OK. However, operations normally should be defined in categories
and users should look for definitions in categories. In some cases
like StreamTaylorSeriesOperations we want to keep functionality
as a package. Technically we can not reuse definition from categories
in such case, but users should understand that the operation are
the same or at least as close as possible do definitions in
categories. The definition above is from MagmaWithUnit, and is
inherited by Ring. StreamTaylorSeriesOperations in general
description says that it 'implements Taylor series arithmetic'.
I am not sure how explict we should be in the description, but
Taylor series _domains_ are declared to be a Ring, and "arithmetic"
operations defined in StreamTaylorSeriesOperations should satisfy
Ring properties. So, any deviations from definition in Ring either
are bugs or rewording wiht the same meaning or maybe a specialization
(ring of power series have some special properties compared to
deneral rings). From my point of view clearest docsting would say
"see corresponding operation in Ring".
To put this differently, at programming language level programmers
can implement quite a different property for operations with the same
name. But in many cases including StreamTaylorSeriesOperations
users should regard operations as "the same" as operations defined
in relevant category, that is Ring in this case.
Should we introduce also a function "\" into StreamTaylorSeriesOperations?
If you wish. For me it is important to avoid breaking things
that used to work, so '/' is important. To say the truth, 'exquo'
is more important as is has uses in other places and is not the
the same as multiplication by reciprocal (which in case of 'exquo'
may be non-existent). Also, if you want tha add new operation
than something like 'leftExquo' probably would fit.
--
You received this message because you are subscribed to the Google Groups "FriCAS -
computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/1f15e0c5-ac4b-4941-8320-be2b6ce39876%40hemmecke.org.
From fdbb9361bab20a9708c351b6096f19d561300b8e Mon Sep 17 00:00:00 2001
From: Ralf Hemmecke <[email protected]>
Date: Sat, 27 Jul 2024 00:06:28 +0200
Subject: division of streams with less memory
---
src/algebra/sttaylor.spad | 162 +++++++++++++++++++++++++++-----------
1 file changed, 117 insertions(+), 45 deletions(-)
diff --git a/src/algebra/sttaylor.spad b/src/algebra/sttaylor.spad
index b805ed07..3880ebaa 100644
--- a/src/algebra/sttaylor.spad
+++ b/src/algebra/sttaylor.spad
@@ -12,7 +12,8 @@
++ References:
++ Description:
++ StreamTaylorSeriesOperations implements Taylor series arithmetic,
-++ where a Taylor series is represented by a stream of its coefficients.
+++ where a Taylor series is represented by a stream of its coefficients,
+++ see corresponding operations in the category Ring.
StreamTaylorSeriesOperations(A) : Exports == Implementation where
A : Ring
RN ==> Fraction Integer
@@ -194,51 +195,122 @@ StreamTaylorSeriesOperations(A) : Exports == Implementation where
zero? s => zro()
map(z +-> z*s, x)
- iDiv : (ST A, ST A, A) -> ST A
- iDiv(x, rst_y, ry0) == delay
- empty? x => empty()
- c0 := frst x * ry0
- concat(c0, iDiv(rst x - c0*rst_y, rst_y, ry0))
-
- i_div1 : (ST A, ST A, A, A) -> ST A
- i_div1(rst_x, rst_y, ry0, c0) == delay
- empty?(rst_y) => map(z +-> z*ry0, rst_x)
- iDiv(rst_x - c0*rst_y, rst_y, ry0)
-
- x exquo y ==
- y0 : A
- for n in 1.. repeat
- n > 1000 => return "failed"
- empty? y => return "failed"
- empty? x => return empty()
- (y0 := frst y) = 0 =>
- frst x = 0 => (x := rst x; y := rst y)
- return "failed"
- -- first entry in y is non-zero
- break
- (ry0u := recip y0) case "failed" => "failed"
- ry0 := ry0u@A
- x0 := frst x
- c0 := x0*ry0
- concat(c0, i_div1(rst x, rst y, ry0, c0))
-
- (x : ST A) / (y : ST A) == delay
- empty? y => error "/: division by zero"
- empty? x => empty()
- (ry0u := recip frst y) case "failed" =>
- error "/: second argument is not invertible"
- ry0 := ry0u@A
- x0 := frst x
- c0 := x0*ry0
- concat(c0, i_div1(rst x, rst y, ry0, c0))
-
- recip x ==
- empty? x => "failed"
- rh1 := recip frst x
- rh1 case "failed" => "failed"
- rh := rh1@A
+-------------------------------------------------------------------
+
+ -- We use
+ -- https://en.wikipedia.org/wiki/Formal_power_series#Multiplicative_inverse.
+ -- For the actual computation formula (non-commutative case) see
+ -- explanation for / below with f=1.
+ -- recip computes the zeroth coefficient and computes arguments
+ -- for restRecip to compute the remaining coefficients.
+ -- restRecip avoids zero multiplication if the tail of the yet computed
+ -- partial stream ends in zeros.
+ -- Assume that x = (a0,a1,a2,...)
+ -- Arguments for restRecip:
+ -- a = (a1, a2, a3,...)
+ -- u = - 1/a0
+ -- n = in this invocation we compute (cn, c(n+1), ...)
+ -- k = the resulting coefficients ck, c(k+1), ..., c(n-1) are equal to 0
+ -- ra = (ak,a(k+1),a(k+2),...) substream of x starting at k
+ -- revc = [c(k-1),c(k-2),...,c(0)] the first k coefficients of the result
+ -- in reverse order, c(k-1)~=0.
+
+ -- local
+ revSum(ra : ST A, rvc : List A) : A ==
+ cc : A := 0
+ for c in rvc while not empty? ra repeat
+ cc := cc + c * frst(ra)
+ ra := rst ra
+ cc
+
+ -- local
+ restRecip(a : ST A, u : A, n : I, k : I, ra : ST A, revc : List A): ST A ==
delay
- concat(rh, iDiv(- rh * rst x, rst x, rh))
+ --assert(k=#revc)
+ --assert(for i in k .. n-1 while not empty? a repeat a := rst a; a=ra)
+ --assert(#revc>0)
+ empty? ra => empty()
+ rvc : List A := revc -- create a local variable inside 'delay'
+ cn : A := u * revSum(ra, rvc) -- u = -1/a(0)
+ zero? cn => concat(cn, restRecip(a, u, n+1, k, rst ra, rvc))
+ for i in k..n-1 repeat rvc := cons(0, rvc) -- fill rvc with zeros
+ concat(cn, restRecip(a, u, n+1, n+1, a, cons(cn, rvc)))
+
+ recip(x : ST A) : UN ==
+ empty? x => "failed"
+ ua : Union(A, "failed") := recip frst x
+ ua case "failed" => "failed"
+ c0 : A := ua @ A
+ concat(c0, restRecip(rst x, -c0, 1, 1, rst x, [c0])) :: UN
+
+-------------------------------------------------------------------
+
+ -- We use
+ -- https://en.wikipedia.org/wiki/Formal_power_series#Division.
+ -- In fact, the Wikipedia formula is for the commutative case.
+ -- We allow / also for the non-commutative case to mean "right-division",
+ -- i.e. if h=f/g then h*g=f. Thus,
+ -- $f=\sum_{n=0}^\infty b_n x^n$,
+ -- $g=\sum_{n=0}^\infty a_n x^n$,
+ -- $h=\sum_{n=0}^\infty c_n x^n$ as in the Wikipedia article, we get:
+ -- $\sum_{n=0}^\infty c_n x^n \sum_{n=0}^\infty a_n x^n
+ -- = \sum_{n-0}^\infty b_n x^n$ and by equating coefficients we get
+ -- for any $n\ge0$ that $b_n = \sum_{i=0}^n c_{n-i} a_i$.
+ -- Therefore $c_0 = b_0 / a_0 = b_0 a_0^{-1}$ and
+ -- $c_n = (b_n - \sum_{i=1}^n c_{n-i} a_i) / a_0$.
+
+ -- / computes the zeroth coefficient and computes arguments
+ -- for restDiv to compute the remaining coefficients.
+ -- restDiv avoids zero multiplication if the tail of the yet computed
+ -- partial stream ends in zeros.
+ -- Assume that x = (b0,b1,b2,...), y = (a0,a1,a2,...).
+ -- Arguments for restDiv:
+ -- a = (a1, a2, a3,...)
+ -- u = - 1/a0
+ -- n = in this invocation we compute (cn, c(n+1), ...)
+ -- k = the resulting coefficients ck, c(k+1), ..., c(n-1) are equal to 0
+ -- ra = (ak,a(k+1),a(k+2),...) substream of x starting at k
+ -- b = (bn, b(n+1), b(n+2),...)
+ -- revc = [c(k-1),c(k-2),...,c(0)] the first k coefficients of the result
+ -- in reverse order, c(k-1)~=0.
+
+ -- local
+ restDiv(a : ST A, u : A, n : I, k : I, ra : ST A, b : ST A,_
+ revc : List A) : ST A == delay
+ --assert(k=#revc)
+ --assert(for i in k .. n-1 while not empty? a repeat a := rst a; a=ra)
+ --assert(#revc>0)
+ empty? b => restRecip(a, u, n, k, ra, revc)
+ rvc : List A := revc -- create a local variable inside 'delay'
+ cn : A := (revSum(ra, rvc) - frst b) * u -- u = -1/a(0)
+ zero? cn => concat(cn, restDiv(a, u, n+1, k, rst ra, rst b, rvc))
+ for i in k..n-1 repeat rvc := cons(0, rvc) -- fill rvc with zeros
+ concat(cn, restDiv(a, u, n+1, n+1, a, rst b, cons(cn, rvc)))
+
+ ((x : ST A) exquo (y : ST A)) : UN ==
+ for n in 1..1000 repeat
+ empty? y => return "failed"
+ empty? x => return empty()
+ not zero? frst y => break
+ not zero? frst x => return "failed"
+ x := rst x
+ y := rst y
+ ua : Union(A, "failed") := recip frst y
+ ua case "failed" => "failed"
+ u : A := ua @ A
+ c0 := frst(x) * u
+ concat(c0, restDiv(rst y, -u, 1, 1, rst y, rst x, [c0]))
+
+ ((x : ST A) / (y : ST A)) : ST A == delay
+ empty? y => error "/: division by zero"
+ empty? x => empty()
+ ua : Union(A, "failed") := recip frst y
+ ua case "failed" => error "/: second argument is not invertible"
+ u : A := ua @ A
+ c0 := frst(x) * u
+ concat(c0, restDiv(rst y, -u, 1, 1, rst y, rst x, [c0]))
+
+-------------------------------------------------------------------
--% coefficients
--
2.34.1
)clear complete
ZZ ==> Integer
QQ ==> Fraction ZZ
S ==> Stream QQ
STT A ==> StreamTaylorSeriesOperations A
INV x ==> ((recip x)$STT(QQ)) :: S
M2 ==> SquareMatrix(2, QQ)
S2 ==> Stream M2
mul2 ==> (*$STT(M2))
div2 ==> (/$STT(M2))
sub2 ==> ((-)$STT(M2))
CHK(x,y) ==> (_
for i in 1..2000 repeat if x.i ~= y.i then (print(i);break);_
print("OK"))
)set stream calc 20
-------------------------------------------------------------------
st1: S := construct [3,0,0,0,2,0,0,0,3,0,4,0,0,5]
st2: S := [i for i in 1..100]
st3: S := [i for i in 1..]
st4: S := [i+1 for i in 1..]
st5: S := [99-i for i in 1..100]
st6: S := [1+i^2 for i in 1..]
)time on
u1 := INV st1; u1.50000; -- 13.5 s
u2 := INV st2; u2.100000; -- 7.2 s
u3 := INV st3; u3.100000; -- 0.3 s
u4 := INV st4; u4.5000; -- 25.6 s
u5 := INV st5; u5.3000; -- 20.0 s
u6 := INV st6; u6.4000; -- 20.6 s
)co sttaylor.spad
v1 := INV st1; v1.50000; -- 13.6 s
v2 := INV st2; v2.100000; -- 6.3 s
v3 := INV st3; v3.100000; -- 0.1
v4 := INV st4; v4.5000; -- 25.9 s
v5 := INV st5; v5.3000; -- 20.4 s
v6 := INV st6; v6.4000; -- 18.4 s
CHK(u1,v1)
CHK(u2,v2)
CHK(u3,v3)
CHK(u4,v4)
CHK(u5,v5)
CHK(u6,v6)
-------------------------------------------------------------------
zz := matrix([[0, 0], [0, 0]])$M2 -- zero
uu := matrix([[1, 0], [0, 1]])$M2 -- unit
ss := matrix([[0,-1], [1, 0]])$M2
tt := matrix([[1, 1], [0, 1]])$M2
m1 := matrix([[0, 1], [0, 0]])$M2
m2 := matrix([[0, 0], [1, 0]])$M2
stx := construct([uu, uu])$S2
str := div2(div2(stx,stx),stx)
one := div2(str,str) -- infinite stream representing one
zero := sub2(str,str) -- infinite stream representing zero
s1 := construct([3::M2, 0, m1, m2, 1])$S2
rs1 := recip(s1)$STT(M2);
CHK(mul2(s1,rs1), one)
CHK(mul2(rs1,s1), one)
s2 := construct([1::M2, 1::M2 + m1, m2, m2])$S2
rs2 := recip(s2)$STT(M2)
CHK(mul2(s2,rs2), one)
CHK(mul2(rs2,s2), one)
s3 := mul2(rs2, s1)
d3 := div2(s3, s1)
CHK(sub2(d3, rs2), zero)
