coercing of log(2)*1.0
Henryk Trappmann
while the general rule of coercing in binary operations seems to be
towards the most unprecise,
for example
RealField(100)(1)*2 == RealField(100)(2) or
RealField(100)(1)*RealField(50)(1) == RealField(50)(1)
this rule seems to be broken when working with log(2)
not RR(1)*log(2) == RR(log(2))
Then I read about the coercing rules and found that
RR.__call__(log(2)) is possible but
RR._coerce_(log(2)) raises a type error
is there reason to do so?
Michel
canonical
morphism SR-->RR. So perhaps coercion is not expected to work?
Michel
William Stein
<bo19...@googlemail.com> wrote:
>
> Hello,
>
> while the general rule of coercing in binary operations seems to be
> towards the most unprecise,
> for example
> RealField(100)(1)*2 == RealField(100)(2) or
> RealField(100)(1)*RealField(50)(1) == RealField(50)(1)
The general rule is not "toward the most unprecise". It's more subtle
than that. At a bare minimum there is never a canonical (automatic)
coercion from elements of R to elements of S unless that coercion
is defined (as a homomorphism) on all of R.
> this rule seems to be broken when working with log(2)
> not RR(1)*log(2) == RR(log(2))
>
> Then I read about the coercing rules and found that
> RR.__call__(log(2)) is possible but
> RR._coerce_(log(2)) raises a type error
>
> is there reason to do so?
As Michel pointed out, log(2) is an element of the symbolic
ring. There is no natural homomorphism from the symbolic
ring to RR. However, there is a natural homomorphism from
RR to the symbolic ring.
sage: parent(log(2))
Symbolic Ring
sage: parent(x)
Symbolic Ring
sage: SR._coerce_(RR(3))
3.00000000000000
sage: RR._coerce_(log(2))
TypeError
You can always force coercion using RR(log(2)), as you know.
-- William
>
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Henryk Trappmann
> However, there is a natural homomorphism from
> RR to the symbolic ring.
to RealField(2) is not valid, because it is no homomorphism at all.
For example let R2 = RealField(2), then
not R2(2.4+1.2)==R2(2.4)+R2(1.2)
Wouldnt it then be more consistent coerce RealFields to higher
precision? There really a homomorphism exists. Then there always would
be a (desirable) difference between rounding and coercing. Rounding
has to be explicit while coercing is automatic.
Of course at this stage I also have to point out that the so called
RealField is no field at all:
not R2(3)+R2(2)-R2(2) == R2(3)
William Stein
<bo19...@googlemail.com> wrote:
>
> On May 31, 3:59 pm, "William Stein" <wst...@gmail.com> wrote:
>> However, there is a natural homomorphism from
>> RR to the symbolic ring.
>
> Hm, if this is the precondition then the coercion of say RealField(52)
> to RealField(2) is not valid, because it is no homomorphism at all.
> For example let R2 = RealField(2), then
> not R2(2.4+1.2)==R2(2.4)+R2(1.2)
>
Of course RealField(53) is not even a
ring -- it fails almost every axiom of a ring.
> Wouldnt it then be more consistent coerce RealFields to higher
> precision?
Suppose you write down an expression involving various digits of precision,
and in order to evaluate it Sage makes a sequence of *automatic*
coercions and outputs the result. Do you want an answer that has
many more digits of precision than have any chance of being
meaningful? E.g., would you really find it natural for
R(2)(0.5) + R(1000)(pi)
to implicitly output an answer with 1000 bits of precision?
> There really a homomorphism exists. Then there always would
> be a (desirable) difference between rounding and coercing. Rounding
> has to be explicit while coercing is automatic.
>
> Of course at this stage I also have to point out that the so called
> RealField is no field at all:
> not R2(3)+R2(2)-R2(2) == R2(3)
Carl Witty
> On May 31, 3:59 pm, "William Stein" <wst...@gmail.com> wrote:
>
> > However, there is a natural homomorphism from
> > RR to the symbolic ring.
>
> Hm, if this is the precondition then the coercion of say RealField(52)
> to RealField(2) is not valid, because it is no homomorphism at all.
> For example let R2 = RealField(2), then
> not R2(2.4+1.2)==R2(2.4)+R2(1.2)
>
> Wouldnt it then be more consistent coerce RealFields to higher
> precision? There really a homomorphism exists. Then there always would
> be a (desirable) difference between rounding and coercing. Rounding
> has to be explicit while coercing is automatic.
RR(R2(2)+R2(3)) != RR(R2(2)) + RR(R2(3))
Floating-point numbers are horrid, of course; trying to use them in a
computer algebra system is just asking for trouble. But you really
can't do without them, so you need to try to make the least bad
choices.
IMHO, giving a+b the precision of the less-precise operand is better
than using the more-precise operand, because the latter has too much
chance of confusing people who don't understand floating point. For
instance, if 1.3+RealField(500)(pi) gave a 500-bit number, I'll bet a
lot of people would assume that this number matched 13/10+pi to almost
500 bits. (Giving a 53-bit number is confusing too, but it's less
dangerous, since you're less likely to believe you have the answer you
wanted.)
Of course, maybe there are other choices that are better than either
of these. We could throw away RealField and always use
RealIntervalField instead; but that's slower, has less special
functions implemented, and has counterintuitive semantics for
comparisons. We could do pseudo-interval arithmetic, and say that
1e6+R2(3) should be a 20-bit number, because we know about 20 bits of
the answer; but to be consistent, we should do similar psuedo-interval
arithmetic even if both operands are the same precision, and then RR
becomes less useful for people who do understand how floating point
works and want to take advantage of it.
Carl
Henryk Trappmann
> > precision?
>
> Suppose you write down an expression involving various digits of precision,
> and in order to evaluate it Sage makes a sequence of *automatic*
> coercions and outputs the result. Do you want an answer that has
> many more digits of precision than have any chance of being
> meaningful? E.g., would you really find it natural for
> R(2)(0.5) + R(1000)(pi)
> to implicitly output an answer with 1000 bits of precision?
everything else would be rounding and I dont want rounding to be
performed automatically. An innocent 1.0*R(1000)(pi) would kill my
laborious obtained precision for pi.
If you would carry over this principle to the integers then 1*RR(pi)
would be 3.
What is the difference? Its just an arbitrariness to coerce to less
precision in RealField.
However it injures two principles that are otherwise (*between*
integers, symbolic ring, real field) always satsified:
1. No automatic rounding.
2. coercing is a homomorphism
Henryk Trappmann
> Actually, there's no homomorphism either way;
> RR(R2(2)+R2(3)) != RR(R2(2)) + RR(R2(3))
homomorphism but perhaps this reasoning has no base.
> IMHO, giving a+b the precision of the less-precise operand is better
> than using the more-precise operand, because the latter has too much
> chance of confusing people who don't understand floating point. For
> instance, if 1.3+RealField(500)(pi) gave a 500-bit number, I'll bet a
> lot of people would assume that this number matched 13/10+pi to almost
> 500 bits.
nobody would assume that 0.3333333 coerced to 500 *decimal* digits
matches 1/3 to 500 digits? I anyway wondered why one can not specify a
base in RealField, or did I merely overlook the opportunity?
> Of course, maybe there are other choices that are better than either
> of these. We could throw away RealField and always use
> RealIntervalField instead; but that's slower, has less special
> functions implemented, and has counterintuitive semantics for
> comparisons. We could do pseudo-interval arithmetic, and say that
> 1e6+R2(3) should be a 20-bit number, because we know about 20 bits of
> the answer; but to be consistent, we should do similar psuedo-interval
> arithmetic even if both operands are the same precision,
> and then RR
> becomes less useful for people who do understand how floating point
> works and want to take advantage of it.
arbitrary to me to coerce down precision, while coercing up precision
with integers. Of course if you consider integer with precision
infinity (as indeed there are no rounding errors and it is an exact
ring) then this makes sense again.
But if so then I want to have something like SymbolicNumber which is
the subset of SymbolicRing that does not contain variables. And that
this SymbolicNumber is coerced automatically down when used with
RealField.
Jason Grout
I hope this is not superfluous, but there is a SymbolicConstant class,
which seems to be automatically used when dealing with numbers and not
variables:
sage: type(SR(3))
<class 'sage.calculus.calculus.SymbolicConstant'>
sage: type(1.0*SR(3))
<class 'sage.calculus.calculus.SymbolicConstant'>
Jason
Gary Furnish
> But if so then I want to have something like SymbolicNumber which is
> the subset of SymbolicRing that does not contain variables. And that
> this SymbolicNumber is coerced automatically down when used with
> RealField.
There are really, really severe issues with coercion out of the
SymbolicRing because of assumptions that are made by the coercion
system. If this was allowed, one could end up with situations where
the code would try to coerce complex numbers into RR, and this should
not be allowed. I'm well aware this is not ideal, it is just that my
hands are largely tied by the coercion system. I'd like to find a
better answer here (and this really becomes less of an issue with the
more advanced features in the new symbolics system I'm working on),
but this is a long term project to get coercion out of SR working.
William Stein
The parent of SymbolicConstant is SymbolicRing. The parent
is what matters for coercion.
sage: parent(SR(3))
Symbolic Ring
It would be conceivable that we could have a SymbolicConstant
ring that is a subring of the full SymbolicRing. I haven't thought
through at all how this would work, but it might address the
OP's confusion. Note that it would problematic in Sage, since
e.g., if x = var('x'), then x + 1 - x is the symbolic 1, so we have
that the sum of two elements of SR would be in SC, i.e., rings
wouldn't be closed under addition. This would take a lot of care
to actually do. I'm not saying any of this should be done.
William
Carl Witty
> On May 31, 10:55 pm, Carl Witty <cwi...@newtonlabs.com> wrote:
> > IMHO, giving a+b the precision of the less-precise operand is better
> > than using the more-precise operand, because the latter has too much
> > chance of confusing people who don't understand floating point. For
> > instance, if 1.3+RealField(500)(pi) gave a 500-bit number, I'll bet a
> > lot of people would assume that this number matched 13/10+pi to almost
> > 500 bits.
>
> Hm, yes, but this binary decimal conversion is another topic, I mean
> nobody would assume that 0.3333333 coerced to 500 *decimal* digits
> matches 1/3 to 500 digits? I anyway wondered why one can not specify a
> base in RealField, or did I merely overlook the opportunity?
Probably the basic reason is that RealField is based on MPFR, and MPFR
operates in binary.
Also, specifying a base, and doing it "right" (so that floating-point
results were correctly rounded to values of the specified number of
digits in that base) would probably be much slower in non-power-of-2
bases (at least I can't think how to do the basic operations as
quickly as in binary).
> > Of course, maybe there are other choices that are better than either
> > of these. We could throw away RealField and always use
> > RealIntervalField instead; but that's slower, has less special
> > functions implemented, and has counterintuitive semantics for
> > comparisons. We could do pseudo-interval arithmetic, and say that
> > 1e6+R2(3) should be a 20-bit number, because we know about 20 bits of
> > the answer; but to be consistent, we should do similar psuedo-interval
> > arithmetic even if both operands are the same precision,
>
> At least RR would then be a ring ;)
what the exact rules might be, but the rules I was envisioning would
not suffice to make RR a ring.
> > and then RR
> > becomes less useful for people who do understand how floating point
> > works and want to take advantage of it.
>
> Ya, I dont want to change floating point, it just seems somewhat
> arbitrary to me to coerce down precision, while coercing up precision
> with integers. Of course if you consider integer with precision
> infinity (as indeed there are no rounding errors and it is an exact
> ring) then this makes sense again.
Consider the sum of the power series (1 + O(x^5)) and (1 + O(x^6)).
Clearly, the sum should be (2 + O(x^5)), not (2+ O(x^6)) -- that is,
we are effectively "rounding" (1 + O(x^6)) to (1 + O(x^5)) before we
perform the addition.
> But if so then I want to have something like SymbolicNumber which is
> the subset of SymbolicRing that does not contain variables. And that
> this SymbolicNumber is coerced automatically down when used with
> RealField.
make Sage work like this.
Carl
Gary Furnish
support something that interacts with symbolics better.
John Cremona
Symbolic constants would surely be a subring of Symbolic Ring, in
which case that would not be a problem.
But coercing symbolic constants into RR or CC is not a simple, (or
even well-defined?) matter. Just think of many-valued nested
radicals; or if a=sqrt(2), b=sqrt(3), c=sqrt(6), would a*b-c
simplify/coerce to 0? This is not stratightforward at all.
John
> William
>
> >
>
Henryk Trappmann
> But coercing symbolic constants into RR or CC is not a simple, (or
> even well-defined?) matter. Just think of many-valued nested
> radicals; or if a=sqrt(2), b=sqrt(3), c=sqrt(6), would a*b-c
> simplify/coerce to 0? This is not stratightforward at all.
I just would evaluate the expression in RR.
And then sqrt(2)*sqrt(3)-sqrt(6) is not 0.
Btw. I just realized that SymbolicRing does assert that
not (sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)).is_zero()
which is wrong. While your example
(sqrt(2)*sqrt(3)-sqrt(6)).is_zero()
is properly recognized.
I guess for radical expressions there is an algorithm that can decide
zeroness. But if we also involve powers of radical expressions and
perhaps logarithms this becomes more difficult. Take for example:
sqrt(2)**log(9,2) = sqrt(2)**(2*log(3,2))=2**log(3,2)=3
which strangely enough is recognized by Sage:
(sqrt(2)**log(9,2)-3).is_zero()
I would be very grateful if anyone could point me to literature about
those decidability problems (whether for radicals, algebraic numbers,
or numbers involving log and exp, or even sin, cos, etc for which I
think there are some undecidability results).
John Cremona
simplication of algebraic expressions, and I don't want to repeat all
that. Briefly, people tend to think this is easy when they look at
examples which only involve square roots of positive reals, since
there is an "obvious" convention that by default we mean the positive
root. It's a lot more complicated when you deal with general
algebraic numbers which have several ways of being embedded into CC.
Even for square roots of negative reals: you might suggest taking the
root with positive imaginary part, but then sqrt(-2)*sqrt(-3) equals
-sqrt(6) and not +sqrt(6).
John
2008/6/1 Henryk Trappmann <bo19...@googlemail.com>:
Henryk Trappmann
> root.
Roots are clearly defined mono-valued functions:
z.nth_root(n)=e^(log(z)/n)
however this function is not continuous in z, as log is not continuous
at the negative real axis. This makes things complicated.
> It's a lot more complicated when you deal with general
> algebraic numbers which have several ways of being embedded into CC.
> Even for square roots of negative reals: you might suggest taking the
> root with positive imaginary part, but then sqrt(-2)*sqrt(-3) equals
> -sqrt(6) and not +sqrt(6).
sqrt(-2)*sqrt(-3)=e^(log(-2)/2+log(-3)/2)=sqrt(6)e^(-pi*i/2-pi*i/
2)=sqrt(6)*(-1)
I never said it is simple but I am sure that there are equality
deciding algorithms.
And I really want to learn about those.
John Cremona
>
>> there is an "obvious" convention that by default we mean the positive
>> root.
>
> We have to distinguish between solutions of polynomials and roots.
> Roots are clearly defined mono-valued functions:
> z.nth_root(n)=e^(log(z)/n)
> however this function is not continuous in z, as log is not continuous
> at the negative real axis. This makes things complicated.
>
>> It's a lot more complicated when you deal with general
>> algebraic numbers which have several ways of being embedded into CC.
>> Even for square roots of negative reals: you might suggest taking the
>> root with positive imaginary part, but then sqrt(-2)*sqrt(-3) equals
>> -sqrt(6) and not +sqrt(6).
>
> by the above definition this can easily be computed:
> sqrt(-2)*sqrt(-3)=e^(log(-2)/2+log(-3)/2)=sqrt(6)e^(-pi*i/2-pi*i/
> 2)=sqrt(6)*(-1)
>
You have only shifted the ambiguity to the multi-values proprty of
log! As I have often expleined to students, the "generic" proprty
that log(ab)=log(a)+log(b) just does not hold as an identity for
arbitrary complex numbers, and there is no branch convention which
will make that true.
> I never said it is simple but I am sure that there are equality
> deciding algorithms.
> And I really want to learn about those.
If they exist, I would like to see them too. But I am not optimistic.
However there must be a vast literature on the subject, which is at
least as old as computer algebra.
John
> >
>
John Cremona
using Google Scholar): "Algebraic simplification a guide for the
perplexed" 1971, has references back to 1960 at least -- and also
mentioned Axiom.
2008/6/1 John Cremona <john.c...@gmail.com>:
Carl Witty
QQbar, the algebraic numbers. For example:
sage: QQbar(sqrt(-2)*sqrt(-3) + sqrt(6)) == 0
True
In QQbar, radicals are defined to mean the principal root (with a
branch cut just below the negative x axis).
This works by (lazily) building number fields (with embeddings into
CC) that contain the numbers in question. Computation in QQbar is
fairly efficient if you only want to prove inequalities (QQbar
automatically attempts to do this using interval arithmetic) or if the
number field is of low degree, but it can be very inefficient to
construct these number fields if they are of moderate or high degree.
Carl
Henryk Trappmann
> Sage has such a decision procedure built in to its implementation of
> QQbar, the algebraic numbers.
Thank you for this hint. I also found AA the algebraic reals if your
computations only involve reals.
> sage: QQbar(sqrt(-2)*sqrt(-3) + sqrt(6)) == 0
> True
QQbar(sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)) == 0
True
AA(sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)) == 0
True
However is it a *bug* that this does not work in the symbolic Ring
with is_zero?
(sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)).is_zero()
False
And good to know that with QQbar or AA the rule of coercing to the
lowest precision works.
AA(sqrt(2))*1.0
1.41421356237309
QQbar(sqrt(2)) * CC(1.0)
1.41421356237309
But back to SymbolicRing and SymbolicConstant.
I have the following improvement
SUGGESTION: when creating sqrt(2) or other roots from integers, then
assign to them the parent AlgebraicReal or AlgebraicNumer accordingly
instead of the too general Symbolic Ring.
Coercing (see above) and calculating with AA and QQbar works well, for
example:
sqrt(AA(sqrt(2))).parent()
Algebraic Real Field
sqrt(-AA(sqrt(2))).parent()
Algebraic Field
For future development maybe one can introduce a RealSymbolicConstant
(and ComplexSymbolicConstant) which allow application of the common
real functions like +,-,*,/,^, exp, log, sin, cos, etc. However the
difficulty I see here is that we must be able to decide whether such
an expression is 0 if we want to reject expr1/expr2 or log(expr2) for
expr2 being 0. We also must be able to decide whether such an
expression is greater 0. When we want to determine the outcome
(RealSymbolicConstant or ComplexSymbolicConstant) of log(expr) or expr
** (1/n). So maybe this isnt possible from a theoretical point of view
as those questions maybe undecidable (opposed to the algebraic case)
which however is not yet clear for me.
William Stein
<bo19...@googlemail.com> wrote:
>
> On 1 Jun., 18:56, Carl Witty <cwi...@newtonlabs.com> wrote:
>> Sage has such a decision procedure built in to its implementation of
>> QQbar, the algebraic numbers.
>
> But I have to say they are quite hidden in the reference manual.
> Thank you for this hint. I also found AA the algebraic reals if your
> computations only involve reals.
>
>> sage: QQbar(sqrt(-2)*sqrt(-3) + sqrt(6)) == 0
>> True
>
> Using QQbar or AA also my previous example works:
> QQbar(sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)) == 0
> True
> AA(sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)) == 0
> True
>
> However is it a *bug* that this does not work in the symbolic Ring
> with is_zero?
> (sqrt(5+2*sqrt(2)*sqrt(3))-sqrt(2)-sqrt(3)).is_zero()
> False
The definition of is_zero in the Sage symbolic ring is currently
*precisely* the definition in Maxima. In fact, most of the functionality
of the Sage symbolic ring is just a light wrapper around Maxima,
though when we found bugs that Maxima couldn't fix (e.g., in
floor/ceiling), we have worked around them.
Work is under way to completely rewrite the symbolic ring.
> And good to know that with QQbar or AA the rule of coercing to the
> lowest precision works.
> AA(sqrt(2))*1.0
> 1.41421356237309
> QQbar(sqrt(2)) * CC(1.0)
> 1.41421356237309
>
> But back to SymbolicRing and SymbolicConstant.
> I have the following improvement
> SUGGESTION: when creating sqrt(2) or other roots from integers, then
> assign to them the parent AlgebraicReal or AlgebraicNumer accordingly
> instead of the too general Symbolic Ring.
>
That's definitely planned.
> Coercing (see above) and calculating with AA and QQbar works well, for
> example:
> sqrt(AA(sqrt(2))).parent()
> Algebraic Real Field
> sqrt(-AA(sqrt(2))).parent()
> Algebraic Field
>
> For future development maybe one can introduce a RealSymbolicConstant
> (and ComplexSymbolicConstant) which allow application of the common
> real functions like +,-,*,/,^, exp, log, sin, cos, etc. However the
> difficulty I see here is that we must be able to decide whether such
> an expression is 0 if we want to reject expr1/expr2 or log(expr2) for
> expr2 being 0. We also must be able to decide whether such an
> expression is greater 0. When we want to determine the outcome
> (RealSymbolicConstant or ComplexSymbolicConstant) of log(expr) or expr
> ** (1/n). So maybe this isnt possible from a theoretical point of view
> as those questions maybe undecidable (opposed to the algebraic case)
> which however is not yet clear for me.
Thanks for your comments.
William
Carl Witty
> On Mon, Jun 2, 2008 at 1:30 AM, Henryk Trappmann
> > I have the following improvement
> > SUGGESTION: when creating sqrt(2) or other roots from integers, then
> > assign to them the parent AlgebraicReal or AlgebraicNumer accordingly
> > instead of the too general Symbolic Ring.
>
> That's definitely planned.
AA(sqrt(2)) is now, I'm not sure that's a good idea, for two reasons.
First, AA and QQbar by design don't maintain enough information to
print nicely. (This could be improved somewhat from the current
state, but not enough to compete with symbolic radical expressions.)
Second, since AA and QQbar incorporate complete decision procedures,
it is easy to construct examples where they are very, very slow; I
think people would often be happier with the less complete but much
faster techniques used in symbolics.
Carl
William Stein
I think the plan is that algebraic elements won't just be generic symbolic
elements, e.g., sqrt(2) would be a generator for ZZ[sqrt(2)]. This has
been discussed a few times. I didn't mean that using AA or QQbar
by default was precisely what is planned.
William
Robert Bradshaw
Yep. Specifically, the plan is for sqrt(2) to become an element of ZZ
[sqrt(2)] *with* and embedding into RR (so stuff like RR(sqrt(2)) or
even 1.123 + sqrt(2) works). We would want to use very nice AA/QQbar
code to compute, say, sqrt(2) + sqrt(3) (the result would live in a
specific number field with embedding). (Nice) number fields with
embedding would coerce into SR.
- Robert
Gary Furnish
using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires going
through the coercion system which was designed to be elegant instead
of fast, so this becomes insanely slow for any serious use. Finally,
this is going to require serious code duplication from symbolics, so
I'm not sure what the big gain is over just using symbolics to do this
in the first place.
William Stein
>
> -1. First, everything cwitty said is correct. Second, if we start
> using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires going
> through the coercion system which was designed to be elegant instead
> of fast, so this becomes insanely slow for any serious use. Finally,
> this is going to require serious code duplication from symbolics, so
> I'm not sure what the big gain is over just using symbolics to do this
> in the first place.
Also, cwitty pointed out that
sage: sum([sqrt(p) for p in prime_range(1000)])
works fine in Sage now, but with your (and my) proposal,
it would be impossible, since it would require constructing
a ring of integers of a number field of degree 2^168..
-- William
John Cremona
Surely that is something which we can train users not to do? I
remember years ago learning the difference in pari/gp between
sum(n=1,1000,1/n)
and
sum(n=1,1000,1.0/n)
and it's a lesson which only needs to be learnt once.
On the other hand pari/gp always has assumed a lot of its users, and
Sage has aspirations to be easy to use too...
John
>
> -- William
>
> >
>
William Stein
>
> 2008/6/2 William Stein <wst...@gmail.com>:
>>
>> On Mon, Jun 2, 2008 at 12:53 PM, Gary Furnish <gfur...@gfurnish.net> wrote:
>>>
>>> -1. First, everything cwitty said is correct. Second, if we start
>>> using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires going
>>> through the coercion system which was designed to be elegant instead
>>> of fast, so this becomes insanely slow for any serious use. Finally,
>>> this is going to require serious code duplication from symbolics, so
>>> I'm not sure what the big gain is over just using symbolics to do this
>>> in the first place.
>>
>> Also, cwitty pointed out that
>>
>> sage: sum([sqrt(p) for p in prime_range(1000)])
>>
>> works fine in Sage now, but with your (and my) proposal,
>> it would be impossible, since it would require constructing
>> a ring of integers of a number field of degree 2^168..
>
> Surely that is something which we can train users not to do?
Why?
Doing that sum works fine now. It's even fun. That
all the following work fine even with our current system
says something...
sage: a = sum([sqrt(p) for p in primes(1000)])
sage: float(a)
3307.7690992952139
sage: RDF(a)
3307.7690993
sage: RealField(1000)(a)
3307.76909929521532337080828563375122663883932558794796752652934831207623440969500409157102585876606649295101934045235119416217303162488345812238151772485033866010535293142064829106852999118135254502207045010375309268967623016528462090375140033057342943102054036545832373683901618507572426685968323860
> I remember years ago learning the difference in pari/gp between
>
> sum(n=1,1000,1/n)
>
> and
>
> sum(n=1,1000,1.0/n)
>
> and it's a lesson which only needs to be learnt once.
>
> On the other hand pari/gp always has assumed a lot of its users, and
> Sage has aspirations to be easy to use too...
>
> John
>
>>
>> -- William
>>
>> >
>>
>
> >
>
--
Gonzalo Tornaria
> Doing that sum works fine now. It's even fun. That
> all the following work fine even with our current system
> says something...
>
> sage: a = sum([sqrt(p) for p in primes(1000)])
> sage: float(a)
> 3307.7690992952139
> sage: RDF(a)
> 3307.7690993
> sage: RealField(1000)(a)
> 3307.76909929521532337080828563375122663883932558794796752652934831207623440969500409157102585876606649295101934045235119416217303162488345812238151772485033866010535293142064829106852999118135254502207045010375309268967623016528462090375140033057342943102054036545832373683901618507572426685968323860
sage: a = e - 1
sage: for i in range(2,100):
....: a = (a-1)*i
....:
sage: float(a)
-1.3491675327339235e+140
sage: RDF(a)
-1.34916753273e+140
sage: RealField(1000)(a)
1.01009999010191094737330783479071750558784854694330666719781458645331617143268163560343706814329184685328514081778158205141322918688406340340355792675106115099850396387113493402064121653102597587035584599866091432567502680436745983344808432577075043013574712846877319486115708287043275014203343194844
sage: a = e - 1
sage: for i in range(2,200):
....: a = (a-1)*i
....:
sage: float(a)
-inf
sage: RealField(100)(a)
-5.6297658631950881996939935452e342
sage: RealField(200)(a)
1.1804891378835970579289419600656355595291746689227892615143e312
sage: RealField(10000)(a)
1.00502499....
--
Gonzalo
Gonzalo Tornaria
'exceptions.RuntimeError'>: maximum recursion depth exceeded
when evaluating float(a)... (sage 3.0)
Gonzalo
William Stein
<joke>
Oh my god, floating point numbers are just broken!
</joke>
Seriously, floating point numbers are what they are.
-- William
Gonzalo Tornaria
> > sage: a = e - 1
> > sage: for i in range(2,200):
> > ....: a = (a-1)*i
> > ....:
> > sage: float(a)
> > -inf
>
> Oh my god, floating point numbers are just broken!
> </joke>
:-)
> Seriously, floating point numbers are what they are.
Fair enough. The point I tried to make (not claiming anything is
broken) is perhaps complemented by the following:
sage: expand(a)
3943289336823952517761816069660925311475679888435866316473712666221797249817016714601521420059923119520886060694598194151288213951213185525309633124764149655567314286353816586186984944719612228107258321201270166459320656137141474266387621212037869516201606287027897843301130159520851620311758504293980894611113948118519486873600000000000000000000000000000000000000000000000*e
- 10718971748644869545882855035558709219326197160960676962072649602582551264719783234505244404668998392076585003319325446678612269966334945842611842879092400407644778161460862429166937616985702801147176754507686032298026763211111234497467428555982733907104528850484797214679822706672415691688065756374587435109932151999457607946738407969860677357914214979720210208723664522819
This is a particular symbolic expression, which has a particular value
in RR, whose best approximation in RDF is far from -inf or
sage: RDF(expand(a))
nan
The "philosophical" question is whether:
sage: RDF(A*e - B)
(where A*e is damn close to B+1)
should be evaluated as
sage: RDF(A)*RDF(e) - RDF(B)
which is definitely nan since RDF(A) == RDF(B) == +inf.
What I mean is, the symbolic expression has (in principle) enough
information to figure out that RDF(A*e-B) is close to 1. Whether it is
reasonable to expect that is the question.
Best, Gonzalo
PS: for the record, this seems to be "broken" in the same way in
mathematica, for one.
Robert Bradshaw
> On Mon, Jun 2, 2008 at 12:53 PM, Gary Furnish
> <gfur...@gfurnish.net> wrote:
>>
>> -1. First, everything cwitty said is correct.
>> Second, if we start
>> using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires
>> going
>> through the coercion system which was designed to be elegant instead
>> of fast, so this becomes insanely slow for any serious use.
The coercion system is designed to be elegant *and* fast. Writing
something like 1 + sqrt(2) is going to require the coercion system no
matter what we do, as is 1 + sqrt(2) + 1/2. Computing QQ(sqrt(2), sqrt
(3)) may take a millisecond or two, but after that coercion into it
from ether ring will be fast.
>> Finally, this is going to require serious code duplication from
>> symbolics, so
>> I'm not sure what the big gain is over just using symbolics to do
>> this
>> in the first place.
Number field element work completely differently than symbolics, I
see little if any code duplication.
> Also, cwitty pointed out that
>
> sage: sum([sqrt(p) for p in prime_range(1000)])
>
> works fine in Sage now, but with your (and my) proposal,
> it would be impossible, since it would require constructing
> a ring of integers of a number field of degree 2^168..
This is the biggest argument against such a proposal, and I'm not
quite sure how best to handle this. One would have to implement large
number fields over sparse polynomials, and only lazily compute the
ring of integers. Or, as John Cremona suggested, "train users." None
of the above are ideal.
I would like to give my motivation for not having sqrt(2) be in SR:
1) Speed. I know you're working very hard to make symbolics much,
much faster than they currently are, but I still don't think it'll be
able to compete in this very specialized domain. Currently:
sage: timeit("((1+sqrt(2))^100).expand()")
5 loops, best of 3: 52.9 ms per loop
sage: timeit("((1+sqrt(2)+sqrt(3))^50).expand()")
5 loops, best of 3: 579 ms per loop
sage: sym_a = sqrt(2) + sqrt(3)
sage: timeit("((1+sym_a)^50).expand()")
5 loops, best of 3: 576 ms per loop
Compared to
sage: K.<a> = NumberField(x^2-2)
sage: timeit("((1+a)^100)")
625 loops, best of 3: 48.4 µs per loop
sage: K.<a> = NumberField(x^4 - 10*x^2 + 1)
sage: timeit("((1+a)^50)")
625 loops, best of 3: 138 µs per loop
That's over three orders of magnitude faster (and it's *not* due to
pexpect/interface overhead as the actual input and output are
relatively small). Making arithmetic involving a couple of radicals
fast should probably be the most important, especially as one starts
making structures over them.
2) I personally don't like having to sprinkle the "expand" and and/or
"simplify" all over the place. Now I don't think symbolics should be
expanded automatically, but stuff like (1+sqrt(2))^2 should be or 1/(1
+i). It's like just getting the question back. (I guess I'm revealing
my bias that I don't think of it as living in SR, but rather a number
field...) On that note, I can't even figure out how to do simplify
"(sqrt(2)-3)/(sqrt(2)-1)" in the symbolics...as opposed to
sage: K.<sqrt2> = NumberField(x^2-2)
sage: (sqrt2-3)/(sqrt2-1)
-2*sqrt2 - 1
sage: 1/(sqrt2-1)
sqrt2 + 1
3) The coercion system works best when things start as high up the
tree as they can, and the Symbolic Ring is like the big black hole at
the bottom that sucks everything in (and makes it slow). There is a
coercion from ZZ[sqrt(2)] (with embedding) to SR, but not the other
way around, and even trying to cast the other way is problematic. I'd
rather that matrix([1, 1/2, 0, sqrt(2)]) land in a matrix space over
the a number field (rather than over SR), and ZZ['x'].gen() + sqrt(2)
be an actual polynomial in x. Also, the SR, though very useful,
somehow seems less rigorous (I'm sure that this is improving).
Of course many of these are a matter of opinion/taste, but I think
they are worth considering.
- Robert
Gary Furnish
<robe...@math.washington.edu> wrote:
>
> On Jun 2, 2008, at 12:55 PM, William Stein wrote:
>
>> On Mon, Jun 2, 2008 at 12:53 PM, Gary Furnish
>> <gfur...@gfurnish.net> wrote:
>>>
>>> -1. First, everything cwitty said is correct.
>
>
>>> Second, if we start
>>> using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires
>>> going
>>> through the coercion system which was designed to be elegant instead
>>> of fast, so this becomes insanely slow for any serious use.
>
> The coercion system is designed to be elegant *and* fast. Writing
> something like 1 + sqrt(2) is going to require the coercion system no
> matter what we do, as is 1 + sqrt(2) + 1/2. Computing QQ(sqrt(2), sqrt
> (3)) may take a millisecond or two, but after that coercion into it
> from ether ring will be fast.
>
critical path that symbolics has to use, the argument that coercion
was written to be fast makes no sense to me.
>>> Finally, this is going to require serious code duplication from
>>> symbolics, so
>>> I'm not sure what the big gain is over just using symbolics to do
>>> this
>>> in the first place.
>
> Number field element work completely differently than symbolics, I
> see little if any code duplication.
>
do much better then maxima, but its not going to be that much better
(maybe if we encapsulate number fields as a special case in SR)
> 2) I personally don't like having to sprinkle the "expand" and and/or
> "simplify" all over the place. Now I don't think symbolics should be
> expanded automatically, but stuff like (1+sqrt(2))^2 should be or 1/(1
> +i). It's like just getting the question back. (I guess I'm revealing
> my bias that I don't think of it as living in SR, but rather a number
> field...) On that note, I can't even figure out how to do simplify
> "(sqrt(2)-3)/(sqrt(2)-1)" in the symbolics...as opposed to
>
> sage: K.<sqrt2> = NumberField(x^2-2)
> sage: (sqrt2-3)/(sqrt2-1)
> -2*sqrt2 - 1
> sage: 1/(sqrt2-1)
> sqrt2 + 1
>
in Mathematica and Maple is wrong. This makes sense for number theory
but not for people using calculus.
> 3) The coercion system works best when things start as high up the
> tree as they can, and the Symbolic Ring is like the big black hole at
> the bottom that sucks everything in (and makes it slow). There is a
> coercion from ZZ[sqrt(2)] (with embedding) to SR, but not the other
> way around, and even trying to cast the other way is problematic. I'd
> rather that matrix([1, 1/2, 0, sqrt(2)]) land in a matrix space over
> the a number field (rather than over SR), and ZZ['x'].gen() + sqrt(2)
> be an actual polynomial in x. Also, the SR, though very useful,
> somehow seems less rigorous (I'm sure that this is improving).
>
When coercion is faster we can consider changing this. My definition
of fast is "<10 cycles if the parents are the same, no dictionary
lookups if one parent is in the other for all common cases, otherwise
reasonablely quick pure Cython code. New and old coercion fail these
tests of sufficiently quick, and I'm not waiting to finish symbolics
until they do pass those tests.
My alternative option is lets throw in a flag, defaults to off
(current behavior) that lets you turn on sqrt/powers as in number
theory by default instead of SR. This makes the code perform as
expected by end users, and advanced users can use number fields if
they know they are appropriate. This is just largely a one if
statement change in the dispatch of sqrt so this should be reasonably
safe.
Robert Bradshaw
>
> On Mon, Jun 2, 2008 at 11:41 PM, Robert Bradshaw
> <robe...@math.washington.edu> wrote:
>>
>> On Jun 2, 2008, at 12:55 PM, William Stein wrote:
>>
>>> On Mon, Jun 2, 2008 at 12:53 PM, Gary Furnish
>>> <gfur...@gfurnish.net> wrote:
>>>>
>>>> -1. First, everything cwitty said is correct.
>>
>> More on this below.
>>
>>>> Second, if we start
>>>> using ZZ[sqrt(2)] and ZZ[sqrt(3)], then sqrt(2)+sqrt(3) requires
>>>> going
>>>> through the coercion system which was designed to be elegant
>>>> instead
>>>> of fast, so this becomes insanely slow for any serious use.
>>
>> The coercion system is designed to be elegant *and* fast. Writing
>> something like 1 + sqrt(2) is going to require the coercion system no
>> matter what we do, as is 1 + sqrt(2) + 1/2. Computing QQ(sqrt(2),
>> sqrt
>> (3)) may take a millisecond or two, but after that coercion into it
>> from ether ring will be fast.
>>
> As long as there are classes in pure python that use MI on the
> critical path that symbolics has to use, the argument that coercion
> was written to be fast makes no sense to me.
Not sure what you mean by "MI" here, could you explain. In any case,
just because coercion isn't as fast as it could be doesn't mean that
it's not written for speed and much faster than it used to be. Of
course there's room for improvement, but right now the focus is
trying to finish the new system (which isn't really that "new"
compared to the change made a year ago) in place.
I'd rather have them be a number field elements (with all the
methods, etc) over providing a wrapping in SR. Otherwise one ends up
with something like pari, where everything just sits in the same parent.
>> 2) I personally don't like having to sprinkle the "expand" and and/or
>> "simplify" all over the place. Now I don't think symbolics should be
>> expanded automatically, but stuff like (1+sqrt(2))^2 should be or
>> 1/(1
>> +i). It's like just getting the question back. (I guess I'm revealing
>> my bias that I don't think of it as living in SR, but rather a number
>> field...) On that note, I can't even figure out how to do simplify
>> "(sqrt(2)-3)/(sqrt(2)-1)" in the symbolics...as opposed to
>>
>> sage: K.<sqrt2> = NumberField(x^2-2)
>> sage: (sqrt2-3)/(sqrt2-1)
>> -2*sqrt2 - 1
>> sage: 1/(sqrt2-1)
>> sqrt2 + 1
>>
> Your going to have a hard time convincing me that the default behavior
> in Mathematica and Maple is wrong. This makes sense for number theory
> but not for people using calculus.
OK, this is a valid point, though the non-calculus portions (and
emphasis) of Sage are (relatively) more significant. Sage is not a
CAS, that is just one (important) piece of it.
Maple does
> 1/(1+I);
1/2 - 1/2 I
at least. Looking to the M's for ideas is good, but they should not
always dictate how we do things--none but Magma has the concept of
parents/elements, and Sage uses a very OO model which differs from
all of them. Why doesn't it make sense for Mathematica/Maple? I think
it's because they view simplification (or even deciding to simplify)
as expensive.
>> 3) The coercion system works best when things start as high up the
>> tree as they can, and the Symbolic Ring is like the big black hole at
>> the bottom that sucks everything in (and makes it slow). There is a
>> coercion from ZZ[sqrt(2)] (with embedding) to SR, but not the other
>> way around, and even trying to cast the other way is problematic. I'd
>> rather that matrix([1, 1/2, 0, sqrt(2)]) land in a matrix space over
>> the a number field (rather than over SR), and ZZ['x'].gen() + sqrt(2)
>> be an actual polynomial in x. Also, the SR, though very useful,
>> somehow seems less rigorous (I'm sure that this is improving).
>>
> When coercion is faster we can consider changing this.
Coercion speed is irrelevant to the issues I mentioned here... and as
coercion+number fields is *currently* faster than what you could hope
to get with SR (the examples above all involve coercion) it wouldn't
help your case either.
> My definition
> of fast is "<10 cycles if the parents are the same,
Python semantics tie our hands a bit here, but I think we're about as
close as we can get.
> no dictionary lookups if one parent is in the other for all common
> cases,
Would this mean hard-coding all common paths? Currently there is a
single dictionary lookup for common cases (and not a Python dict).
> otherwise reasonablely quick pure Cython code.
Yes, it should be fast, but only has to be done once and then it's
cached. Of course the code specific to the ring/elements is as fast
or slow as whoever wrote it.
> New and old coercion fail these
> tests of sufficiently quick, and I'm not waiting to finish symbolics
> until they do pass those tests.
Thanks for your concise definition--if you have any input of how to
make things faster without hard-coding tons of special cases I'd be
very glad to hear (though the current focus is getting things merged
back into the main sage branch before we focus on optimization again).
> My alternative option is lets throw in a flag, defaults to off
> (current behavior) that lets you turn on sqrt/powers as in number
> theory by default instead of SR. This makes the code perform as
> expected by end users, and advanced users can use number fields if
> they know they are appropriate. This is just largely a one if
> statement change in the dispatch of sqrt so this should be reasonably
> safe.
"Perform as expected" is what we disagree on, though I'd call us both
end- and advanced users. I generally dislike the idea of flags the
user sets that change behavior, but it's an option to consider.
- Robert
John Cremona
2008/6/3 Robert Bradshaw <robe...@math.washington.edu>:
Robert Bradshaw
> MI = Multiple Inheritance?
That's the only thing I can think of, but there isn't any of that in
the coercion model...
- Robert
root
>
> > 1/(1+I);
> 1/2 - 1/2 I
Axiom does
1/(1+%i)
1
------
1 + %i
which is of type Fraction Complex Integer, that is a fraction
whose numerator and denominator are of type Complex(Integer).
You can ask Axiom to place the result in a different type:
1/(1+%i)::Complex Fraction Integer
1 1
- - - %i
2 2
and the type is Complex Fraction Integer, that is a Complex
number whose real and imaginary components are Fraction(Integer).
The "simple form" depends on what is canonical for the type.
Tim
Gary Furnish
get used in coercion, and use MI, making them impossible to cythonize.
simplification), but (x+1)^2 shouldn't get simplified under any
circumstances. This has (little) do with speed (for this small of
exponent) and everything to do with being consistent with the high
degree cases and keeping expressions uncluttered.
> at least. Looking to the M's for ideas is good, but they should not
> always dictate how we do things--none but Magma has the concept of
> parents/elements, and Sage uses a very OO model which differs from
> all of them. Why doesn't it make sense for Mathematica/Maple? I think
> it's because they view simplification (or even deciding to simplify)
> as expensive.
>
>>> 3) The coercion system works best when things start as high up the
>>> tree as they can, and the Symbolic Ring is like the big black hole at
>>> the bottom that sucks everything in (and makes it slow). There is a
>>> coercion from ZZ[sqrt(2)] (with embedding) to SR, but not the other
>>> way around, and even trying to cast the other way is problematic. I'd
>>> rather that matrix([1, 1/2, 0, sqrt(2)]) land in a matrix space over
>>> the a number field (rather than over SR), and ZZ['x'].gen() + sqrt(2)
>>> be an actual polynomial in x. Also, the SR, though very useful,
>>> somehow seems less rigorous (I'm sure that this is improving).
>>>
>> When coercion is faster we can consider changing this.
>
> Coercion speed is irrelevant to the issues I mentioned here... and as
> coercion+number fields is *currently* faster than what you could hope
> to get with SR (the examples above all involve coercion) it wouldn't
> help your case either.
that for endusers, sum(sqrt(p)) behaves as expected.
>
>> My definition
>> of fast is "<10 cycles if the parents are the same,
>
> Python semantics tie our hands a bit here, but I think we're about as
> close as we can get.
>
>> no dictionary lookups if one parent is in the other for all common
>> cases,
>
> Would this mean hard-coding all common paths? Currently there is a
> single dictionary lookup for common cases (and not a Python dict).
>
statements away (and not a bunch of virtual calls either. They cost
performance too. No more then one should be necessary for the common
case (the code to handle this can probably go in the
addition/multiplication handlers)). Then if that fails we can take
the cached dict lookup route. Make the common case fast at the
expense of the uncommon case.
>> otherwise reasonablely quick pure Cython code.
>
> Yes, it should be fast, but only has to be done once and then it's
> cached. Of course the code specific to the ring/elements is as fast
> or slow as whoever wrote it.
>
>> New and old coercion fail these
>> tests of sufficiently quick, and I'm not waiting to finish symbolics
>> until they do pass those tests.
>
> Thanks for your concise definition--if you have any input of how to
> make things faster without hard-coding tons of special cases I'd be
> very glad to hear (though the current focus is getting things merged
> back into the main sage branch before we focus on optimization again).
>
It is more important for coercion to be fast (especially if we are
using it internally in algorithms) then for it to be pretty (although
it can still be designed in a mathematically rigorous way, the code
that actually implements it may not be pretty)
>> My alternative option is lets throw in a flag, defaults to off
>> (current behavior) that lets you turn on sqrt/powers as in number
>> theory by default instead of SR. This makes the code perform as
>> expected by end users, and advanced users can use number fields if
>> they know they are appropriate. This is just largely a one if
>> statement change in the dispatch of sqrt so this should be reasonably
>> safe.
>
> "Perform as expected" is what we disagree on, though I'd call us both
> end- and advanced users. I generally dislike the idea of flags the
> user sets that change behavior, but it's an option to consider.
>
understand why he can't perform a sum of the sqrt of some primes. If
we are to be a viable alternative for non-research mathematicians we
can't run off and implement things that drastically change the
complexity of simple operations.
> - Robert
>
>
> >
>
William Stein
> The average calculus student coming from maple is not going to
> understand why he can't perform a sum of the sqrt of some primes. If
> we are to be a viable alternative for non-research mathematicians we
> can't run off and implement things that drastically change the
> complexity of simple operations.
Let me take a moment to totally hijack this thread and say that
I think we should drastically change the complexity of simple
operations... for the better :-) There's just a *shocking* number
of "simple operations" out there in Maple/Mathematica, etc.,
which we can and should do much faster in Sage.
OK, back to the thread...
-- William
Soroosh Yazdani
<snip>
Your going to have a hard time convincing me that the default behavior
in Mathematica and Maple is wrong. This makes sense for number theory
but not for people using calculus.
Hmm, I must say from using maple on expressions like this, I found the answers that it gave were completely pointless at times, and I was forced to run expand and simplify many times, and use many tricks to get the answer in the most simplified form. At times, the easiest way was to evaluate the expressions, and then use LLL to get the expression back. Admittedly I am a number theorist, although at the time when I was using maple extensively, I was still in undergrad. When I started using MAGMA, the learning curve seemed considerably higher, but completely worth the trouble for algebraic expressions. The current proposal by Robert seems to be the best of both world for my taste. (And I agree with him that a lot of it is a matter of taste.)
As an aside, one of my gripes teaching calculus is that the students don't simplify expressions like sqrt(2)sqrt(3) or 1/(1+i), and I would prefer if SAGE doesn't contribute to such atrocity. :)
Soroosh
Robert Bradshaw
>>> As long as there are classes in pure python that use MI on the
>>> critical path that symbolics has to use, the argument that coercion
>>> was written to be fast makes no sense to me.
>>
>> Not sure what you mean by "MI" here, could you explain. In any case,
>> just because coercion isn't as fast as it could be doesn't mean that
>> it's not written for speed and much faster than it used to be. Of
>> course there's room for improvement, but right now the focus is
>> trying to finish the new system (which isn't really that "new"
>> compared to the change made a year ago) in place.
>>
> Sets, and in particular a bunch of the category functionality (homset)
> get used in coercion, and use MI, making them impossible to cythonize.
Ah, yes. Homsets. They're not used anywhere in the critical path
though. (If so, that should be fixed.)
>>
I agree that (x+1)^2 shouldn't get simplified, but for me this has a
very different feel than (1+I)^2 or (1+sqrt(2))^2.
n-th root would have a similar speed increase, but other than those
two cases I don't see one wanting algebraic extensions (short of
explicitly making a number field).
>>
>>> My definition
>>> of fast is "<10 cycles if the parents are the same,
>>
>> Python semantics tie our hands a bit here, but I think we're about as
>> close as we can get.
>>
>>> no dictionary lookups if one parent is in the other for all common
>>> cases,
>>
>> Would this mean hard-coding all common paths? Currently there is a
>> single dictionary lookup for common cases (and not a Python dict).
>>
> Common cases should be no more then a virtual call and a few if
> statements away (and not a bunch of virtual calls either. They cost
> performance too. No more then one should be necessary for the common
> case (the code to handle this can probably go in the
> addition/multiplication handlers)). Then if that fails we can take
> the cached dict lookup route. Make the common case fast at the
> expense of the uncommon case.
I am -1 for hard-coding knowledge and logic about ZZ, QQ, RR, RDF,
CC, CDF, ... into the coercion model.
If we can change the complexity for the better we should. It's a
tradeoff of speed for code with a small number of radicals vs. speed
with a large number of radicals.
- Robert
Gary Furnish
for determining the result of something in parent1 op something in
parent2 (maybe it was with someone else?) I'm also -1 for hard-coding
knowledge and logic about ZZ,QQ, etc into the coercion model. I am +1
for hardcoding it into the elements of say, ZZ,QQ,RR and then having
them call the coercion model only if those hardcodes can't figure the
situation out.
Robert Bradshaw
> I think we had a discussion on irc about how homsets still got used
> for determining the result of something in parent1 op something in
> parent2 (maybe it was with someone else?)
I sure hope not. If so, that needs to change (but I'm pretty sure it
isn't).
> I'm also -1 for hard-coding
> knowledge and logic about ZZ,QQ, etc into the coercion model. I am +1
> for hardcoding it into the elements of say, ZZ,QQ,RR and then having
> them call the coercion model only if those hardcodes can't figure the
> situation out.
That sounds much better, though I'm still not a fan.
Gary Furnish
expressions x^20+20*x^19+......... doesn't tell me much about where an
expression comes from. (x+5)^20 tells me a bunch. Expanding
expressions generally causes information loss for many calculus and
physics problems and going overboard can be bad (although this isn't
an issue for number theory). Furthermore sqrt(2)sqrt(3) is not
necessarily equal to sqrt(6), so not simplifying there is appropriate
in most circumstances.
Gary Furnish
<robe...@math.washington.edu> wrote:
>
> On Jun 3, 2008, at 11:06 AM, Gary Furnish wrote:
>
>> I think we had a discussion on irc about how homsets still got used
>> for determining the result of something in parent1 op something in
>> parent2 (maybe it was with someone else?)
>
> I sure hope not. If so, that needs to change (but I'm pretty sure it
> isn't).
>
something is in if you have parent1 op parentt2 that is new in this
system (but I can't find said function). Allegedly said function has
to use Homsets, so performance is going to be horrible(and this makes
sense, because its what I have to use now). I consider homsets to be
a gigantic flaw in coercion that absolutely have to be fixed for me to
consider using more of the coercion system in symbolics.
>> I'm also -1 for hard-coding
>> knowledge and logic about ZZ,QQ, etc into the coercion model. I am +1
>> for hardcoding it into the elements of say, ZZ,QQ,RR and then having
>> them call the coercion model only if those hardcodes can't figure the
>> situation out.
>
> That sounds much better, though I'm still not a fan.
>
Sure, its kindof ugly, but it will give us the performance we need,
and I don't see a better way to allow us to use coercions all over the
place without having performance drop to 0.
Henryk Trappmann
> At a bare minimum there is never a canonical (automatic)
> coercion from elements of R to elements of S unless that coercion
> is defined (as a homomorphism) on all of R.
I dont want to be heretical by why is it so important that coercion is
totally defined?
Considering all the current discussion, I get the impression that
allowing partial coercion would solve much of the arising problems
nearly without effort.
SR can be easily partially coerced to any number system. Just try to
evaluate the SR term in that number system. Either it succeeds for
example
log(2) or sqrt(2) directly translates to RR by log(RR(2)) and
sqrt(RR(2)) or it fails, for example RR(x) or RR(1)/RR(0).
While sqrt(RR(-2)) is somehow mixed, depending on the sign of the
argument sqrt returns either something of a complex or a real field.
Which is quite desirable. So that means if one would try to coerce
from SR to RR it can happen that the result is a CC, what however
would perfectly fit, i.e. would perfectly be what one expects when one
writes: sqrt(-2)*1.0.
I hope everyone agrees that the above two (but at least the above
first) paragraph(s) is what one would expect/desire.
Then for coercion would only matter the exactness of a structure and
the subset relation. The general rule would be: If there one structure
is less exact then coerce to this structure. If we have two structures
with equal exactness coerce to the substructure if possible, in doubt
coerce to the superstructure.
This would also solve the performance and representation issue. As
long as no coercing takes places just compute in whatever Ring you
want. But if it comes to coercion/evaluation/decision it may take
somewhat longer (what also would be what I expect).
Robert Bradshaw
> On Tue, Jun 3, 2008 at 12:11 PM, Robert Bradshaw
> <robe...@math.washington.edu> wrote:
>>
>> On Jun 3, 2008, at 11:06 AM, Gary Furnish wrote:
>>
>>> I think we had a discussion on irc about how homsets still got used
>>> for determining the result of something in parent1 op something in
>>> parent2 (maybe it was with someone else?)
>>
>> I sure hope not. If so, that needs to change (but I'm pretty sure it
>> isn't).
>>
> Well allegedly there is some function that computes what parent
> something is in if you have parent1 op parentt2 that is new in this
> system (but I can't find said function). Allegedly said function has
> to use Homsets, so performance is going to be horrible(and this makes
> sense, because its what I have to use now).
When a new morphism is created it needs a parent, which is a Homset
that may be looked up/created at that time. This is probably what you
are talking about. However, this is a single (tiny) constant overhead
over the entire runtime of Sage. I'm sure this could be further
optimized, but creating all the homsets ZZ -> QQ -> RR -> CC, etc.
will be done long before any symbolic code gets hit.
> I consider homsets to be
> a gigantic flaw in coercion that absolutely have to be fixed for me to
> consider using more of the coercion system in symbolics.
Ironically, other people see it as a plus that coercion has been
given a more categorical founding.
>
>>> I'm also -1 for hard-coding
>>> knowledge and logic about ZZ,QQ, etc into the coercion model. I
>>> am +1
>>> for hardcoding it into the elements of say, ZZ,QQ,RR and then having
>>> them call the coercion model only if those hardcodes can't figure
>>> the
>>> situation out.
>>
>> That sounds much better, though I'm still not a fan.
>>
> Sure, its kindof ugly, but it will give us the performance we need,
> and I don't see a better way to allow us to use coercions all over the
> place without having performance drop to 0.
One may be able to eek out a bit more performance by doing this, but
it's not as if performance is awful in the current model.
- Robert
Robert Bradshaw
This is a good question, and has been extensively discussed before,
but I will try and justify it a bit here. In Sage there are two kinds
of conversions from one ring to another: casting and coercion. If I
have an element a in S, and I want to view it as an element in R, I
can write R(a). This should work whenever it's possible to make sense
of a in R. Coercion, on the other hand, is implicit and thus should
be fairly conservative. To reduce ambiguity, it should also go
exactly one direction (unless, perhaps, there is a canonical
isomorphism). It makes more sense for there to be a coercion RR -> SR
than the other way around, so that when one writes "sqrt(x) * 2.4" it
tries to give a result in SR rather than an error due to failure to
give a result in RR.
I'm not sure how allowing partial coercions would help the situation.
The problems are (1) given a and b in different rings, how to quickly
find a' and b' so that a' + b' makes sense and dispatch to that
operation and (2) whether or not sqrt(2) should start out as a SR
object, or something more specific the same way integers and
rationals do, and then coerced down to SR if need be.
- Robert
root
>The problems are (1) given a and b in different rings, how to quickly
>find a' and b' so that a' + b' makes sense and dispatch to that
>operation and (2) whether or not sqrt(2) should start out as a SR
>object, or something more specific the same way integers and
>rationals do, and then coerced down to SR if need be.
Axiom does this in two different ways, depending on whether the
code is interpreted or compiled.
Interpreted code does the lookup dynamically and looks at the domains.
Domains (types) can export a coerce operation that will coerce things
from other types into the current type, e.g.
Polynomial(Integer) exports two coerce operations:
coerce : Symbol -> Polynomial(Integer)
coerce : Integer -> Polynomial(Integer)
and the interpreter can automatically use these to lift the input
types to the Polynomial type. Searching for and using these coerce
operations take time in the interpreter. Thus,
1 + x
finds 9 plus operations that do not apply and one that can apply which
allows the interpreter to automatically coerce
1 -> Polynomial(Integer)
x -> Polynomial(Integer)
and Polynomial(Integer) has a plus operation
+(POLY(INT),POLY(INT)) -> POLY(INT)
so the result gets automatically coerced to Polynomial(Integer).
The other way occurs during compiles. The compiler requires you
to specify the types explicitly so it can lay down fast code.
At runtime, compiled code does a hardcoded index into the vector
that implements the Domain type to invoke the proper function.
Of course, Sage uses Python so compiling for speed isn't possible.
So allowing classes to export coerce operations that can be
automatically applied by the interpreter is effective, if not
the most efficient method. Classes are the only place that can
encapsulate the "special case knowledge" necessary to do the
coercions. If the coercions are all buried in some external,
general purpose coercion routine with special case knowledge then,
by design, this routine needs to know everything about everything.
This makes it impossible to add new classes that "just work".
Tim
Bill Page
>
> On Jun 3, 2008, at 11:17 AM, Gary Furnish wrote:
>> I consider homsets to be a gigantic flaw in coercion that
>> absolutely have to be fixed for me to consider using more
>> of the coercion system in symbolics.
>
> Ironically, other people see it as a plus that coercion has
> been given a more categorical founding.
>
Absolutely! :-)
BTW, where can I read more about these categorical concepts that are
currently built-in or planned for Sage?
Regards,
Bill Page.
William Stein
Gary Furnish
<robe...@math.washington.edu> wrote:
>
> On Jun 3, 2008, at 11:17 AM, Gary Furnish wrote:
>
>> On Tue, Jun 3, 2008 at 12:11 PM, Robert Bradshaw
>> <robe...@math.washington.edu> wrote:
>>>
>>> On Jun 3, 2008, at 11:06 AM, Gary Furnish wrote:
>>>
>>>> I think we had a discussion on irc about how homsets still got used
>>>> for determining the result of something in parent1 op something in
>>>> parent2 (maybe it was with someone else?)
>>>
>>> I sure hope not. If so, that needs to change (but I'm pretty sure it
>>> isn't).
>>>
>> Well allegedly there is some function that computes what parent
>> something is in if you have parent1 op parentt2 that is new in this
>> system (but I can't find said function). Allegedly said function has
>> to use Homsets, so performance is going to be horrible(and this makes
>> sense, because its what I have to use now).
>
> When a new morphism is created it needs a parent, which is a Homset
> that may be looked up/created at that time. This is probably what you
> are talking about. However, this is a single (tiny) constant overhead
> over the entire runtime of Sage. I'm sure this could be further
> optimized, but creating all the homsets ZZ -> QQ -> RR -> CC, etc.
> will be done long before any symbolic code gets hit.
>
can't access members of homset without leaving cython and using
python.
>> I consider homsets to be
>> a gigantic flaw in coercion that absolutely have to be fixed for me to
>> consider using more of the coercion system in symbolics.
>
> Ironically, other people see it as a plus that coercion has been
> given a more categorical founding.
>
of categories to be bad. Implementations that make extensive use of
MI and are thus impossible to cythonize or access without py dict
lookups are not good implementations of categories. If coercion was
implemented with 100% pure Cython code (with an eye for speed where it
is needed), I would be significantly less upset with it then I am now
where people tell me that if I need more speed I'm out of luck.
>>
>>>> I'm also -1 for hard-coding
>>>> knowledge and logic about ZZ,QQ, etc into the coercion model. I
>>>> am +1
>>>> for hardcoding it into the elements of say, ZZ,QQ,RR and then having
>>>> them call the coercion model only if those hardcodes can't figure
>>>> the
>>>> situation out.
>>>
>>> That sounds much better, though I'm still not a fan.
>>>
>> Sure, its kindof ugly, but it will give us the performance we need,
>> and I don't see a better way to allow us to use coercions all over the
>> place without having performance drop to 0.
>
> One may be able to eek out a bit more performance by doing this, but
> it's not as if performance is awful in the current model.
>
system anywhere near as much as symbolics in Sage does.
> - Robert
>
>
> >
>
Robert Bradshaw
>
> On Tue, Jun 3, 2008 at 2:48 PM, Robert Bradshaw
> <robe...@math.washington.edu> wrote:
>>
>> When a new morphism is created it needs a parent, which is a Homset
>> that may be looked up/created at that time. This is probably what you
>> are talking about. However, this is a single (tiny) constant overhead
>> over the entire runtime of Sage. I'm sure this could be further
>> optimized, but creating all the homsets ZZ -> QQ -> RR -> CC, etc.
>> will be done long before any symbolic code gets hit.
>>
> This is not what I'm talking about. What I'm talking about is you
> can't access members of homset without leaving cython and using
> python.
Of course homsets could be optimized, but once created they aren't
used in the coercion itself. I don't know why you'd need to access
the members of homset unless you were doing some very high-level
programming. The domain and codomain are cdef attributes of
Morphisms, which are in Cython (and that took a fair amount of work
as they was a lot of MI going on with them until a year ago).
>>> I consider homsets to be
>>> a gigantic flaw in coercion that absolutely have to be fixed for
>>> me to
>>> consider using more of the coercion system in symbolics.
>>
>> Ironically, other people see it as a plus that coercion has been
>> given a more categorical founding.
>>
> No, I consider categories to be good. I consider bad implementations
> of categories to be bad. Implementations that make extensive use of
> MI and are thus impossible to cythonize or access without py dict
> lookups are not good implementations of categories.
That depends on what they are being used for, but categories lend
themselves very naturally to multiple inheritance because of what
they are mathematically. I understand wanting, .e.g., arithmetic to
be super fast, but I don't see the gain in hacks/code duplication to
strip out multiple inheritance from and cythonize categories. Python
is a beautiful language, but it comes with a cost (e.g. dict
lookups). Other than initial homset creation, what would be faster
(for you) if homsets and categories were written in Cython?
> If coercion was implemented with 100% pure Cython code (with an eye
> for speed where it is needed),
The critical path for doing arithmetic between elements is 100% pure
Cython code. The path for discovering coercions has Python in it, but
until (if ever) all Parents are re-written in Cython this isn't going
to be fully optimal anyways. And it only happens once (compared to
the old model where it happened every single time a coercion was
needed).
Of course, much of the discover part could and should be optimized.
But right now we've already got enough on our plate trying to get the
changes we have pushed through. (Any help on this front would be
greatly appreciated.)
> I would be significantly less upset with it then I am now
> where people tell me that if I need more speed I'm out of luck.
That's not the message I'm trying to send--I'm sure there's room for
improvement (though I have a strong distaste for hard-coding special
cases in all over the place). I don't think "make everything Cython"
is going to solve the problem...
Nick Alexander
>> it's not as if performance is awful in the current model.
>>
> For the things you do. There is no code that pushes the coercion
> system anywhere near as much as symbolics in Sage does.
Please explain why the symbolics code is any more taxing on the
coercion system than the existing code.
Nick
Gary Furnish
Maybe the solution here is to run an_element on each parent and then
multiply them (and in fact I do this as a last case resort) and then
get the parent of the result, but this is a pretty bad way to do
things as it requires construction of elements during coercion.
Furthermore, this isn't even implemented in some classes, the default
is written in python, etc. Even if those issues were dealt with,
having to multiply two elements to figure out where the result is a
bad design.
> Of course, much of the discover part could and should be optimized.
> But right now we've already got enough on our plate trying to get the
> changes we have pushed through. (Any help on this front would be
> greatly appreciated.)
>
>> I would be significantly less upset with it then I am now
>> where people tell me that if I need more speed I'm out of luck.
>
> That's not the message I'm trying to send--I'm sure there's room for
> improvement (though I have a strong distaste for hard-coding special
> cases in all over the place). I don't think "make everything Cython"
> is going to solve the problem...
No but special cases/writing my own fast coercion code seems to work
significantly better then trying to use your system. This is
definitely a bad situation, because we shouldn't need another set of
coercion codes that deal with the cases where the main coercion code
is too slow. Either coercion has to be designed to be fast, or
symbolics is going to have to reach into the innards of the coercion
framework to make it possible to deal with quickly. It would be even
better if we could have symbolicring over RR or symbolicring over ZZ
or what not, but this is 100% impossible unless the coercion framework
is done 100% in cython, with special cases, with speed as a top
priority instead of beauty.
mabshoff
wrote:
<SNIP>
Hi,
> > If coercion was implemented with 100% pure Cython code (with an eye
> > for speed where it is needed),
>
> The critical path for doing arithmetic between elements is 100% pure
> Cython code. The path for discovering coercions has Python in it, but
> until (if ever) all Parents are re-written in Cython this isn't going
> to be fully optimal anyways. And it only happens once (compared to
> the old model where it happened every single time a coercion was
> needed).
>
> Of course, much of the discover part could and should be optimized.
> But right now we've already got enough on our plate trying to get the
> changes we have pushed through. (Any help on this front would be
> greatly appreciated.)
before the coercion rewrite since the coercion rewrite has been going
on much longer and there is substantial effort in there that is also
highly desired by the sage-combinat people for example. Once coercion
is in [hopefully around Dev1] it is desired to improve the existing
code [I am not qualified to speculate here what can and will be done
in that regard] and I am sure that Robert, Gary and everybody else
interested will find some time during Dev1 to sit down and look at the
code together to come up with a workable solution. Sage is about
incremental improvements and we should remember: "The perfect is the
enemy of the good" ;)
I know very, very well that Gary is very serious about performance and
that what he considers fast is way beyond the normal expectation. His
effort to achieve this is beyond what most people would consider
reasonable and I am very glad he is working on Symbolics ;) So let's
not get into an intense discussion here and now since that will not
resolve any problem.
<SNIP>
Cheers,
Michael
Bill Page
>
> On Tue, Jun 3, 2008 at 2:45 PM, Bill Page wrote:
>>
>> On Tue, Jun 3, 2008 at 4:48 PM, Robert Bradshaw wrote:
>>>
>>> On Jun 3, 2008, at 11:17 AM, Gary Furnish wrote:
>>> ...
>>>> I consider homsets to be a gigantic flaw in coercion that
>>>> absolutely have to be fixed for me to consider using more
>>>> of the coercion system in symbolics.
>>>
>>> Ironically, other people see it as a plus that coercion has
>>> been given a more categorical founding.
>>>
>>
>> Absolutely! :-)
>>
>> BTW, where can I read more about these categorical concepts
>> that are currently built-in or planned for Sage?
>>
>
> This is very relevant:
>
> http://wiki.sagemath.org/days7/coercion
>
Thanks for the reference. No mention of "homsets" here. :-( Only one
mention in /days7/coercion/todo ...
But reading this makes me feel a little, well ah, "light-headed". At
best it seems like something very hastily grafted-on to the design. Is
this part of Sage about which you would like comments and more
discussion? Or is there more information somewhere that I am missing?
I am afraid that there is not much here that category theorists are
likely to find interesting. I have many many questions, but mostly I
wonder what the overall intention of this construction really is in
Sage? Is it only relevant to coercion? How does it related to the
concept of "parent" - which seems equally ill-defined to me? What is
the relationship to inheritance in Python? Is the intention to give
all mathematical objects defined in Sage some categorical structure?
What about categories themselves as mathematical structures - e.g. a
category is a kind of directed graph with additional structure?
Regards,
Bill Page.
Robert Bradshaw
Ah, yes. We talked about this before, and I implemented the analyse
and explain functions which carefully avoid all Python:
http://cython.org/coercion/hgwebdir.cgi/sage-coerce/rev/1a5c8ccfd0df
>
>> Of course, much of the discover part could and should be optimized.
>> But right now we've already got enough on our plate trying to get the
>> changes we have pushed through. (Any help on this front would be
>> greatly appreciated.)
>>
>>> I would be significantly less upset with it then I am now
>>> where people tell me that if I need more speed I'm out of luck.
>>
>> That's not the message I'm trying to send--I'm sure there's room for
>> improvement (though I have a strong distaste for hard-coding special
>> cases in all over the place). I don't think "make everything Cython"
>> is going to solve the problem...
> No but special cases/writing my own fast coercion code seems to work
> significantly better then trying to use your system.
Writing a huge number of special cases is almost always going to be
faster, not doubt about it. The Sage coercion code needs to be
extremely generic as it handles all kinds of objects (and I'm excited
to have symbolics that respect this).
> This is
> definitely a bad situation, because we shouldn't need another set of
> coercion codes that deal with the cases where the main coercion code
> is too slow.
Yes.
> Either coercion has to be designed to be fast, or
> symbolics is going to have to reach into the innards of the coercion
> framework to make it possible to deal with quickly. It would be even
> better if we could have symbolicring over RR or symbolicring over ZZ
> or what not, but this is 100% impossible unless the coercion framework
> is done 100% in cython, with special cases, with speed as a top
> priority instead of beauty.
Would it be 95% possible if it's 95% written in Cython, with only a
5% performance hit :-).
The best balance I see is you can hard code ZZ/RR/... for speed if
you want (you're worried about virtual function calls anyways), and
then call off to coercion in the generic cases. You're going to have
to do this a bit anyways, as the coercion model doesn't handle stuff
like "where is the sin(x) if x is in R?" which is handled via the
normal OO sense based on the type of x (and may depend on x).
To help with collaboration, what you want out of the coercion model
is, given R and S (and an operation op), what will be the result of
op between elements of R and S, and you want this to be as fast as
possible (e.g. 100% Cython, no homsets, ...).
- Robert
Gary Furnish
> - Robert
>
>
> >
>
Robert Bradshaw
>
> On Tue, Jun 3, 2008 at 5:48 PM, William Stein wrote:
>>
>> On Tue, Jun 3, 2008 at 2:45 PM, Bill Page wrote:
>>>
>>> On Tue, Jun 3, 2008 at 4:48 PM, Robert Bradshaw wrote:
>>>>
>>>> On Jun 3, 2008, at 11:17 AM, Gary Furnish wrote:
>>>> ...
>>>>> I consider homsets to be a gigantic flaw in coercion that
>>>>> absolutely have to be fixed for me to consider using more
>>>>> of the coercion system in symbolics.
>>>>
>>>> Ironically, other people see it as a plus that coercion has
>>>> been given a more categorical founding.
>>>>
>>>
>>> Absolutely! :-)
>>>
>>> BTW, where can I read more about these categorical concepts
>>> that are currently built-in or planned for Sage?
>>>
>>
>> This is very relevant:
>>
>> http://wiki.sagemath.org/days7/coercion
>>
>
> Thanks for the reference. No mention of "homsets" here. :-( Only one
> mention in /days7/coercion/todo ...
>
> But reading this makes me feel a little, well ah, "light-headed". At
> best it seems like something very hastily grafted-on to the design. Is
> this part of Sage about which you would like comments and more
> discussion? Or is there more information somewhere that I am missing?
The discussion on that page is (unfortunately) best understood by
someone who already is fairly familiar with Sage. It will also
probably be re-written/cleaned up now that it has been implemented...
> I am afraid that there is not much here that category theorists are
> likely to find interesting. I have many many questions, but mostly I
> wonder what the overall intention of this construction really is in
> Sage? Is it only relevant to coercion?
The interesting coercion stuff (in my opinion) was all done last
June, and much of the work outlined here is trying to unify the API
and make things more consistent across all of Sage (though there is
more than just that).
> How does it related to the
> concept of "parent" - which seems equally ill-defined to me?
A Parent is an Object in the category of Sets, though in the context
of coercion one usually assumes one can some kind of arithmetic
(binary operations) on their elements.
> What is
> the relationship to inheritance in Python? Is the intention to give
> all mathematical objects defined in Sage some categorical structure?
> What about categories themselves as mathematical structures - e.g. a
> category is a kind of directed graph with additional structure?
A big push of this model is to divorce the relationship between
inheritance (which primarily has to do with implementation) and
mathematical categorization. This will allow categories to play a
larger role (though work in that area can be done after the merge as
it is not disruptive). For example, Z/nZ[x] is implemented the same
whether or not p is a prime, but the category it lives in (and hence
the methods one can do on it) vary. As another example, matrix spaces
are algebras iff they are square, but one doesn't want to have a
special class for square matrices over <insert optimized type here>.
Generic algorithms can be put on the categories such that if g is in
category C, and C has a method x, then g.x() will work.
- Robert
David Harvey
On Jun 3, 2008, at 10:04 PM, Robert Bradshaw wrote:
>> How does it related to the
>> concept of "parent" - which seems equally ill-defined to me?
>
> A Parent is an Object in the category of Sets,
huh? Don't you mean to say something more like "a parent is an object
of a concrete category", i.e. a category C with a faithful functor
f : C -> Set, such that the "elements" (as understood by Sage) of the
parent P are exactly the elements of f(P)?
> though in the context
> of coercion one usually assumes one can some kind of arithmetic
> (binary operations) on their elements.
david
Robert Bradshaw
> On Jun 3, 2008, at 10:04 PM, Robert Bradshaw wrote:
>
>>> How does it related to the
>>> concept of "parent" - which seems equally ill-defined to me?
>>
>> A Parent is an Object in the category of Sets,
>
> huh? Don't you mean to say something more like "a parent is an object
> of a concrete category", i.e. a category C with a faithful functor
> f : C -> Set, such that the "elements" (as understood by Sage) of the
> parent P are exactly the elements of f(P)?
Yes, that is a much more precise (and correct) explanation of what I
meant.
- Robert
John Cremona
> But right now we've already got enough on our plate trying to get the
> changes we have pushed through. (Any help on this front would be
> greatly appreciated.)
Robert,
Can you be more specific about how others can help with this?
John
> >
>
Bill Page
Thanks for the elaboration. I hope that this follow-up is not too far
off topic and I don't want to distract you (or any of the Sage
developers) from your main tasks. I wrote this mostly as notes and
questions to myself - but if you or anyone else have some time or the
inclination to correct or comment as time permits, that would be
great.
On Tue, Jun 3, 2008 at 10:04 PM, Robert Bradshaw wrote:
>
> On Jun 3, 2008, at 4:50 PM, Bill Page wrote:
>
>> How does it [categories in Sage] relate to the
>> concept of "parent" - which seems equally ill-defined to me?
>
> A Parent is an Object in the category of Sets, though in the context
> of coercion one usually assumes one can some kind of arithmetic
> (binary operations) on their elements.
> ...
> On Jun 3, 2008, at 7:11 PM, David Harvey wrote:
> ...
>> Don't you mean to say something more like "a parent is an object
>> of a concrete category", i.e. a category C with a faithful functor
>> f : C -> Set, such that the "elements" (as understood by Sage) of
>> the parent P are exactly the elements of f(P)?
Ok (and thanks also for the clarification, David). There are of course
two different uses of "object" here: 1) object of some category, 2)
Python object. All Python objects have a 'type', i.e. belong to some
Python class.
So in Sage 3.0.2 I observe for example:
sage: type(1)
<type 'sage.rings.integer.Integer'>
sage: parent(1)
Integer Ring
sage: type(parent(1))
<type 'sage.rings.integer_ring.IntegerRing_class'>
sage: category(parent(1))
Category of rings
sage: type(category(parent(1)))
<class 'sage.categories.category_types.Rings'>
and
sage: type(1.0)
<type 'sage.rings.real_mpfr.RealNumber'>
sage: parent(1.0)
Real Field with 53 bits of precision
sage: type(parent(1.0))
<type 'sage.rings.real_mpfr.RealField'>
sage: category(parent(1.0))
Category of fields
sage: type(category(parent(1.0)))
<class 'sage.categories.category_types.Fields'>
These seem consistent to me, albeit rather complex. However I am not
sure I understand the following:
sage: parent(IntegerRing())
<type 'sage.rings.integer_ring.IntegerRing_class'>
sage: parent(RealField())
<type 'sage.rings.real_mpfr.RealField'>
Could you explain why the parent of an object of some category is a type?
>
>> What is the relationship to inheritance in Python? Is the intention
>> to give all mathematical objects defined in Sage some categorical
>> structure? What about categories themselves as mathematical
>> structures - e.g. a category is a kind of directed graph with additional
>> structure?
>
> A big push of this model is to divorce the relationship between
> inheritance (which primarily has to do with implementation) and
> mathematical categorization.
At first blush it still seems strange to me to distinguish between the
parent of an object and the type of an object. The instances of a type
are often considered elements. Is the invention of 'parent' in Sage
due to some inherent limitation of Python classes? Do you really want
to say more abstractly that a certain type *represents* a certain
parent? E.g. Does 'sage.rings.real_mpfr.RealNumber' represent
'RealField()'? Can there be other representations of 'RealField()' in
Sage? How can I find out for example if 'sage.rings.integer.Integer'
represents 'IntegerRing()'?
From my experience with Axiom I am sympathetic to the separation of
implementation and specification but in Axiom "domains" are both
parents and types. A domain can inherit it's implementation from
another domain (via add). A domain can also serve as the
representation for some other domain (via Rep). For example
IntegerMod(5) is the representation of PrimeField(5). IntegerMod(5) in
turn is represented by SingleInteger, etc. In Axiom, at the deepest
level all types are ultimately provided by Lisp. These types are not
really Axiom domains and so I suppose the situation is analogous to
Python types and Sage parents.
Axiom categories on the other hand form a lattice. Domains belong to
categories by assertion or by virtue of the subdomain relationship. In
Axiom there is an operator named 'has' which provides the same
information as 'category' in Sage so that:
x has y
is just 'y == category(x)'. For example 'Float has Field' is equivalent to:
sage: Fields() == category(RealField())
True
Although in Axiom a domain may satisfy more than one category. Will
this be possible in Sage in the future?
What information does Sage actually have about a given category? For
example, what does Sage know about the relationship between Rings()
and Fields()? E.g. functors?
None of this has much to do with "category theory" as such. Categories
in Axiom are only vaguely related to categories in Mathematics. So I
suppose from your description that categories in Sage will play a
similar role but will be more strongly related to the mathematical
concept?
> This will allow categories to play a larger role (though work in that
> area can be done after the merge as it is not disruptive). For
> example, Z/nZ[x] is implemented the same whether or not p is
> a prime, but the category it lives in (and hence the methods one
> can do on it) vary.
Are you suggesting that Z/pZ should be considered a field just because
p is prime?
sage: category(Zmod(7))
Category of rings
sage: Zmod(7).is_field()
True
Perhaps an alternative would be to provide a coercion to some object
whose parent is in Fields()?
> As another example, matrix spaces are algebras iff they are square,
> but one doesn't want to have a special class for square matrices over
> <insert optimized type here>.
I am not sure what you mean by "matrix spaces" as such but the
category of matrices does make sense to me. See for example:
http://unapologetic.wordpress.com/2008/06/02/the-category-of-matrices-i
Maybe it is interesting because in this category the morphisms are
matrices and the objects are just the row and column dimensions -
non-negative integers. How would this fit into the category model in
Sage?
There is a monoidal sub-category of square matrices.
>Generic algorithms can be put on the categories such
> that if g is in category C, and C has a method x, then g.x() will
> work.
Do you mean 'category(parent(g))==C'? How will such a generic method x
of C operate on g? Do you have some example coding that implements
this idea?
Regards,
Bill Page.
David Harvey
I'm not an expert on these things any more, but I can tell you what
happened back in the day (when the sage version number was slightly
smaller):
Only "elements" have real parents. Since ZZ = IntegerRing() is not an
"element", it doesn't have a parent. If X does not have a parent,
then parent(X) returns the type of X instead. I'm not sure this is a
good idea, but there it is.
Of course then you're wondering why ZZ is not an element of the
category of the rings. That I do not know. But you can try:
sage: IntegerRing().parent()
to see that IntegerRing() simply doesn't have a parent() method,
which is why the global version parent(IntegerRing()) returns the
type of IntegerRing() instead.
One difficulty with ZZ being both a parent and an element is that
Cython does not support multiple inheritance.
BTW I should mention that the whole idea of "parents" vs "types" is
one of the main conceptual things inherited from Magma.
sorry gotta go cannot answer rest of email, hopefully someone else
can....
david
Robert Bradshaw
Sure. Currently we've made all the underlying changes, and we're in
the process of fixing doctests (mostly they're things like typos/
unimplemented functionality rather than deep mathematical
peculiarities). The way to get started is:
- Download and build Sage 2.10.1 http://sagemath.org/dist/src/
sage-2.10.1.tar
- Install the latest Cython http://sage.math.washington.edu/home/
robertwb/cython/cython-0.9.6.14.p1.spkg
- Pull from the coercion repo and build http://cython.org/coercion/
hgwebdir.cgi/sage-coerce/
- Choose a (set of) file(s) from the todo list at http://
wiki.sagemath.org/days7/coercion/todo and try to get all doctests to
pass there. We're using the wiki to keep track of who's doing what.
Some hints:
- __cmp__ has been changed to _cmp_ in SageObject, and __cmp__ no
longer works due to a default __richcmp__ (that's just how Python
works). In retrospect, this probably should not have been done on top
of everything else, 'cause it's somewhat independent of coercion and
the source of seemingly most failures.
- A useful command is coercion_traceback(). This will give the
traceback of all exceptions caught during the coercion process.
Hopefully this is enough to get started, and feel free to ask questions.
- Robert
Robert Bradshaw
> Robert,
>
> Thanks for the elaboration. I hope that this follow-up is not too far
> off topic and I don't want to distract you (or any of the Sage
> developers) from your main tasks. I wrote this mostly as notes and
> questions to myself - but if you or anyone else have some time or the
> inclination to correct or comment as time permits, that would be
> great.
I'm sure a lot of other people have similar questions, and it will be
good to put it down in writing.
David's explanation of this is right on. We need parent() to work in
some sensible way on non-Elements (e.g. Python ints, objects from
outside Sage) to be able do some kind of coercion reasoning on them.
Python is a dynamically typed language, ever object has a type.
>
>>
>>> What is the relationship to inheritance in Python? Is the intention
>>> to give all mathematical objects defined in Sage some categorical
>>> structure? What about categories themselves as mathematical
>>> structures - e.g. a category is a kind of directed graph with
>>> additional
>>> structure?
>>
>> A big push of this model is to divorce the relationship between
>> inheritance (which primarily has to do with implementation) and
>> mathematical categorization.
>
> At first blush it still seems strange to me to distinguish between the
> parent of an object and the type of an object. The instances of a type
> are often considered elements. Is the invention of 'parent' in Sage
> due to some inherent limitation of Python classes? Do you really want
> to say more abstractly that a certain type *represents* a certain
> parent? E.g. Does 'sage.rings.real_mpfr.RealNumber' represent
> 'RealField()'? Can there be other representations of 'RealField()' in
> Sage? How can I find out for example if 'sage.rings.integer.Integer'
> represents 'IntegerRing()'?
One should think of the type of an object as what specifies its
implementation (including its internal representation). In this sense
it is like a class in Java, C++, or any other programming language.
(In Python there is a subtle (and mostly historical) difference
between types and classes depending on whether or not it's written in
C, but for the most part they can be used interchangeably). Python
has multiple inheritance, but Cython does not (the C format of a
compiled Python class is a C struct, so this is not so much a Cython
limitation as a Python/C API limitation.)
Most of the time there is a 1-to-1 correspondence between Parents and
the types of its Elements. For example, and element of the Parent RR
(= RealField(53) is represented by sage.rings.real_mpfr.RealNumber.
If I do RR(3) I get back an object of type
sage.rings.real_mpfr.RealNumber representing the floating point
number 3 to 53 bits of precision, and it's Parent is RR. RealField
(100) is a new Parent, whose elements are again of type
sage.rings.real_mpfr.RealNumber but to 100 bits of precision.
Sometimes, however, the same type is used for multiple parents. There
are several integer mod types (depending on whether or not the
modulus fits in a single word) and they are used for both generic Z/
nZ and prime-order finite fields. Thus we have
sage: [type(mod(-1, 100^k)) for k in [2..5]]
[<type 'sage.rings.integer_mod.IntegerMod_int'>,
<type 'sage.rings.integer_mod.IntegerMod_int64'>,
<type 'sage.rings.integer_mod.IntegerMod_int64'>,
<type 'sage.rings.integer_mod.IntegerMod_gmp'>]
sage: [parent(mod(-1, 100^k)) for k in [2..5]]
[Ring of integers modulo 10000,
Ring of integers modulo 1000000,
Ring of integers modulo 100000000,
Ring of integers modulo 10000000000]
On the other hand, some parents may have elements of several types.
The "SymbolicRing" is one such example.
sage: [type(a) for a in v]
[<class 'sage.calculus.calculus.SymbolicConstant'>,
<class 'sage.functions.constants.Pi'>,
<class 'sage.calculus.calculus.SymbolicArithmetic'>,
<class 'sage.calculus.calculus.SymbolicComposition'>]
sage: [parent(a) for a in v]
[Symbolic Ring, Symbolic Ring, Symbolic Ring, Symbolic Ring]
> From my experience with Axiom I am sympathetic to the separation of
> implementation and specification but in Axiom "domains" are both
> parents and types. A domain can inherit it's implementation from
> another domain (via add). A domain can also serve as the
> representation for some other domain (via Rep). For example
> IntegerMod(5) is the representation of PrimeField(5). IntegerMod(5) in
> turn is represented by SingleInteger, etc. In Axiom, at the deepest
> level all types are ultimately provided by Lisp. These types are not
> really Axiom domains and so I suppose the situation is analogous to
> Python types and Sage parents.
This sounds like the correct analogy. In Sage as well elements of
prime finite fields are the same type as elements of Z/nZ for
composite n, but they have different parents.
> Axiom categories on the other hand form a lattice. Domains belong to
> categories by assertion or by virtue of the subdomain relationship. In
> Axiom there is an operator named 'has' which provides the same
> information as 'category' in Sage so that:
>
> x has y
>
> is just 'y == category(x)'. For example 'Float has Field' is
> equivalent to:
>
> sage: Fields() == category(RealField())
> True
>
> Although in Axiom a domain may satisfy more than one category. Will
> this be possible in Sage in the future?
Categories in Sage aren't as fully developed as they could be, but
they do form a lattice.
sage: C = category(FiniteField(13)); C
Category of fields
sage: C.is_subcategory(Rings())
True
One of the changes in this latest coercion push is to allow Parents
to belong to several categories.
> What information does Sage actually have about a given category? For
> example, what does Sage know about the relationship between Rings()
> and Fields()? E.g. functors?
Only the bare minimum at this point.
> None of this has much to do with "category theory" as such. Categories
> in Axiom are only vaguely related to categories in Mathematics. So I
> suppose from your description that categories in Sage will play a
> similar role but will be more strongly related to the mathematical
> concept?
Yes. This will also be a place to put generic implementations of
algorithms. For example, one has the category of Euclidian domains,
where one can define a generic gcd algorithm that will then get
attached (as a method) to every element of that category.
>> This will allow categories to play a larger role (though work in that
>> area can be done after the merge as it is not disruptive). For
>> example, Z/nZ[x] is implemented the same whether or not p is
>> a prime, but the category it lives in (and hence the methods one
>> can do on it) vary.
>
> Are you suggesting that Z/pZ should be considered a field just because
> p is prime?
No, it might be expensive to test for primality.
> sage: category(Zmod(7))
> Category of rings
>
> sage: Zmod(7).is_field()
> True
>
> Perhaps an alternative would be to provide a coercion to some object
> whose parent is in Fields()?
The is_field method tests whether or not this specific instance
happens to be a field, which is true in this case. But we do provide
a coercion
sage: F = GF(7)
sage: F.coerce_map_from(Zmod(7))
Natural morphism:
From: Ring of integers modulo 7
To: Finite Field of size 7
>
>> As another example, matrix spaces are algebras iff they are square,
>> but one doesn't want to have a special class for square matrices over
>> <insert optimized type here>.
>
> I am not sure what you mean by "matrix spaces" as such but the
> category of matrices does make sense to me. See for example:
>
> http://unapologetic.wordpress.com/2008/06/02/the-category-of-
> matrices-i
>
> Maybe it is interesting because in this category the morphisms are
> matrices and the objects are just the row and column dimensions -
> non-negative integers. How would this fit into the category model in
> Sage?
>
> There is a monoidal sub-category of square matrices.
A MatrixSpace is a parent
sage: M = matrix(QQ, 2); M
[0 0]
[0 0]
sage: parent(M)
Full MatrixSpace of 2 by 2 dense matrices over Rational Field
It is a module (via scalar multiplication) over its basering (in the
case above, the integers) and so (should) lie in the category of Q-
modules. In the case of square matrices, it is actually an algebra
over its basering, and its category is the category of Q-algebras.
>> Generic algorithms can be put on the categories such
>> that if g is in category C, and C has a method x, then g.x() will
>> work.
>
> Do you mean 'category(parent(g))==C'?
Yes.
> How will such a generic method x
> of C operate on g? Do you have some example coding that implements
> this idea?
Not at the moment. Say a and b have parent P in category C. One types
a.gcd(b)
Then, a does not have a method "gcd" it will try
category(parent(a)).element_gcd(a,b)
which will be a generic implementation of the euclidean algorithm
(assuming C is a category supporting this).
- Robert
Bill Page
>
> On Jun 4, 2008, at 4:07 AM, Bill Page wrote:
>>
>> These seem consistent to me, albeit rather complex. However I am
>> not sure I understand the following:
>>
>> sage: parent(IntegerRing())
>> <type 'sage.rings.integer_ring.IntegerRing_class'>
>> sage: parent(RealField())
>> <type 'sage.rings.real_mpfr.RealField'>
>>
>> Could you explain why the parent of an object of some category is a
>> type?
>
> David's explanation of this is right on. We need parent() to work in
> some sensible way on non-Elements (e.g. Python ints, objects from
> outside Sage) to be able do some kind of coercion reasoning on them.
> Python is a dynamically typed language, ever object has a type.
>
So is this essentially an admission that the concept of 'parent' is
superfluous and could have in fact been implemented just at Python
'type'? I am not sure that I would actually advocate this now, but for
argument's sake: Why not do it this way all of the time? Python
classes all end up as (potential) parents. I think that this more or
less is what at least some Magma developers had in mind. E.g. We read
in the preface to the Magma HTML Help Document:
http://magma.maths.usyd.edu.au/magma/htmlhelp/preface.htm
"Every object created during the course of a computation is associated
with a unique parent algebraic structure. The type of an object is
then simply its parent structure."
> ...
> One should think of the type of an object as what specifies its
> implementation (including its internal representation). In this sense
> it is like a class in Java, C++, or any other programming language.
> (In Python there is a subtle (and mostly historical) difference
> between types and classes depending on whether or not it's written
> in C, but for the most part they can be used interchangeably).
> Python has multiple inheritance, but Cython does not (the C
> format of a compiled Python class is a C struct, so this is not so
> much a Cython limitation as a Python/C API limitation.)
>
The lack of multiple inheritance in Cython seems to be a major design
obstacle. Although a C struct is not a real class, it is not clear to
me why this is not adequate in compiled code. Wouldn't it be
sufficient to merge (flatten) the components of the multiple ancesters
into one struct? If not (e.g. if it is necessary to retain some
knowledge of their origins), why not provide nested structs? But I
suppose that this is quite off-topic...
> Most of the time there is a 1-to-1 correspondence between Parents
> and the types of its Elements. For example, and element of the Parent
> RR (= RealField(53) is represented by sage.rings.real_mpfr.RealNumber.
> If I do RR(3) I get back an object of type sage.rings.real_mpfr.RealNumber
> representing the floating point number 3 to 53 bits of precision, and it's
> Parent is RR. RealField (100) is a new Parent, whose elements are
> again of type sage.rings.real_mpfr.RealNumber but to 100 bits of
> precision.
>
> Sometimes, however, the same type is used for multiple parents.
> There are several integer mod types (depending on whether or not
> the modulus fits in a single word) and they are used for both generic
> Z/nZ and prime-order finite fields. Thus we have
>
> sage: [type(mod(-1, 100^k)) for k in [2..5]]
> [<type 'sage.rings.integer_mod.IntegerMod_int'>,
> <type 'sage.rings.integer_mod.IntegerMod_int64'>,
> <type 'sage.rings.integer_mod.IntegerMod_int64'>,
> <type 'sage.rings.integer_mod.IntegerMod_gmp'>]
>
> sage: [parent(mod(-1, 100^k)) for k in [2..5]]
> [Ring of integers modulo 10000,
> Ring of integers modulo 1000000,
> Ring of integers modulo 100000000,
> Ring of integers modulo 10000000000]
>
Hmmm... what I think I see here is multiple Python types (_int,
_int64, _gmp) being used for the *same* parent, i.e. "Ring of integers
modulo n" (for some parameter n). Simple inheritance should make this
easy, no?
> On the other hand, some parents may have elements of several
> types. The "SymbolicRing" is one such example.
>
> sage: [type(a) for a in v]
> [<class 'sage.calculus.calculus.SymbolicConstant'>,
> <class 'sage.functions.constants.Pi'>,
> <class 'sage.calculus.calculus.SymbolicArithmetic'>,
> <class 'sage.calculus.calculus.SymbolicComposition'>]
>
> sage: [parent(a) for a in v]
> [Symbolic Ring, Symbolic Ring, Symbolic Ring, Symbolic Ring]
>
"SymbolicRing" is another thing that worries me a little. Exactly what
(mathematically) is a symbolic ring? E.g. Why is it a ring and not a
field? Operations like 'sin(sin(sin(x))).parent()' doesn't look much
like a ring to me. Probably most other computer algebra systems would
just call this an "expression tree" or AST or something like that. Why
do we need different types (sub-classes)?
>> From my experience with Axiom I am sympathetic to the separation of
>> implementation and specification but in Axiom "domains" are both
>> parents and types. A domain can inherit it's implementation from
>> another domain (via add). A domain can also serve as the
>> representation for some other domain (via Rep). For example
>> IntegerMod(5) is the representation of PrimeField(5). IntegerMod(5) in
>> turn is represented by SingleInteger, etc.
> ...
> This sounds like the correct analogy. In Sage as well elements of
> prime finite fields are the same type as elements of Z/nZ for
> composite n, but they have different parents.
>
To the best of my knowledge, the formal concept of "representation" in
a computer algebra systems is unique to Axiom. Each domain has a
specified 'Rep'consisting of some domain expression. E.g.
Rep == Integer
or
Rep == Union(Record(car:%, cdr:%), "nil")
% stands for "this domain", so this definition amounts to some
recursive structure. Of course recursively defined classes are also
possible in Python but they require much more discipline on the part
of the programmer.
Now there are two operators 'rep' and 'per' that provide coercion to
and from Rep. So we might say for example that some operation * in
this domain is implemented by * from Rep:
x * y == per( rep(x) * rep(y) )
The trouble is: I don't see anything here that could be called a
"parent". Or perhaps, they are all "parents" except for those initial
domains for which we never explicitly specify a Rep.
>
> Categories in Sage aren't as fully developed as they could be,
> but they do form a lattice.
>
> sage: C = category(FiniteField(13)); C
> Category of fields
> sage: C.is_subcategory(Rings())
> True
>
Ah, I think that is very useful. Does the "Category of fields" inherit
any generic operations from "Rings()"?
> One of the changes in this latest coercion push is to allow Parents
> to belong to several categories.
>
That sounds good. So does this mean that the 'category' method now
returns a list?
> ... This will also be a place to put generic implementations of
> algorithms. For example, one has the category of Euclidian
> domains, where one can define a generic gcd algorithm that will
> then get attached (as a method) to every element of that category.
> ...
>>
>> sage: Zmod(7).is_field()
>> True
>>
>> Perhaps an alternative would be to provide a coercion to some
>> object whose parent is in Fields()?
>
> The is_field method tests whether or not this specific instance
> happens to be a field, which is true in this case. But we do provide
> a coercion
>
> sage: F = GF(7)
> sage: F.coerce_map_from(Zmod(7))
>
> Natural morphism:
> From: Ring of integers modulo 7
> To: Finite Field of size 7
>
Excellent. But, hmm... this is supposed to be a morphism in some
category. What category?
sage: category(Zmod(7))
Category of rings
sage: category(FiniteField(7))
Category of fields
sage: category(FiniteField(7)).is_subcategory(category(Zmod(7)))
True
------
Ok, good. We just need the injection.
>> ...
>> I am not sure what you mean by "matrix spaces" as such but
>> the category of matrices does make sense to me. See for
>> example:
>>
>> http://unapologetic.wordpress.com/2008/06/02/the-category-of-matrices-i
>>
>> Maybe it is interesting because in this category the morphisms
>> are matrices and the objects are just the row and column
>> dimensions - non-negative integers. How would this fit into the
>> category model in Sage?
>>...
>
> A MatrixSpace is a parent
>
> sage: M = matrix(QQ, 2); M
> [0 0]
> [0 0]
> sage: parent(M)
> Full MatrixSpace of 2 by 2 dense matrices over Rational Field
>
> It is a module (via scalar multiplication) over its basering (in the
> case above, the integers) and so (should) lie in the category of
> Q-modules.
Ok, it is an object in some category. Maybe a good name would be the
category "VectorSpace"? (Q-modules is ok). What makes them a
MatrixSpace?
> In the case of square matrices, it is actually an algebra
> over its basering, and its category is the category of Q-algebras.
>
Anyway, I suppose this taxonomy is not such a big deal. One only needs
to learn the language - of course it has turned out in Axiom that this
is not that simple a problem - hence interactive tools for querying
and navigating the domain and category hierarchies (hyperdoc).
> ...
Thanks for continuing the discussion.
Regards,
Bill Page.
William Stein
>
> On Wed, Jun 4, 2008 at 1:54 PM, Robert Bradshaw wrote:
>>
>> On Jun 4, 2008, at 4:07 AM, Bill Page wrote:
>> ...
>>>
>>> These seem consistent to me, albeit rather complex. However I am
>>> not sure I understand the following:
>>>
>>> sage: parent(IntegerRing())
>>> <type 'sage.rings.integer_ring.IntegerRing_class'>
>>> sage: parent(RealField())
>>> <type 'sage.rings.real_mpfr.RealField'>
>>>
>>> Could you explain why the parent of an object of some category is a
>>> type?
>>
>> David's explanation of this is right on. We need parent() to work in
>> some sensible way on non-Elements (e.g. Python ints, objects from
>> outside Sage) to be able do some kind of coercion reasoning on them.
>> Python is a dynamically typed language, ever object has a type.
>>
>
> So is this essentially an admission that the concept of 'parent' is
> superfluous and could have in fact been implemented just at Python
> 'type'?
That is not the case, except that of course *anything* that can be
implemented, can be implemented in any turing complete language.
> I am not sure that I would actually advocate this now, but for
> argument's sake: Why not do it this way all of the time? Python
> classes all end up as (potential) parents. I think that this more or
> less is what at least some Magma developers had in mind. E.g. We read
> in the preface to the Magma HTML Help Document:
>
> http://magma.maths.usyd.edu.au/magma/htmlhelp/preface.htm
>
> "Every object created during the course of a computation is associated
> with a unique parent algebraic structure. The type of an object is
> then simply its parent structure."
>
You took that quote out of context. The full quote is: "The
theoretical basis for the design
of Magma is founded on the concepts and methodology of modern algebra.
The central notion is that of an algebraic structure. Every object
created during
the course of a computation is associated with a unique parent
algebraic structure.
The type of an object is then simply its parent structure."
The word "type" is being used in a formal sense as defined in the excellent
paper "The Magma Algebra System I: The User Language" by Bosma, Cannon,
and Playoust.
There are infinitely many different parents, one of each value
of n. For a given n, the ring Z/nZ is a specific mathematical
object, which is itself for many many reasons worth modeling
in a computer. If nothing else, it is the place on which to hang
functions like "multiplicative_generator()". Any computer algebra
system that doesn't have first class actual objects that represent
basic mathematical structures such as rings, fields, etc., feels
extremely barren if one has used Magma for a sufficient amount
of time.
Specific rings, fields, vector spaces, etc., are all basic objects of
mathematics, that are just as important to fully implement and model
in a CAS as polynomials, rational numbers, vectors, etc. I firmly
believe Magma was the first ever system to really get this right,
and Sage merely copies this. I really greatly appreciate what
Cannon et al. did...
>
>> On the other hand, some parents may have elements of several
>> types. The "SymbolicRing" is one such example.
>>
>> sage: [type(a) for a in v]
>> [<class 'sage.calculus.calculus.SymbolicConstant'>,
>> <class 'sage.functions.constants.Pi'>,
>> <class 'sage.calculus.calculus.SymbolicArithmetic'>,
>> <class 'sage.calculus.calculus.SymbolicComposition'>]
>>
>> sage: [parent(a) for a in v]
>> [Symbolic Ring, Symbolic Ring, Symbolic Ring, Symbolic Ring]
>>
>
> "SymbolicRing" is another thing that worries me a little. Exactly what
> (mathematically) is a symbolic ring? E.g. Why is it a ring and not a
> field? Operations like 'sin(sin(sin(x))).parent()' doesn't look much
> like a ring to me. Probably most other computer algebra systems would
> just call this an "expression tree" or AST or something like that. Why
> do we need different types (sub-classes)?
To me, the SymbolicRing is the mathematical structure in which
calculus teachers, engineers, etc., are kept happy. In Sage
that was 100% its purpose. Ask Marshall Hampton, who is teaching
calculus right now using Sage, about what he needs to make
teaching calculus easy, and the answer is "The Symbolic Ring",
plus graphics and interact.
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Bill Page
>
> On Wed, Jun 4, 2008 at 7:35 PM, Bill Page wrote:
>>
>> On Wed, Jun 4, 2008 at 1:54 PM, Robert Bradshaw wrote:
>>>
>>> David's explanation of this is right on. We need parent() to work
>>> in some sensible way on non-Elements (e.g. Python ints, objects
>>> from outside Sage) to be able do some kind of coercion reasoning
>>> on them. Python is a dynamically typed language, ever object has
>>> a type.
>>
>> So is this essentially an admission that the concept of 'parent' is
>> Python 'type'?
>
> That is not the case, except that of course *anything* that can
> be implemented, can be implemented in any turing complete
> language.
>
I did not mean to pose this question in a trivial way. I meant to
suggest that the concept of type in the Python programming language,
i.e. classes, could directly implement the concept of "parent" as used
in Sage without abusing either concept at all. Can you suggest an
example that demonstrates that this is not the case?
> ...
> There are infinitely many different parents, one of each value
> of n. For a given n, the ring Z/nZ is a specific mathematical
> object, which is itself for many many reasons worth modeling
> in a computer. If nothing else, it is the place on which to hang
> functions like "multiplicative_generator()".
Python classes can also take parameters.
> Any computer algebra system that doesn't have first class
> actual objects that represent basic mathematical structures
> such as rings, fields, etc., feels extremely barren if one has
> used Magma for a sufficient amount of time.
>
Sure. Rings and fields are categories. Magma was not the first system
to represent these mathematical structures in a computer.
> Specific rings, fields, vector spaces, etc., are all basic objects of
> mathematics, that are just as important to fully implement and
> model in a CAS as polynomials, rational numbers, vectors, etc.
> I firmly believe Magma was the first ever system to really get this
> right, and Sage merely copies this. I really greatly appreciate
> what Cannon et al. did...
>
I am sure Magma is an interesting system - everything I have been
reading about it confirms that - but I do not think it was the first.
I am also not yet convinced that it has implemented these general
mathematically concepts in any substantially better way then, for
example Axiom. Some things it does seem peculiar and ad hoc to me.
These are the things that standout in my mind when I think about in
what ways Sage resembles Magma. Of course there are some specific
computations that Magma does much more efficiently than any other
system. But that is not the point. We are talking here about general
design issues.
> ...
>> "SymbolicRing" is another thing that worries me a little. Exactly
>> what (mathematically) is a symbolic ring? E.g. Why is it a ring
>> and not a field? Operations like 'sin(sin(sin(x))).parent()' doesn't
>> look much like a ring to me. Probably most other computer algebra
>> systems would just call this an "expression tree" or AST or
>> something like that.
> ...
> To me, the SymbolicRing is the mathematical structure in which
> calculus teachers, engineers, etc., are kept happy. In Sage
> that was 100% its purpose. Ask Marshall Hampton, who is teaching
> calculus right now using Sage, about what he needs to make
> teaching calculus easy, and the answer is "The Symbolic Ring",
> plus graphics and interact.
>
I don't doubt that doing symbolic mathematics by computer is very
useful for teaching calculus - after all that is what is normally
meant by the term "calculus", i.e. a method of calculation based on
the symbolic manipulation of expressions. The name "Symbolic Ring"
however makes it sound more "algebraic",i.e. concerned with structure,
relation and quantity. So I was looking for a more formal mathematical
definition of this sort. I did not expect it was a name intended just
to keep certain people happy. :-(
> ...
Regards,
Bill Page.
William Stein
>
> On Wed, Jun 4, 2008 at 11:06 PM, William Stein wrote:
>>
>> On Wed, Jun 4, 2008 at 7:35 PM, Bill Page wrote:
>>>
>>> On Wed, Jun 4, 2008 at 1:54 PM, Robert Bradshaw wrote:
>>>>
>>>> David's explanation of this is right on. We need parent() to work
>>>> in some sensible way on non-Elements (e.g. Python ints, objects
>>>> from outside Sage) to be able do some kind of coercion reasoning
>>>> on them. Python is a dynamically typed language, ever object has
>>>> a type.
>>>
>>> So is this essentially an admission that the concept of 'parent' is
>>> superfluous and could have in fact been implemented just as
>>> Python 'type'?
>>
>> That is not the case, except that of course *anything* that can
>> be implemented, can be implemented in any turing complete
>> language.
>>
>
> I did not mean to pose this question in a trivial way. I meant to
> suggest that the concept of type in the Python programming language,
> i.e. classes, could directly implement the concept of "parent" as used
> in Sage without abusing either concept at all. Can you suggest an
> example that demonstrates that this is not the case?
I think Robert's example of the integers modulo n is an excellent example
of precisely this. I'm confused about why you don't see this.
>
>> ...
>> There are infinitely many different parents, one of each value
>> of n. For a given n, the ring Z/nZ is a specific mathematical
>> object, which is itself for many many reasons worth modeling
>> in a computer. If nothing else, it is the place on which to hang
>> functions like "multiplicative_generator()".
>
> Python classes can also take parameters.
I didn't know that. I thought the only way to create a Python class
is for the Python interpreter to execute Python code that looks like this:
class Foo(...):
...
That makes a new class called Foo. How are you going to make, at
runtime, new classes for each of Z/nZ say, for 1 <= n < 10^5, i.e.,
something like this in Sage:
v = [Integers(n) for n in range(1,10^5)]
I do not think it is possible to concisely create a hundred thousand
separate Python classes like that.
>> Any computer algebra system that doesn't have first class
>> actual objects that represent basic mathematical structures
>> such as rings, fields, etc., feels extremely barren if one has
>> used Magma for a sufficient amount of time.
>>
>
> Sure. Rings and fields are categories.
I meant "Rings" and "fields" not as categories but specific rings
and specified fields as mathematical objects in their own right.
In Magma and Sage there is virtually no difference between
"a polynomial ring over the rational numbers" and a specific
polynomial over the rationals -- both are first class instances
in exactly the same way.
> Magma was not the first system
> to represent these mathematical structures in a computer.
I believe Magma was the first to very systematically represent
the objects of the categories as first class objects systematically
throughout the system. I know of no other system that does this
(except Sage).
>
>> Specific rings, fields, vector spaces, etc., are all basic objects of
>> mathematics, that are just as important to fully implement and
>> model in a CAS as polynomials, rational numbers, vectors, etc.
>> I firmly believe Magma was the first ever system to really get this
>> right, and Sage merely copies this. I really greatly appreciate
>> what Cannon et al. did...
>>
>
> I am sure Magma is an interesting system - everything I have been
> reading about it confirms that - but I do not think it was the first.
> I am also not yet convinced that it has implemented these general
> mathematically concepts in any substantially better way then, for
> example Axiom.
I think Magma is brilliant. Disclaimer: I used to be one of the
biggest Magma zealots on the planet, so take that with a grain of salt!
It is also not my intention here to make any negative comments
at all about Axiom (say) -- only positive comments about Magma.
> Some things it does seem peculiar and ad hoc to me.
> These are the things that standout in my mind when I think about in
> what ways Sage resembles Magma. Of course there are some specific
> computations that Magma does much more efficiently than any other
> system. But that is not the point. We are talking here about general
> design issues.
Yet efficiency and design are closely linked. When I first started using
Magma after years and years of C++, I was amazed how they were able
to unify so many ideas and bring such a vast amount of quality code
together in a way that worked pretty well. I soon started writing, in the
Magma interpreter, vastly more efficient code than I could ever write in C++,
and code that accomplished a lot more. Much of this was for me directly
a result of the excellent design ideas that went into Magma. This was
a real eye opener for me.
>> ...
>>> "SymbolicRing" is another thing that worries me a little. Exactly
>>> what (mathematically) is a symbolic ring? E.g. Why is it a ring
>>> and not a field? Operations like 'sin(sin(sin(x))).parent()' doesn't
>>> look much like a ring to me. Probably most other computer algebra
>>> systems would just call this an "expression tree" or AST or
>>> something like that.
>> ...
>> To me, the SymbolicRing is the mathematical structure in which
>> calculus teachers, engineers, etc., are kept happy. In Sage
>> that was 100% its purpose. Ask Marshall Hampton, who is teaching
>> calculus right now using Sage, about what he needs to make
>> teaching calculus easy, and the answer is "The Symbolic Ring",
>> plus graphics and interact.
>>
>
> I don't doubt that doing symbolic mathematics by computer is very
> useful for teaching calculus - after all that is what is normally
> meant by the term "calculus", i.e. a method of calculation based on
> the symbolic manipulation of expressions. The name "Symbolic Ring"
> however makes it sound more "algebraic",i.e. concerned with structure,
> relation and quantity. So I was looking for a more formal mathematical
> definition of this sort. I did not expect it was a name intended just
> to keep certain people happy. :-(
It's "symbolic" in the sense of "symbolic calculus" as opposed to
"numerical" analysis. The name wasn't chosen to keep people happy,
but the object itself was created to do so. I think it might have been
called the SymbolicCalculusRing or SymbolicExpressionRing
at one point but that was too long.
>
>> ...
Robert Bradshaw
> On Wed, Jun 4, 2008 at 1:54 PM, Robert Bradshaw wrote:
>>
>> On Jun 4, 2008, at 4:07 AM, Bill Page wrote:
>> ...
>>>
>>> These seem consistent to me, albeit rather complex. However I am
>>> not sure I understand the following:
>>>
>>> sage: parent(IntegerRing())
>>> <type 'sage.rings.integer_ring.IntegerRing_class'>
>>> sage: parent(RealField())
>>> <type 'sage.rings.real_mpfr.RealField'>
>>>
>>> Could you explain why the parent of an object of some category is a
>>> type?
>>
>> David's explanation of this is right on. We need parent() to work in
>> some sensible way on non-Elements (e.g. Python ints, objects from
>> outside Sage) to be able do some kind of coercion reasoning on them.
>> Python is a dynamically typed language, ever object has a type.
>
> So is this essentially an admission that the concept of 'parent' is
> superfluous and could have in fact been implemented just at Python
> 'type'?
Not at all. It is the realization that not all Python objects will
have a Parent (e.g. the builtin objects, objects from other
libraries, Sage Parents, etc.) but they will have a type, which is a
much weaker notion but can be reasoned with a bit (e.g. the builtin
Python 'int' type coerces nicely into the IntegerRing).
> I am not sure that I would actually advocate this now, but for
> argument's sake: Why not do it this way all of the time? Python
> classes all end up as (potential) parents.
Types are all about the implementations of things, they synonymous
with the "classes" of Object Oriented programming, and are
insufficient (and the wrong vehicle) to carry deeper mathematical
structure. For example, an element of Z/5Z and an element of Z/7Z
have the same type (i.e. they're instances of the same class, with
the same methods, internal representation, etc), but not the same
parent (which is data).
> I think that this more or
> less is what at least some Magma developers had in mind. E.g. We read
> in the preface to the Magma HTML Help Document:
>
> http://magma.maths.usyd.edu.au/magma/htmlhelp/preface.htm
>
> "Every object created during the course of a computation is associated
> with a unique parent algebraic structure. The type of an object is
> then simply its parent structure."
William Stein adequately addressed this point.
>
>> ...
>> One should think of the type of an object as what specifies its
>> implementation (including its internal representation). In this sense
>> it is like a class in Java, C++, or any other programming language.
>> (In Python there is a subtle (and mostly historical) difference
>> between types and classes depending on whether or not it's written
>> in C, but for the most part they can be used interchangeably).
>> Python has multiple inheritance, but Cython does not (the C
>> format of a compiled Python class is a C struct, so this is not so
>> much a Cython limitation as a Python/C API limitation.)
>>
>
> The lack of multiple inheritance in Cython seems to be a major design
> obstacle. Although a C struct is not a real class, it is not clear to
> me why this is not adequate in compiled code. Wouldn't it be
> sufficient to merge (flatten) the components of the multiple ancesters
> into one struct? If not (e.g. if it is necessary to retain some
> knowledge of their origins), why not provide nested structs? But I
> suppose that this is quite off-topic...
The exact details of why this isn't as easy as it seems are quite
technical, and off-topic. The short answer is that with a major re-
design of Cython multiple inheritance could be supported at the cost
of an extra indirection (even on non-multiply inherited types) on
ever C member/method lookup. (This is done in C++, but it would be
even more complicated to deal with the way Python handles extension
types). One of the primary reasons for Cython is speed, which is not
worth sacrificing for this feature.
As stated before, there are 4 distinct parents here (incidentally
instances of the "IntegerMod_Ring" class).
> Simple inheritance should make this easy, no?
What can be inherited is inherited from a common IntegerMod_abstract
class. Subclassing is used where methods can be overridden to be
faster (e.g. arithmetic can be made much faster if the data is a C
int, rather than an arbitrary-length mpz_t integer type).
>
>> On the other hand, some parents may have elements of several
>> types. The "SymbolicRing" is one such example.
sage: v = [SR(1), pi, x+1, sin(x)] # I forgot this important input
line in my example.
>> sage: [type(a) for a in v]
>> [<class 'sage.calculus.calculus.SymbolicConstant'>,
>> <class 'sage.functions.constants.Pi'>,
>> <class 'sage.calculus.calculus.SymbolicArithmetic'>,
>> <class 'sage.calculus.calculus.SymbolicComposition'>]
>>
>> sage: [parent(a) for a in v]
>> [Symbolic Ring, Symbolic Ring, Symbolic Ring, Symbolic Ring]
>>
>
> "SymbolicRing" is another thing that worries me a little. Exactly what
> (mathematically) is a symbolic ring? E.g. Why is it a ring and not a
> field? Operations like 'sin(sin(sin(x))).parent()' doesn't look much
> like a ring to me. Probably most other computer algebra systems would
> just call this an "expression tree" or AST or something like that.
Most symbolic systems (Magma and Axiom excluded) have very little
notion of algebras, rings, modules, etc.--virtually everything is an
expression or list of expressions. SymbolicRing is a way to fit these
abstract "expression trees" into a larger framework of Sage. It not
only keeps engineers and calculus students happy, it keeps me happy
when just want to sit down and "solve for x" or compute a symbolic
integral. It's also a statement that rigor should not get in the way
of usability (an important component of being a viable alternative to
the M's).
Nope, nothing I see here could reasonably be called a Parent. Is
there an object that formally represents "the ring of integers" that
one could pass to a function?
>
>>
>> Categories in Sage aren't as fully developed as they could be,
>> but they do form a lattice.
>>
>> sage: C = category(FiniteField(13)); C
>> Category of fields
>> sage: C.is_subcategory(Rings())
>> True
>>
>
> Ah, I think that is very useful. Does the "Category of fields" inherit
> any generic operations from "Rings()"?
Very little, because up 'till now there wasn't an easy way to use
them. But this is the plan.
>
>> One of the changes in this latest coercion push is to allow Parents
>> to belong to several categories.
>>
>
> That sounds good. So does this mean that the 'category' method now
> returns a list?
There are two methods, category (returning the first cat in the list)
and categories (returning the whole list). The set of categories is
an ordered list, and things are attempted in the first category (and
its supercategories) before moving on to the next. Most objects, of
course, will belong to one category (not including its super
categories).
>
>> ... This will also be a place to put generic implementations of
>> algorithms. For example, one has the category of Euclidian
>> domains, where one can define a generic gcd algorithm that will
>> then get attached (as a method) to every element of that category.
>> ...
>>>
>>> sage: Zmod(7).is_field()
>>> True
>>>
>>> Perhaps an alternative would be to provide a coercion to some
>>> object whose parent is in Fields()?
>>
>> The is_field method tests whether or not this specific instance
>> happens to be a field, which is true in this case. But we do provide
>> a coercion
>>
>> sage: F = GF(7)
>> sage: F.coerce_map_from(Zmod(7))
>>
>> Natural morphism:
>> From: Ring of integers modulo 7
>> To: Finite Field of size 7
>>
>
> Excellent. But, hmm... this is supposed to be a morphism in some
> category. What category?
sage: mor = F.coerce_map_from(Zmod(7))
sage: parent(mor)
Set of Homomorphisms from Ring of integers modulo 7 to Finite Field
of size 7
sage: parent(mor).homset_category()
Category of rings
(The reason we have to ask for "homset_category" is because the
category of a homset is (generically at least) Sets.
Yeah, Q is a field after all.
> What makes them a MatrixSpace?
It's just a term used for the set of all matrices of a specified
dimension and basering.
>
>> In the case of square matrices, it is actually an algebra
>> over its basering, and its category is the category of Q-algebras.
>>
>
> Anyway, I suppose this taxonomy is not such a big deal. One only needs
> to learn the language - of course it has turned out in Axiom that this
> is not that simple a problem - hence interactive tools for querying
> and navigating the domain and category hierarchies (hyperdoc).
Yeah. Introspection is important for something like this.
>
>> ...
>
> Thanks for continuing the discussion.
No problem. I am hopeful that this discussion is useful for more than
just you and I.
- Robert
Gonzalo Tornaria
>
>
> I didn't know that. I thought the only way to create a Python class
> is for the Python interpreter to execute Python code that looks like this:
>
> class Foo(...):
> ...
>
> That makes a new class called Foo. How are you going to make, at
> runtime, new classes for each of Z/nZ say, for 1 <= n < 10^5, i.e.,
> something like this in Sage:
>
> v = [Integers(n) for n in range(1,10^5)]
>
> I do not think it is possible to concisely create a hundred thousand
> separate Python classes like that.
FWIW, I think classes in python are "just" instances of class
classobj. This latter class classobj (from new import classobj) is
what implements the semantics of python's "new-style classes".
Also classes can be constructed at runtime, as instances of class
classobj. Moreover, one can also subclass "classobj" or "type" to
produce metaclasses with different semantics (or whatever).
Search for "python metaclass", e.g.
http://www.onlamp.com/pub/a/python/2003/04/17/metaclasses.html
WRT your example:
sage: def classinteger(m):
... class A:
... def __init__(self, n):
... self.__n = n % m
... def __repr__(self):
... return "%d mod %d" % (self.__n, m)
... A.__name__ = "classintegers mod %d" % m
... return A
sage: classinteger(5)(3)
3 mod 5
sage: time v = [classinteger(n) for n in range(1,10^5)]
CPU time: 3.64 s, Wall time: 4.14 s
sage: time w = [Integers(n) for n in range(1,10^5)]
CPU time: 15.75 s, Wall time: 16.21 s
sage: v[1256](34251235)
499 mod 1257
sage: w[1256](34251235)
499
Note that classinteger(n) is an instance of class 'classobj' in this
simple example.
sage: type(classinteger(n))
<type 'classobj'>
sage: classinteger(7)
<class __main__.classintegers mod 7 at 0x3697770>
But one could define a metaclass so this class works fine as a
first-class object (e.g. it prints nicely, can define members acting
on the class itself, etc):
I was going to say that using classes built at runtime like this was
fine but not a good option for performance reasons, but the example
above shows the overhead of python class creation time is not so
significant.
I guess the part of creating a class which will take time is the
creation of the dictionary; but then one could have the dictionary
pre-loaded somehow.
If, on the other hand, one doesn't care about parent creation time but
worries about element operation time, it might be possible to, at
class creation time, compile the class in cython so that e.g. the
modulus is a hard-coded compile constant... (don't know if this is of
any advantage, it is just for the sake of an example).
Best, Gonzalo
Bill Page
>
> On Wed, Jun 4, 2008 at 10:16 PM, Bill Page wrote:
>> ...
>> Python classes can also take parameters.
>
> I didn't know that. I thought the only way to create a Python class
> is for the Python interpreter to execute Python code that looks like this:
>
> class Foo(...):
> ...
>
> That makes a new class called Foo. How are you going to make, at
> runtime, new classes for each of Z/nZ say, for 1 <= n < 10^5, i.e.,
> something like this in Sage:
>
> v = [Integers(n) for n in range(1,10^5)]
>
> I do not think it is possible to concisely create a hundred thousand
> separate Python classes like that.
>
See Python MetaClasses. E.g.
http://www.onlamp.com/pub/a/python/2003/04/17/metaclasses.html
sage: X = type('X', (object,), dict(a=1))
sage: x=X()
sage: x.a
1
Regards,
Bill Page.
Bill Page
> On Thu, Jun 5, 2008 at 1:47 AM, William Stein wrote:
>> runtime, new classes for each of Z/nZ say, for 1 <= n < 10^5,
>> i.e., something like this in Sage:
>>
>> v = [Integers(n) for n in range(1,10^5)]
>>
>> I do not think it is possible to concisely create a hundred thousand
>> separate Python classes like that.
> WRT your example:
>
> sage: def classinteger(m):
> ... class A:
> ... def __init__(self, n):
> ... self.__n = n % m
> ... def __repr__(self):
> ... return "%d mod %d" % (self.__n, m)
> ... A.__name__ = "classintegers mod %d" % m
> ... return A
> sage: classinteger(5)(3)
> 3 mod 5
> sage: time v = [classinteger(n) for n in range(1,10^5)]
> CPU time: 3.64 s, Wall time: 4.14 s
> sage: time w = [Integers(n) for n in range(1,10^5)]
> CPU time: 15.75 s, Wall time: 16.21 s
> sage: v[1256](34251235)
> 499 mod 1257
Excellent example. Thank you!
> I was going to say that using classes built at runtime like this
> was fine but not a good option for performance reasons, but the
> example above shows the overhead of python class creation time
> is not so significant.
>
Agreed.
Python can treat types as fully first class objects. I think this
strongly supports William's initial decision to choose Python as the
basis for a new computer algebra system comparable to Magma and Axiom.
Regards,
Bill Page.
William Stein
Cool. I stand corrected. I'm glad this is fully supported in Python,
at least for new style classes.
OK, I have to come clean and admit that already actually knew all
about metaclasses, but considered using them this way so
unnatural and ugly that I would not recommend it.
> I was going to say that using classes built at runtime like this was
> fine but not a good option for performance reasons, but the example
> above shows the overhead of python class creation time is not so
> significant.
It is very impressive how fast Python classobj is,
but it is still an order of magnitude slower than object instantiation,
though your flawed benchmark incorrectly suggests the opposite.
Two points:
(1) It makes no sense to compare with Integers(n), which is
doing all kinds of other things (it could for all we know for now
be factoring n), etc. You should compare with a minimal
class MyIntegers. Creating Python class instances can be massively
optimized by doing the class creation in Cython, like we do in integer.pyx.
(2) It is important to consider the performance implications if the class you're
constructing has lots of methods. I tried the above benchmark but
with 8 trivial methods added to the class, and it took 5 times as long
to run! Of course, this could easily be mitigated with inheritance,
but is worth
pointing out. Also, every single classobj takes *memory* for all methods
that are special for it. Again inheritance could mitigate this problem.
Here's an example benchmark but using MyIntegers instead of Integer.
In it, creating a class takes four times as long as instantiating an instance:
{{{id=54|
def classinteger(m):
class A:
def __init__(self, n):
self.__n = n % m
def __repr__(self):
return "%d mod %d" % (self.__n, m)
A.__name__ = "classintegers mod %d" % m
return A
///
}}}
{{{id=57|
class MyIntegers:
def __init__(self, m):
self.__m = m
///
}}}
{{{id=55|
time v = [classinteger(n) for n in range(1,10^5)]
///
CPU time: 2.89 s, Wall time: 2.98 s
}}}
{{{id=53|
time w = [MyIntegers(n) for n in range(1,10^5)]
///
CPU time: 0.63 s, Wall time: 0.63 s
}}}
{{{id=52|
2.89 / 0.63
///
4.58730158730159
}}}
Bill Page
> ...
> On Thu, Jun 5, 2008 at 12:23 PM, Gonzalo Tornaria wrote:
>> sage: def classinteger(m):
>> ... class A:
>> ... def __init__(self, n):
>> ... self.__n = n % m
>> ... def __repr__(self):
>> ... return "%d mod %d" % (self.__n, m)
>> ... A.__name__ = "classintegers mod %d" % m
>> ... return A
>> sage: classinteger(5)(3)
> OK, I have to come clean and admit that already actually knew
> all about metaclasses, but considered using them this way so
> unnatural and ugly that I would not recommend it.
>
Do you still consider the example code like that given above by
Gonzalo "unnatural and ugly"? It seems like pretty standard Python to
me. There of course various ways of packaging the metaclass machinery
to provide an even cleaner user interface.
The point I am making is that Python already has everything you need
to implement the concept of the "parent" of an element as a type, i.e.
directly as an instance of the class that creates it. This could
result in considerable simplification compared to the current
situation. Right now a Sage user has to deal with three separate
type-related concepts: type (Python class), parent (an object in some
category) and category. Metaclasses are a good match for Categories.
Thus you can have everything in one package and remain closer to
conventional Python programming in spirit.
Regards,
Bill Page.
Mike Hansen
>>> ... class A:
>>> ... def __init__(self, n):
>>> ... self.__n = n % m
>>> ... def __repr__(self):
>>> ... return "%d mod %d" % (self.__n, m)
>>> ... A.__name__ = "classintegers mod %d" % m
>>> ... return A
>>> sage: classinteger(5)(3)
>> ...
>> OK, I have to come clean and admit that already actually knew
>> all about metaclasses, but considered using them this way so
>> unnatural and ugly that I would not recommend it.
>>
>
> Do you still consider the example code like that given above by
> Gonzalo "unnatural and ugly"? It seems like pretty standard Python to
> me. There of course various ways of packaging the metaclass machinery
> to provide an even cleaner user interface.
To me, that feels much more "unnatural and ugly" than the current
situation. Perhaps that's just because I'm used to things the way
they are currently.
> The point I am making is that Python already has everything you need
> to implement the concept of the "parent" of an element as a type, i.e.
> directly as an instance of the class that creates it. This could
> result in considerable simplification compared to the current
> situation. Right now a Sage user has to deal with three separate
> type-related concepts: type (Python class), parent (an object in some
> category) and category. Metaclasses are a good match for Categories.
> Thus you can have everything in one package and remain closer to
> conventional Python programming in spirit.
For me, a class represents an implementation / model of a mathematical
object whether it be a multivariate polynomial or a multivariate
polynomial ring. I don't see why one would want to tie the
implementation of the parent structure together with that of its
elements.
In your setup, I'm still not sure where place things like
.multiplicative_generators(). Sure, you could make it a class method
so that you don't need an instance to access it, but then all of that
class's instance have that method which was not the intention.
--Mike
Robert Bradshaw
Parents are much, much more than types, and I see little to gain by
trying to unify the two. To me, "Python class" and "object in some
category" are very different concepts.
There are, however, many drawbacks. First, there is far from a one-
two-one correspondence between types an parents. Some parents have
elements of varying type, and some types are used for a variety of
parents. Though metaclasses can be used to mitigate this somewhat, I
think it would add much more complexity than it would remove. Also,
the __call__ method often does much more than just construct a new
object.
It would also introduce enormous performance overhead--the cost of
calling classinteger above increases linearly with the number of
methods the class has. Perhaps this could be limited by pre-creating
the methods and adding them to the dict all at once, but but this
would be even more unnatural. More importantly, this would mean that
cdef'd (Cython) elements/parents are out. (Well, one could generate a
cython file at run time, compile it, load the .so and then use that
"dynamically generated type," but that would be so insanely slow). Is
this a deficiency of cdef'd types? Yes. But this is what makes them
so fast--members and methods can be statically resolved at compile time.
- Robert
Bill Page
>
> On Jun 4, 2008, at 7:35 PM, Bill Page wrote:
> Types are all about the implementations of things, they
> synonymous with the "classes" of Object Oriented programming,
> and are insufficient (and the wrong vehicle) to carry deeper
> mathematical structure.
That is precisely the claim that I am trying to dispute. Axiom for
example specifically takes the point of view that mathematical
structure should be implemented as types in a sufficiently rich
programming language, i.e. one where types are first order objects. I
think that Axiom (and later the programming language Aldor that was
originally implemented as the "next generation" library compiler for
Axiom) has demonstrated that this view is viable.
> For example, an element of Z/5Z and an element of Z/7Z have
> the same type (i.e. they're instances of the same class, with
> the same methods, internal representation, etc), but not the same
> parent (which is data).
>
I think Gonsalo Tornaria in this thread has demonstrated that this
need not be the case.
> ...
>> "SymbolicRing" is another thing that worries me a little. Exactly what
>> (mathematically) is a symbolic ring? E.g. Why is it a ring and not a
>> field? Operations like 'sin(sin(sin(x))).parent()' doesn't look much
>> like a ring to me. Probably most other computer algebra systems
>> would just call this an "expression tree" or AST or something like
>> that.
>
> Most symbolic systems (Magma and Axiom excluded) have very
> little notion of algebras, rings, modules, etc.--virtually everything is
> an expression or list of expressions. SymbolicRing is a way to fit
> these abstract "expression trees" into a larger framework of Sage.
> It not only keeps engineers and calculus students happy, it keeps
> me happy when just want to sit down and "solve for x" or compute
> a symbolic integral. It's also a statement that rigor should not get
> in the way of usability (an important component of being a viable
> alternative to the M's).
>
It seems to me that symbolic computation need not lack rigor. And in a
computer algebra system rigor should not compromise usability. These
are design challenges that we face. This is one of the things that
strongly-typed systems like Axiom and Magma (and now Sage) offer over
these other systems.
> ...
>>
>> To the best of my knowledge, the formal concept of "representation"
>> in a computer algebra systems is unique to Axiom. Each domain
>> has a specified 'Rep'consisting of some domain expression. E.g.
>>
>> Rep == Integer
>>
>> or
>>
>> Rep == Union(Record(car:%, cdr:%), "nil")
>>
>> % stands for "this domain", so this definition amounts to some
>> recursive structure. Of course recursively defined classes are also
>> possible in Python but they require much more discipline on the
>> part of the programmer.
>>
>> Now there are two operators 'rep' and 'per' that provide coercion to
>> and from Rep. So we might say for example that some operation *
>> in this domain is implemented by * from Rep:
>>
>> x * y == per( rep(x) * rep(y) )
>>
>> The trouble is: I don't see anything here that could be called a
>> "parent". Or perhaps, they are all "parents" except for those
>> initial domains for which we never explicitly specify a Rep.
>
> Nope, nothing I see here could reasonably be called a Parent.
> Is there an object that formally represents "the ring of integers"
> that one could pass to a function?
>
Sure. For example in OpenAxiom one can write:
(1) -> Integer has Ring
(1) true
Type: Boolean
(2) -> Pick(x,A,B)==if x=0 then A else B
Type: Void
(3) -> i:=1
(3) 1
Type: PositiveInteger
(4) -> n:Pick(i,Integer,Float):=1
Compiling function Pick with type (PositiveInteger,Domain,Domain)
-> Domain
(4) 1.0
Type: Float
And of course even easier in Sage :-)
sage: category(IntegerRing())
Category of rings
sage: def Pick(x,A,B):
if x==0:
return A
else:
return B
....:
sage: i=1
sage: n=Pick(i,IntegerRing(),RealField())(1)
sage: n; parent(n)
1.00000000000000
Real Field with 53 bits of precision
> ...
>>
>>> One of the changes in this latest coercion push is to allow
>>> Parents to belong to several categories.
>>>
>>
>> That sounds good. So does this mean that the 'category' method
>> now returns a list?
>
> There are two methods, category (returning the first cat in the
> list) and categories (returning the whole list). The set of categories
> is an ordered list, and things are attempted in the first category
> (and its supercategories) before moving on to the next. Most
> objects, of course, will belong to one category (not including its
> super categories).
>
Here is something to consider: Instead of having two methods, one that
return a list, you could implement some operations on categories that
yield new categories. One such operation (because the categories form
a lattice) might be 'Join'. Then if a parent belongs to more than one
category you could just return their Join. One of the methods provided
with category could be one that lists it's component categories, or
more generally the expression (if any) that created it.
From category theory it is easy to imagine other operations that one
might want to perform on categories to create new categories.
> ...
Regards,
Bill Page.
Robert Bradshaw
> On Thu, Jun 5, 2008 at 4:16 AM, Robert Bradshaw wrote:
>>
>> On Jun 4, 2008, at 7:35 PM, Bill Page wrote:
>> ...
>> Types are all about the implementations of things, they
>> synonymous with the "classes" of Object Oriented programming,
>> and are insufficient (and the wrong vehicle) to carry deeper
>> mathematical structure.
>
> That is precisely the claim that I am trying to dispute. Axiom for
> example specifically takes the point of view that mathematical
> structure should be implemented as types in a sufficiently rich
> programming language, i.e. one where types are first order objects. I
> think that Axiom (and later the programming language Aldor that was
> originally implemented as the "next generation" library compiler for
> Axiom) has demonstrated that this view is viable.
Perhaps Sage is demonstrating that not taking this view is viable
too :-).
>
>> For example, an element of Z/5Z and an element of Z/7Z have
>> the same type (i.e. they're instances of the same class, with
>> the same methods, internal representation, etc), but not the same
>> parent (which is data).
>>
>
> I think Gonsalo Tornaria in this thread has demonstrated that this
> need not be the case.
By making elements of Z/5Z and Z/7Z be instances of different
(dynamically generated) classes. This has drawbacks too, also
mentioned in this thread. Personally, I find it conceptually easier
to think of Z/5Z and Z/7Z as distinct instances of the same class.
>
>> ...
>>> "SymbolicRing" is another thing that worries me a little. Exactly
>>> what
>>> (mathematically) is a symbolic ring? E.g. Why is it a ring and not a
>>> field? Operations like 'sin(sin(sin(x))).parent()' doesn't look much
>>> like a ring to me. Probably most other computer algebra systems
>>> would just call this an "expression tree" or AST or something like
>>> that.
>>
>> Most symbolic systems (Magma and Axiom excluded) have very
>> little notion of algebras, rings, modules, etc.--virtually
>> everything is
>> an expression or list of expressions. SymbolicRing is a way to fit
>> these abstract "expression trees" into a larger framework of Sage.
>> It not only keeps engineers and calculus students happy, it keeps
>> me happy when just want to sit down and "solve for x" or compute
>> a symbolic integral. It's also a statement that rigor should not get
>> in the way of usability (an important component of being a viable
>> alternative to the M's).
>>
>
> It seems to me that symbolic computation need not lack rigor. And in a
> computer algebra system rigor should not compromise usability. These
> are design challenges that we face. This is one of the things that
> strongly-typed systems like Axiom and Magma (and now Sage) offer over
> these other systems.
Yes, I'm not trying to throw away rigor. I'm just saying it's a place
where I can write "a+b" without having to think about what the "+"
really means, it's just formal. Perhaps it could equally be called
"SetOfUnevaluatedExpressions" but that's a bit verbose.
Your Integer object also stands for the ring of Integers. So can you
do stuff like
sage: ZZ.is_commutative()
True
sage: ZZ.krull_dimension()
1
sage: ZZ.is_finite()
False
>
>> ...
>>>
>>>> One of the changes in this latest coercion push is to allow
>>>> Parents to belong to several categories.
>>>>
>>>
>>> That sounds good. So does this mean that the 'category' method
>>> now returns a list?
>>
>> There are two methods, category (returning the first cat in the
>> list) and categories (returning the whole list). The set of
>> categories
>> is an ordered list, and things are attempted in the first category
>> (and its supercategories) before moving on to the next. Most
>> objects, of course, will belong to one category (not including its
>> super categories).
>>
>
> Here is something to consider: Instead of having two methods, one that
> return a list, you could implement some operations on categories that
> yield new categories. One such operation (because the categories form
> a lattice) might be 'Join'. Then if a parent belongs to more than one
> category you could just return their Join. One of the methods provided
> with category could be one that lists it's component categories, or
> more generally the expression (if any) that created it.
That is a very good idea, I think we should do this.
- Robert
Bill Page
>
> Your Integer object also stands for the ring of Integers. So can
> you do stuff like
>
> sage: ZZ.is_commutative()
> True
> sage: ZZ.krull_dimension()
> 1
> sage: ZZ.is_finite()
> False
>
(1) -> Integer has Finite
(1) false
Type: Boolean
(2) -> Integer has CommutativeRing
(2) true
Type: Boolean
(3) -> characteristic()$Integer
(3) 0
Type: NonNegativeInteger
Finite and CommutativeRing are categories. characteristic is a method
of Integer.
Regards,
Bill Page.
Carl Witty
> The point I am making is that Python already has everything you need
> to implement the concept of the "parent" of an element as a type, i.e.
> directly as an instance of the class that creates it. This could
> result in considerable simplification compared to the current
> situation. Right now a Sage user has to deal with three separate
> type-related concepts: type (Python class), parent (an object in some
> category) and category. Metaclasses are a good match for Categories.
> Thus you can have everything in one package and remain closer to
> conventional Python programming in spirit.
very well in Aldor, and that it's possible in Python. However, I
don't think it's a good idea for Sage, because of pragmatic
limitations due to Python/Cython. These points were all made by other
people in this thread, but I'm combining them into one convenient list
here for debating purposes.
1. Creating classes at run-time in Python is slow and memory-hungry;
the more methods you add to the class, the slower and more memory-
hungry it gets.
2. Creating classes at run-time in Cython is impossible. If Cython
was extended to support creating classes at run-time (which it may be,
at some point, since there are Cython developers who want it to
compile the entire Python language) there would need to be significant
language extensions (and corresponding compiler extensions) to let
these run-time be classes be as fast as the former compile-time
classes.
3. While it is possible in Python to add methods to types, these
methods are inherited by values of that type. Currently, we have
methods on parents that are not inherited by methods on the values; we
don't want to lose that.
Here are some other points that were made, but seem to me less
important than the above:
4. Currently, if P is a parent and T is a type, P(...) and T(...) both
have meanings which are not the same. If we merged parents and types,
we would have to leave the function-call syntax for constructors, and
come up with a new syntax for conversion (coercion). This would be
annoying and non-backward-compatible, but is not a fatal flaw.
5. People have mentioned that we have cases where values of the same
type can have different parents, like IntegerMod_int which is the
implementation for both GF(7) and Integers(7). This particular case
is not a problem at all... we could just have GF(7) inherit from
Integers(7).
6. People have mentioned that we have cases where values of different
types share the same parent, such as symbolic constants and symbolic
expressions both being in the symbolic ring. This may sometimes not
be a problem; for instance, parent(x)==SR could be replaced by
isinstance(x, SR), even though type(x)==SR would not work.
In summary, I think we could unify parents and types with a large
amount of work and a significant performance loss; or with a very
large amount of work (including a lot of work on Cython) and a minor
performance loss. In either case, we would have a functionality loss
(the ability to put methods on parents that are not inherited by
values).
I don't think it's worth it.
Carl
Gonzalo Tornaria
>
> On Thu, Jun 5, 2008 at 1:18 PM, William Stein wrote:
> > OK, I have to come clean and admit that already actually knew
> > all about metaclasses, but considered using them this way so
> > unnatural and ugly that I would not recommend it.
> >
>
>
> Do you still consider the example code like that given above by
> Gonzalo "unnatural and ugly"? It seems like pretty standard Python to
> me. There of course various ways of packaging the metaclass machinery
> to provide an even cleaner user interface.
I do think my code (as posted) is "unnatural and ugly". Just an
example, not even a class. I'd really love to figure out a proper
model of parent-element relationship using metaclasses, etc.
I tried very hard (several times) to find such a model, or at least "a
cleaner user interface", and I played with some interesting ideas but
cooked no cake.
Maybe you have some cool ideas on how one can make this work (and
appealing to others).
> category) and category. Metaclasses are a good match for Categories.
Also cool. I remember I was hit by the python version of a few of the
paradoxes in early days set thery. YMMV.
Best, Gonzalo
Robert Bradshaw
I think that sums it up well (though personally I think types and
parents are distinct enough concepts that keeping them separate is a
good idea anyways). On (2), creating classes at runtime will
certainly be supported, that's how "normal" classes are done now, but
creating cdef classes at runtime without invoking gcc (or some other
sort of compiler) at runtime is perhaps not even possible. For (6) we
do a lot of "parent(x) is SR", especially in the fast paths of the
coercion model (due to caching, parents are usually unique, and if
not we follow it up by the (relatively expensive) ==, so even with a
very large amount of work, that would still slow things down.
This thread got slightly off topic to the original question though:
where should sqrt(2) belong? Z[sqrt(2)], SR, or somewhere else?
- Robert
Nick Alexander
> where should sqrt(2) belong? Z[sqrt(2)], SR, or somewhere else?
I always want my data to start as close to the initial object as
possible. In this case, Z[sqrt(2)] \into SR and not vice versa -- so
sqrt(2) should be in Z[sqrt(2)].
Nick
Bill Page
>
> OK, I'm convinced that unifying parents and types is elegant and
> works very well in Aldor, and that it's possible in Python. However,
> I don't think it's a good idea for Sage, because of pragmatic
> limitations due to Python/Cython. These points were all made by
> other people in this thread, but I'm combining them into one
> convenient list here for debating purposes.
Thanks for summarizing the arguments against unifying parents and
types in Sage. As I stated near the start of this thread, it was not
my intention to recommend that Sage actually take this approach. Sage
is 3.0 and the current approach is established and working well for
the end user. My regret however is that it seems to make new
development in Sage more complicated, with a steeper learning curve
for new developers than it really needs to be.
On Thu, Jun 5, 2008 at 9:09 PM, Gonzalo Tornaria
<torn...@math.utexas.edu> wrote:
>
>>
>> On Thu, Jun 5, 2008 at 1:18 PM, William Stein wrote:
>> > OK, I have to come clean and admit that already actually knew
>> > all about metaclasses, but considered using them this way so
>> > unnatural and ugly that I would not recommend it.
>> >
>>
> On 6/5/08, Bill Page <bill...@newsynthesis.org> wrote:
>> Do you still consider the example code like that given above by
>> Gonzalo "unnatural and ugly"? It seems like pretty standard
>> Python to me. There of course various ways of packaging the
>> metaclass machinery to provide an even cleaner user interface.
>
> I do think my code (as posted) is "unnatural and ugly". Just an
> example, not even a class. I'd really love to figure out a proper
> model of parent-element relationship using metaclasses, etc.
>
It seems to me that the parent-element relationship is already there.
It is inherent in the Python concept of an instance of a class. In
this model classes *are* parents. Instances are elements. Your example
code presents a function that generates classes based on a parameter -
a standard "class factory". As you say, what might be more interesting
is a class that behaves like a category, i.e. whose instances
themselves are classes. Python's type module is an example of such a
metaclass.
I think this is all very well and concisely described here:
http://www.onlamp.com/pub/a/python/2003/04/17/metaclasses.html
A Primer on Python Metaclass Programming
by David Mertz
> I tried very hard (several times) to find such a model, or at least
> "a cleaner user interface", and I played with some interesting ideas
> but cooked no cake.
>
> Maybe you have some cool ideas on how one can make this work
> (and appealing to others).
>
Have you looked at the 'type' class? I think that is pretty cool.
>> category) and category. Metaclasses are a good match for
>> Categories.
>
> Also cool. I remember I was hit by the python version of a few of
> the paradoxes in early days set theory. YMMV.
>
I don't think you should let apparent paradoxes in set theory worry
you. OpenAxiom for example has a domain called Domain. Domain is an
instance of Domain so we can write:
x:Domain:=Domain
and the world does not end..
And would I really would like to get as much mileage as possible from
native Python constructions.
In the Axiom world we have a bit of a problem: Unlike Axiom itself,
the new generation Axiom library compiler - Aldor - is not free
(although just last year it finally became at least "semi-open"). In
my opinion, to progress rapidly and as intended when Axiom was still a
commercial product Axiom badly needs Aldor. The existing compiler in
Axiom - SPAD - is significantly harder to use and to extend (although
Gabriel Dos Reis has started to do some wonderful new things with it
in the OpenAxiom fork of the original open source Axiom project.). I
think this is one of the main reasons why Axiom has had a relatively
difficult time attracting new uses and developers.
Even Aldor however is already a rather "old" language. There has been
a lot of progress in programming language design since the mid 1990's.
Python is one result. It seems to me that Python has many of the
attributes of the kind of language that the Axiom developers were
trying to achieve. This is one of the primary reasons why I was
interested in the Sage project in the first place. Sage is not Axiom,
but it is (more or less) based on Magma and Magma has at least some of
the features that I think make Axiom rather special as a computer
algebra system. So like several other people before me, e.g. the
developers of Sage and the developers of SymPy, I am now thinking
quite seriously about how much of the Axiom design could actually be
accommodated in Python. Implementing domains and categories as Python
classes is an essential part of that strategy.
Regards,
Bill Page.
John Cremona
I agree: for exactly the same reason that I expect 2 to be in ZZ and not in SR!
What about sqrt(8)? Would you put that in the maximal order
Z[sqrt(2)] or in its own order Z[sqrt(8)]? I would recommenfd the
former, otherwise every time you type sqrt(n) for some integer n then
you would be forcing the factorization of n (on the grounds that there
is no known algorithm for finding the square-free part of an integer
which is faster than factorization).
Similarly for other algebraic integers alpha (at least those which
have a "symbolic" definition like sqrt(n)).
What about sqrt(1/2)? (or anyother non-integral algebraic number)?
This does not belong to any finitely-generated ring (I mean f.g. as
Z-module of course) so I would put it straight into the field
Q(sqrt(2)).
John
Robert Bradshaw
Yes, this is what I was thinking to--I think I've mentioned it
before, but I think there should be a coercion parent(a) -> parent
(sqrt(a)).
- Robert