On scientific computing, Python and Julia

[Luca De Feo]

I guess these two posts by G. Hoare (creator of the Rust language) are
of interest to the community:

http://graydon2.dreamwidth.org/3186.html
http://graydon2.dreamwidth.org/189377.html

See also the presentation linked therein:

http://nbviewer.ipython.org/gist/StefanKarpinski/b8fe9dbb36c1427b9f22

I am very impressed by the way Julia's generic functions (aka
multimethods) solve many of the problems we constatly have to face
when designing Sage's type system (expression problem).

Luca

[Robert Bradshaw]

Julia is certainly very impressive.

I'm not sure multimethods alone are enough to solve issues with Sage's
type system (e.g. coercion, where the result of a+b may happen in a
completely new domain) but they could certainly help.

- Robert

[Nils Bruin]

On Thursday, July 17, 2014 10:35:32 PM UTC-7, Robert Bradshaw wrote:

I'm not sure multimethods alone are enough to solve issues with Sage's
type system (e.g. coercion, where the result of a+b may happen in a
completely new domain) but they could certainly help.

Indeed. If you want multiple dispatch in python you can have it:

 class multi_dispatch(object):
    def __init__(self):
        self.dispatch = []
    def register_function(self,signature,function):
        self.dispatch.append( (signature,function) )
    def register(self,*signature):
        def register_decorator(function):
            self.register_function(signature,function)
            return function
        return register_decorator
    def __call__(self,*args):
        for sig in self.dispatch:
            if len(args) == len(sig) and all(isinstance(a,b) for a,b in zip(args,sig[0])):
                return sig[1](*args)
        raise TypeError("Signature not implemented")

With this it is now quite straightforward to write multiple dispatch at the cost of using a decorator:

T = multi_dispatch()

@T.register(int,int)
def int_int(a,b):
    return "called with integers (%d,%d)"%(a,b)

@T.register(int,float)
def int_float(a,b):
    return "called with integer and float (%d,%f)"%(a,b)

@T.register(float,float)
def float_float(a,b):
    return "called with floats (%f,%f)"%(a,b)

You can then just call T:

>>> T(1,2)
'called with integers (1,2)'
>>> T(1,2.0)
'called with integer and float (1,2.000000)'
>>> T(1.0,2.0)
'called with floats (1.000000,2.000000)'
>>> T(1.0,2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "dispatch.py", line 15, in __call__
    raise TypeError("Signature not implemented")
TypeError: Signature not implemented

The key is in optimizing the signature lookup (and the kind of lookup can definitely be refined a *lot* in the above implementation--linear search by registration order is probably not optimal) and finding a useful way of describing signatures and properly prioritizing which signature should get precedence (if you think C3 is complicated for multiple inheritance, brace yourself for multiple dispatch).

In a dynamically typed language, multiple dispatch has a rather high cost. LISP MOP implementations have a lot of experience. In the end it's just a little helper tool, though. It doesn't really do things you can't otherwise do.

[William Stein]

Magma uses "multiple dispatch" extensively to deal with what our coercion model code addresses, and it's definitely a step back in my experience. 

> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To post to this group, send email to sage-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-devel.
> For more options, visit https://groups.google.com/d/optout.

[rjf]

As Nils says, Lisp has dealt with this in CLOS  (Common Lisp Object System.

Among other advantages, it is possible to properly compile such calls when the

type discrimination can be done at compile time.

   If you have enough 

mathematical categories, the discrimination based on implementation structural

types leaves much to be decided by other mechanisms which may or may not

be convenient to the user.

perhaps the simplest is adding floats to exact rationals.  You could muddy up

the result by converting a rational to the nearest float, or you could convert the

float to a rational  (exactly, since every float is integer X 2^ integer, and hence

an exact rational.

It gets much worse if you toss in complex numbers, say.

People have thought about such things in the computer algebra community for

maybe 4 or 5 decades. (eg. see Axiom)

In lisp I think this would be

  (defmethod T((a integer) (b integer))  

       (format nil "called with integers ~s, ~s"  a b))

[Nils Bruin]

On Friday, July 18, 2014 7:26:00 AM UTC-7, William wrote:

Magma uses "multiple dispatch" extensively to deal with what our coercion model code addresses, and it's definitely a step back in my experience. 

I agree that in scenarios for which coercion was designed, it works better than multiple-dispatch based solutions. We still have serious problems, but they seem to be inherent to the task at hand.

Multiple dispatch does something different, however: it provides a decentralized way of extending functionality.

Suppose for instance that I'm implementing a class for genus 1 models and I want to make a routine that constructs an elliptic curve out of it when a point is given. Pythonic would be

<genus_one_model>.elliptic_curve(<point>)

However, this is at odds with how elliptic curves get constructed otherwise:

EllipticCurve(<elliptic curve data>)

so it would be nice if this were to work:

EllipticCurve(<genus_one_model>,<point>)

but for that to work one would need to *edit* EllipticCurve to include yet another case (together with code to detect when the new constructor should be called). If we'd find this happens a lot we'd probably end up writing something like

EllipticCurve(first_data, *rest):
    if has_attribute(first_data,"elliptic_curve"):
        return first_data.elliptic_curve(*rest)
    ...

but that is just an ad-hoc partial implementation of multiple dispatch. It's also not very self-documenting, because when we're implementing the "elliptic_curve" method it isn't immediately clear that EllipticCurve will be picking this up.
 
It's conceivable that especially constructors/factories would benefit from some uniform dispatch system, since these tend to be the routines that accept the widest variety of input. That said, it does look like single dispatch solves about 80%, and python has a nice low cost for that (by spelling it as a method call).

[Incidentally, when you google "python multimethod" or "python multiple dispatch" you'll find a dozen or so other python implementations, although most of those run into namespace issues that the above doesn't suffer from]

[Robert Bradshaw]

Yes, but here T must explicitly implement to multiple dispatch. One
can't take an existing object and "add a method" based on the type of
the second argument.

We also have @coerce_binop
https://github.com/sagemath/sage/blob/6.2/src/sage/structure/element.pyx#L3260
which is nice. Here we're conflating the issue of types and Parents
(where the latter can be viewed, somewhat, as dynamically created
parameterized types), and it's often the latter that affects the
choice of implementation (or even semantic meaning) of a method.

- Robert

[Stefan Karpinski]

Hi Sage folks! I gave that presentation on multiple dispatch at Strange Loop last year. Since someone pointed this thread out to me and the topic of multiple dispatch and type promotion is near and dear to my heart, I thought I'd chime in. Maybe my comments will be of some use. (Btw, since some people prefer their notebooks with some explanation, if anyone would care to see a video of a similar but not quite identical talk, this one from CodeMesh last year is pretty close, although it spends less time on multiple dispatch and type promotion.)

On Thursday, July 17, 2014 10:35:32 PM UTC-7, Robert Bradshaw wrote:

I'm not sure multimethods alone are enough to solve issues with Sage's type system (e.g. coercion, where the result of a+b may happen in a completely new domain) but they could certainly help.

I'm not familiar with Sage's issues with coercion (or even what you mean exactly by coercion), so take this with a grain of salt. I'll hazard a guess that you mean the problem of how to let the user specify what to do when it encounters something like a + b where a and b are different types, one or both of which are user defined, and which may or may not know anything about each other. When developing Julia, initially we didn't realize that multiple dispatch could handle this kind of problem so elegantly, but it turns out it can, which I attempted to explain in my talk and which I will say a bit about below.

On Thursday, July 17, 2014 11:56:43 PM UTC-7, Nils Bruin wrote:

In a dynamically typed language, multiple dispatch has a rather high cost. LISP MOP implementations have a lot of experience.

If Julia has shown anything, it's that you *can* have ubiquitous multiple dispatch in a dynamic language and have it be very fast – it is commonly a zero-cost abstraction in Julia.

First let me explain how multiple dispatch allows a promotion system that users can hook into for their own types by explaining how 1 + 2.5 works in Julia. 1 and 2.5 are of type Int and Float64 respectively, both of which are defined in Julia just like any user-defined type would be. When you do 1 + 2.5, there is no specific method for +(::Int, ::Float64), so a very generic fallback method is called:

+(x::Number, y::Number) = +(promote(x,y)...)

This tries to add a pair of generic numbers by calling promote(1,2.5) to convert 1 and 2.5 to a common type, and then calling the + function on the resulting like-typed pair. The promote function itself is just a generic function, defined in Julia, that calls the promote_rule function, another generic function, which acts as a user-extensible lookup table, mapping pairs of argument types the desired promotion type. In this case, there's a promote_rule method indicating that promote(1,2.5) should convert its arguments to Float64 (using the convert function – yet another user-extensible generic function). The resulting Float64 pair (1.0,2.5) is added in the usual way.

That is nice and all, but how can it possibly be fast? The trick is that Julia's type inference knows about generic functions and multiple dispatch and can actually analyze the promote function's behavior, determine that all it does for (Int,Float64) arguments is convert the Int to Float64, then inline that behavior into code for the expression +(promote(x,y)...). As a result, the LLVM code that does 1 + 2.5 is just this:

julia> @code_llvm 1 + 2.5

define double @"julia_+;19364"(i64, double) {

top:

  %2 = sitofp i64 %0 to double, !dbg !1885

  %3 = fadd double %2, %1, !dbg !1885

  ret double %3, !dbg !1885

}

I.e. just a signed int to floating-point conversion of the first argument and a floating-point add of the result of that with the second argument. Or in native x86_64 code:

julia> @code_native 1 + 2.5

.section __TEXT,__text,regular,pure_instructions

Filename: promotion.jl

Source line: 158

push RBP

mov RBP, RSP

Source line: 158

vcvtsi2sd XMM1, XMM0, RDI

vaddsd XMM0, XMM1, XMM0

pop RBP

ret

Despite the layers of abstraction, the meat of the addition method is just two instructions: vcvtsi2sd and vaddsd (the rest is function prologue and epilogue).

One can be forgiven for suspecting that this might only work for simple built-in types like Int and Float64, but that's not the case – it's completely general. Among my favorite examples is the SIUnits package. This package defines an SIQuantity wrapper type that tracks SI units in the type system. To teach Julia to do 1s + 2.5s (where `s` is a constant exported by SIUnits bound to the SIQuantity representing one SI second), all one needs is to provide a method for generically adding two SIQuantity values with the same units:

function +{T,S,m,kg,s,A,K,mol,cd}(

    x::SIQuantity{T,m,kg,s,A,K,mol,cd},y::SIQuantity{S,m,kg,s,A,K,mol,cd})

    val = x.val+y.val

    SIQuantity{typeof(val),m,kg,s,A,K,mol,cd}(val)

end

It's a bit verbose, I'll admit, but it buys you an awful lot. For example:

julia> using SIUnits, SIUnits.ShortUnits # allow `s` for `seconds`

julia> 1s + 2.5s

3.5 s

That's kind of cool, but what's really cool is that even with the additional layers of abstraction, this still boils down to just two CPU instructions:

julia> @code_native 1s + 2.5s

.section __TEXT,__text,regular,pure_instructions

Filename: /Users/stefan/.julia/SIUnits/src/SIUnits.jl

Source line: 139

push RBP

mov RBP, RSP

Source line: 139

vcvtsi2sd XMM1, XMM0, RDI

vaddsd XMM0, XMM1, XMM0

Source line: 140

pop RBP

ret

With no further definitions you also get much more complex behavior, like, for example, this:

julia> (7//5)m/s + (12873918273918273918273981273981273981273980 + 2im)m/s

64369591369591369591369906369906369906369907//5 + 2//1*im m s⁻¹

julia> typeof(ans)

SIQuantity{Complex{Rational{BigInt}},1,0,-1,0,0,0,0} (constructor with 1 method)

The SIUnits package doesn't mention Complex, Rational or BigInt, yet this works.

On Thursday, July 17, 2014 11:56:43 PM UTC-7, Nils Bruin wrote:

In the end it's just a little helper tool, though. It doesn't really do things you can't otherwise do.

Of course, you can do anything in one turing complete language that can be done in any other, but that's clearly not a good argument, since if it was, we'd all be programming in machine code. If multiple dispatch is "just a little helper tool" then so are object orientation and functional programming. In a sense, both are subsumed by multiple dispatch / generic functions – multiple dispatch is the object-oriented side of the coin while generic functions are the functional side of the coin.

So is this at all helpful for Sage? I'm not sure. The things that you need to make multiple dispatch really useful for something as basic as integer-float addition are:

1. The + operator as a generic function – otherwise the multiple dispatch can't come into play.
   
2. Enough type information and inlining to allow + to still be sufficiently fast.

I'm not sure if either of these are possible in Sage. The speed requirement can be alleviated by only doing multiple dispatch on user-defined types where + doesn't need to be so fast. Alternatively, you can limit operations that are less performance critical than +, such as the EllipticCurve constructor in Nils' example. However, multiple dispatch should not be dismissed out of hand since it can provide elegant solutions to type promotion and mixed-type operations (and many other issues), and it can be made very fast, although this is admittedly not trivial to do.

[Nils Bruin]

On Friday, July 18, 2014 12:59:37 PM UTC-7, Robert Bradshaw wrote:

Yes, but here T must explicitly implement to multiple dispatch. One
can't take an existing object and "add a method" based on the type of
the second argument.

Do you mean that

a.bar(b)
a.bar(c)

would dispatch to different methods

a.bar_B(b)
a.bar_C(c)

based on properties (probably type) of b,c?
 
You could do that via exactly the same mechanism. You may have to jump through some hoops to get the right binding behaviour. And of course, inheritance of "method" on subclasses of type(a) might get messy.

Use could be spelled similar to

class foo(...):
    bar = binding_multimethod()

    @bar.register(int,float)
    def i_f(self, n, epsilon):
        ...
    @bar.register(int,int)
    def i_i(self,n,m):
        ...
 
and monkeypatching would be straightforward too:

@foo.bar.register(float,float)
def f_f(self,epsilon,delta):
    ...

and could be desirable if for some reason the code for f_f is best grouped according to another parameter.

From a language design point of view, you'd get a strange mix, because you'd be using the "a.bar" spelling, which is purpose-made for single dispatch, and mix it with another dispatch mechanism: if you have a more general dispatch mechanism, why not use it for the first argument as well? For integration with the rest of python it might be appropriate to go with the hybrid, though. It would be the kind of thing you can just add to an existing class.

We also have @coerce_binop
https://github.com/sagemath/sage/blob/6.2/src/sage/structure/element.pyx#L3260
which is nice.

Yes, that does a good job for what it's for. However, not all multiple dispatch problems are coercion problems, and they're not all for binops.
 

Here we're conflating the issue of types and Parents
(where the latter can be viewed, somewhat, as dynamically created
parameterized types),

You don't have to do multiple dispatch on type. In a statically typed language you would want to do that, though, because that's the only case you may be able to optimize at compile-time. All the other ones can only try to do some runtime work to try the most likely candidates first.

The most obvious place where I think multimethods might have a place in sage at the moment is in constructors/factories and conversion calls (which are almost spelled as constructors anyway).

[Luca De Feo]

> I'm not sure multimethods alone are enough to solve issues with Sage's
> type system (e.g. coercion, where the result of a+b may happen in a
> completely new domain) but they could certainly help.

I purposedly said "many", and not "all", nor even "most". The most
important feature to a CAS missing in Julia is multiple inheritance,
in my opinion.

However, Julia multimethods are backed up by a powerful coercion
system, so I do not understand the "step back" criticism.

Luca

[Nils Bruin]

On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:

However, Julia multimethods are backed up by a powerful coercion
system, so I do not understand the "step back" criticism.

That comment wasn't made with respect to Julia, because that would be comparing the coercion facilities of a CAS to those of a programming language. Coercion in a CAS tends to be a *lot* more complicated than what programming languages are designed for. As an example:

Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in F_q[[x]] (finite field with q elements).

Do you expect your CAS to make sense of that addition? Sage does. It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power series in x over F_q) . You can argue whether it's desirable for a system to try to be that smart, but all computer algebra systems I know are a "step back" relative to this. Programming languages do not tend to have type models that would even allow you to try and make sense of this kind of question.
 

[Stefan Karpinski]

On Sat, Jul 19, 2014 at 8:22 AM, Nils Bruin <nbr...@sfu.ca> wrote:

Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in F_q[[x]] (finite field with q elements).

Do you expect your CAS to make sense of that addition? Sage does. It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power series in x over F_q) . You can argue whether it's desirable for a system to try to be that smart, but all computer algebra systems I know are a "step back" relative to this. Programming languages do not tend to have type models that would even allow you to try and make sense of this kind of question.

This can be pretty straightforwardly handled with multiple dispatch promotion in a parametric type system. Similar things are done in real code regularly in Julia. Although I'm not going to try to define these particular types, the polynomial, power series, and finite field example is quite possible given appropriate type definitions – people have done similar things. An example that already works and is similar in spirit, is figuring out that a Complex{BigInt} plus a Vector{Rational{Int}} should be a Vector{Complex{Rational{BigInt}}}. This falls out of a surprisingly small set of generic function methods and promotion rules (which are just methods):

julia> big(3)^100*im + [1, 1//2, 2]

3-element Array{Complex{Rational{BigInt}},1}:

 1//1+515377520732011331036461129765621272702107522001//1*im

 1//2+515377520732011331036461129765621272702107522001//1*im

 2//1+515377520732011331036461129765621272702107522001//1*im

It doesn't quite get you to a CAS (nor is it meant to), but it's getting close. I suspect this sort of thing can also be done in Haskell (although the lack of dependent types may make defining the desired types tough) and can almost certainly be done in a system like Idris. Julia kind of cheats here by being dynamic – having runtime values as type parameters is easy if you don't try to do type checking in advance.

In any case, I'm mainly just advocating that you don't dismiss multiple dispatch out of hand as a technique for expressing complex, composable, user-extensible promotion rules. You can do this sort of thing quite nicely with multiple dispatch and it might be helpful in Sage. I'm quite ignorant on this point, however, since I don't understand how Sage currently does things and what the pain points are.

[Francesco Biscani]

Hello Stefan,

thanks for the explanations, this all looks really interesting to me.

On 19 July 2014 22:10, Stefan Karpinski <ste...@karpinski.org> wrote:

This can be pretty straightforwardly handled with multiple dispatch promotion in a parametric type system. Similar things are done in real code regularly in Julia. Although I'm not going to try to define these particular types, the polynomial, power series, and finite field example is quite possible given appropriate type definitions – people have done similar things. An example that already works and is similar in spirit, is figuring out that a Complex{BigInt} plus a Vector{Rational{Int}} should be a Vector{Complex{Rational{BigInt}}}. This falls out of a surprisingly small set of generic function methods and promotion rules (which are just methods):

julia> big(3)^100*im + [1, 1//2, 2]

3-element Array{Complex{Rational{BigInt}},1}:

 1//1+515377520732011331036461129765621272702107522001//1*im

 1//2+515377520732011331036461129765621272702107522001//1*im

 2//1+515377520732011331036461129765621272702107522001//1*im

It's funny because I have been thinking about this for quite a long time from a modern C++ perspective, and I had come to the conclusion that with constructs like std::enable_if you can get to these results with template metaprogramming in C++11. The tricky bit is to write arithmetic operators that are generic but do not "swallow up" every type. In this specific case, you would have the Vector's operator+ kick in instead of the Complex's one (somewhere the Vector's operator+ should be defined to have a "higher rank" than the operator+ of Complex), then you would need to deduce the common type from Complex{BigInt} + Complex{Rational{...}} (possibly via decltype()), and proceed to a recursion down in the type hierarchy. Of course as you point out you would need all the ancillary functions to actually perform the conversions from one type to the other.

All this would entail quite a bit of template trickery, but I find it interesting that it could be implemented all at compile time.

(I have a primitive prototype version of something like this implemented in a project of mine, hopefully eventually I will get around to do it properly :)

Anyway, cheers for the interesting discussion.

  Francesco.

[Stefan Karpinski]

Yes, you can definitely do most of this with template metaprogramming, although, as you say it's pretty tricky. I've found that runtime generics + dispatch are easier to work with, but I don't think there's much difference in power – except in cases where you want to write code that needs to actually do the dispatch at runtime, in which case you need to resort to double dispatch or something like that in C++.

--

[Nils Bruin]

On Friday, July 18, 2014 4:01:10 PM UTC-7, Stefan Karpinski wrote:

If Julia has shown anything, it's that you *can* have ubiquitous multiple dispatch in a dynamic language and have it be very fast – it is commonly a zero-cost abstraction in Julia.

I find this an interesting claim. Do you have some easy, high level pointers to what ideas are used to get to zero-cost? (I realize you have a "commonly" disclaimer here, so perhaps I'm hoping for too much).

Consider the following example:

F(a,b,c)

I could see how one can dispatch *very* quickly based on exact types: make a hash of ("F",type(a),type(b),type(c)) and look up in a hash-table. Other than the increased cost of creating the hash key, the work required here is about the same as what python has to do every time anyway (look up the binding of "F" which, depending on the context, may just be as bad as doing a hash table lookup)

However, dispatch on exact types is virtually useless in a setting that uses inheritance to any significant degree. You need to take inheritance into account for type dispatch. That means that signature lookup has to happen with respect to a partial ordering: We want the most specific signature that generalizes the types of the arguments--if that's not unique we have an ambiguity. It seems to me the complexity of this will always depend on the relevant part of the type inheritance tree. How do you get from there to zero cost, or at least cost independent of the complexity of the type tree? By hoping that type inference gets you down to a small enough corner of the type tree?

[Stefan Karpinski]

On Mon, Jul 21, 2014 at 11:50 AM, Nils Bruin <nbr...@sfu.ca> wrote:

On Friday, July 18, 2014 4:01:10 PM UTC-7, Stefan Karpinski wrote:

If Julia has shown anything, it's that you *can* have ubiquitous multiple dispatch in a dynamic language and have it be very fast – it is commonly a zero-cost abstraction in Julia.

I find this an interesting claim. Do you have some easy, high level pointers to what ideas are used to get to zero-cost? (I realize you have a "commonly" disclaimer here, so perhaps I'm hoping for too much).

 

Yes, that's definitely *the* question. You certainly can't do dispatch at runtime – even a hash table lookup is far too slow to implement something as basic as +, and as you point out, hashing on exact argument types doesn't solve the problem anyway. The answer is a combination of techniques, which may or may not be applicable to Sage – please bear with me if they aren't. The keys techniques are:

1. code specialization
2. type inference
   
3. inlining
   

If you're evaluating F(a,b,c), you can generate a version of F that is specialized to the exact types of a, b and c and then cache that. When generating the specialized version of F for those types, you can use the known types of its arguments to infer as much as possible about the types of values in the code (more on that below). When you know exactly what method to call, you can potentially inline that call, often through many layers of calls, which is how all the abstraction overhead is eliminated. There is still one slow lookup at the top level that either generates specialized code or finds it pre-generated in a method cache, but after that, there should be few slow lookups as long as the type inference works.

The algorithm Julia uses for type inference works by walking through a program, starting with the types of its input values, and abstractly interpreting the code. Instead of applying the code to values, it applies the code to types, following all branches concurrently and tracking a superposition of all the possible states the program could be in, including all the types each expression could assume. (There's a certain similarity to the execution of a non-deterministic finite automaton or the process of quantum computation.) Since programs can iterate arbitrarily long and recurse to arbitrary depths, this process must be iterated until it reaches a fixed point. Once that point is reached, the result is that the program is annotated with upper bounds on the types that each variable and expression can assume. You can often get really good type information – we get concrete types for a majority of expressions, but that Julia's standard library is carefully designed to help here, and expressing polymorphism with dispatch on types is also helpful. Even when you can't get precise type information, you can sometimes cut down on the number of things that need to be checked to determine which method ought to be called at each call site, potentially reducing it to a few fast checks. You can find more details here: http://arxiv.org/abs/1209.5145.

One thing that's crucial for type promotion in such a system is that the type inference algorithm be able to fully analyze the promotion code. Otherwise when you define +(x::Number, y::Number) = +(promote(x,y)...) there's no hope that the promote call can get inlined, and without that call getting inlined, this whole thing can't boil down to just a few instructions. But as long as type inference can analyze the promote function, then you can get it to inline – not just for + but for any operation that is simple enough, even if it's user-defined. There's a certain brittleness here since a minor change to the promotion system can suddenly make everything very slow. But this brittleness isn't exposed to users: they can safely use the promotion system any way they care to; only *changing* the promotion system (or type inference) is fragile – and we don't do that very often.

I hope that all makes some sense. I'm not sure if doing this kind of code specialization in Sage is realistic, but type inference on Python code is certainly possible. As you've pointed out, there is no obvious way to make a completely generic version of code using multiple dispatch really fast (although polymorphic inline caches might help a lot).

[Nils Bruin]

On Monday, July 21, 2014 1:51:37 PM UTC-7, Stefan Karpinski wrote:

I hope that all makes some sense. I'm not sure if doing this kind of code specialization in Sage is realistic, but type inference on Python code is certainly possible. As you've pointed out, there is no obvious way to make a completely generic version of code using multiple dispatch really fast (although polymorphic inline caches might help a lot).

Thanks! Yes, it does.  It seems to me that the sage coercion system is a little too squishy (I suspect this can be made into a technical term, but I don't have a definition handy right now) to be made to fit. For one thing, it's presently not mainly driven by python "types" but by "parents". This allows a separation between "mathematical meaning" and "implementation" that is sometimes convenient. I'm not sure how easy it would be to fully map "parents" onto a (parametrized) type hierarchy. Axiom definitely tried to do that, but I haven't see Axiom applied in areas where I have seen Sage and Magma applied, so I don't know how effective it would be.

Type "promotion" code in sage can be pretty much anything and can change fairly dynamically. So any optimizations along the lines you're suggesting would only be possible at a JIT compiler level, where any compiled code caches would probably have to be invalidated when the type/parent/coercion graph gets modified. It certainly seems unlikely we'd be able to utilize these ideas without significantly altering python/cython. I guess "Sulia" or "Juxiom" might be the next generation computer algebra system.

[Stefan Karpinski]

On Mon, Jul 21, 2014 at 2:12 PM, Nils Bruin <nbr...@sfu.ca> wrote:

It certainly seems unlikely we'd be able to utilize these ideas without significantly altering python/cython. I guess "Sulia" or "Juxiom" might be the next generation computer algebra system.

Glad the explanation was helpful, even if immediately not applicable. The idea of a CAS based on Julia is pretty interesting and might be really fruitful. We've got plenty to chew on with numerical work, so that will have to remain in the future for now.

So any optimizations along the lines you're suggesting would only be possible at a JIT compiler level, where any compiled code caches would probably have to be invalidated when the type/parent/coercion graph gets modified.

We actually still have this problem in Julia – cached generated code can get out of sync when new method definitions are introduced. It's not usually a big problem, but it can be annoying in interactive work.

[Fernando Perez]

On Sun, Jul 20, 2014 at 8:43 PM, Stefan Karpinski <ste...@karpinski.org> wrote:

Yes, you can definitely do most of this with template metaprogramming, although, as you say it's pretty tricky.

Some people, when confronted with a problem, think "I know, I'll use template metaprogramming." Now they have a type-parametric problem factory...

Sorry, couldn't resist ;) Too many memories of C++ template error messages...

Awesome discussion, BTW! Thanks Stefan for the careful explanations, much appreciated.

--
Fernando Perez (@fperez_org; http://fperez.org)
fperez.net-at-gmail: mailing lists only (I ignore this when swamped!)
fernando.perez-at-berkeley: contact me here for any direct mail

[rjf]

On Saturday, July 19, 2014 8:22:39 AM UTC-7, Nils Bruin wrote:

On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:

However, Julia multimethods are backed up by a powerful coercion
system, so I do not understand the "step back" criticism.

That comment wasn't made with respect to Julia, because that would be comparing the coercion facilities of a CAS to those of a programming language. Coercion in a CAS tends to be a *lot* more complicated than what programming languages are designed for. As an example:

Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in F_q[[x]] (finite field with q elements).

Do you expect your CAS to make sense of that addition? Sage does.

A  CAS that has representations for those two objects will very likely make
sense of that addition, so Sage is hardly unique.

It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power series in x over F_q) . You can argue whether it's desirable for a system to try to be that smart, but all computer algebra systems I know are a "step back" relative to this. Programming languages do not tend to have type models that would even allow you to try and make sense of this kind of question.

A harder question is when the coercion is not so obvious,  As a simple example, is
the polynomial x^0 coerced to the integer 1?   How do you add two bigfloats of different precision?  Or add a float to an exact rational?  Do you allow 1/(1+i) or do you coerce by rationalizing the denominator?
The speed of dispatch and specialization has probably been explored in the context of compiling Lisp, ad nauseum.  I don't know if Julia adds something to the mix, but I am skeptical that it has so much to do with the title of this thread  (scientific computing).  While Julia might be a superior vehicle for writing scientifc computing programs for any number of reasons, the issues it addresses may be irrelevant to the applications of computers in scientific disciplines, which tend to be written in stuff like Matlab, (or clones), C, or FORTRAN.  It is not necessarily a technology problem as much as a marketing problem.

From personal experience I can attest to the fact that showing someone  P who is relatively happy with programming language X that there is a "better" language Y for what he is doing, is  not sufficient to convince P to rewrite all his programs and those of his friends/students in language Y.  Usually.   And if you DO make such a case, you find that the new programs in language Y use almost none of the features that make Y better than X.  They are just programs transliterated into Y.   Like FORTRAN programs, GOTOs and all, written in Lisp...

But have fun, anyway.

 

[Robert Bradshaw]

On Wed, Jul 23, 2014 at 5:47 PM, rjf <fat...@gmail.com> wrote:
>
>
> On Saturday, July 19, 2014 8:22:39 AM UTC-7, Nils Bruin wrote:
>>
>> On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:
>>>
>>> However, Julia multimethods are backed up by a powerful coercion
>>> system, so I do not understand the "step back" criticism.
>>>
>> That comment wasn't made with respect to Julia, because that would be
>> comparing the coercion facilities of a CAS to those of a programming
>> language. Coercion in a CAS tends to be a *lot* more complicated than what
>> programming languages are designed for. As an example:
>>
>> Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in
>> F_q[[x]] (finite field with q elements).
>>
>> Do you expect your CAS to make sense of that addition? Sage does.
>
> A CAS that has representations for those two objects will very likely make
> sense of that addition, so Sage is hardly unique.

Show me one.

>> It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power
>> series in x over F_q) . You can argue whether it's desirable for a system to
>> try to be that smart, but all computer algebra systems I know are a "step
>> back" relative to this. Programming languages do not tend to have type
>> models that would even allow you to try and make sense of this kind of
>> question.
>
>
> A harder question is when the coercion is not so obvious, As a simple
> example, is
> the polynomial x^0 coerced to the integer 1?

The *polynomial* x^0, yes.

> How do you add two bigfloats of different precision?

Return the result in lower precision.

> Or add a float to an exact rational?

Coerce to a rational (same rule as above).

> Do you allow 1/(1+i) or do you coerce by rationalizing the denominator?

That depends on what you're going to do with it. Viewed as elements of
C, or Q[i], yes, rationalize the denominator.

> The speed of dispatch and specialization has probably been explored in the
> context of compiling Lisp, ad nauseum.

Which is one of the reasons it didn't take decades to get to where it is.

> I don't know if Julia adds something
> to the mix, but I am skeptical that it has so much to do with the title of
> this thread (scientific computing). While Julia might be a superior
> vehicle for writing scientifc computing programs for any number of reasons,
> the issues it addresses may be irrelevant to the applications of computers
> in scientific disciplines, which tend to be written in stuff like Matlab,
> (or clones), C, or FORTRAN. It is not necessarily a technology problem as
> much as a marketing problem.
>
> From personal experience I can attest to the fact that showing someone P
> who is relatively happy with programming language X that there is a "better"
> language Y for what he is doing, is not sufficient to convince P to rewrite
> all his programs and those of his friends/students in language Y. Usually.

Often P isn't happy with language X. Or, in my personal case, P is
interested in languages in general (especially so called "goldilocks
languages").

> And if you DO make such a case, you find that the new programs in language Y
> use almost none of the features that make Y better than X. They are just
> programs transliterated into Y. Like FORTRAN programs, GOTOs and all,
> written in Lisp...

Though I've seen this, that hasn't been my experience in general.

> But have fun, anyway.

Oh, we will :).

- Robert

[John Cremona]

On 24 July 2014 04:22, Robert Bradshaw <robe...@math.washington.edu> wrote:
> On Wed, Jul 23, 2014 at 5:47 PM, rjf <fat...@gmail.com> wrote:
>>
>>
>> On Saturday, July 19, 2014 8:22:39 AM UTC-7, Nils Bruin wrote:
>>>
>>> On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:
>>>>
>>>> However, Julia multimethods are backed up by a powerful coercion
>>>> system, so I do not understand the "step back" criticism.
>>>>
>>> That comment wasn't made with respect to Julia, because that would be
>>> comparing the coercion facilities of a CAS to those of a programming
>>> language. Coercion in a CAS tends to be a *lot* more complicated than what
>>> programming languages are designed for. As an example:
>>>
>>> Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in
>>> F_q[[x]] (finite field with q elements).
>>>
>>> Do you expect your CAS to make sense of that addition? Sage does.
>>
>> A CAS that has representations for those two objects will very likely make
>> sense of that addition, so Sage is hardly unique.
>
> Show me one.

Not Magma, anyway:

Magma V2.20-6 Thu Jul 24 2014 09:18:47 on hilbert [Seed = 4286344687]
Type ? for help. Type <Ctrl>-D to quit.
> R1<x,y> := PolynomialRing(Integers(),2);
> f := x+y^2+5;
> f;
x + y^2 + 5
> q := 125;
> Fq := GF(q);
> R2<x> := PowerSeriesRing(Fq);
> g := R2!(3+x+x^4);
> g;
3 + x + x^4
> f+g;

>> f+g;
^
Runtime error in '+': Bad argument types
Argument types given: RngMPolElt, RngSerPowElt[FldFin]

>
>>> It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power
>>> series in x over F_q) . You can argue whether it's desirable for a system to
>>> try to be that smart, but all computer algebra systems I know are a "step
>>> back" relative to this. Programming languages do not tend to have type
>>> models that would even allow you to try and make sense of this kind of
>>> question.
>>
>>

[rjf]

On Wednesday, July 23, 2014 8:22:39 PM UTC-7, Robert Bradshaw wrote:

On Wed, Jul 23, 2014 at 5:47 PM, rjf <fat...@gmail.com> wrote:
>
>
> On Saturday, July 19, 2014 8:22:39 AM UTC-7, Nils Bruin wrote:
>>
>> On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:
>>>
>>> However, Julia multimethods are backed up by a powerful coercion
>>> system, so I do not understand the "step back" criticism.
>>>
>> That comment wasn't made with respect to Julia, because that would be
>> comparing the coercion facilities of a CAS to those of a programming
>> language. Coercion in a CAS tends to be a *lot* more complicated than what
>> programming languages are designed for. As an example:
>>
>> Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series in
>> F_q[[x]] (finite field with q elements).
>>
>> Do you expect your CAS to make sense of that addition? Sage does.
>
> A  CAS that has representations for those two objects will very likely make
> sense of that addition, so Sage is hardly unique.

Show me one.

Certainly Macsyma has been able to add taylor series and polynomials for, oh, 40 years.

If it had taylor series over F_q, it would have to make sense of that addition as well.

Perhaps Axiom and Fricas have such odd Taylor series objects;  they don't seem to be

something that comes up much, and I would expect they have some uncomfortable

properties.

Note that even adding  5 mod 13   and the integer 10  is potentially uncomfortable,

and the rather common operation of Hensel lifting requires doing arithmetic in a combination

of fields (or rings) of different related sizes.

 

>> It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power
>> series in x over F_q) . You can argue whether it's desirable for a system to
>> try to be that smart, but all computer algebra systems I know are a "step
>> back" relative to this. Programming languages do not tend to have type
>> models that would even allow you to try and make sense of this kind of
>> question.
>
>
> A harder question is when the coercion is not so obvious,  As a simple
> example, is
> the polynomial x^0 coerced to the integer 1?

The *polynomial* x^0, yes.

That means that 1 is not a polynomial. Adds to some checking, if you assume that

a polynomial has a main variable, but 1  does not.

As does the representation of zero as a polynomial with no terms, or a polynomial

with one term, say 0*x^0.

 

> How do you add two bigfloats of different precision?

Return the result in lower precision.

I would say that is wrong.  Macsyma converts both to the globally set precision.

There is no reason to believe that a "low precision" representation is

inaccurate, and that one can raise the precision by adding (binary) zeros on the right.

Or not.

 

> Or add a float to an exact rational?

Coerce to a rational (same rule as above).

That's not the rule used in most programming languages, where the result is

returned as a float.  Like FORTRAN did, I guess.  Any float can be converted

to an exact rational, but do you want to do that?  Sometimes.  That's

why it is sticky.

 

> Do you allow 1/(1+i) or do you coerce by rationalizing the denominator?

That depends on what you're going to do with it. Viewed as elements of
C, or Q[i], yes, rationalize the denominator.

So now you are going to require users to have a background in, what,

modern algebra? 

> The speed of dispatch and specialization has probably been explored in the
> context of compiling Lisp, ad nauseum.

Which is one of the reasons it didn't take decades to get to where it is.

I can't parse this.  Are you being sarcastic and saying it took Common Lisp

decades to compile object-oriented code?   I don't think this is true, generally.

But maybe I don't understand. 

[William Stein]

On Wed, Jul 30, 2014 at 9:11 PM, rjf <fat...@gmail.com> wrote:
>
>
> On Wednesday, July 23, 2014 8:22:39 PM UTC-7, Robert Bradshaw wrote:
>>
>> On Wed, Jul 23, 2014 at 5:47 PM, rjf <fat...@gmail.com> wrote:
>> >
>> >
>> > On Saturday, July 19, 2014 8:22:39 AM UTC-7, Nils Bruin wrote:
>> >>
>> >> On Saturday, July 19, 2014 5:43:57 AM UTC-7, defeo wrote:
>> >>>
>> >>> However, Julia multimethods are backed up by a powerful coercion
>> >>> system, so I do not understand the "step back" criticism.
>> >>>
>> >> That comment wasn't made with respect to Julia, because that would be
>> >> comparing the coercion facilities of a CAS to those of a programming
>> >> language. Coercion in a CAS tends to be a *lot* more complicated than
>> >> what
>> >> programming languages are designed for. As an example:
>> >>
>> >> Consider A+B where A is a polynomial in ZZ[x,y] and B is a power series
>> >> in
>> >> F_q[[x]] (finite field with q elements).
>> >>
>> >> Do you expect your CAS to make sense of that addition? Sage does.
>> >
>> > A CAS that has representations for those two objects will very likely
>> > make
>> > sense of that addition, so Sage is hardly unique.
>>
>> Show me one.
>
>
> Certainly Macsyma has been able to add taylor series and polynomials for,
> oh, 40 years.
> If it had taylor series over F_q, it would have to make sense of that
> addition as well.

How do you define an element of F_q in Maxima? I couldn't find
anything after a cursory look at [1]. I'm trying to figure out if
Maxima counts as "a CAS that has representations for those two
objects".

[1] http://maxima.sourceforge.net/docs/manual/en/maxima_29.html

> Perhaps Axiom and Fricas have such odd Taylor series objects; they don't
> seem to be
> something that comes up much, and I would expect they have some
> uncomfortable
> properties.

They do come up frequently in algebraic geometry, number theory,
representation theory, combinatorics, coding theory and many other
areas of pure and applied mathematics.

> Note that even adding 5 mod 13 and the integer 10 is potentially
> uncomfortable,

sage: a = Mod(5,13); a
5
sage: parent(a)
Ring of integers modulo 13
sage: type(a)
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
sage: a + 10
2
sage: parent(a+10)
Ring of integers modulo 13

That was really comfortable. It's basically the same in Magma and
GAP. In PARI it's equally comfortable.

? Mod(5,13)+10
%2 = Mod(2, 13)

> and the rather common operation of Hensel lifting requires doing arithmetic
> in a combination
> of fields (or rings) of different related sizes.
>
>>
>>
>> >> It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power
>> >> series in x over F_q) . You can argue whether it's desirable for a
>> >> system to
>> >> try to be that smart, but all computer algebra systems I know are a
>> >> "step
>> >> back" relative to this. Programming languages do not tend to have type
>> >> models that would even allow you to try and make sense of this kind of
>> >> question.
>> >
>> >
>> > A harder question is when the coercion is not so obvious, As a simple
>> > example, is
>> > the polynomial x^0 coerced to the integer 1?
>>
>> The *polynomial* x^0, yes.
>
>
> That means that 1 is not a polynomial. Adds to some checking, if you assume
> that
> a polynomial has a main variable, but 1 does not.
> As does the representation of zero as a polynomial with no terms, or a
> polynomial
> with one term, say 0*x^0.

There are many 1's.

>> > How do you add two bigfloats of different precision?
>>
>> Return the result in lower precision.
>
>
> I would say that is wrong. Macsyma converts both to the globally set
> precision.

You can't be serious -- that's a total a recipe for disaster in a
complicated system. Unbelievable.

> There is no reason to believe that a "low precision" representation is
> inaccurate, and that one can raise the precision by adding (binary) zeros on
> the right.
>
> Or not.
>
>>
>>
>> > Or add a float to an exact rational?
>>
>> Coerce to a rational (same rule as above).
>
>
> That's not the rule used in most programming languages, where the result is
> returned as a float. Like FORTRAN did, I guess. Any float can be converted
> to an exact rational, but do you want to do that? Sometimes. That's
> why it is sticky.

It's not what Sage does. Robert wrote the wrong thing (though his
"see above" means he just got his words confused).
In Sage, one converts the rational to the parent ring of the float,
since that has lower precision.

sage: 1.3 + 2/3
1.96666666666667

Systematically, throughout the system we coerce from higher precision
(more info) naturally to lower precision. Another example is

Z ---> Z/13Z

If you add an exact integer ("high precision") to a number mod 13
("less information"), you get a number mod 13 (as above).

>
>>
>>
>> > Do you allow 1/(1+i) or do you coerce by rationalizing the denominator?
>>
>> That depends on what you're going to do with it. Viewed as elements of
>> C, or Q[i], yes, rationalize the denominator.
>
>
> So now you are going to require users to have a background in, what,
> modern algebra?

Here's what happens:

sage: 1/(1+i)
-1/2*I + 1/2

If a user does have a background in "modern" algebra (which was
developed in the 1800's by some French kid?), then they have more
options available. At least way, way more options than they will
ever have in your CAS (Macsyma).

Just because Sage is capable of doing X, and X is only something that
a user with a background in quasi-functional triads can understand,
does not mean that ever user of Sage is assumed to have a background
in quasi-functional triads.

William

--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

[rjf]

I think you have to look at the documentation for tellrat,  and maybe algebraic:true

There is a fairly serious questions as to whether one should distinguish between

3 mod 13  and  the integer 3,   or should one make a distinction between

addition mod 13  and  addition in Z.

 

Maxima tends toward not marking the objects, and puts the burden on the operation.

Axiom, and apparently Sage, requires that you mark the objects and define the operations

and coercions, and then let the type system do its job, or not.

[1] http://maxima.sourceforge.net/docs/manual/en/maxima_29.html

> Perhaps Axiom and Fricas have such odd Taylor series objects;  they don't
> seem to be
> something that comes up much, and I would expect they have some
> uncomfortable
> properties.

They do come up frequently in algebraic geometry, number theory,
representation theory, combinatorics, coding theory and many other
areas of pure and applied mathematics.

   You can say

whatever you want about some made-up computational problems

in  "pure mathematics" but I think you are just bluffing regarding applications.

> Note that even adding  5 mod 13   and the integer 10  is potentially
> uncomfortable,

sage: a = Mod(5,13); a
5
sage: parent(a)
Ring of integers modulo 13
sage: type(a)
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
sage: a + 10
2
sage: parent(a+10)
Ring of integers modulo 13

That was really comfortable.  

Unfortunately, it 

(a) looks uncomfortable, and (

b) got the answer wrong.

The sum of 5 mod13 and 10  is 15, not 2.

The calculation  (( 5 mod 13)+10) mod 13 is 2.

Note that the use of mod in the line above, twice, means different things.

One is a tag on the number 5 and the second is an operation akin to remainder.

 

  It's basically the same in Magma and
GAP.   In PARI it's equally comfortable.

? Mod(5,13)+10
%2 = Mod(2, 13)


> and the rather common operation of Hensel lifting requires doing arithmetic
> in a combination
> of fields (or rings) of different related sizes.

If you knew about Hensel lifting, I would think that you would comment here, and that

you also would know why 15 rather than 2 might be the useful answer. 

coercing the integer 10  into   something like 10 mod 13   or  perhaps -3 mod 13 is

a choice, probably the wrong one.

presumably one could promote your variable a to have a value in Z and do

the addition again. Whether comfortably or not.

 

>
>>
>>
>> >> It returns an answer in F_q[[x]][y] (i.e., a polynomial in y over power
>> >> series in x over F_q) . You can argue whether it's desirable for a
>> >> system to
>> >> try to be that smart, but all computer algebra systems I know are a
>> >> "step
>> >> back" relative to this. Programming languages do not tend to have type
>> >> models that would even allow you to try and make sense of this kind of
>> >> question.
>> >
>> >
>> > A harder question is when the coercion is not so obvious,  As a simple
>> > example, is
>> > the polynomial x^0 coerced to the integer 1?
>>
>> The *polynomial* x^0, yes.
>
>
> That means that 1 is not a polynomial. Adds to some checking, if you assume
> that
> a polynomial has a main variable, but 1  does not.
> As does the representation of zero as a polynomial with no terms, or a
> polynomial
> with one term, say 0*x^0.

There are many 1's.

that's a choice;  it may make good sense if you are dealing with a raft of algebraic structures

with various identities.

I'm not convinced it is a good choice if you are dealing with the bulk of applied analysis and

scientific computing applications amenable to symbolic manipulation.

Where one is one and all alone and ever more shall be so.

 

>> > How do you add two bigfloats of different precision?
>>
>> Return the result in lower precision.
>
>
> I would say that is wrong.  Macsyma converts both to the globally set
> precision.

You can't be serious -- that's a total a recipe for disaster in a
complicated system.  Unbelievable.

I am serious, and this method has been used since, oh, 1974 I would guess.

Of course you can easily rebind that precision, e.g.

do_careful_stuff(a,b):= block([fpprec:fpprec*2 ],    .... compute with a,b, at doubled precision)

I don't know why you think it is a recipe for disaster. Perhaps you can explain?

> There is no reason to believe that a "low precision" representation is
> inaccurate, and that one can raise the precision by adding (binary) zeros on
> the right.
>
> Or not.
>
>>
>>
>> > Or add a float to an exact rational?
>>
>> Coerce to a rational (same rule as above).
>
>
> That's not the rule used in most programming languages, where the result is
> returned as a float.  Like FORTRAN did, I guess.  Any float can be converted
> to an exact rational, but do you want to do that?  Sometimes.  That's
> why it is sticky.

It's not what Sage does.  Robert wrote the wrong thing (though his
"see above" means he just got his words confused).
In Sage, one converts the rational to the parent ring of the float,
since that has lower precision.

It seems like your understanding of "precision" is more than a little fuzzy.

You would not be alone in that, but writing programs for other people to

use sometimes means you need to know enough more to keep from

offering them holes to fall into.

Also, as you know (machine)  floats don't form a ring. 

Also, as you know some rationals cannot be converted to machine floats.

( or do you mean software floats with an unbounded exponent?  

Since these are roughly the same as rationals with power-of-two denominators / multipliers).

What is the precision of the parent ring of a float??  Rings with precision???

sage: 1.3 + 2/3
1.96666666666667

 

arguably, this is wrong.

1.3, if we presume this is a machine double-float in IEEE form, is

about

1.300000000000000044408920985006261616945266723633....

or EXACTLY

5854679515581645 / 4503599627370496

note that the denominator is 2^52.

adding 2/3 to that gives

26571237801485927 / 13510798882111488

EXACTLY.

which can be written

1.9666666666666667110755876516729282836119333903...

So you can either do "mathematical addition" by default or do

"addition yada yada".

Systematically, throughout the system we coerce from higher precision
(more info) naturally to lower precision.  Another example is

     Z   ---> Z/13Z

If you add an exact integer ("high precision") to a number mod 13
("less information"), you get a number mod 13 (as above).

This is a choice but it is hardly defensible.   

Here is an extremely accurate number 0.5

even though it can be represented in low precision, if I tell you it is

accurate to 100 decimal digits, that is independent of its representation.

If I write only 0.5,  does that mean that 0.5+0.04   = 0.5?  by your rule of

precision, the answer can only have 2 digits, the precision of 0.5, and so correctly rounded,

the answer is 0.5.  Tada.  [Seriously, some people do something very close to this.

Wolfram's Mathematica, using software floats. One of the easy ways of mocking it.]

I think that numbers mod 13  are perfectly precise, and have full information.

Now if you were representing integers modulo some huge modulus as

nearby floating-point numbers, I guess you would lose some information.

There is an excellent survey on What Every Computer Scientist Should Know About Floating Point...

by David Goldberg.    easily found via Google.

I recommend it, often.

>
>>
>>
>> > Do you allow 1/(1+i) or do you coerce by rationalizing the denominator?
>>
>> That depends on what you're going to do with it. Viewed as elements of
>> C, or Q[i], yes, rationalize the denominator.

In other words, your coercion fails you and you have to (indirectly) specify the

operation by reconsidering 

(a) the type of the numerator

(b) the type of the denominator

(c) the type of the target 

Is this kind of multi-inheritance supported by Python?  Sage?

>
>
> So now you are going to require users to have a background in, what,
> modern algebra?

Here's what happens:

sage: 1/(1+i)
-1/2*I + 1/2

If a user does have a background in "modern" algebra (which was
developed in the 1800's by some French kid?), then they have more
options available.   At least way, way more options than they will
ever have in your CAS (Macsyma).

I think that this approach was explored (and convincingly demonstrated to

be unworkable) in Axiom.  It was

also explored in this work..

http://www.eecs.berkeley.edu/Pubs/TechRpts/1983/CSD-83-160.pdf

When I said "what, modern algebra"  I was referring to the term used by

Birkhoff and MacLane, Survey of Modern Algebra. 

If you have a preferred name for this subject, please tell me.

 I think you may find

that parts of Maxima explicitly refer to it  or the earlier and in some ways

more constructive book by van der Waerden.  

Certainly the Berkeley project described in the report above

delves into the consequences of building software this way.

I think there is a choice to be made between making everything possible

but nothing easy,  and making a few common operations 

relatively easy, automatically  (but the less common -- more difficult)

Just because Sage is capable of doing X, and X is only something that
a user with a background in quasi-functional triads can understand,
does not mean that ever user of Sage is assumed to have a background
in quasi-functional triads.

No, but if you require that someone be facile with the construction of 

rings, fields, (etc) in order to do polynomial arithmetic, you have lost

a lot of users, and you also have a potentially exponential growth in coercions.

RJF

[William]

On Sat, Aug 2, 2014 at 11:16 PM, rjf <fat...@gmail.com> wrote:
> On Wednesday, July 30, 2014 10:23:03 PM UTC-7, William wrote:
>> [1] http://maxima.sourceforge.net/docs/manual/en/maxima_29.html
>>
>> > Perhaps Axiom and Fricas have such odd Taylor series objects;  they
>> > don't
>> > seem to be
>> > something that comes up much, and I would expect they have some
>> > uncomfortable
>> > properties.
>>
>> They do come up frequently in algebraic geometry, number theory,
>> representation theory, combinatorics, coding theory and many other
>> areas of pure and applied mathematics.
>
>
>    You can say
> whatever you want about some made-up computational problems
> in  "pure mathematics" but I think you are just bluffing regarding
> applications.

Let me get this straight -- you honestly thinking I'm bluffing regarding applications of:

     "power series in one variable over a finite field"

Are you seriously claiming that "power series over a finite field" have no applications?

>> > Note that even adding  5 mod 13   and the integer 10  is potentially
>> > uncomfortable,
>>
>> sage: a = Mod(5,13); a
>> 5
>> sage: parent(a)
>> Ring of integers modulo 13
>> sage: type(a)
>> <type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
>> sage: a + 10
>> 2
>> sage: parent(a+10)
>> Ring of integers modulo 13
>>
>> That was really comfortable.  
>
>

> Unfortunately, it   (a) looks uncomfortable, and (b) got the answer wrong.

>
> The sum of 5 mod13 and 10  is 15, not 2.
>
> The calculation  (( 5 mod 13)+10) mod 13 is 2.
>
> Note that the use of mod in the line above, twice, means different things.
> One is a tag on the number 5 and the second is an operation akin to
> remainder.

The answer from Sage is correct, and your remarks to the contrary are just typical of rather shallow understanding of mathematics.

>>   It's basically the same in Magma and
>> GAP.   In PARI it's equally comfortable.
>>
>> ? Mod(5,13)+10
>> %2 = Mod(2, 13)
>>
>>
>> > and the rather common operation of Hensel lifting requires doing
>> > arithmetic
>> > in a combination
>> > of fields (or rings) of different related sizes.
>
>
> If you knew about Hensel lifting, I would think that you would comment here,
> and that
> you also would know why 15 rather than 2 might be the useful answer.

The answer output by Pari is neither 15 nor 2.  It is "Mod(2, 13)".

> coercing the integer 10  into   something like 10 mod 13   or  perhaps -3
> mod 13 is
> a choice, probably the wrong one.

It is a completely canonical choice, and the only possible sensible one to make in this situation.  That's why every single math software system with any real algebraic sophistication has made this (same) choice.

> presumably one could promote your variable a to have a value in Z and do
> the addition again. Whether comfortably or not.

sage: a = Mod(5,13); a
5

sage: a + 10
2

sage: b = a.lift(); parent(b)
Integer Ring
sage: b + 10
15

>> There are many 1's.
>
>
> that's a choice;  it may make good sense if you are dealing with a raft of
> algebraic structures
> with various identities.

We are.

> I'm not convinced it is a good choice if you are dealing with the bulk of
> applied analysis and
> scientific computing applications amenable to symbolic manipulation.

Sage is definitely not only dealing with the that.  

>> It's not what Sage does.  Robert wrote the wrong thing (though his
>> "see above" means he just got his words confused).
>> In Sage, one converts the rational to the parent ring of the float,
>> since that has lower precision.
>
>
> It seems like your understanding of "precision" is more than a little fuzzy.
> You would not be alone in that, but writing programs for other people to
> use sometimes means you need to know enough more to keep from
> offering them holes to fall into.
>
> Also, as you know (machine)  floats don't form a ring.
>
> Also, as you know some rationals cannot be converted to machine floats.
>
> ( or do you mean software floats with an unbounded exponent?  
> Since these are roughly the same as rationals with power-of-two denominators
> / multipliers).

Yes, Sage uses http://www.mpfr.org/

> What is the precision of the parent ring of a float??  Rings with
> precision???

sage: R = parent(1.399390283048203480923490234092043820348); R
Real Field with 133 bits of precision
sage: R.precision()
133

>
>>
>> sage: 1.3 + 2/3
>> 1.96666666666667
>
>  
> arguably, this is wrong.
>
> 1.3, if we presume this is a machine double-float in IEEE form, is
> about
> 1.300000000000000044408920985006261616945266723633....
> or EXACTLY
> 5854679515581645 / 4503599627370496
>
> note that the denominator is 2^52.
>
> adding 2/3 to that gives
>
> 26571237801485927 / 13510798882111488
>
> EXACTLY.
>
>
> which can be written
> 1.9666666666666667110755876516729282836119333903...
>
> So you can either do "mathematical addition" by default or do
> "addition yada yada".
>

In Sage both operands are first converted to a common structure in which the operation can be performed, and then the operation is performed.

>
>>
>> Systematically, throughout the system we coerce from higher precision
>> (more info) naturally to lower precision.  Another example is
>>
>>      Z   ---> Z/13Z
>>
>> If you add an exact integer ("high precision") to a number mod 13
>> ("less information"), you get a number mod 13 (as above).
>
>
> This is a choice but it is hardly defensible.  
> Here is an extremely accurate number 0.5
> even though it can be represented in low precision, if I tell you it is
> accurate to 100 decimal digits, that is independent of its representation.
>
> If I write only 0.5,  does that mean that 0.5+0.04   = 0.5?  by your rule of
> precision, the answer can only have 2 digits, the precision of 0.5, and so
> correctly rounded,
> the answer is 0.5.  Tada.  [Seriously, some people do something very close
> to this.
> Wolfram's Mathematica, using software floats. One of the easy ways of
> mocking it.]

The string representation of a number in Sage is not the same thing as that number.    Here's an example of how to make the image of 1/2, but stored with very few bits of precision, then add 0.04 to it:

sage: RealField(2)(0.5) + 0.04
0.50

If one needs more precision tracking regarding exactly what happens when doing arithmetic (and using functions) with real or complex numbers, Sage also have efficient interval arithmetic, based on, e.g., http://gforge.inria.fr/projects/mpfi/. 

>
> I think that numbers mod 13  are perfectly precise, and have full
> information.
>
> Now if you were representing integers modulo some huge modulus as
> nearby floating-point numbers, I guess you would lose some information.
>
> There is an excellent survey on What Every Computer Scientist Should Know
> About Floating Point...
> by David Goldberg.    easily found via Google.
>
> I recommend it, often.

Paul Zimmerman writes many useful things about what mathematicians should know about floating point -- http://www.loria.fr/~zimmerma/papers/. 

William

[rjf]

On Sunday, August 3, 2014 10:11:57 PM UTC-7, William wrote:

On Sat, Aug 2, 2014 at 11:16 PM, rjf <fat...@gmail.com> wrote:
> On Wednesday, July 30, 2014 10:23:03 PM UTC-7, William wrote:
>> [1] http://maxima.sourceforge.net/docs/manual/en/maxima_29.html
>>
>> > Perhaps Axiom and Fricas have such odd Taylor series objects;  they
>> > don't
>> > seem to be
>> > something that comes up much, and I would expect they have some
>> > uncomfortable
>> > properties.
>>
>> They do come up frequently in algebraic geometry, number theory,
>> representation theory, combinatorics, coding theory and many other
>> areas of pure and applied mathematics.
>
>
>    You can say
> whatever you want about some made-up computational problems
> in  "pure mathematics" but I think you are just bluffing regarding
> applications.

Let me get this straight -- you honestly thinking I'm bluffing regarding applications of:

     "power series in one variable over a finite field"

Yes.

I read your statement above to say, among other things, that

{such odd power series objects} come up frequently in ... many area of ... applied mathematics.

 

Are you seriously claiming that "power series over a finite field" have no applications?

I can imagine you could invent an "application", so I would not say  these objects 

have "no applications".    I can even invent one application for you.  The application is

"show a structure that can be created by a computer program that consists of a power

series etc etc."   This sort of smacks of the self referentiality that comes from answering

the command:

Use X in a sentence.

by the response

"Use X in a sentence. 

>> > Note that even adding  5 mod 13   and the integer 10  is potentially
>> > uncomfortable,
>>
>> sage: a = Mod(5,13); a
>> 5
>> sage: parent(a)
>> Ring of integers modulo 13
>> sage: type(a)
>> <type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
>> sage: a + 10
>> 2
>> sage: parent(a+10)
>> Ring of integers modulo 13
>>
>> That was really comfortable.  
>
>
> Unfortunately, it   (a) looks uncomfortable, and (b) got the answer wrong.
>
> The sum of 5 mod13 and 10  is 15, not 2.
>
> The calculation  (( 5 mod 13)+10) mod 13 is 2.
>
> Note that the use of mod in the line above, twice, means different things.
> One is a tag on the number 5 and the second is an operation akin to
> remainder.

The answer from Sage is correct, and your remarks to the contrary are just typical of rather shallow understanding of mathematics.

Are you familiar with TeX's   \bmod ?

It is possible to define anything that comes from Sage as correct.

Mathematica does that kind of thing often, declaring its reported bugs to be features.  Sometimes

I agree, sometimes not.

>>   It's basically the same in Magma and
>> GAP.   In PARI it's equally comfortable.
>>
>> ? Mod(5,13)+10
>> %2 = Mod(2, 13)
>>
>>
>> > and the rather common operation of Hensel lifting requires doing
>> > arithmetic
>> > in a combination
>> > of fields (or rings) of different related sizes.
>
>
> If you knew about Hensel lifting, I would think that you would comment here,
> and that
> you also would know why 15 rather than 2 might be the useful answer.

The answer output by Pari is neither 15 nor 2.  It is "Mod(2, 13)".

Well, if you simply insist that arithmetic requires the conversion of 10 to a number in Z_13,

then that answer is OK.

 

> coercing the integer 10  into   something like 10 mod 13   or  perhaps -3
> mod 13 is
> a choice, probably the wrong one.

It is a completely canonical choice, and the only possible sensible one to make in this situation.

Which one, -3 or 10?

They can't both be canonical.

 

 That's why every single math software system with any real algebraic sophistication has made this (same) choice.

I think the issue  comes up as to the definition of a system with "real algebraic sophistication"  and for that

matter, whether you are familiar with all of them.

While you are probably aware that I am not a big fan of Mathematica, it is quite popular.  This is what it

has in its documentation.

Modulus->n

is an option that can be given in certain algebraic functions to specify that integers should be treated modulo n. 

It seems that Mathematica does not have built-in functionality to manipulate objects such as Mod(2,13) .

In fact, your whole premise seems to be that mathematically sophisticated <whatever>  must be based on

some notion of strongly-typed hierarchical  mathematical categories, an idea that has proven unpopular

with users, even if it has seemed to system builders, repeatedly and apparently incorrectly, as the golden fleece,

the universal solvent, etc.   

I think the historical record shows not only the marketing failure of Axiom, but also mod Reduce, and some subsystem

written in Maple.  Also the Newspeak system of a PhD student of mine, previously mentioned.  

Now you may think that (for example) Magma is just what you? the world?  needs (except only that it is not free).

That does not mean it should be used for scientific computing.

> presumably one could promote your variable a to have a value in Z and do
> the addition again. Whether comfortably or not.

sage: a = Mod(5,13); a
5
sage: a + 10
2
sage: b = a.lift(); parent(b)
Integer Ring
sage: b + 10
15

Thanks.  Can the lifting be done to another computational structure, say  modulo 13^(2^n)?

>> There are many 1's.
>
>
> that's a choice;  it may make good sense if you are dealing with a raft of
> algebraic structures
> with various identities. 

We are.

> I'm not convinced it is a good choice if you are dealing with the bulk of
> applied analysis and
> scientific computing applications amenable to symbolic manipulation.

Sage is definitely not only dealing with the that.  

Actually, I expect that Sage, as a front end to several systems, does not adequately

"deal" with ideas like "1".  Which might be used in any number of the subsystems

for concepts like matrix identity under multiplication.  Or a power series, or a polynomial...

 

>> It's not what Sage does.  Robert wrote the wrong thing (though his
>> "see above" means he just got his words confused).
>> In Sage, one converts the rational to the parent ring of the float,
>> since that has lower precision.
>
>
> It seems like your understanding of "precision" is more than a little fuzzy.
> You would not be alone in that, but writing programs for other people to
> use sometimes means you need to know enough more to keep from
> offering them holes to fall into.
>
> Also, as you know (machine)  floats don't form a ring.
>
> Also, as you know some rationals cannot be converted to machine floats.
>
> ( or do you mean software floats with an unbounded exponent?  
> Since these are roughly the same as rationals with power-of-two denominators
> / multipliers).

Yes, Sage uses http://www.mpfr.org/

While that solves the problem of overflow or underflow, it still does not allow

the representation of 1/3  as a (binary radix) float of some finite precision.

 

> What is the precision of the parent ring of a float??  Rings with
> precision???

sage: R = parent(1.399390283048203480923490234092043820348); R
Real Field with 133 bits of precision
sage: R.precision()
133

I'm confused here.  Are you now saying R is a Field?  Is a Real Field a kind of

Field?  Is there a multiplicative inverse in R for the object 3.0 ?  Or is a Real Field

like a "Dutch Oven"   which is actually not an oven, but a pot?

>
>>
>> sage: 1.3 + 2/3
>> 1.96666666666667
>
>  
> arguably, this is wrong.
>
> 1.3, if we presume this is a machine double-float in IEEE form, is
> about
> 1.300000000000000044408920985006261616945266723633....
> or EXACTLY
> 5854679515581645 / 4503599627370496
>
> note that the denominator is 2^52.
>
> adding 2/3 to that gives
>
> 26571237801485927 / 13510798882111488
>
> EXACTLY.
>
>
> which can be written
> 1.9666666666666667110755876516729282836119333903...
>
> So you can either do "mathematical addition" by default or do
> "addition yada yada".
>

In Sage both operands are first converted to a common structure in which the operation can be performed, and then the operation is performed.

Well, there are a large number of common structures. 

One that gets the exact answer is  to convert to arbitrary-precision rationals.

One that does not always get the exact answer is some kind of floating-point representation.

It is of course possible to convert these object to (say) polynomials or Taylor series..

>
>>
>> Systematically, throughout the system we coerce from higher precision
>> (more info) naturally to lower precision.  Another example is
>>
>>      Z   ---> Z/13Z
>>
>> If you add an exact integer ("high precision") to a number mod 13
>> ("less information"), you get a number mod 13 (as above).
>
>
> This is a choice but it is hardly defensible.  
> Here is an extremely accurate number 0.5
> even though it can be represented in low precision, if I tell you it is
> accurate to 100 decimal digits, that is independent of its representation.
>
> If I write only 0.5,  does that mean that 0.5+0.04   = 0.5?  by your rule of
> precision, the answer can only have 2 digits, the precision of 0.5, and so
> correctly rounded,
> the answer is 0.5.  Tada.  [Seriously, some people do something very close
> to this.
> Wolfram's Mathematica, using software floats. One of the easy ways of
> mocking it.]

The string representation of a number in Sage is not the same thing as that number.

I think that is kind of taken for granted as the relationship between input/output forms

and whatever is in the computer.  Electrons.  And even those are not the abstraction

"number".

 

   Here's an example of how to make the image of 1/2, but stored with very few bits of precision, then add 0.04 to it:

sage: RealField(2)(0.5) + 0.04
0.50

Hm.  Can one change the meaning of  "+"   so that returns  0.54?

What is the abstraction behind 0.04 ?   RealField(n)(0.04)   for some integer n?   

And is there a difference between 0.10  and 0.1000000000000000000000000000000000

(There is a difference in Mathematica.  Mathematica claims to be the best system for scientific computing. :)

 

If one needs more precision tracking regarding exactly what happens when doing arithmetic (and using functions) with real or complex numbers, Sage also have efficient interval arithmetic, based on, e.g., http://gforge.inria.fr/projects/mpfi/. 

I think that mpfi is probably a fine interval arithmetic package, but I think it

is not doing precision tracking, really. It's doing reliable (if pessimistic)

arithmetic.

 

>
> I think that numbers mod 13  are perfectly precise, and have full
> information.
>
> Now if you were representing integers modulo some huge modulus as
> nearby floating-point numbers, I guess you would lose some information.
>
> There is an excellent survey on What Every Computer Scientist Should Know
> About Floating Point...
> by David Goldberg.    easily found via Google.
>
> I recommend it, often.

Paul Zimmerman writes many useful things about what mathematicians should know about floating point -- http://www.loria.fr/~zimmerma/papers/. 

I have enjoyed reading some of Zimmerman's papers, and he and his authors have

said many interesting things. The Goldberg article tries to dispel commonly held

myths about machine floating-point arithmetic.  Myths often held by computer

programmers and of course mathematicians. Zimmerman has addressed many

fun problems in clever ways.

Inventing fast methods is fun, although (my opinion)  multiplying integers of astronomical size is hardly

mainstream scientific computing.  

Not to say that someone might claim that

this problem occurs frequently in many computations in pure and applied

mathematics...

RJF

[Robert Bradshaw]

Yes, they can "both" be canonical: they are equal (as elements of
Z/13Z). There is the choice of internal representation, and what
string to use when displaying them to the user, but that doesn't
affect its canonicity as a mathematical object.

>> That's why every single math software system with any real algebraic
>> sophistication has made this (same) choice.
>
>
> I think the issue comes up as to the definition of a system with "real
> algebraic sophistication" and for that
> matter, whether you are familiar with all of them.
>
> While you are probably aware that I am not a big fan of Mathematica, it is
> quite popular. This is what it
> has in its documentation.
>
> Modulus->n
> is an option that can be given in certain algebraic functions to specify
> that integers should be treated modulo n.
>
> It seems that Mathematica does not have built-in functionality to manipulate
> objects such as Mod(2,13) .

Nope. Or function fields. Or p-adics. Let alone more complicated
objects like points on algebraic varieties.

> In fact, your whole premise seems to be that mathematically sophisticated
> <whatever> must be based on
> some notion of strongly-typed hierarchical mathematical categories, an idea
> that has proven unpopular
> with users, even if it has seemed to system builders, repeatedly and
> apparently incorrectly, as the golden fleece,
> the universal solvent, etc.
>
> I think the historical record shows not only the marketing failure of Axiom,
> but also mod Reduce, and some subsystem
> written in Maple. Also the Newspeak system of a PhD student of mine,
> previously mentioned.
> Now you may think that (for example) Magma is just what you? the world?
> needs (except only that it is not free).
> That does not mean it should be used for scientific computing.

I can just see you sitting there wondering why Sage hasn't died the
tragic death your crystal ball said it would ages ago...

>> > presumably one could promote your variable a to have a value in Z and do
>> > the addition again. Whether comfortably or not.
>>
>> sage: a = Mod(5,13); a
>> 5
>> sage: a + 10
>> 2
>> sage: b = a.lift(); parent(b)
>> Integer Ring
>> sage: b + 10
>> 15
>
>
> Thanks. Can the lifting be done to another computational structure, say
> modulo 13^(2^n)?

Sure.

sage: b = Integers(13^1024)(a); b
5
sage: b^(4*13^1023)
1

R is a (very well studied) approximation to a field.

>> >> sage: 1.3 + 2/3
>> >> 1.96666666666667
>> >
>> >
>> > arguably, this is wrong.
>> >
>> > 1.3, if we presume this is a machine double-float in IEEE form, is
>> > about
>> > 1.300000000000000044408920985006261616945266723633....
>> > or EXACTLY
>> > 5854679515581645 / 4503599627370496
>> >
>> > note that the denominator is 2^52.
>> >
>> > adding 2/3 to that gives
>> >
>> > 26571237801485927 / 13510798882111488
>> >
>> > EXACTLY.
>> >
>> >
>> > which can be written
>> > 1.9666666666666667110755876516729282836119333903...
>> >
>> > So you can either do "mathematical addition" by default or do
>> > "addition yada yada".
>> >
>>
>> In Sage both operands are first converted to a common structure in which
>> the operation can be performed, and then the operation is performed.
>
>
> Well, there are a large number of common structures.
> One that gets the exact answer is to convert to arbitrary-precision
> rationals.
> One that does not always get the exact answer is some kind of floating-point
> representation.

It is more conservative to convert operands to the domain with less
precision. We consider the rationals to have infinite precision, our
real "fields" a specified have finite precision. This lower precision
is represented in the output, similar to how significant figures are
used in any other scientific endeavor.

This is similar to 10 + (5 mod 13), where the right hand side has
"less precision" (in particular there's a canonical map one way, many
choices of lift the other.

Also, when a user writes 1.3, they are more likely to mean 13/10 than
5854679515581645 / 4503599627370496, but by expressing it in decimal
form they are asking for a floating point approximation. Note that
real literals are handled specially to allow immediate conversion to a
higher-precision parent.

sage: QQ(1.3)
13/10
sage: (1.3).exact_rational()
5854679515581645/4503599627370496
sage: a = RealField(200)(1.3); a
1.3000000000000000000000000000000000000000000000000000000000
sage: a.exact_rational()
522254864384171839551137680010877845819715972979407671472947/401734511064747568885490523085290650630550748445698208825344

Yep, RealField(53)(0.04), which is a bit arbitrary.

> And is there a difference between 0.10 and
> 0.1000000000000000000000000000000000

Yep.

sage: a = 0.1000000000000000000000000000000000
sage: a.precision()
113

> (There is a difference in Mathematica. Mathematica claims to be the best
> system for scientific computing. :)

Few systems don't make that claim about themselves :).

I've personally "applied" multiplying astronomically sized before
(thought the result itself is squarely in the domain of pure math):
http://www.nsf.gov/news/news_summ.jsp?cntn_id=115646/

- Robert

[John Cremona]

>>
>> Are you seriously claiming that "power series over a finite field" have no
>> applications?
>
>
> I can imagine you could invent an "application", so I would not say these
> objects
> have "no applications". I can even invent one application for you. The
> application is
> "show a structure that can be created by a computer program that consists of
> a power
> series etc etc." This sort of smacks of the self referentiality that comes
> from answering
> the command:
>
> Use X in a sentence.
>
> by the response
> "Use X in a sentence.
>>

Power series over finite fields play a very important part in the
algorithms used in elliptic curve cryptography. That is a real world
application, and not one which has been invented for the sake of it.

John

[rjf]

I don't understand this.

I use the word canonical to mean  unique distinguished exemplar.

So there can't be two.  If there are two distinguishable items, then one or the other or neither might be canonical. Not both.

 

>>  That's why every single math software system with any real algebraic
>> sophistication has made this (same) choice.
>
>
> I think the issue  comes up as to the definition of a system with "real
> algebraic sophistication"  and for that
> matter, whether you are familiar with all of them.
>
> While you are probably aware that I am not a big fan of Mathematica, it is
> quite popular.  This is what it
> has in its documentation.
>
> Modulus->n
> is an option that can be given in certain algebraic functions to specify
> that integers should be treated modulo n.
>
> It seems that Mathematica does not have built-in functionality to manipulate
> objects such as Mod(2,13) .

Nope. Or function fields. Or p-adics. Let alone more complicated
objects like points on algebraic varieties.

OK, but of course any of these concepts could be added  -- or built-in in

a future version.   If someone cared.

 

> In fact, your whole premise seems to be that mathematically sophisticated
> <whatever>  must be based on
> some notion of strongly-typed hierarchical  mathematical categories, an idea
> that has proven unpopular
> with users, even if it has seemed to system builders, repeatedly and
> apparently incorrectly, as the golden fleece,
> the universal solvent, etc.
>
> I think the historical record shows not only the marketing failure of Axiom,
> but also mod Reduce, and some subsystem
> written in Maple.  Also the Newspeak system of a PhD student of mine,
> previously mentioned.
> Now you may think that (for example) Magma is just what you? the world?
> needs (except only that it is not free).
> That does not mean it should be used for scientific computing.

I can just see you sitting there wondering why Sage hasn't died the
tragic death your crystal ball said it would ages ago...

No, I think that Sage is fairly moribund.  Look around at the postings in sage-devel.

 

>> > presumably one could promote your variable a to have a value in Z and do
>> > the addition again. Whether comfortably or not.
>>
>> sage: a = Mod(5,13); a
>> 5
>> sage: a + 10
>> 2
>> sage: b = a.lift(); parent(b)
>> Integer Ring
>> sage: b + 10
>> 15
>
>
> Thanks.  Can the lifting be done to another computational structure, say
> modulo 13^(2^n)?

Sure.

sage: b = Integers(13^1024)(a); b
5
sage: b^(4*13^1023)
1

And what structure is that?  Does Sage know about   Z_{nonprime} ?

I'm still confused.   Is the term "Real Field" in Sage  the (or some)  real field?

If it is an approximation to a field, but not a field, why are you calling it a field?

If that doesn't get you in trouble, why doesn't it?  Does Real Field inherit from

Field?  Does every non-zero element have an inverse?  

Does Sage have other um, approximations, in its nomenclature?

Why do you say that?  You can always exactly convert a float number  in radix b to

an equal number of higher precision in radix b by appending zeros.

So it is more conserving (of values) to do so, rather than clipping off

bits from the other.

We consider the rationals to have infinite precision, our
real "fields" a specified have finite precision. This lower precision
is represented in the output, similar to how significant figures are
used in any other scientific endeavor.

Thanks for distinguishing between "field" and field.  You don't seem

to understand the concept of precision though. You seem to think

that a number of low precision has some inaccuracy or uncertainty.

Which it doesn't.   0.5 is the same number as 0.500000.

Unless you believe that 0.5  is  the interval [0.45000000000....01,  0.549999..........9]

which you COULD believe  -- some people do believe this.

But you probably don't want to ask them about scientific computing.

 

This is similar to 10 + (5 mod 13), where the right hand side has
"less precision" (in particular there's a canonical map one way, many
choices of lift the other.

Also, when a user writes 1.3, they are more likely to mean 13/10 than
5854679515581645 / 4503599627370496, but by expressing it in decimal
form they are asking for a floating point approximation. Note that
real literals are handled specially to allow immediate conversion to a
higher-precision parent.

What do you know when a user writes 1.3, really?  You want the user

to believe that Sage uses decimal arithmetic?  Seriously?  How far

are you going to try to carry that illusion?  

You imagine a user who understands rings and fields, and knows that

Real Field is not a field, but knows so little about computers that he

thinks 1.3  is 13/10 exactly?   (By the way, I have no problem  if

1.3 actually produces 13/10,  and to get a float, you have to try to

convert it explicitly, and the conversion might even come up with an interval

or an error bound or something that leads to "reliable" computing.

Rather than some guessy-how-many-bits stab in the dark thing that

prints as 1.3

sage: QQ(1.3)
13/10
sage: (1.3).exact_rational()
5854679515581645/4503599627370496
sage: a = RealField(200)(1.3); a
1.3000000000000000000000000000000000000000000000000000000000
sage: a.exact_rational()
522254864384171839551137680010877845819715972979407671472947/401734511064747568885490523085290650630550748445698208825344

I assume it is possible to calculate all kinds of things if you carefully specify them,

in Sage.  After all, it has all those programs to call, including sympy.  The

issue iwe have been discussing is really what does it do "automatically". 

So 0.1- 0.10000000000000000000000000000....

is 0.0?????????????????   where ? = undetermined?

and anyone who writes x+0.1  is relegating that sum to 1 digit "precision"? 

> (There is a difference in Mathematica.  Mathematica claims to be the best
> system for scientific computing. :)

Few systems don't make that claim about themselves :).

I've noticed.

 Assume there is an application that involves multiplication of polynomials.

You can multiply polynomials by encoding them as large integers,

multiplying, and decoding.  Sometimes called Kronecker's Trick.

So there are lots of applications.  Are they stupid tricks? Probably.

RJF

 

- Robert

[Robert Bradshaw]

On Wed, Aug 6, 2014 at 7:00 PM, rjf <fat...@gmail.com> wrote:
>> >> > coercing the integer 10 into something like 10 mod 13 or
>> >> > perhaps
>> >> > -3
>> >> > mod 13 is
>> >> > a choice, probably the wrong one.
>> >>
>> >> It is a completely canonical choice, and the only possible sensible one
>> >> to
>> >> make in this situation.
>> >
>> > Which one, -3 or 10?
>> > They can't both be canonical.
>>
>> Yes, they can "both" be canonical: they are equal (as elements of
>> Z/13Z). There is the choice of internal representation, and what
>> string to use when displaying them to the user, but that doesn't
>> affect its canonicity as a mathematical object.
>>
> I don't understand this.
> I use the word canonical to mean unique distinguished exemplar.
> So there can't be two. If there are two distinguishable items, then one or
> the other or neither might be canonical. Not both.

The are two representations of the same canonical object.

>> >> > presumably one could promote your variable a to have a value in Z and
>> >> > do
>> >> > the addition again. Whether comfortably or not.
>> >>
>> >> sage: a = Mod(5,13); a
>> >> 5
>> >> sage: a + 10
>> >> 2
>> >> sage: b = a.lift(); parent(b)
>> >> Integer Ring
>> >> sage: b + 10
>> >> 15
>> >
>> >
>> > Thanks. Can the lifting be done to another computational structure, say
>> > modulo 13^(2^n)?
>>
>> Sure.
>>
>> sage: b = Integers(13^1024)(a); b
>> 5
>> sage: b^(4*13^1023)
>> 1
>>
> And what structure is that? Does Sage know about Z_{nonprime} ?

Of course, as illustrated.

sage: Integers(13^1024)
Ring of integers modulo 4764...1

>> >> > What is the precision of the parent ring of a float?? Rings with
>> >> > precision???
>> >>
>> >> sage: R = parent(1.399390283048203480923490234092043820348); R
>> >> Real Field with 133 bits of precision
>> >> sage: R.precision()
>> >> 133
>> >>
>> > I'm confused here. Are you now saying R is a Field? Is a Real Field a
>> > kind
>> > of
>> > Field? Is there a multiplicative inverse in R for the object 3.0 ? Or
>> > is a
>> > Real Field
>> > like a "Dutch Oven" which is actually not an oven, but a pot?
>>
>> R is a (very well studied) approximation to a field.
>
>
> I'm still confused. Is the term "Real Field" in Sage the (or some) real
> field?
>
> If it is an approximation to a field, but not a field, why are you calling
> it a field?

Because it's shorter to type and easier to find/discover than
ApproximateRealField or something like that.

> If that doesn't get you in trouble, why doesn't it? Does Real Field inherit
> from
> Field? Does every non-zero element have an inverse?

Of course it suffers from the same issues that (standard) floating
point numbers do in any language, user be aware (and we at least
document that).

> Does Sage have other um, approximations, in its nomenclature?

Sure. RealField(123)[x]. Power series rings. P-adics.

Clipping bits (or digits) is exactly how one is taught to deal with
significant figures in grade school, and follows the principle of
least surprise (though floating point numbers like to play surprises
on you no matter what). It's also what floating point arithmetic does
when the exponent is different.

>> We consider the rationals to have infinite precision, our
>> real "fields" a specified have finite precision. This lower precision
>> is represented in the output, similar to how significant figures are
>> used in any other scientific endeavor.
>
> Thanks for distinguishing between "field" and field. You don't seem
> to understand the concept of precision though.

That's a bold claim. My Ph.D. thesis depended on understanding issues
of precision. I'll admit explaining it to a layman can be difficult.

> You seem to think
> that a number of low precision has some inaccuracy or uncertainty.
> Which it doesn't. 0.5 is the same number as 0.500000.
> Unless you believe that 0.5 is the interval [0.45000000000....01,
> 0.549999..........9]
> which you COULD believe -- some people do believe this.
> But you probably don't want to ask them about scientific computing.

No, I don't think that at all. Sage also has the concept of real
intervals distinct from real numbers.

There's 0.5 the real number equal to one divided by two. There's also
0.5 the IEEE floating point number, which is a representative for an
infinite number of real numbers in a small interval.

>> This is similar to 10 + (5 mod 13), where the right hand side has
>> "less precision" (in particular there's a canonical map one way, many
>> choices of lift the other.
>>
>> Also, when a user writes 1.3, they are more likely to mean 13/10 than
>> 5854679515581645 / 4503599627370496, but by expressing it in decimal
>> form they are asking for a floating point approximation. Note that
>> real literals are handled specially to allow immediate conversion to a
>> higher-precision parent.
>
> What do you know when a user writes 1.3, really?

Quite a bit. They probably didn't mean pi. If they really cared, they
could have been more specific. At least we recognize this ambiguity
but we can't let it paralyze us.

> You want the user
> to believe that Sage uses decimal arithmetic? Seriously? How far
> are you going to try to carry that illusion?

If they don't immediately specify a new domain, we'll treat it as
having 53 bits. It's syntactic sugar.

> You imagine a user who understands rings and fields, and knows that
> Real Field is not a field, but knows so little about computers that he
> thinks 1.3 is 13/10 exactly? (By the way, I have no problem if
> 1.3 actually produces 13/10, and to get a float, you have to try to
> convert it explicitly, and the conversion might even come up with an
> interval
> or an error bound or something that leads to "reliable" computing.
> Rather than some guessy-how-many-bits stab in the dark thing that
> prints as 1.3

Again, it's syntactic sugar.

>> sage: QQ(1.3)
>> 13/10
>> sage: (1.3).exact_rational()
>> 5854679515581645/4503599627370496
>> sage: a = RealField(200)(1.3); a
>> 1.3000000000000000000000000000000000000000000000000000000000
>> sage: a.exact_rational()
>> 522254864384171839551137680010877845819715972979407671472947/401734511064747568885490523085290650630550748445698208825344
>
> I assume it is possible to calculate all kinds of things if you carefully
> specify them,
> in Sage. After all, it has all those programs to call, including sympy.
> The
> issue iwe have been discussing is really what does it do "automatically".

I don't think any answer is going to be right for everyone, given the
diversity of users we have. I personally think we've found a nice
happy medium, but you're free to disagree.

>> > Hm. Can one change the meaning of "+" so that returns 0.54?
>> > What is the abstraction behind 0.04 ? RealField(n)(0.04) for some
>> > integer n?
>>
>> Yep, RealField(53)(0.04), which is a bit arbitrary.
>>
>> > And is there a difference between 0.10 and
>> > 0.1000000000000000000000000000000000
>>
>> Yep.
>>
>> sage: a = 0.1000000000000000000000000000000000
>> sage: a.precision()
>> 113
>
>
> So 0.1- 0.10000000000000000000000000000....
> is 0.0????????????????? where ? = undetermined?

sage: .1- 0.10000000000000000000000000000
0.000000000000000
sage: parent(_)
Real Field with 53 bits of precision

> and anyone who writes x+0.1 is relegating that sum to 1 digit "precision"?

As I mentioned previously, we default to the (somewhat arbitrary)
minimum of 53 bits of precision.

>> > Inventing fast methods is fun, although (my opinion) multiplying
>> > integers
>> > of astronomical size is hardly
>> > mainstream scientific computing.
>> >
>> > Not to say that someone might claim that
>> > this problem occurs frequently in many computations in pure and applied
>> > mathematics...
>>
>> I've personally "applied" multiplying astronomically sized before
>> (thought the result itself is squarely in the domain of pure math):
>> http://www.nsf.gov/news/news_summ.jsp?cntn_id=115646/
>
> Assume there is an application that involves multiplication of polynomials.
> You can multiply polynomials by encoding them as large integers,
> multiplying, and decoding. Sometimes called Kronecker's Trick.
>
> So there are lots of applications. Are they stupid tricks? Probably.

You claim "this idea has practical implications for efficient
programs" in http://www.cs.berkeley.edu/~fateman/papers/polysbyGMP.pdf
Now you claim it's stupid. Maybe it's both.

- Robert

[rjf]

On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw wrote:

The are two representations of the same canonical object.

The (computer algebra) use of the term, as in "simplified to a canonical form"  means

the representation is canonical.  It doesn't make much sense to claim that all these

are canonical:   1+1, 2,  2*x^0,  sin(x)^2+cos(x)^2 + exp(0).     

>>
> And what structure is that?  Does Sage know about   Z_{nonprime} ?

Of course, as illustrated.

sage: Integers(13^1024)
Ring of integers modulo 4764...1

How much does it know? Does it know that it is not a field, but that Integers(13) is a field?

> I'm still confused.   Is the term "Real Field" in Sage  the (or some)  real
> field?
>
> If it is an approximation to a field, but not a field, why are you calling
> it a field?

Because it's shorter to type and easier to find/discover than
ApproximateRealField or something like that.

> If that doesn't get you in trouble, why doesn't it?  Does Real Field inherit
> from
> Field?  Does every non-zero element have an inverse?

Of course it suffers from the same issues that (standard) floating
point numbers do in any language, user be aware (and we at least
document that).

And you know that everyone reads the documentation?

No, it doesn't suffer from the same issues as in other languages, because

those other languages probably don't refer to it as a field.

 

> Does Sage have other um, approximations, in its nomenclature?

Sure. RealField(123)[x]. Power series rings. P-adics.

These approximations are approximations by their nature.  If you are

computing with a power series, the concept inherently includes an error term

which you are aware of.  Real Field is (so far as I know) a concept that

should have the properties of a field.  The version in Sage does not.

It's like saying someone isn't pregnant.  well only a little pregnant. 

>

.... snip.... 

>> It is more conservative to convert operands to the domain with less
>> precision.
>
> Why do you say that?  You can always exactly convert a float number  in
> radix b to
> an equal number of higher precision in radix b by appending zeros.
> So it is more conserving (of values) to do so, rather than clipping off
> bits from the other.

Clipping bits (or digits) is exactly how one is taught to deal with
significant figures in grade school, and follows the principle of
least surprise (though floating point numbers like to play surprises
on you no matter what). It's also what floating point arithmetic does
when the exponent is different.

It is of course also taught in physics and chemistry labs, and I used this

myself in the days when slide-rules were used and you could read only

3 or so significant figures.  That doesn't make it suitable for a computer

system.  There are many things you learn along the way that are simplified

versions of the more fully elaborated systems of higher math.

What did you know about the branch cuts in the complex logarithm

or  log(-1)  when you were first introduced to log?

 

>> We consider the rationals to have infinite precision, our
>> real "fields" a specified have finite precision. This lower precision
>> is represented in the output, similar to how significant figures are
>> used in any other scientific endeavor.
>
> Thanks for distinguishing between "field" and field.  You don't seem
> to understand the concept of precision though.

That's a bold claim. My Ph.D. thesis depended on understanding issues
of precision. I'll admit explaining it to a layman can be difficult.

Is your thesis available online?  I would certainly look at it and see

how you define precision.

 

> You seem to think
> that a number of low precision has some inaccuracy or uncertainty.
> Which it doesn't.   0.5 is the same number as 0.500000.
> Unless you believe that 0.5  is  the interval [0.45000000000....01,
> 0.549999..........9]
> which you COULD believe  -- some people do believe this.
> But you probably don't want to ask them about scientific computing.

No, I don't think that at all. Sage also has the concept of real
intervals distinct from real numbers.

I suppose that the issue here is I don't know what you mean by "real number".

Sage has something in it called "real".  Mathematics uses that

term too (e.g. Dedekind cut).  They appear to be different. 

There's 0.5 the real number equal to one divided by two. There's also
0.5 the IEEE floating point number, which is a representative for an
infinite number of real numbers in a small interval.

Can you cite a source for that last statement?  While I suppose you

can decide anything you wish about your computer, the (canonical?)

explanation is that the IEEE ordinary floats correspond to a subset

of the EXACT rational numbers equal to <some integer> X 2^<some integer>

which is not an interval at all.

(there are also inf and nan things which are not ordinary in that sense)

So, citation?  (And I don't mean citing a physicist, or someone who

learned his arithmetic in 4th grade and hasn't re-evaluated it since.

A legitimate numerical analyst.) 

>> This is similar to 10 + (5 mod 13), where the right hand side has
>> "less precision" (in particular there's a canonical map one way, many
>> choices of lift the other.
>>
>> Also, when a user writes 1.3, they are more likely to mean 13/10 than
>> 5854679515581645 / 4503599627370496, but by expressing it in decimal
>> form they are asking for a floating point approximation. Note that
>> real literals are handled specially to allow immediate conversion to a
>> higher-precision parent.
>
> What do you know when a user writes 1.3, really?

Quite a bit. They probably didn't mean pi. If they really cared, they
could have been more specific. At least we recognize this ambiguity
but we can't let it paralyze us.

I think there is, in the community that consumes "scientific computing"

since 1960 or so, a set of expectations about "1.3".    You can use this,

or try to alter the expectations for your system.  For example, if you

displayed it as 13/10, that would be a change.  Or if you displayed it

as [1.25,1.35].   But then you are using a notation 1.25, and what

does that mean,  [[1.245,1.255], ....  ]  etc. 

> You want the user
> to believe that Sage uses decimal arithmetic?  Seriously?  How far
> are you going to try to carry that illusion?

If they don't immediately specify a new domain, we'll treat it as
having 53 bits. It's syntactic sugar.

So it sounds like you actually read the input as  13/10, because only then can

you  approximate it to higher precision than 53 bits or whatever.   Why not just admit this instead of talking

about 1.3.

 

> You imagine a user who understands rings and fields, and knows that
> Real Field is not a field, but knows so little about computers that he
> thinks 1.3  is 13/10 exactly?    (By the way, I have no problem  if
> 1.3 actually produces 13/10,  and to get a float, you have to try to
> convert it explicitly, and the conversion might even come up with an
> interval
> or an error bound or something that leads to "reliable" computing.
> Rather than some guessy-how-many-bits stab in the dark thing that
> prints as 1.3

Again, it's syntactic sugar.

For something to be syntactic sugar, you have to specify what it means underneath.

It seems that 1.3  without further context is syntactic sugar for 13/10 .

sage:1.3   is syntactic sugar for

   

sage: 13/10.RealField(53)

Or something like that.   

>> sage: QQ(1.3)
>> 13/10
>> sage: (1.3).exact_rational()
>> 5854679515581645/4503599627370496
>> sage: a = RealField(200)(1.3); a
>> 1.3000000000000000000000000000000000000000000000000000000000
>> sage: a.exact_rational()
>> 522254864384171839551137680010877845819715972979407671472947/401734511064747568885490523085290650630550748445698208825344
>
> I assume it is possible to calculate all kinds of things if you carefully
> specify them,
> in Sage.  After all, it has all those programs to call, including sympy.
> The
> issue iwe have been discussing is really what does it do "automatically".

I don't think any answer is going to be right for everyone, given the
diversity of users we have. I personally think we've found a nice
happy medium, but you're free to disagree.

Yes.   

... 

>> sage: a = 0.1000000000000000000000000000000000
>> sage: a.precision()
>> 113
>
>
> So 0.1- 0.10000000000000000000000000000....
> is 0.0?????????????????   where ? = undetermined?

sage: .1- 0.10000000000000000000000000000
0.000000000000000
sage: parent(_)
Real Field with 53 bits of precision

> and anyone who writes x+0.1  is relegating that sum to 1 digit "precision"?

As I mentioned previously, we default to the (somewhat arbitrary)
minimum of 53 bits of precision.

OK,  so you are saying that .1   has 53 bits of precision even though it appears to have, oh, about 3 bits.

Are you familiar with the problems such a design decision causes in Mathematica? 

>> > Inventing fast methods is fun, although (my opinion)  multiplying
>> > integers
>> > of astronomical size is hardly
>> > mainstream scientific computing.
>> >
>> > Not to say that someone might claim that
>> > this problem occurs frequently in many computations in pure and applied
>> > mathematics...
>>
>> I've personally "applied" multiplying astronomically sized before
>> (thought the result itself is squarely in the domain of pure math):
>> http://www.nsf.gov/news/news_summ.jsp?cntn_id=115646/
>
>  Assume there is an application that involves multiplication of polynomials.
> You can multiply polynomials by encoding them as large integers,
> multiplying, and decoding.  Sometimes called Kronecker's Trick.
>
> So there are lots of applications.  Are they stupid tricks? Probably.

You claim "this idea has practical implications for efficient
programs" in http://www.cs.berkeley.edu/~fateman/papers/polysbyGMP.pdf
Now you claim it's stupid. Maybe it's both.

Yes it is both.  A stupid trick because (assuming you are grabbing someone

else's library code for huge integer multiplication) you might as well grab

someone else's library code for fast polynomial multiplication.  And that

would generally be faster.

The practical implication is,  if your system constrains you to using a naive polynomial

multiplication program,  but it also has available a super-fast integer multiplication program,

as well as  [generally very important] fast encoding and decoding, you have a choice of

how to do polynomial multiplication, notably for high degree polynomials with small, esp. finite

field, coefficients.  (This is not usually a choice a CAS designer has to make -- it is

easier to program a non-naive polynomial multiplication program.)

RJF

[Robert Bradshaw]

On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
>
>
> On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw wrote:
>>
>>
>>
>> The are two representations of the same canonical object.
>
>
> The (computer algebra) use of the term, as in "simplified to a canonical
> form" means
> the representation is canonical. It doesn't make much sense to claim that
> all these
> are canonical: 1+1, 2, 2*x^0, sin(x)^2+cos(x)^2 + exp(0).

The point was that there's a canonical domain in which to do the computation.

>> > And what structure is that? Does Sage know about Z_{nonprime} ?
>>
>> Of course, as illustrated.
>>
>> sage: Integers(13^1024)
>> Ring of integers modulo 4764...1
>
>
> How much does it know? Does it know that it is not a field, but that
> Integers(13) is a field?

sage: Integers(13).is_field()
True
sage: Integers(13^1024).is_field()
False

>> > I'm still confused. Is the term "Real Field" in Sage the (or some)
>> > real
>> > field?
>> >
>> > If it is an approximation to a field, but not a field, why are you
>> > calling
>> > it a field?
>>
>> Because it's shorter to type and easier to find/discover than
>> ApproximateRealField or something like that.
>>
>> > If that doesn't get you in trouble, why doesn't it? Does Real Field
>> > inherit
>> > from
>> > Field? Does every non-zero element have an inverse?
>>
>> Of course it suffers from the same issues that (standard) floating
>> point numbers do in any language, user be aware (and we at least
>> document that).
>
> And you know that everyone reads the documentation?
> No, it doesn't suffer from the same issues as in other languages, because
> those other languages probably don't refer to it as a field.

The issues of floating point errors and rounding are much larger than
the question of whether every element has an inverse. You seem very
fixated on the name.

We also have an object called the ring of integers, but really it's
the ring of integers that fits into the memory of your computer.
Should we not call it a Ring?

>> > Does Sage have other um, approximations, in its nomenclature?
>>
>> Sure. RealField(123)[x]. Power series rings. P-adics.
>
> These approximations are approximations by their nature. If you are
> computing with a power series, the concept inherently includes an error term
> which you are aware of. Real Field is (so far as I know) a concept that
> should have the properties of a field. The version in Sage does not.
> It's like saying someone isn't pregnant. well only a little pregnant.

They're no more approximate by nature than the real numbers.

The p-adic numbers form a field. For any choice of representation some
of them can be represented exactly on a computer, most can't. When
doing computations with p-adic numbers one is typically chooses a
precision (e.g. how many digits, not unlike a choice of number of
bits) to use.

Power series (to make things concrete, say the power series in one
variable over the integers) form a ring. For any choice of
representation some of them can be represented exactly on a computer,
most can't. When doing computations with power series one is typically
chooses a precision (e.g. how many terms, not unlike a choice of
number of bits) to use.

Real numbers form a field. For any choice of representation some of
them can be represented exactly on a computer, most can't. When doing
computations with real numbers...

>> >> It is more conservative to convert operands to the domain with less
>> >> precision.
>> >
>> > Why do you say that? You can always exactly convert a float number in
>> > radix b to
>> > an equal number of higher precision in radix b by appending zeros.
>> > So it is more conserving (of values) to do so, rather than clipping off
>> > bits from the other.
>>
>> Clipping bits (or digits) is exactly how one is taught to deal with
>> significant figures in grade school, and follows the principle of
>> least surprise (though floating point numbers like to play surprises
>> on you no matter what). It's also what floating point arithmetic does
>> when the exponent is different.
>
>
> It is of course also taught in physics and chemistry labs, and I used this
> myself in the days when slide-rules were used and you could read only
> 3 or so significant figures. That doesn't make it suitable for a computer
> system. There are many things you learn along the way that are simplified
> versions of the more fully elaborated systems of higher math.
> What did you know about the branch cuts in the complex logarithm
> or log(-1) when you were first introduced to log?

Only being able to store 53 significant bits is completely analogous
to only being able to read 3 significant (decimal) figures. I think
the analogy is very suitable for a computer system. It can clearly be
made much more rigorous and precise.

Or are you seriously proposing when adding 3.14159 and 1e-100 it makes
more sense, by default, to pad the left hand side with zeros (whether
in binary or decimal) and return 3.1415900000...0001 as the result?

>> >> We consider the rationals to have infinite precision, our
>> >> real "fields" a specified have finite precision. This lower precision
>> >> is represented in the output, similar to how significant figures are
>> >> used in any other scientific endeavor.
>> >
>> > Thanks for distinguishing between "field" and field. You don't seem
>> > to understand the concept of precision though.
>>
>> That's a bold claim. My Ph.D. thesis depended on understanding issues
>> of precision. I'll admit explaining it to a layman can be difficult.
>
> Is your thesis available online? I would certainly look at it and see
> how you define precision.

Yes, it's online.

I define precision (e.g. computing a value to a given precision) as
the negated log of the absolute (respectively relative) error, but
most importantly it's something you can have more or less of, and lose
due to rounding, etc. Perhaps I could have used the term "accuracy"
instead. I use the actual errors themselves in my analysis.

Generally, however, precision is something attached to a read "field"
and specifies the number of bits used to represent its elements (and,
consequently, exactly what rounding errors are introduced when doing
arithmetic).

How about the docs on error analysis for LAPACK, which are presumably
written by an expert:
https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-51614DEE-9DF8-4D45-80F5-3B25CB7FF748.htm

There it says "you often need to solve a system Ax = b, where the data
(the elements of A and b) are not known exactly." So what are the IEEE
floats sitting in your computer when you're doing this computation?
They're representatives for the actual (either unknown or impossible
to represent) elements of your data. And, yes, they're also exact
rational numbers--once again they can be both--but the point is that
the input and output floating point numbers are viewed as
perturbations of the (actual) real field elements you care about.

Here they also talk about "loss of precision" similar in spirit to the
precision (of a value) as above.

In this case the user gives us a decimal literal. Yes, this literal is
equal to 13/10. We defer interpreting this as a 53-bit binary floating
point number long enough for the user to tell us to interpret it
differently. This prevents surprises like

sage: RealField(100)(float(1.3))
1.3000000000000000444089209850

or, more subtly

sage: sqrt(RealField(100)(float(1.3)))
1.1401754250991379986106491649

instead of

sage: sqrt(RealField(100)(1.3))
1.1401754250991379791360490256

When you write 1.3, do you really think 5854679515581645 /
4503599627370496, or is your head really thinking "the closest thing
to 13/10 that I can get given my choice of floating point
representation?" I bet it's the latter, which is why we do what we do.

- Robert

[rjf]

On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:

On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
>
>
> On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw wrote:
>>
>>
>>
>> The are two representations of the same canonical object.
>
>
> The (computer algebra) use of the term, as in "simplified to a canonical
> form"  means
> the representation is canonical.  It doesn't make much sense to claim that
> all these
> are canonical:   1+1, 2,  2*x^0,  sin(x)^2+cos(x)^2 + exp(0).

The point was that there's a canonical domain in which to do the computation.

I have not previously encountered the term "canonical domain".  There is
a CAS literature which includes the concept of simplification to a canonical form.
There is also a useful concept of a zero-equivalence test, whereby E1-E2
can be shown to be zero, although there is not necessarily a simplification
routine that will "canonically simplify"  E1  to E3 and also E2 to E3.
 

Just that for a mathematician to call something a field when it isn't
would, I think, be a serious mistake.  You seem to think that it is
ok for a computer system to do so.    I certainly agree that there are
more important issues with floating point computation than the
fact that these numbers do not constitute a real closed field.
 

We also have an object called the ring of integers, but really it's
the ring of integers that fits into the memory of your computer.
Should we not call it a Ring?

The domain of arbitrary-precision integers is an excellent model of the
ring of integers.  It is true that one can specify a computation that would
fill up the memory of all the computers in existence. or even all the atoms
in the (known?) universe.  Presumably a well-constructed support system
will give an error message on much smaller examples.   I assume
that your Real Field  operation of   division would give an error if the
result is inexact.   

>> > Does Sage have other um, approximations, in its nomenclature?
>>
>> Sure. RealField(123)[x]. Power series rings. P-adics.
>
> These approximations are approximations by their nature.  If you are
> computing with a power series, the concept inherently includes an error term
> which you are aware of.  Real Field is (so far as I know) a concept that
> should have the properties of a field.  The version in Sage does not.
> It's like saying someone isn't pregnant.  well only a little pregnant.

They're no more approximate by nature than the real numbers.

Huh?  How so? I was not aware that real numbers (at least the ones
that you can construct) are approximate by nature.  But maybe you
can direct me to some reference on this.

The p-adic numbers form a field. For any choice of representation some
of them can be represented exactly on a computer, most can't. When
doing computations with p-adic numbers one is typically chooses a
precision (e.g. how many digits, not unlike a choice of number of
bits) to use.

Power series (to make things concrete, say the power series in one
variable over the integers) form a ring. For any choice of
representation some of them can be represented exactly on a computer,
most can't. When doing computations with power series one is typically
chooses a precision (e.g. how many terms, not unlike a choice of
number of bits) to use.

Truncated power series with coefficients that are typically exact rational
numbers  (but could be, you say , elements of a finite field)  form a
computation structure that has exact operations.   Just because you
associate some other concept that is outside the computer with that
TPS  doesn't mean the TPS is approximate.

 

Real numbers form a field. For any choice of representation some of
them can be represented exactly on a computer, most can't. When doing
computations with real numbers...

So you would say that 3  is approximate, because maybe it is pi, but pi
cannot be represented as an integer.  This point of view seems to me
to be peculiar.

>> >> It is more conservative to convert operands to the domain with less
>> >> precision.
>> >
>> > Why do you say that?  You can always exactly convert a float number  in
>> > radix b to
>> > an equal number of higher precision in radix b by appending zeros.
>> > So it is more conserving (of values) to do so, rather than clipping off
>> > bits from the other.
>>
>> Clipping bits (or digits) is exactly how one is taught to deal with
>> significant figures in grade school, and follows the principle of
>> least surprise (though floating point numbers like to play surprises
>> on you no matter what). It's also what floating point arithmetic does
>> when the exponent is different.
>
>
> It is of course also taught in physics and chemistry labs, and I used this
> myself in the days when slide-rules were used and you could read only
> 3 or so significant figures.  That doesn't make it suitable for a computer
> system.  There are many things you learn along the way that are simplified
> versions of the more fully elaborated systems of higher math.
> What did you know about the branch cuts in the complex logarithm
> or  log(-1)  when you were first introduced to log?

Only being able to store 53 significant bits is completely analogous
to only being able to read 3 significant (decimal) figures.

Actually this analogy is false.  The 3 digits (sometimes 4) from a slide
rule are the best that can be read out because of the inherent uncertainty
in the rulings and construction of the slide rule, the human eye reading
the lines, etc.   So if I read my slide rule and say 0.25  it is because I  think
it is closer to 0.25  than 0.24 or 0.26   There is uncertainty there.
If a floating point number is computed as 0.25, there is no uncertainty in
the representation per se.  It is 1/4, exactly a binary fraction, etc.
Now you could use this representation in various ways, e.g.
0.25+-0.01    storing 2 numbers representing a center and a "radius"
or an interval or ..../   But the floating point number itself is simply
a computer representation of a particular rational number    aaa x 2^bbb
Nothing more, nothing less.  And in particular it does NOT mean
that bits 54,55,56... are uncertain.  Those bits do not exist in the representation
and are  irrelevant for ascertaining the value of the number aaa x 2^bbb.

So the analogy is false.



On the other hand, the
 

I think
the analogy is very suitable for a computer system. It can clearly be
made much more rigorous and precise.

What you are referring to is sometimes called significance arithmetic,
and it has been thoroughly discredited. 
Sadly, Wolfram the physicist put it in Mathematica.


 

Or are you seriously proposing when adding 3.14159 and 1e-100 it makes
more sense, by default, to pad the left hand side with zeros (whether
in binary or decimal) and return 3.1415900000...0001 as the result?

If you did so, you would preserve the  identity  (a+b)-a   =  b
 
If you round to some number of bits, say 53,  with a=3.14159  and b=1e-100,
the left side is 0, and the right side  is 1e-100.  The relative error in the answer
is, um, infinite.

Now if the user specified the kind of arithmetic explicitly, or even implicitly
by saying "use IEEE754 binary floating point arithmetic everywhere" then
I could go along with that.

>> >> We consider the rationals to have infinite precision, our
>> >> real "fields" a specified have finite precision. This lower precision
>> >> is represented in the output, similar to how significant figures are
>> >> used in any other scientific endeavor.
>> >
>> > Thanks for distinguishing between "field" and field.  You don't seem
>> > to understand the concept of precision though.
>>
>> That's a bold claim. My Ph.D. thesis depended on understanding issues
>> of precision. I'll admit explaining it to a layman can be difficult.
>
> Is your thesis available online?  I would certainly look at it and see
> how you define precision.

Yes, it's online.

I define precision (e.g. computing a value to a given precision) as
the negated log of the absolute (respectively relative) error, but
most importantly it's something you can have more or less of, and lose
due to rounding, etc. Perhaps I could have used the term "accuracy"
instead. I use the actual errors themselves in my analysis.

I looked at your thesis (we exchanged mail on this.)
Yes you could have used  absolute error or relative error.
Taking a log is optional.
I think that if you said you increase the error by cancellation or rounding,
that would be OK.

 

Generally, however, precision is something attached to a read "field"
and specifies the number of bits used to represent its elements (and,
consequently, exactly what rounding errors are introduced when doing
arithmetic).

The number of bits used to represent a number n  is one of the parameters
that determines the result of rounding, yes.  Other parameters are the presence
or absence of guard digits, and rules for resolving ties.

It does not say whether a number is accurate or not.

 

What this means is that the numbers in the system A and b often come
 from the real world and you don't know their true values, just the values
that you measured.
When you enter the data into the computer you put the values for A and b
that are the closest numbers to what you believe.  A major question in
computational linear algebra is how to measure the sensitivity of the computed
solution to slight changes in the initial data.  The numbers used in solving
the system are exactly the represented floating-point numbers.

They're representatives for the actual (either unknown or impossible
to represent) elements of your data. And, yes, they're also exact
rational numbers--once again they can be both--but the point is that
the input and output floating point numbers are viewed as
perturbations of the (actual) real field elements you care about.

That does not give you or anyone else permission to say  "oh, this
data might not be right so instead of solving the system as best as
I can, I will use some off-hand rule to clip bits off, since it's just a
perturbation of something".   No, you compute the best you can
with the data you have.
 

Here they also talk about "loss of precision" similar in spirit to the
precision (of a value) as above.

eh, precision occurs once in the link you cite...

"..... then the error in the solution x is also large; you may even encounter a complete loss of precision.  ...."

The relevant term is actually "complete loss of precision"   which does have a useful meaning, which
is that the relative error is infinite. 

I suspect it is not what python does. 
It is what Macsyma does if you write 1.3b0   to indicate "bigfloat".
 

- Robert

[Erik Massop]

On Wed, 6 Aug 2014 19:00:22 -0700 (PDT)
rjf <fat...@gmail.com> wrote:

> On Tuesday, August 5, 2014 7:59:00 PM UTC-7, Robert Bradshaw wrote:
> >
> > On Tue, Aug 5, 2014 at 6:36 PM, rjf <fat...@gmail.com>
> > wrote:

...

> > > Which one, -3 or 10?
> > > They can't both be canonical.
> >
> > Yes, they can "both" be canonical: they are equal (as elements of
> > Z/13Z). There is the choice of internal representation, and what
> > string to use when displaying them to the user, but that doesn't
> > affect its canonicity as a mathematical object.
> >
> I don't understand this.
> I use the word canonical to mean unique distinguished exemplar.
> So there can't be two. If there are two distinguishable items, then one or
> the other or neither might be canonical. Not both.

It seems to me that both of you mean "unique distinguished exemplar"
when speaking of canonical. The problem seems to be that for "unique"
to make sense you should agree on an ambient set. For RJF the ambient
set seems to be ZZ, the set of integers, while for Robert the ambient
set seems to be ZZ/13ZZ, the set of integers modulo 13.

Let's do a more elaborate example: Let ZZ be the ring of integers and
let 3ZZ be its ideal {3k : k in ZZ}. Take ZZ/3ZZ to be the quotient
ring (that is, the ring of integers modulo 3, which happens to be a
field), and take b to be the element represented by the integer 2.

Then we have:
I) The integer 2 is not a canonical representation of b.
II) The element 2 of {0,1,2} is a canonical representation of b.
III) The element b of ZZ/3ZZ is a canonical representation of b.
IV) The element 2 of ZZ/3ZZ is a canonical representation of b.
V) The element 5 of ZZ/3ZZ is a canonical representation of b.
VI) The elements 2 and 5 of ZZ/3ZZ are both canonical representations
of b.
These are the proofs:
I) The integer 5, which is not 2, also represents b.
II) The elements 0 and 1 of {0,1,2} represent elements of ZZ/3ZZ that
are not b.
III) This is a trivial: If c is an element of ZZ/3ZZ that is b, then it
is b.
IV) The element 2 of ZZ/3ZZ - that is to say, the element represented
by the integer 2 - is b. The result follows by assertion III.
V) Analogous to IV.
VI) Immediate from IV and V.

The importance of agreeing on the ambient set becomes apparent when
comparing assertions I and IV.

The words "of ZZ/3ZZ" are important indeed in assertion VI: If we had
read "of ZZ", then this assertion would be pointing directly at its own
counterexample.


Regards,

Erik Massop

[Erik Massop]

On Fri, 8 Aug 2014 18:57:21 -0700 (PDT), rjf <fat...@gmail.com> wrote:

> On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:
>
> > On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
> >
> > > On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw
> > > wrote:

...

> Just that for a mathematician to call something a field when it isn't
> would, I think, be a serious mistake. You seem to think that it is
> ok for a computer system to do so. I certainly agree that there are
> more important issues with floating point computation than the
> fact that these numbers do not constitute a real closed field.

In my experience computer scientists are often more precise than
mathematicians.

It seems quite common to me to (intentionally) confuse mathematical
objects with their imperfect realizations. For instance when I use my
compasses to draw something that looks a lot like a circle, I would call
what I drew a circle. Yet I am sure that the points that I drew do not
have equal distance to the center.

Similarly I would call RealField(prec) a field, because it represents
the field of real numbers (e.g. Dedekind cuts). At the same time I am
sure that some mathematical properties that I know from the mathematical
field of real numbers will not hold in RealField(prec). Luckily Sage
readily admits that RealField(prec) is not an exact representation of
the mathematical field of real numbers:
sage: F = RealField(100)
sage: F.is_exact()
False

I agree that it is imprecise that when using F.is_field(), we are
asking whether the mathematical object of real numbers is a field, while
when using F.is_exact() we are asking whether the representation is
exact.

By the way, there also is AA or AlgebraicRealField() in Sage, which is
an exact representation of the real closed field of algebraic real
numbers.

...

> The domain of arbitrary-precision integers is an excellent model of
> the ring of integers. It is true that one can specify a computation
> that would fill up the memory of all the computers in existence. or
> even all the atoms in the (known?) universe. Presumably a
> well-constructed support system will give an error message on much
> smaller examples. I assume that your Real Field operation of
> division would give an error if the result is inexact.

RealField(prec) is an inexact approximation of a field (as witnessed by
is_exact), so I never expect the division to be exact. In fact I don't
expect addition, negation, multiplication, or subtraction to be exact
either. Indeed, they are not exact:
sage: a = 1e-58
sage: a
1.00000000000000e-58
sage: (1+a)-1
0.000000000000000

> > > > > Does Sage have other um, approximations, in its nomenclature?
> > > >
> > > > Sure. RealField(123)[x]. Power series rings. P-adics.
> > >
> > > These approximations are approximations by their nature. If you
> > > are computing with a power series, the concept inherently includes
> > > an error term which you are aware of. Real Field is (so far as I
> > > know) a concept that should have the properties of a field. The
> > > version in Sage does not. It's like saying someone isn't pregnant.
> > > well only a little pregnant.
> >
> > They're no more approximate by nature than the real numbers.
>
> Huh? How so? I was not aware that real numbers (at least the ones
> that you can construct) are approximate by nature. But maybe you
> can direct me to some reference on this.

What do you mean by "can construct"? The computable reals perhaps? Those
would indeed form an actual field. I am very unsure about how practical
they are though.

Might we then not also have a ring of computable power series? That is,
those power series given by an algorithm that given n returns the n'th
term in finite time?

...

> > Power series (to make things concrete, say the power series in one
> > variable over the integers) form a ring. For any choice of
> > representation some of them can be represented exactly on a
> > computer, most can't. When doing computations with power series one
> > is typically chooses a precision (e.g. how many terms, not unlike a
> > choice of number of bits) to use.
>
> Truncated power series with coefficients that are typically exact
> rational numbers (but could be, you say , elements of a finite field)
> form a computation structure that has exact operations. Just because
> you associate some other concept that is outside the computer with
> that TPS doesn't mean the TPS is approximate.

I would argue that a ring of truncated power series with rational
coefficients is an inexact representation of the ring of (untruncated)
power series with rational coefficients, just as RealField(prec) is an
inexact representation of the field of real numbers.

It is important to realize that my use of "inexact" here refers to
TPS/RealField(prec) as a representation of the ring of (untruncated)
power series with rational coefficients/the field of real numbers, and
not to TPS/RealField(prec) itself.

The object PowerSeriesRing(QQ, 'x') in Sage is not a ring of truncated
power series QQ[x]/x^n. In it not a ring at all. In fact == is not even
transitive:
sage: R.<x> = PowerSeriesRing(QQ, 'x')
sage: a = O(x^20)
sage: a in R
True
sage: (0 == x^30 + a) and (x^30 + a == x^30)
True
sage: 0 == x^30
False

> > Real numbers form a field. For any choice of representation some of
> > them can be represented exactly on a computer, most can't. When
> > doing computations with real numbers...
>
> So you would say that 3 is approximate, because maybe it is pi, but pi
> cannot be represented as an integer. This point of view seems to me
> to be peculiar.

Indeed, I would argue that 3 in RealField(100) is approximate. For
instance the (computable) real number 3.00...001, where the ellipsis is
1000 zeroes, gets mapped to 3 in RealField(100), but is not equal to 3
as an actual real number.

When speaking of exact/inexact I think it is important, to keep track of
what is to be represented. Hence the element 3 of RealField(100) is an
inexact approximation of a real number since there are many reals that
become 3 in RealField(100). The element 3 of IntegerRing() is an exact
approximation of an integer, since there is only one integer that
becomes 3 in IntegerRing().

I think it is again important to remember what is to be represented. If
you use floats to express numbers of the form aaa x 2^bbb (with
appropriate conditions on aaa and bbb), then sure, floats are an exact
representation. However there are many real numbers that are not of
this form. So, if you use floats for real numbers, then they are an
inexact representation. In that case the bits 54, 55, 56, ... of the
real number are lost in the 53 bit floating point representation. Hence
given the 53 bit floating point representation, the bits 54, 55,
56, ... of the real number that it is supposed to represent, are indeed
uncertain. (This is simplified since the bit-numbering depends on the
exponent.)

If you happen to know that the real number that is represented is of the
form aaa x 2^bbb (with appropriate conditions on aaa and bbb), then you
can reconstruct all those extra bits, and the representation is arguably
exact. But how should Sage know that the represented real number is of
this form?

...

> > Or are you seriously proposing when adding 3.14159 and 1e-100 it
> > makes more sense, by default, to pad the left hand side with zeros
> > (whether in binary or decimal) and return 3.1415900000...0001 as the
> > result?
>
> If you did so, you would preserve the identity (a+b)-a = b

How would the real number pi be represented in this system? Do you have
to explicitly say that you want the first n digits? Or is the padding
scheme smart enough to add pi's digits? In the latter case, how do you
keep checking (in)equality efficient?

How does (1/3)*3 give? Does it compare equal to 1? I mean, I would
imagine 1/3 to be 0.333..., with more digits 3 as I need them, coming
from larger and larger 0 paddings of 1 and 3. Supposing this, (1/3)*3
would be 0.999..., with more digits 9 as I need them. Will comparing
with 1 ever finish?

...

> Now if the user specified the kind of arithmetic explicitly, or even
> implicitly by saying "use IEEE754 binary floating point arithmetic
> everywhere" then I could go along with that.

I think you can specify this to some extend. At least
sage: RealField?
gives mentions opportunities to tweak precision and rounding. You can
use AA if you want exact arithmetic (of real algebraic numbers). I don't
think there is a field of computable reals in Sage.

...

> > > > > You seem to think that a number of low precision has some
> > > > > inaccuracy or uncertainty. Which it doesn't. 0.5 is the same
> > > > > number as 0.500000. Unless you believe that 0.5 is the
> > > > > interval [0.45000000000....01, 0.549999..........9] which you
> > > > > COULD believe -- some people do believe this.
> > > > > But you probably don't want to ask them about scientific
> > > > > computing.

Again, I think you should keep track of ambient sets here. For sure, the
real numbers 0.5 and 0.5000000 are the same, and real numbers have no
uncertainty (when representing real numbers). However, users do not type
real numbers. They type finite strings of characters that they probably
mean to represent some real number. The _strings_ 0.5 and 0.5000000 are
not the same. Hence it might be reasonable to assign different meaning
to these strings, as long as the user also has a way to precisely
express what is meant. Yet, I'm no expert on scientific computing, so I
will leave this be.

...

> > > > There's 0.5 the real number equal to one divided by two. There's
> > > > also 0.5 the IEEE floating point number, which is a
> > > > representative for an infinite number of real numbers in a small
> > > > interval.
> > >
> > > Can you cite a source for that last statement? While I suppose
> > > you can decide anything you wish about your computer, the
> > > (canonical?) explanation is that the IEEE ordinary floats
> > > correspond to a subset of the EXACT rational numbers equal to
> > > <some integer> X 2^<some integer> which is not an interval at all.
> > > (there are also inf and nan things which are not ordinary in that
> > > sense)
> > >
> > > So, citation? (And I don't mean citing a physicist, or someone
> > > who learned his arithmetic in 4th grade and hasn't re-evaluated it
> > > since. A legitimate numerical analyst.)

I am not a numerical analyst, so perhaps the following is naive, but I
would really like know where it is flawed or why my assumptions are
unreasonable.

I call a real number exactly representable when it occurs as an exact
value of a float with some bounded number of bits for its exponent.
Then when using floats to represent real numbers, you have to decide
for each real number x by which exactly representable real number f(x)
it should be represented. I can think of two desirable properties:
* exactly representable real numbers should be represented by
themselves,
* the relation <= should be preserved.
These properties yield that for each exactly representable real number y
the set {x real : f(x) = y} is an interval around y. Moreover this
interval is contained in every interval (a, b) where a and b are exactly
representable real numbers with a < y and y < b.

Okay, we need some condition on the mantissa and exponent, for if the
mantissa and exponent are allowed to be arbitrary integers, a function
with these properties fails to exists: Let f be a partial function with
the required properties, then the intervals {x real : f(x) = y}, being
contained in the intervals (a,b) with a and b exactly representable
with a < y < b, necessarily contain only y. Since there are only
countable many exactly representable real numbers, each with preimage
containing only itself, the domain of definition of f is countable. In
particular f is not a total function.

...

> What this means is that the numbers in the system A and b often come
> from the real world and you don't know their true values, just the
> values that you measured.
> When you enter the data into the computer you put the values for A and
> b that are the closest numbers to what you believe. A major question
> in computational linear algebra is how to measure the sensitivity of
> the computed solution to slight changes in the initial data. The
> numbers used in solving the system are exactly the represented
> floating-point numbers.
>
> > They're representatives for the actual (either unknown or impossible
> > to represent) elements of your data. And, yes, they're also exact
> > rational numbers--once again they can be both--but the point is that
> > the input and output floating point numbers are viewed as
> > perturbations of the (actual) real field elements you care about.
>
> That does not give you or anyone else permission to say "oh, this
> data might not be right so instead of solving the system as best as
> I can, I will use some off-hand rule to clip bits off, since it's just
> a perturbation of something". No, you compute the best you can
> with the data you have.

You cannot just assume either that the data is exact.

Also it makes more sense to me not to do the best you can, but the best
you can do efficiently. Trying to do the best you can is of no use if it
means running out of memory or taking a year. If doing things
efficiently involves some clipping, then that's okay with me,
especially if the clipping is configurable.


Regards,

Erik Massop

[Ralf Stephan]

I would like to thank Eric Massop for his excellent and enlightening posts.

[Dr. David Kirkby (Kirkby Microwave Ltd)]

On 9 Aug 2014 02:57, "rjf" <fat...@gmail.com> wrote:
> On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:

>> The point was that there's a canonical domain in which to do the computation.
>
>
> I have not previously encountered the term "canonical domain".  There is
> a CAS literature which includes the concept of simplification to a canonical form.

Richard,
Is it really necessary or desirable to quote such large sections of a message?

I have not counted the number of words in your email, but it fills 18 screen fills on my mobile phone. I am sure that you could trim a large amount out, whilst preserving sufficient that someone can follow the context.

Dave.

[Bill Hart]

On Thursday, 31 July 2014 06:11:28 UTC+2, rjf wrote:

Note that even adding  5 mod 13   and the integer 10  is potentially uncomfortable,

and the rather common operation of Hensel lifting requires doing arithmetic in a combination

of fields (or rings) of different related sizes.

 

Oh what a shame I missed the party. And I don't have time to read all the wonderfully eloquent nonsense in this thread. But this above really wasn't a good example.

The new flint interpreter, written in Julia, did this in week 1 I think. It is not even a coercion from my perspective (I call it ad hoc polymorphism):

julia> R = ResidueRing(ZZ, 13)

Residue ring of ZZ modulo 13 (constructor with 4 methods)

julia> a = R(5) + 10

2

julia> typeof(a)

Residue ring of ZZ modulo 13 (constructor with 4 methods)

Now, Nemo (that's the name of the flint interpreter) doesn't handle "A+B where A is a polynomial in ZZ[x,y] and B is a power series in F_q[[x]] (finite field with q elements)", for two reasons:

1) I haven't implemented power series yet (the other objects are there, and give me another week on the power series).

2) I wouldn't implement such a coercion. The problem is that, at best, it will create a ring that the user did not define. So I am firmly in the camp that believes the user should create the ring they want these added in, then convert at least one of the objects to that ring. Then the coercion system can take over from there. Whether Nemo will even go that far, I'm not sure. It's just a little flint interpreter and it certainly doesn't need to be overengineered! (Note to self: the above example is an example of overengineering.)

In fact, I wouldn't want that example to even make sense. The variable x is a generator of a certain ring, with type information attached to it, and shouldn't be used in two different ways. You are effectively redefining x. When I type parent(x), what will it tell me? Or does x also get lifted to the common ring before I even try to add objects from the two different rings? Suddenly the example makes me feel like I don't understand any more what Sage is even going to do when I type stuff into it.

Magma does silly things like this. You define a univariate polynomial ring over x and another over y. It happily adds the two together, forgetting that x is not y, so long as it can coerce the coefficient rings to a common ring (or something like that -- it depends how people felt on a certain day), which is really not useful when that gives you the wrong answer entirely, which it does (I have nontrivial examples where this really matters).

Two important things to note:

1) I did not give a reason "this would be hard to do in Julia". Because it wouldn't. It's so trivial as to be embarrassing. And yes it is zero cost. No looking up things in hash tables or whatever. Julia can even programmatically generate coercions on-the-fly. And I do this already in Nemo!! Powerful metaprogramming + dependent and parametric types + multimethods + type inference + jit + conversion functions + promotion rules + bare metal types + efficient C FFI = something Python can never do and will never do.

The proof is in the pudding. Just generic multivariate polynomial multiplication is already 3 times faster than Sage. And I didn't even try yet. Believe it or not, you can in theory actually resoundingly beat C at generic algorithms with Julia. And that is not an empty claim. We are already at exactly the same speed as a fast hand written C generic rings implementation. And I didn't even try yet! I've literally applied one trivial optimisation, which Sage certainly uses. In theory, a language like Julia can do *much* better.

2) I'm not defending Sage here, even though Richard's example was completely trivial. Obviously I check the answers from the flint interpreter against what comes out of Sage, Magma and Pari. And Sage leaves me scratching my head pretty often. Magma errs on the side of not defining questionable things at all, or if it does, it gives answers you really don't expect, because it loses type information when doing coercion. Sage errs on the side of defining the balmy and bizarre. Pari errs on the side of not being a serious programming language.

Writing a little interpreter has certainly opened my eyes to the vagaries of computer algebra and mathematics in general. I'm a little more sympathetic as a whole to the difficulties of writing computer algebra systems. And I now see mathematics itself as a patchwork of competing definitions and nomenclature and a set of noncompatible axiom systems, rather than a unified whole.

But I think Julia is good for a few years yet. It's one of the few truly 21st century programming languages available. I don't know if Julia is going to be the language to displace Python or not. It has some seriously difficult problems to solve ahead of it (fortunately they don't include parallel programming, a package manager which works, an IPython interface, efficient C, C++ or Python interfaces, an extremely modern syntax or language features, a well-fleshed-out standard library, scientific computing capabilities, high performance, profiling and timing, disassembly features, etc.). I'm sure you can sense from Stefan's posts here that the Julia developers are no idiots. A language like Julia has the potential to be a game changer for computer algebra. Python doesn't. I've been saying that for years, and the reasons haven't changed one bit. You can check my git repository and see that for years I've been working on toy programming languages based on LLVM, with type inference, parameterised type systems, a REPL, Jit, all the things Julia is. I never got around to implementing serious metaprogramming, but there are blog and even sage-devel posts going back years where I've advocated that.

Computer scientists have a lot more to offer the mathematics community than it is currently letting them do. I'm glad for the existence of Sage. It's a great community of great people, and a useful tool. It's an outstanding accomplishment. But technology is moving at a far faster pace than the past. Sage is not keeping up. Not even nearly. It will be superceded, and I want to be on record as having said so (again).

[William A Stein]

 If so, we'll just include whatever supersedes sage in Sage...   

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
To post to this group, send email to sage-...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.

--
William Stein
Professor of Mathematics
University of Washington
http://wstein.org

wst...@uw.edu

[Bill Hart]

Too late. Whatever succeeded Sage already included Sage in it.

Anyway, that is a brilliant idea. But you should start now, not when it is too late. And I mean that with all sincerity.

On Wednesday, 20 August 2014 06:19:21 UTC+2, wstein wrote:

To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscribe@googlegroups.com.

To post to this group, send email to sage-...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.

[rjf]

Sorry I've been away for a week, and sorry for quoting so much of 

Erik's note.  It is unusual for me to disagree with virtually EVERY statement in such a large collection.  Obviously we have different contexts for almost everything.  Now I would not say "violent disagreement" because almost each point is trivial, yet wrong from my perspective. 

On Saturday, August 9, 2014 7:35:13 PM UTC-7, Erik Massop wrote:

On Fri, 8 Aug 2014 18:57:21 -0700 (PDT), rjf <fat...@gmail.com> wrote:

> On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:
>
> > On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
> >
> > > On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw
> > > wrote:
...
> Just that for a mathematician to call something a field when it isn't
> would, I think, be a serious mistake.  You seem to think that it is
> ok for a computer system to do so.    I certainly agree that there are
> more important issues with floating point computation than the
> fact that these numbers do not constitute a real closed field.

In my experience computer scientists are often more precise than

mathematicians. )

Logicians,  a subclass of mathematicians, tend to be precise. Computer scientists have to be precise when writing programs (or correct programs, anyway .)

It seems quite common to me to (intentionally) confuse mathematical
objects with their imperfect realizations. For instance when I use my
compasses to draw something that looks a lot like a circle, I would call
what I drew a circle. Yet I am sure that the points that I drew do not
have equal distance to the center.

It is obvious that a picture of something resembling a circle is not a circle in the abstract, and there is no confusion.   There IS a confusion if you have a bit arrangement in a computer word that you say is a canonical 2, and another bit arrangement in a computer word that differs and you say that is a canonical 2.

More later.

 

[rjf]

On Saturday, August 9, 2014 7:35:13 PM UTC-7, Erik Massop wrote:

....

 

Similarly I would call RealField(prec) a field, because it represents
the field of real numbers (e.g. Dedekind cuts). At the same time I am
sure that some mathematical properties that I know from the mathematical
field of real numbers will not hold in RealField(prec).

So it seems to me that you should not call it a Field.

Here's the inverse.  Or maybe contrapositive.  As a riddle:

How many legs does a sheep have if you call its tail a leg?

(Answer later...)

 

Luckily Sage
readily admits that RealField(prec) is not an exact representation of
the mathematical field of real numbers:
  sage: F = RealField(100)
  sage: F.is_exact()
  False

Hardly luck, but random insertion of a fact in a way similar to the riddle.

(Answer later)

 

I agree that it is imprecise that when using F.is_field(), we are
asking whether the mathematical object of real numbers is a field, while
when using F.is_exact() we are asking whether the representation is
exact.

By the way, there also is AA or AlgebraicRealField() in Sage, which is
an exact representation of the real closed field of algebraic real
numbers.

...
> The domain of arbitrary-precision integers is an excellent model of
> the ring of integers.  It is true that one can specify a computation
> that would fill up the memory of all the computers in existence. or
> even all the atoms in the (known?) universe.  Presumably a
> well-constructed support system will give an error message on much
> smaller examples.   I assume that your Real Field  operation of
> division would give an error if the result is inexact.

RealField(prec) is an inexact approximation of a field (as witnessed by
is_exact), so I never expect the division to be exact.

Oh, so it just gives wrong answers without a clue.

 

In fact I don't
expect addition, negation, multiplication, or subtraction to be exact
either. Indeed, they are not exact:
  sage: a = 1e-58
  sage: a
  1.00000000000000e-58
  sage: (1+a)-1
  0.000000000000000

That's because it is not a field, and you know it.  So why call it a field?

 

> > > > > Does Sage have other um, approximations, in its nomenclature?
> > > >
> > > > Sure. RealField(123)[x]. Power series rings. P-adics.
> > >
> > > These approximations are approximations by their nature.  If you
> > > are computing with a power series, the concept inherently includes
> > > an error term which you are aware of.  Real Field is (so far as I
> > > know) a concept that should have the properties of a field.  The
> > > version in Sage does not. It's like saying someone isn't pregnant.
> > > well only a little pregnant.
> >
> > They're no more approximate by nature than the real numbers.
>
> Huh?  How so? I was not aware that real numbers (at least the ones
> that you can construct) are approximate by nature.  But maybe you
> can direct me to some reference on this.

What do you mean by "can construct"? The computable reals perhaps? Those
would indeed form an actual field. I am very unsure about how practical
they are though.

How about the rational numbers?  Those you can easily construct.

The computable reals are not fun to compute with in most "application" circumstances.

 

Might we then not also have a ring of computable power series? That is,
those power series given by an algorithm that given n returns the n'th
term in finite time?

It is not usually phrased that way, but given a collection of base functions as explicit power series

and operations like composition, there is a computational domain of truncated power series.

It is usually set up to NOT be a ring because of truncation issues.

...
> > Power series (to make things concrete, say the power series in one
> > variable over the integers) form a ring. 

For any choice of
> > representation some of them can be represented exactly on a
> > computer, most can't. When doing computations with power series one
> > is typically chooses a precision (e.g. how many terms, not unlike a
> > choice of number of bits) to use.
>
> Truncated power series with coefficients that are typically exact
> rational numbers  (but could be, you say , elements of a finite field)
> form a computation structure that has exact operations.   Just because
> you associate some other concept that is outside the computer with
> that TPS  doesn't mean the TPS is approximate.

I would argue that a ring of truncated power series with rational
coefficients is an inexact representation of the ring of (untruncated)
power series with rational coefficients, just as RealField(prec) is an
inexact representation of the field of real numbers.

Well then, I would argue that it is not.  It is a computational domain

with its own useful properties that can be nicely characterized, and

happens not to be a ring.

 

It is important to realize that my use of "inexact" here refers to
TPS/RealField(prec) as a representation of the ring of (untruncated)
power series with rational coefficients/the field of real numbers, and
not to TPS/RealField(prec) itself.

Basically, this is analogous to saying that floating point numbers

form a field, which they don't.  Let the coefficients be 0 or 1, and the

base of the power series x is 1/2.

or maybe p-adics.

 

The object PowerSeriesRing(QQ, 'x') in Sage is not a ring of truncated
power series QQ[x]/x^n. In it not a ring at all. In fact == is not even
transitive:
  sage: R.<x> = PowerSeriesRing(QQ, 'x')
  sage: a = O(x^20)
  sage: a in R
  True
  sage: (0 == x^30 + a) and (x^30 + a == x^30)
  True
  sage: 0 == x^30
  False

As is the case for every computer algebra system I know about that implements

truncated power series.   Incidentally, Macsyma implements a different, additional, notion

of  power series wehre they are represented by a formula for the arbitrary nth term.

Since we agree that TPS do not form a ring,  why does Sage lie about it?

 

> > Real numbers form a field. For any choice of representation some of
> > them can be represented exactly on a computer, most can't. When
> > doing computations with real numbers...
>
> So you would say that 3  is approximate, because maybe it is pi, but pi
> cannot be represented as an integer.  This point of view seems to me
> to be peculiar.

Indeed, I would argue that 3 in RealField(100) is approximate.

I can't see why you would think so.

 

For
instance the (computable) real number 3.00...001, where the ellipsis is
1000 zeroes, gets mapped to 3 in RealField(100), but is not equal to 3
as an actual real number.

So the number k = 3.000....1 is not in RealField(100).  Who cares?

That does not justify treating 3 as a huge fuzzball.

After all if you believe that all the numbers in RealField are fuzzballs,

why not represent k as the number 5 in RealField(100)? 

Do you have a theory of error in RealField(100)?  I suspect not.

When speaking of exact/inexact I think it is important, to keep track of
what is to be represented.

Nope.

Absolute nonsense.

Arithmetic is done on what you are given, not on what someone in the back room is thinking.

If you want to compute with numbers that have error bounds or error distributions, there

are ways of representing them and doing arithmetic on them.

 

Hence the element 3 of RealField(100) is an
inexact approximation of a real number since there are many reals that
become 3 in RealField(100).

Nonsense.

 

The element 3 of IntegerRing() is an exact
approximation of an integer, since there is only one integer that
becomes 3 in IntegerRing().

Why do you know that?  What if the integer is read off a (weight) scale that

has gradations at (say) 150, 155, 160, ...   and so the integers 153,154,155,156,157   all become 155.

 

.... 

> > to only being able to read 3 significant (decimal) figures.
>

(RJF said this....) 

> Actually this analogy is false.  The 3 digits (sometimes 4) from a
> slide rule are the best that can be read out because of the inherent
> uncertainty in the rulings and construction of the slide rule, the
> human eye reading the lines, etc.   So if I read my slide rule and say
> 0.25  it is because I think it is closer to 0.25  than 0.24 or 0.26
> There is uncertainty there. If a floating point number is computed as
> 0.25, there is no uncertainty in the representation per se.  It is
> 1/4, exactly a binary fraction, etc. Now you could use this
> representation in various ways, e.g. 0.25+-0.01    storing 2 numbers
> representing a center and a "radius" or an interval or ..../   But the
> floating point number itself is simply a computer representation of a
> particular rational number    aaa x 2^bbb  Nothing more, nothing less.
> And in particular it does NOT mean that bits 54,55,56... are
> uncertain.  Those bits do not exist in the representation and are
> irrelevant for ascertaining the value of the number aaa x 2^bbb.
>
> So the analogy is false.

I think it is again important to remember what is to be represented. If
you use floats to express numbers of the form aaa x 2^bbb (with
appropriate conditions on aaa and bbb), then sure, floats are an exact
representation. However there are many real numbers that are not of
this form.

So you cannot represent them, and it is mostly pointless to discuss what

you would do with them arithmetically.   If I said that I can't represent a

shoe in floating-point, and then proceeded to describe how to do arithmetic

on shoes, socks, etc,  you would probably find this dubious.

 

So, if you use floats for real numbers, then they are an
inexact representation. In that case the bits 54, 55, 56, ... of the
real number are lost in the 53 bit floating point representation. Hence
given the 53 bit floating point representation, the bits 54, 55,
56, ... of the real number that it is supposed to represent, are indeed
uncertain. (This is simplified since the bit-numbering depends on the
exponent.)

No they are not uncertain.  They are not represented at all.  They are

not part of the number.  If I told you that  3 was of uncertain color --

it might be red, blue, green or some combination,  would that enhance

your arithmetic?

 

If you happen to know that the real number that is represented is of the
form aaa x 2^bbb (with appropriate conditions on aaa and bbb), then you
can reconstruct all those extra bits, and the representation is arguably
exact.

No, the number in the computer is exactly the number that is the number.  It is not

some other number. 

But how should Sage know that the represented real number is of
this form?

Because that is the number in the computer.  If you want to represent a range,

you might use two numbers.   If the computer numbers are fuzzballs , you cannot

represent a range because the endpoints are fuzzballs,  so you need a range for

them.  Infinite recursion to represent one "real" number?   Eh, not such a great idea.

...
> > Or are you seriously proposing when adding 3.14159 and 1e-100 it
> > makes more sense, by default, to pad the left hand side with zeros
> > (whether in binary or decimal) and return 3.1415900000...0001 as the
> > result?
>
> If you did so, you would preserve the  identity  (a+b)-a   =  b

How would the real number pi be represented in this system? Do you have
to explicitly say that you want the first n digits? Or is the padding
scheme smart enough to add pi's digits? In the latter case, how do you
keep checking (in)equality efficient?

Well, there is a large literature on how to represent large ranges with acceptable

efficiency.  Papers by Douglas Priest, and more recently, Jonathan Shewchuk.

The real number pi is not a computer-representable floating-point number in base 2.

There are other numbers that can be represented, like 22/7. 

How does (1/3)*3 give? Does it compare equal to 1? I mean, I would
imagine 1/3 to be 0.333..., with more digits 3 as I need them, coming
from larger and larger 0 paddings of 1 and 3. Supposing this, (1/3)*3
would be 0.999..., with more digits 9 as I need them. Will comparing
with 1 ever finish?

...
> Now if the user specified the kind of arithmetic explicitly, or even
> implicitly by saying "use IEEE754 binary floating point arithmetic
> everywhere" then I could go along with that.

I think you can specify this to some extend. At least
  sage: RealField?
gives mentions opportunities to tweak precision and rounding.

There have been various efforts to axiomatize floating-point numbers.

I don't know if the Sage people refer to any of this, but calling something

RealField does not inspire confidence.

 

You can
use AA if you want exact arithmetic (of real algebraic numbers). I don't
think there is a field of computable reals in Sage.

I am not proposing to use computable reals in any practical sense.

Sort of like Sashimi.   I am happy to look at it, but not eat it.

 

...
> > > > > You seem to think that a number of low precision has some
> > > > > inaccuracy or uncertainty. Which it doesn't.   0.5 is the same
> > > > > number as 0.500000. Unless you believe that 0.5  is  the
> > > > > interval [0.45000000000....01, 0.549999..........9] which you
> > > > > COULD believe  -- some people do believe this.
> > > > > But you probably don't want to ask them about scientific
> > > > > computing.

Again, I think you should keep track of ambient sets here. For sure, the
real numbers 0.5 and 0.5000000 are the same, and real numbers have no
uncertainty (when representing real numbers). However, users do not type
real numbers. They type finite strings of characters that they probably
mean to represent some real number. The _strings_ 0.5 and 0.5000000 are
not the same. Hence it might be reasonable to assign different meaning
to these strings, as long as the user also has a way to precisely
express what is meant. Yet, I'm no expert on scientific computing, so I
will leave this be.

I think I finally agree with you on that last sentence. 

...
> > > > There's 0.5 the real number equal to one divided by two. There's
> > > > also 0.5 the IEEE floating point number, which is a
> > > > representative for an infinite number of real numbers in a small
> > > > interval.
> > >
> > > Can you cite a source for that last statement?  While I suppose
> > > you can decide anything you wish about your computer, the
> > > (canonical?) explanation is that the IEEE ordinary floats
> > > correspond to a subset of the EXACT rational numbers equal to
> > > <some integer> X 2^<some integer> which is not an interval at all.
> > > (there are also inf and nan things which are not ordinary in that
> > > sense)
> > >
> > > So, citation?  (And I don't mean citing a physicist, or someone
> > > who learned his arithmetic in 4th grade and hasn't re-evaluated it
> > > since. A legitimate numerical analyst.)

I am not a numerical analyst, so perhaps the following is naive, but I
would really like know where it is flawed or why my assumptions are
unreasonable.

I call a real number exactly representable when it occurs as an exact
value of a float with some bounded number of bits for its exponent.
Then when using floats to represent real numbers, you have to decide
for each real number x by which exactly representable real number f(x)
it should be represented.

OK, you are positing   (a) not representable real number x and (b) representable number y.

You cannot compute representable e = x-y    because then  x is representable by y+e.

So whatever you mean by x  you cannot represent.  It because pointless to talk about

it.

 

I can think of two desirable properties:
 * exactly representable real numbers should be represented by
   themselves,

All floating point numbers exactly represent themselves.

 

 * the relation <= should be preserved.
These properties yield that for each exactly representable real number y
the set {x real : f(x) = y} is an interval around y.

What's f?

If you are doing interval arithmetic, please read up on it.  the literature is

for the most part quite accessible.

 

Moreover this
interval is contained in every interval (a, b) where a and b are exactly
representable real numbers with a < y and y < b.

Okay, we need some condition on the mantissa and exponent, for if the
mantissa and exponent are allowed to be arbitrary integers, a function
with these properties fails to exists: Let f be a partial function with
the required properties, then the intervals {x real : f(x) = y}, being
contained in the intervals (a,b) with a and b exactly representable
with a < y < b, necessarily contain only y. Since there are only
countable many exactly representable real numbers, each with preimage
containing only itself, the domain of definition of f is countable. In
particular f is not a total function.

huh?  I read this twice, and refuse to read it again.

 

...
> What this means is that the numbers in the system A and b often come
> from the real world and you don't know their true values, just the
> values that you measured.

 

> When you enter the data into the computer you put the values for A and
> b that are the closest numbers to what you believe.  A major question
> in computational linear algebra is how to measure the sensitivity of
> the computed solution to slight changes in the initial data.  The
> numbers used in solving the system are exactly the represented
> floating-point numbers.
>
> > They're representatives for the actual (either unknown or impossible
> > to represent) elements of your data. And, yes, they're also exact
> > rational numbers--once again they can be both--but the point is that
> > the input and output floating point numbers are viewed as
> > perturbations of the (actual) real field elements you care about.
>
> That does not give you or anyone else permission to say  "oh, this
> data might not be right so instead of solving the system as best as
> I can, I will use some off-hand rule to clip bits off, since it's just
> a perturbation of something".   No, you compute the best you can
> with the data you have.

You cannot just assume either that the data is exact.

If you think that there is a better assumption that gives you more information than

computing the best you can, plus condition numbers, you are free to offer it

to the numerical linear algebra specialists who have been laboring under this

model for, oh, 60 years.

 

Also it makes more sense to me not to do the best you can, but the best
you can do efficiently. Trying to do the best you can is of no use if it
means running out of memory or taking a year. If doing things
efficiently involves some clipping, then that's okay with me,
especially if the clipping is configurable.

Since the cost of a floating-point operation is generally constant

these days, you do not get higher efficiency by zeroing-out trailing bits.

Finding results to particular accuracy,  e.g. finding a bounding interval for the root of a

polynomial,  can be done by specifying the polynomial and a width.  In this case

"the best" makes no sense.

Assuming that floating-point numbers are dirty out of ignorance, without context,

is fundamentally a bad idea.

RJF

Regards,

Erik Massop

[Bill Hart]

You didn't tell us the answers to the riddles. Or did I miss them.

[rjf]

On Wednesday, August 20, 2014 7:08:24 PM UTC-7, Bill Hart wrote:

You didn't tell us the answers to the riddles. Or did I miss them.

oops/   4 legs.  calling the tail a leg doesn't make it a leg.

 

[Bill Hart]

By the way, when reading the names sage gives to various implemented things, you should always realise that they are abbreviations.

PowerSeriesRing is an abbreviation for ComputableModelOfAPowerSeriesRingCheckDocumentationForDetailsAndOptions.

The lack of transitivity of equality in sage does bother me at times. But I'm not sure that's one of the cases that bothers me. The semantics are quite useful in practice.

[Robert Bradshaw]

On Fri, Aug 8, 2014 at 6:57 PM, rjf <fat...@gmail.com> wrote:
>
>
> On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:
>>
>> On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
>> >
>> >
>> > On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw wrote:
>> >>
>> >>
>> >>
>> >> The are two representations of the same canonical object.
>> >
>> >
>> > The (computer algebra) use of the term, as in "simplified to a canonical
>> > form" means
>> > the representation is canonical. It doesn't make much sense to claim
>> > that
>> > all these
>> > are canonical: 1+1, 2, 2*x^0, sin(x)^2+cos(x)^2 + exp(0).
>>
>> The point was that there's a canonical domain in which to do the
>> computation.
>
> I have not previously encountered the term "canonical domain". There is
> a CAS literature which includes the concept of simplification to a canonical
> form.
> There is also a useful concept of a zero-equivalence test, whereby E1-E2
> can be shown to be zero, although there is not necessarily a simplification
> routine that will "canonically simplify" E1 to E3 and also E2 to E3.

You have to think beyond just the limited domain of a computer *algebra* system.

If I want to do arithmetic between a \in Z and b \in Z/nZ, I could
either lift b to Z or push a down to Z/nZ. Only one of these maps is
canonical.

>> We also have an object called the ring of integers, but really it's
>> the ring of integers that fits into the memory of your computer.
>> Should we not call it a Ring?
>
> The domain of arbitrary-precision integers is an excellent model of the
> ring of integers. It is true that one can specify a computation that would
> fill up the memory of all the computers in existence. or even all the atoms
> in the (known?) universe. Presumably a well-constructed support system
> will give an error message on much smaller examples. I assume
> that your Real Field operation of division would give an error if the
> result is inexact.

Such a system would be pedantic to the point of being unuseful.

I would argue that most floating point numbers are either (1)
real-world measurements or (2) intermediate results, both of which are
(again, likely) approximations to the value they're representing. When
they are measured/stored, they are truncated due to the "construction
of the [machine], the [sensor] reading the [values], etc." Thus the
analogy.

> On the other hand, the
>
>>
>> I think
>> the analogy is very suitable for a computer system. It can clearly be
>> made much more rigorous and precise.
>
> What you are referring to is sometimes called significance arithmetic,
> and it has been thoroughly discredited.
> Sadly, Wolfram the physicist put it in Mathematica.

Nope, that's not what I'm referring to.

>> Or are you seriously proposing when adding 3.14159 and 1e-100 it makes
>> more sense, by default, to pad the left hand side with zeros (whether
>> in binary or decimal) and return 3.1415900000...0001 as the result?
>
>
> If you did so, you would preserve the identity (a+b)-a = b
>
> If you round to some number of bits, say 53, with a=3.14159 and b=1e-100,
> the left side is 0, and the right side is 1e-100. The relative error in
> the answer
> is, um, infinite.
>
> Now if the user specified the kind of arithmetic explicitly, or even
> implicitly
> by saying "use IEEE754 binary floating point arithmetic everywhere" then
> I could go along with that.

You would suggest that this be IEEE754 be requested by the user
(perhaps globally) before using it? Is that how maxima works? (I think
not.)

>> > So it sounds like you actually read the input as 13/10, because only
>> > then
>> > can
>> > you approximate it to higher precision than 53 bits or whatever. Why
>> > not
>> > just admit this instead of talking
>> > about 1.3.
>>
>> In this case the user gives us a decimal literal. Yes, this literal is
>> equal to 13/10. We defer interpreting this as a 53-bit binary floating
>> point number long enough for the user to tell us to interpret it
>> differently. This prevents surprises like
>>
>> sage: RealField(100)(float(1.3))
>> 1.3000000000000000444089209850
>>
>> or, more subtly
>>
>> sage: sqrt(RealField(100)(float(1.3)))
>> 1.1401754250991379986106491649
>>
>> instead of
>>
>> sage: sqrt(RealField(100)(1.3))
>> 1.1401754250991379791360490256
>>
>> When you write 1.3, do you really think 5854679515581645 /
>> 4503599627370496, or is your head really thinking "the closest thing
>> to 13/10 that I can get given my choice of floating point
>> representation?" I bet it's the latter, which is why we do what we do.
>
> I suspect it is not what python does.

It is, in the degenerate case that Python only has one native choice
of floating point representation. It's also what (to bring things full
circle), Julia does too.

> It is what Macsyma does if you write 1.3b0 to indicate "bigfloat".

You're still skirting the question of whether *you* mean when you write 1.3.

- Robert

[Robert Bradshaw]

On Wed, Aug 20, 2014 at 4:47 PM, rjf <fat...@gmail.com> wrote:
>
> On Saturday, August 9, 2014 7:35:13 PM UTC-7, Erik Massop wrote:
>
>> In fact I don't
>> expect addition, negation, multiplication, or subtraction to be exact
>> either. Indeed, they are not exact:
>> sage: a = 1e-58
>> sage: a
>> 1.00000000000000e-58
>> sage: (1+a)-1
>> 0.000000000000000
>
>
> That's because it is not a field, and you know it. So why call it a field?

Because it is an approximation to a field, and that's just too much of
a mouthful to say everywhere (though
ComputableModelOfAPowerSeriesRingCheckDocumentationForDetailsAndOptions
has a nice ring to it :-).

As Erik said:

> It seems quite common to me to (intentionally) confuse mathematical
> objects with their imperfect realizations.

The actual reason for RealField is probably because that's what Magma
called it, and there wasn't a strong motivation to arbitrarily
disagree with that. (Most other systems don't even have this as an
object to call out or construct.)

You're (intentionally?) missing the point.

Given the linear system Ax = b, where A and b are given in terms of
floating point numbers, one could

(1) Return x' that is the closest (according to some chosen rounding
scheme) to the x exact solution of Ax = b, treating the entries of A
and b as exact rational numbers.

(2) Solve for x approximately, in such a way that x is the exact
solution to a perturbed problem, and one is able to bound the error
wrt the exact solution in terms of condition numbers, etc.

Approach (2) is saying, to quote you, "oh, this data might not be

right so instead of solving the system as best as I can, I will use

some off-hand rule to clip bits off" where the "clipping" here is not
treating the unrepresented lower order bits as all 0. The input is not
taken to be exact (meaning exactly representing the system you really
want to solve).

All (floating point) linear algebra software that I'm aware of take
approach (2).

- Robert

[Erik Massop]

Dear RJF,
Dear list,


This was an interesting read and I have learnt some things:

* I wrongly assumed floating point in computer context to mean sign
bit, fixed range for mantissa, fixed range for exponent, special
cases for -inf, +inf, +0, -0, nan.

* I think Sage's RealField(prec) is more accurately called
fixed-precision floating-point than arbitrary-precision floating-point,
since range of mantissa and exponent are fixed once prec is known.

* Priest's Algorithms for Arbitrary Precision Floating Point Arithmetic
is an interesting read, and made me realize the above.

* I wrongly assumed the "truncated power series" to be QQ[[t]]/t^n for
some fixed non-negative integer n.

* I should not be participating in threads like these. They take too
much time that should be spent on work.


I also have a question:

* Is there a definition of computational domain? If so, do you have some
pointers? That'd be interesting.


Regards,

Erik Massop

> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.

[Harald Schilly]

On Wednesday, August 20, 2014 4:54:35 AM UTC+2, Bill Hart wrote:

Oh what a shame I missed the party. And I don't have time to read all the wonderfully eloquent nonsense in this thread.

So, since we both are the only remaining hangover crowd, I've one last question about Julia:
 

But I think Julia is good for a few years yet. It's one of the few truly 21st century programming languages available. 

All your arguments center around technical arguments: speed, efficiently maintenance of "exponentially" growing complexities in codebases, etc. Those are without question desirable properties. What's missing for me is the other side for the actual humble user of this.

Python's roots are in ABC (which does not stand for abstract base classes, but for the effort to produce a human friendly programming language as a basic successor). Does Julia have such goals?
Second question is about multi-dispatch, and goes in the same direction: Given I have an object behind x, is there a x.[TAB] which lists me what I can "do" with it? Or does this playful way of discovering functionalities work differently?

-- H

[kcrisman]

* I think Sage's RealField(prec) is more accurately called
fixed-precision floating-point than arbitrary-precision floating-point,
since range of mantissa and exponent are fixed once prec is known.

Yes, and although I am not an expert in such things, I always point this out in tutorials.  I don't know how often people make this mistake, though.

* I should not be participating in threads like these. They take too
much time that should be spent on work.

Yup, correct again!  But you had some valuable contributions and learned something, so it wasn't a total waste.

- kcrisman

[Bill Hart]

On Thursday, 21 August 2014 10:58:36 UTC+2, Harald Schilly wrote:

On Wednesday, August 20, 2014 4:54:35 AM UTC+2, Bill Hart wrote:

Oh what a shame I missed the party. And I don't have time to read all the wonderfully eloquent nonsense in this thread.

So, since we both are the only remaining hangover crowd, I've one last question about Julia:
 

But I think Julia is good for a few years yet. It's one of the few truly 21st century programming languages available. 

All your arguments center around technical arguments: speed, efficiently maintenance of "exponentially" growing complexities in codebases, etc. Those are without question desirable properties. What's missing for me is the other side for the actual humble user of this.

Python's roots are in ABC (which does not stand for abstract base classes, but for the effort to produce a human friendly programming language as a basic successor). Does Julia have such goals?

I think it manifestly  does. Julia was inspired by Lisp, Python and a third language that I forgot. The syntax is very pythonic, where the language doesn't fundamentally deviate from the paradigm python belongs to. I would say I'm far more productive in Julia than I could be in Python. It makes easy things easy, intermediate things intermediate and difficult things feasible. Python does not do this. Many difficult things are nearly impossible in Python and many intermediate things are very hard.

Second question is about multi-dispatch, and goes in the same direction: Given I have an object behind x, is there a x.[TAB] which lists me what I can "do" with it?

Yes. Note the dot notation is not used in Julia. It uses a different paradigm (multi-methods, not classes). The Julia approach is much easier and far less restrictive than the usual OO paradigm. 

Given a type A I can do methodswith(A) and it will list all methods that use type A *anywhere* in their prototype.

But the focus in Julia is usually on functions, not objects. So it is more usual in practice to be writing methods(A) where A is the name of a *constructor* for a type (usually also called A). The constructor for A is a generic function with numerous methods, and methods(A) will give you a short list of methods that can be used to construct A (which usually won't have A at all in their prototype -- unless you count the return value -- yes, it would appear explicitly in the prototype in Python and implicitly in the prototype in C++, but not in Julia which doesn't overload on return type).

 

Or does this playful way of discovering functionalities work differently?

It does. But it is extremely practical.

There's also an online help system baked right into Julia. So you can get help on any function, method, type, etc. And this is generated automatically, though you can add your own documentation to that system easily.

It's also trivial to extend Julia's assistance to the user in whatever way you like. One of my favourite things is the ability to be able to overload the show method for any given type and any given value of any type. Doing the former will change the way the type is displayed. E.g.

julia> R, x = PolynomialRing(ZZ, "x")

(Univariate polynomial ring in x over ZZ,x)

Note that the type R is really Poly{ZZ, :x} in Julia's notation. But we overload the show method to display ("Univariate polynomial ring in ", string(:x), " over ", ZZ).

And the polynomial x would usually be displayed as Poly(2, 2, 0xa56c4856727f0a40). We overload the show method to display "x".

I personally think Julia dominates most other languages entirely. I would certainly include Python in that. I can't think of a single feature of Python I would miss in Julia, including easy to use syntax. 

Some languages I think that aren't completely dominated by Julia are Rust, D and Nimrod, and to some degree C++, though it certainly completely dominates more areas of C++ than those other languages.

When I say "some languages" I'm only referring to languages with imperative syntax. There are functional languages that aren't strictly comparable.

 

Note, Julia is not some obscure unheard of language. It gets 150,000 visits a month to its website. And it is gaining popularity rapidly. There are already hundreds of user contributed packages for Julia (which by the way you can download into any Julia session connected to the internet by typing Pkg.add("name of package").

I doubt it would be this popular if it had awkward syntax.

Take a look at my Nemo documentation and see if you feel comfortable with the syntax:

http://www.flintlib.org/nemo.pdf

Note that Nemo is Julia, in the same way as Sage is Python. There is no preparser used in Nemo.

Julia is not the best possible language for computer algebra. I list some of the issues I have with it in my nemo.pdf above. But the Julia developers are actively thinking about some of those. It's even possible that Julia will always suffice and I'll not need to finish writing my own language. There are a couple of things which I consider to be fundamental design decisions in Julia which I think make it impossible to get the best possible performance. But the Julia language designers have rebuffed me on that point, stating that they don't consider them to be fundamental design decisions at all. They are indeed looking at ways of at least mitigating, if not removing those performance deficits in later versions of the language. We'll see. I know how hard it is, because I've done some of those things in my own toy language. Writing languages is very hard and can take decades to get right.

Bill.

-- H

[Bill Hart]

Sorry, my answer to your question about tab completion was incomplete. I posted it before rereading it.

You can define the A.b syntax in Julia if you should so desire. It's essentially just another kind of method overload in Julia.

And then, it supports A.<tab> giving a list of all the things that could follow the dot, just as Python would.

For example Pkg.<tab> yields:

Pkg.Cache         Pkg.Resolve       Pkg.edit          Pkg.register

Pkg.DEFAULT_META  Pkg.Types         Pkg.eval          Pkg.release

Pkg.Dir           Pkg.Write         Pkg.fix           Pkg.resolve

Pkg.Entry         Pkg.add           Pkg.fixup         Pkg.rm

Pkg.Generate      Pkg.available     Pkg.free          Pkg.status

Pkg.META_BRANCH   Pkg.build         Pkg.generate      Pkg.tag

Pkg.Pkg           Pkg.cd            Pkg.installed     Pkg.update

Pkg.Query         Pkg.checkout      Pkg.init

Pkg.Read          Pkg.clone         Pkg.pin

Pkg.Reqs          Pkg.dir           Pkg.publish

But you have to ask yourself how useful that actually is. If you want to write Python in Julia, I suppose it is useful. But the Julia paradigm is different to the Python paradigm. You might actually find that Julia is not generically used that way very much.

[Bill Hart]

I should probably add, since it might excite you Harald, that most syntax in Julia does not actually exist in the language proper. Much of the syntax is defined in library code, e.g. array access, the dot notation and all sorts of other things.

But because of the way Julia operates, all of that does not come with a runtime penalty. You can see from Stefan's examples above that many abstractions in Julia come with zero "runtime" cost. So an array access, or a range syntax, or a dot notation (well, in the next version anyway), etc, can all be overloaded with user defined functions for any type, or typeclass, or union of types or parameterised class of types, or whatever. But those function calls don't result in a runtime penalty when used.

Thus the language is formally extensible at no cost and in an extremely simple way that the user can easily do themselves. Again, total domination of Python.

This is similar to the way Lisp was the first language to have OO, since the language was formally extensible by the user, and so people just wrote macros which did the OO stuff. The only difference is, that was hard to do unless you understood a lot about metaprogramming in Lisp. In Julia, anyone can do it.

julia> getindex(a::Array{Int, 1}, b::Int) = print("test")

getindex (generic function with 140 methods)

julia> a = [1, 2]

2-element Array{Int,1}:

 1

 2

julia> a[2]

test

Bill.

[Volker Braun]

Incidentally, the whole argument type-dependent function dispatch is at the core of the GAP language (though you have to recompile the GAP kernel to change it). And GAP predates K&R C. 

On Thursday, August 21, 2014 4:56:12 PM UTC+1, Bill Hart wrote:

Thus the language is formally extensible at no cost and in an extremely simple way that the user can easily do themselves. 

 

Depends on your point of view. If you want to tie together different third-party libraries then it can quickly become a liability that you can just redefine everything. 

julia> getindex(a::Array{Int, 1}, b::Int) = print("test")

getindex (generic function with 140 methods)

Type this into the REPL on your colleague's computer who just left to get a coffee and see if he ever figures that one out. My bet would be that he's not going to get any work done until the next restart ;-)

[Bill Hart]

On Thursday, 21 August 2014 18:23:48 UTC+2, Volker Braun wrote:

Incidentally, the whole argument type-dependent function dispatch is at the core of the GAP language (though you have to recompile the GAP kernel to change it). And GAP predates K&R C. 

There are other languages which have used it. That's not particularly novel.

It's just the "other" way of doing object orientation.

There is one significant issue with it. To do arithmetic efficiently on composite types, you need to be able to replace c = a*b with mul(c, a, b).

Using the class based approach, this is possible using metaprogramming. We do this in flint. Unfortunately, the flint C++ module takes minutes per file to compile and blows up many version of gcc and all versions of MSVC (internal errors).

So it's not such a huge loss that you can't use template class metaprogramming in a language like Julia (template class metaprogramming has other major issues which make it essentially unusable in practice anyway, such as needing to recompile to change anything).

I don't know a solution to the c = a*b => mul(c, a, b) problem. Nimrod claims to have found some sort of solution, called term rewriting macros. But I really don't know enough about them to evaluate their effectiveness as a solution to this problem.

Julia will have to solve the problem some day, though it will be hard.

In my language I (theoretically) get around it by defining two kinds of type. One is an arithmetic type, which has the above semantics. The other is a usual type which has the usual semantics. But I've not fully implemented the solution, so I don't know to what extent it can really be a solution to the problem. Either way, it will be really, really hard to implement in the "compiler".

Bill.

[Bill Hart]

Oh, I should mention, the c = a*b => mul(c, a, b) problem is probably soluble in D with a much cleaner and easier to compile solution than C++ offers. D metaprogramming is strictly more powerful than that of C++.

[Bill Hart]

Let me be even more specific. If you could come up with a workable solution to the c = a*b => mul(c, a, b) problem, you could write a paper on it and get a lot of kudos for doing so. Even after 50 years, it's still essentially an open problem in modern compiler design.

Move semantics and rvalue assignment in C++ gets you partially the way there. That's a significant advance, technologically, based around the observation that previously, languages were conflating two concepts which should be distinguished. There is probably another equally simple observation that will lead to a solution to this problem. But currently, I don't think it has been solved fully. 

I could be wrong. I'm a mathematician after all, not a computer scientist.

Bill.

[mmarco]

So, would it be thinkable of to move sage from Python to Julia? Sounds like a titanic task, but sounds like if there are so many advantages in Julia with respect to Python, it could be worth it.

[Bill Hart]

Julia can already call Python functions (and I don't mean in some theoretical, technical sense, I mean very practically via an interface designed explicitly for such). So it's not necessary to "move" Sage from Python to Julia. Other Scientific Python projects haven't done this.

There are other reasons why this is not necessary, because of other technical advances that will become available in the coming few years (they are working in labs right now).

Instead, the Sage project should do two things in my opinion:

1) Make Julia part of the Sage ecosystem (it's already had a huge investment from the scientific and statistical communities, so this is a no-brainer if you want to embrace that community)

2) Invest in the technologies that are making Julia successful (jit, dependent typing, metaprogramming, type checking and inference, etc.)

Whether 2 involves rewriting some functionality in Julia, or simply finding ways of adding such functionality to Python is really neither here nor there.

What Sage can't do is just stagnate and ignore progress. If it does, it will be brushed aside as if it wasn't even there, as has happened over and over again in the history of computer algebra! And it's happening to Sage. A few years ago, people at conferences were excitedly demonstrating stuff in Sage. This year, they've moved back to Magma.

As I have already mentioned, the technology behind Julia can in theory do generic programming even faster than C or C++. So it is a strictly superior technology to what has been available in the past. It's a genuine innovation, whether the Julia designers see it that way or not.

The Julia developers are not like me. They like Python and embrace it (many of them are pythonistas). They would never go around claiming Julia is in any way superior to Python. But I am not one of them and can go around making such claims, because in the sense that I mean it, it is true.

Bill.

[D. S. McNeil]

Having been dabbling in Julia myself, I can agree it's definitely
worth a look for people interested in numerical programming and lispy
metaprogramming facilities. For people coming from other languages it
can take a while to figure out where exactly you're supposed to *put*
things, and the type system takes some getting used to, but there are
some really neat features once you push through that.

I don't know how typical my experience is, but their "benchmark" page
put me off the language for quite a while. (Explaining exactly why
would be even more off-topic than this thread already is, but this SO
link http://stackoverflow.com/questions/9968578/speeding-up-julias-poorly-written-r-examples
has arguments on both sides, including from Stefan). But one
commercial with an annoyingly stupid jingle doesn't mean the product
isn't worth buying, and I'm glad I finally decided to push mute.

It's still a very small community -- I woke up one day and was in the
top ten on StackOverflow in terms of answering Julia questions, and
it's not because I've spent a lot of time doing it -- and the
documentation "ecology" is as yet underdeveloped as a result.

There are also a number of things they were surprisingly slow to
embrace (modules, for one; the counterarguments were pretty silly),
but they seem to come around eventually, so I'm optimistic about its
future.

Still not sure I'd recommend it as a language for beginners, or as an
every-day multitool, but what it's good at it's very good at.


Doug

[William A Stein]

On Thu, Aug 21, 2014 at 9:45 PM, Bill Hart <goodwi...@googlemail.com> wrote:
> Julia can already call Python functions (and I don't mean in some
> theoretical, technical sense, I mean very practically via an interface
> designed explicitly for such). So it's not necessary to "move" Sage from
> Python to Julia. Other Scientific Python projects haven't done this.
>
> There are other reasons why this is not necessary, because of other
> technical advances that will become available in the coming few years (they
> are working in labs right now).
>
> Instead, the Sage project should do two things in my opinion:
>
> 1) Make Julia part of the Sage ecosystem (it's already had a huge investment
> from the scientific and statistical communities, so this is a no-brainer if
> you want to embrace that community)
>
> 2) Invest in the technologies that are making Julia successful (jit,
> dependent typing, metaprogramming, type checking and inference, etc.)
>
> Whether 2 involves rewriting some functionality in Julia, or simply finding
> ways of adding such functionality to Python is really neither here nor
> there.
>
> What Sage can't do is just stagnate and ignore progress. If it does, it will
> be brushed aside as if it wasn't even there, as has happened over and over
> again in the history of computer algebra! And it's happening to Sage. A few
> years ago, people at conferences were excitedly demonstrating stuff in Sage.
> This year, they've moved back to Magma.

To fill out this point a little, in my experience (recently, and in
the past), such people choose to use Magma instead of Sage because
Magma is capable of doing X, where X is something important they need
for their work. In some cases, Magma does X slowly, with frequent
crashes, etc.. But it does X nonetheless; whereas, programming Sage
to do X from scratch would take the right expert about 2 years. This
illustrates how Sage is failing to be a viable alternative to Magma.
If Magma does X, and Sage doesn't, and people care about X, then Sage
is not a viable alternative.

Unfortunately, I don't think see "Sage switching to Julia", or Sage
having better support for "(jit, dependent typing, metaprogramming,
type checking and inference, etc.)" in itself would be enough of a
reason for people to have the option to choose Sage over Magma.
Such advantages don't address the **real difficulty** in implementing
a first version of X that works, since implementing X fundamentally
takes reading and understanding papers, working out the right
algorithms, etc. -- it's almost entirely a human expert and design
problem. At least that's my experience watching/making this happen
many times over the last decade. Please keep in mind that much of
your (Bill Hart's) implementation work has been looking at
implementations of various algorithms and coming up with a new
approach that is better in certain ways -- that's certainly very
important and useful, but it is a different problem than doing the
first general implementation of code for computing with some class of
mathematical objects. It's like the difference between building the
first ever car and inventing a better wheel.

Here are two concrete examples of X's where I've seen this exact
situation in the last few months, and for which there is no open
source implementation in any language, despite a demand for one over
the last 10 years:

- (a) quaternion algebras over number fields (topic of Sage Days
61 next week -- https://sites.google.com/site/sagedays61/)

- (b) functional fields, mainly Florian Hess's Riemann-Roch spaces
algorithm; this supports basic algebraic geometry constructions, which
Magma has and Sage doesn't.

There are dozens of other similar examples. They are hard to get
funded with grants, because in all cases the proposal basically says:
"replicate something that's already been done... years ago". Also,
the type of funding one wants is for a person to work full time on the
problem for 1-2 years. A grant proposal that says "We plan to do Y,
and [lots about Y], but we need X [a little about X]," has a chance to
get funded (depending on Y). However, there are many X's (like a and
b above) where I certainly can't make a compelling case for any Y
except "make sage a viable alternative".

If we actually had "real money" to support open source math software
development, I'm certain there could be useful open source
implementations of (a) and (b) above within 2 years. I can think of
specific people who would love to carry out these core implementation
problems, but instead will be teaching a lot of Calculus instead, and
probably not enjoying it.

Anyway, a lack of funding is a big obstruction to creating a viable
open source free alternative to the Ma's for all.

-- William

[Dima Pasechnik]

On 2014-08-21, Bill Hart <goodwi...@googlemail.com> wrote:
> ------=_Part_5_2037022158.1408650350021
> Content-Type: text/plain; charset=UTF-8

>
> Julia can already call Python functions (and I don't mean in some
> theoretical, technical sense, I mean very practically via an interface
> designed explicitly for such). So it's not necessary to "move" Sage from
> Python to Julia. Other Scientific Python projects haven't done this.
>
> There are other reasons why this is not necessary, because of other
> technical advances that will become available in the coming few years (they
> are working in labs right now).
>
> Instead, the Sage project should do two things in my opinion:
>
> 1) Make Julia part of the Sage ecosystem (it's already had a huge
> investment from the scientific and statistical communities, so this is a
> no-brainer if you want to embrace that community)
>
> 2) Invest in the technologies that are making Julia successful (jit,
> dependent typing, metaprogramming, type checking and inference, etc.)
>
> Whether 2 involves rewriting some functionality in Julia, or simply finding
> ways of adding such functionality to Python is really neither here nor
> there.
>
> What Sage can't do is just stagnate and ignore progress. If it does, it
> will be brushed aside as if it wasn't even there, as has happened over and
> over again in the history of computer algebra! And it's happening to Sage.
> A few years ago, people at conferences were excitedly demonstrating stuff
> in Sage. This year, they've moved back to Magma.

Perhaps, it's merely thanks to Simons Foundation making Magma free for all
universities in USA?

Oh, perhaps you mean that Magma has switched over to Julia? :-)

>
> As I have already mentioned, the technology behind Julia can in theory do
> generic programming even faster than C or C++. So it is a strictly superior
> technology to what has been available in the past. It's a genuine
> innovation, whether the Julia designers see it that way or not.

Wait, are you saying there is an optimising Julia compiler available?
Without it, it sounds a bit like claims by Java poeple saying that their
JIT compilers are so great that they beat C compiled with gcc...
My cursory look into Julia only poped out JIT stuff, not a real complier...

>
> The Julia developers are not like me. They like Python and embrace it (many
> of them are pythonistas). They would never go around claiming Julia is in
> any way superior to Python. But I am not one of them and can go around
> making such claims, because in the sense that I mean it, it is true.

What personally put me off Julia when it was announced was that
"We are power Matlab users" in the top of the blog post announcing
Julia project back in Feb 2012. I thought "oh yes, swell, we do need a better
Matlab nchoosek() and 'everything is a matrix' greatness :-)"...

Dima

[Bill Hart]

Conversely, Magma often does things much faster than Sage and without the crashes. I don't think the two projects are strictly comparable on such a metric. It depends on the domain you are talking about. By and large, Magma is faster than Sage for the things I am interested in. 

But you know I'm no Magma fan. Its closed source model is indefensible, even if people believe it is responsible for its success (its not).

 

But it does X nonetheless; whereas, programming Sage
to do X from scratch would take the right expert about 2 years.   This
illustrates how Sage is failing to be a viable alternative to Magma.
If Magma does X, and Sage doesn't, and people care about X, then Sage
is not a viable alternative.

You are absolutely right. And I fully support what you are trying to do to fix this problem (i.e. raise funds).

But let me divert the discussion slightly. It isn't the only issue. Once those experts find a way to develop quaternion algebras or whatever else you want (there was some talk of doing this in Antic for a while, but it died down), the question is, which system will they choose to invest their time in developing such code.

We already have the answer to that question, and unfortunately it is a closed source system.

Sage needs to be attractive to experts as a platform for doing such work. And there are a number of key things they will care about. Ease of development is the main one. Performance is another. I think Julia has it over Sage in both cases. Whilst you do need to be sufficiently smart to understand what a type is. And some of the mathematicians I have met don't meet this standard. But those that do, recognise that the type checker can save an enormous amount of implementation work.

 

Unfortunately, I don't think see "Sage switching to Julia", or Sage
having better support for "(jit, dependent typing, metaprogramming,
type checking and inference, etc.)" in itself would be enough of a
reason for people to have the option to choose Sage over  Magma.

No, it wouldn't. But those things would enable what is currently impossible in Sage.

As I said, I don't think Sage switching to Julia is something people should be talking about. An independent project which uses Julia instead of Python is much more practical.

I already made a choice along those lines. I could have chosen Sage as a higher level language for flint and spent a lot of time speeding up Cython code which does generic arithmetic in Sage. I didn't, because it would not meet any of the criteria I would set for such a project.

Fundamentally, the Cython and Python people aren't willing to adopt technological advances fast enough. And that is not a new problem.

 

Such advantages don't address the **real difficulty** in implementing
a first version of X that works,

Yes, they do. And that's the fundamental point you are missing. To the extent that it is easier to do job X in project Y than project Z, and to the extent that the result will be more impressive in project Y than project Z, people will choose project Y instead of project Z. That's a fact of life. We don't use Algol any more for a reason. It's not that Algol wasn't designed for this sort of thing. It was. We don't use it because it was superceded.

 

since implementing X fundamentally
takes reading and understanding papers, working out the right
algorithms, etc. -- it's almost entirely a human expert and design
problem.    

I do understand your point. I've experienced that myself. 

One data point you should consider is that there is a reason people here in Europe can't get grant funding for working on implementing project X in Sage. It's because Sage is not a European project. I shouldn't say too much about this for obvious reasons. But talk to me about it some time. (If you don't think I've tried to get people here in Europe to use Sage, you are wrong.)

 

The man hour problem, however, is greatly diminished by the tools and infrastructure we use.

If you think that somehow the problem with Sage is that there isn't enough maths research being paid for, you are wrong there too. Plenty of people are getting grants to do the research just fine. It's the implementation of that research (coding) that is the problem. That's the problem that needs addressing. And it is a multifaceted problem. It's not just a money problem.

At least that's my experience watching/making this happen
many times over the last decade.    

I don't discount your experience at all.

 

Please keep in mind that much of
your (Bill Hart's) implementation work has been looking at
implementations of various algorithms and coming up with a new
approach that is better in certain ways -- that's certainly very
important and useful, but it is a different problem than doing the
first general implementation of code for computing with some class of
mathematical objects.    It's like the difference between building the
first ever car and inventing a better wheel.

I'd love to be doing the latter, not the former. But the state of affairs in the computer algebra industry has been such a shambles for such a long time that it would be an utter waste of my time to even begin. People consistently cut off the easy problems, solve them again (often in an inferior way to the way they've been solved in the past), publish a paper on it and call it progress. Once they hit truly difficult problems, they start moaning about how hard it is, and most significantly how hard it is *in their system of choice*. And they give up, and move on to other low hanging fruit.

The *bulk* of implementation that is happening is not useful. That's not to say people don't exist who are capable and have the will to tackle the hard problems. But their tools and experience make the job too hard for them. So they pay students to do it for them, often over, and over, and over again.

 

Here are two concrete examples of X's where I've seen this exact
situation in the last few months, and for which there is no open
source implementation in any language, despite a demand for one over
the last 10 years:

    - (a) quaternion algebras over number fields (topic of Sage Days
61 next week -- https://sites.google.com/site/sagedays61/)

Yeah, I thought of going to that. But I've decided against. My time is better spent on improving Antic. And I have so little time these days.

 

    - (b) functional fields, mainly Florian Hess's Riemann-Roch spaces
algorithm; this supports basic algebraic geometry constructions, which
Magma has and Sage doesn't.

I don't know anything about this, but might have heard it referred to once or twice in the corridors here. 

There are dozens of other similar examples.   They are hard to get
funded with grants, because in all cases the proposal basically says:
"replicate something that's already been done... years ago".  Also,
the type of funding one wants is for a person to work full time on the
problem for 1-2 years.    A grant proposal that says "We plan to do Y,
and [lots about Y], but we need X [a little about X]," has a chance to
get funded (depending on Y).  However, there are many X's (like a and
b above) where I certainly can't make a compelling case for any Y
except "make sage a viable alternative".

That won't work here in Europe. Sage did the right thing in embracing various European projects, like Pari, Gap, Singular, etc. But I can tell you as someone with boots on the ground, that is definitely not leading to a solution to very real problems we have here in Europe.

Sage actually, in the final analysis, risks being completely displaced by a European funded software project.

 

If we actually had "real money" to support open source math software
development, I'm certain there could be useful open source
implementations of (a) and (b) above within 2 years.   I can think of
specific people who would love to carry out these core implementation
problems, but instead will be teaching a lot of Calculus instead, and
probably not enjoying it.

Well, we agree on that much at least.

 

Anyway, a lack of funding is a big obstruction to creating a viable
open source free alternative to the Ma's for all.

It is. As I've said, I fully support your strategy in that regard. Actually, I've already talked to people here in Europe, trying to get people to think about using such technology.

Here are some of the complaints I hear:

* I don't want my code running on machines under the control of someone else's lab

* Why should I pay for software that is (or actually isn't) open source. What's the difference between doing that and downloading some other open source software myself and putting it on machines I control.

* Why would I want to invest in a technology where the money I'm paying is being used for Sage development and not for providing support and services relevant to the technology I paid for. I don't use Sage, and nor can I, so why would I want to invest in that.

Trust me, these are very real issues that affect me here and now, today! And I'm desperate for a solution that will save me having to go down another, technologically inferior route. Out of pure pragmatism, I've completely abandoned projects I was working on that would eventually give me a platform I can use, because I need one quickly. A month ago, if possible!

 

 -- William

[Bill Hart]

Definitely a factor. But a minor one.

William is pretty much right about why people are going back to Magma. But a lot of people are missing the other reasons....

 

Oh, perhaps you mean that Magma has switched over to Julia? :-)

Noooo. Julia is not the issue here. I honestly couldn't give a toss about Julia. I've been very open about the fact that I think there are things wrong with Julia and that I will likely move my work from Julia to some other language in the future (when such a language becomes available). But it is definitely the best language around for this sort of thing.

The issue here is the technology Julia is based on. Without it, you are going to very rapidly fall behind. And I said so very publicly, years before Julia was even announced (I'd have to check, but possibly even before private development of Julia began). At the time my "criticism" was posted to sage-flame and ridiculed.

Anyway, to address your point. Magma is technologically superior to Sage in many ways (and in other ways is completely backward).

Magma has been actively investing in parallel and GPU technology. And I happen to know they employed at least one person working on LLVM and web-based technologies.

Magma has its own language, and doesn't need Julia. They can (just as Sage can), embrace the new technological innovations, without throwing away everything they have.

Magma isn't a static target for which development stopped when Sage began. It is a moving target, that in some areas will always outstrip Sage. The only way to beat it, is to make sure Sage can move even faster than they can in embracing new technology. It needs to be possible to get a Magma beating solution off the ground in half the time Magma can do it. And that is nowhere near the case.

 

>
> As I have already mentioned, the technology behind Julia can in theory do
> generic programming even faster than C or C++. So it is a strictly superior
> technology to what has been available in the past. It's a genuine
> innovation, whether the Julia designers see it that way or not.

Wait, are you saying there is an optimising Julia compiler available?

You don't need an optimising compiler for Julia. Precisely the same backend that is used for the Clang C/C++ compiler underlies Julia. It's over 1 million lines of code. The difference is, it is jit compiled.

 

Without it, it sounds a bit like claims by Java poeple saying that their
JIT compilers are so great that they beat C compiled with gcc...

No, I said generic programming can be faster in Julia than in C. Not that generic programming can be faster than some particular C compiler. The technology is fundamentally superior to C. That's a concept people are going to take a long time to understand and accept, but it's true.

In theory, I can do generic programming in Julia that is faster (at runtime) than you can do in *any* C compiler. And I don't just mean a little bit. I mean a lot. I don't even think the Julia people fully realise this yet (I might be wrong about that, I don't know). And it's not been realised in practice for many cases. But it is true.

On the other hand, because the Julia Jit is based on the technology behind Clang, it can sometimes outperform gcc, just because Clang can. GCC isn't working on their own Jit technology for nothing!

 

The point is, Julia has a massive optimising compiler behind it. People don't realise this. LLVM is NOT a virtual machine, nor is it a peephole optimiser. It is an optimising compiler. It performs most of the optimisations *at the back end* that would normally be done by the compiler *front end*, because it carries "rich type information" right through. Type erasure still occurs at execution time, obviously. But LLVM can optimise what would otherwise be done by the front end of a compiler because it retains enough semantic information to do so.

And LLVM can be used to automatically optimise any program language that uses it for GPU's and, and, and.

However, none of that is what I am referring to. I'm referring to Julia's ability to do zero cost abstractions, meta programming all on top of a Jit. You cannot do any of those things in C. It's strictly more powerful.

My cursory look into Julia only poped out JIT stuff, not a real complier...

Static compilation is coming to Julia. But not because it will speed things up. It already has all the speed of an optimising compiler at the console, as you type!

>
> The Julia developers are not like me. They like Python and embrace it (many
> of them are pythonistas). They would never go around claiming Julia is in
> any way superior to Python. But I am not one of them and can go around
> making such claims, because in the sense that I mean it, it is true.

What personally put me off Julia when it was announced was that
"We are power Matlab users" in the top of the blog post announcing
Julia project back in Feb 2012. I thought "oh yes, swell, we do need a better
Matlab nchoosek() and 'everything is a matrix' greatness :-)"...

They actually wanted to be the next R. But lots of Matlab users saw potential.

In my opinion, Julia is Ma* agnostic. They have a heavy bias towards numerical stuff at present, and they need an algebraic stack. But I'm having no problem writing one.

You can literally load Julia in bare-module mode. It throws all the numeric stuff in the toilet, and you can define your own semantics on top of Julia. Kind of like the prelude in Haskell.

This needs to be done by some serious algebraists. But it's not hard work.

[Bill Hart]

On Friday, 22 August 2014 00:08:20 UTC+2, Bill Hart wrote:

My cursory look into Julia only poped out JIT stuff, not a real complier...

Static compilation is coming to Julia. But not because it will speed things up. It already has all the speed of an optimising compiler at the console, as you type!

I should qualify that. Technically, jit compilation and type inference begin at the function boundary in Julia. So stuff you type at the top level is not optimised. I think that is a mistake, as it's unnecessary. In the language I was working on (using LLVM) type inference and optimisation begin at the top level (which for me is a function). There is a small runtime hit with global variables, but that is all.

In particular, in the language I was working on, loops at the top level are jit compiled. In Julia they are not. But it's almost irrelevant, except when doing timing experiments to see how fast Julia really goes. You might easily think to type a loop into the top level and time it and wonder what all the fuss is about. Put that loop inside a function and you'll soon see.

julia> function sumit(n::Int)

          s = 0

          for i = 1:n

            s += i

          end

          return s

       end

sumit (generic function with 1 method)

julia> @time sumit(1000000000)

elapsed time: 0.913325036 seconds (78060 bytes allocated)

500000000500000000

Bill.

[Bill Hart]

Another thing that I think is terribly important is that the entire technology stack, LLVM, all its libraries and packages, its console, the IJulia web interface and Julia itself, all work out of the box on Windows 64, natively, not as a Cygwin app (they use MinGW64/MSYS2 to build Julia).

And I've personally verified this. After some initial troubles due to the fact that flint.dll is not built correctly on Windows 64, Nemo works perfectly on Windows and from the IJulia interface running on my local Windows 64 machine.

This doesn't solve any problems for Sage, as many mathematical libraries still don't natively support Windows 64 despite it being a decade old. But it is a significant feature of Julia (and of the Nemo project I am writing with Julia).

Bill.

[Bill Hart]

I just found time to actually read the links posted by the OP on Graydon Hoare's history of programming language development, culminating in Julia.

Oh I wish I could write like that! The guy is clearly a genius of immense proportions. 

I truly, truly wish that the computer algebra and number theoretic communities were engaging people like this and leveraging their insight and vision to the fullest.

Bill.

[Bill Hart]

As evidence of my claims that Julia is technically capable of doing generic programming faster than you can do it in C, I offer Fateman benchmark times for various systems.

Julia 0.3 was released today and gives a nice speedup over 0.2, hence the better times for Julia.

Fateman benchmark

================

Magma: 247s

Pari/GP: 161s

Sage: 135 s

Bland (generic C + flint): 42s

Nemo (Julia 0.3 + flint): 36s

The Fateman benchmark for us is to precompute f = 1 + x + y + z + t; p = f^30; and then to time p*(p + 1);

Whilst you can compute the Fateman benchmark faster than this with a special bignum structure and sparse distributed format for the polynomials, the point here is that the generic C version in Bland and the Julia version in Nemo both use flint's fmpz integer type, recursive dense format (!) and the same algorithms.

Another benchmark, called the Pearce benchmark, which computes f = (x + y + 2z^2 + 3t^3 + 5u^5 + 1)^16; g = (u + t + 2z^2 + 3y^3 + 5x^5 + 1)^16; and then times f*g; dropped from 90s down to 72s. This is a sparse benchmark, but we are just using recursive dense representation.

I claim it's possible for Julia to do better still. And I'm sure we'll see that in a later release.

Bill.

[Pierre]

I thought we would have another argument about real numbers this morning, but I see the topic has changed back to Julia. And I have to say that I find the Nemo project very interesting, in spite of the over-aggressive publicity it is receiving here.

You don't often get to see a new piece of software in its infancy, with the opportunity to request features from the developpers so early on. So here's a suggestion for Bill Hart; you write:

>Nemo is Julia, in the same way as Sage is Python

I'm not sure if I fully understand this, but I'm afraid I do. And I don't think it's a good decision. Let me explain.

I'm very jealous of some colleagues in numerical analysis here at my university whose ecosystem is ipython+numpy+matplotlib (pretty much, i may forget something). All three can be installed with easy_install at the prompt; it's modular (don't load matplotlib if you don't plan to plot); it's lightweight.

why shouldn't nemo by a Julia module that people load explicitly?

There are many good reasons why Sage cannot just be a module for Python. However I remember conversations in this forum years ago when people were musing about how great it would be if, by some miracle, Sage could indeed be a (collection of) module(s): it would be so simple and healthy. (Wasn't the Purple Sage project started with this "simplicity" in mind?) And it would be much easier to explain to people what Sage is. This may be a minor point, but if you've followed the conversation on reddit a few hours ago:

http://www.reddit.com/r/math/comments/2e3qla/william_stein_says_sage_has_overall_failed/

then you will see how many misconceptions people have about Sage, and how scary the technology appears to some.

Again there are very good reasons why Sage is what it is, but I see no reason why Nemo should be compelled to take over Julia like Sage has had to take over Python. You even say that no preparsing is needed! I'm not sure if, for example, you want to insist that integer constants be interpreted automatically as elements of ZZ. It seems to me a minor nuisance to have to write a few extra ZZ's, if in exchange the system can be that much simpler. I'm wondering what other people think.

Oh, and I ain't trying Nemo myself until installation is along the lines of "julia --easy_install nemo", sorry... I would love to try it on SageMathCloud if it's available.

best,
Pierre

[Bill Hart]

On Friday, 22 August 2014 11:49:20 UTC+2, Pierre wrote:

I thought we would have another argument about real numbers this morning, but I see the topic has changed back to Julia. And I have to say that I find the Nemo project very interesting, in spite of the over-aggressive publicity it is receiving here.

I wasn't going to mention it here at all. But I decided that on the tail end of a very long thread it wouldn't hurt much. It seems like a lot more people were reading the thread than I thought.

I don't intend to plug Nemo on sage-devel!

 

You don't often get to see a new piece of software in its infancy, with the opportunity to request features from the developpers so early on. So here's a suggestion for Bill Hart; you write:

>Nemo is Julia, in the same way as Sage is Python

I'm not sure if I fully understand this, but I'm afraid I do. And I don't think it's a good decision. Let me explain.

I'm very jealous of some colleagues in numerical analysis here at my university whose ecosystem is ipython+numpy+matplotlib (pretty much, i may forget something). All three can be installed with easy_install at the prompt; it's modular (don't load matplotlib if you don't plan to plot); it's lightweight.

why shouldn't nemo by a Julia module that people load explicitly?

At the moment this is possible. But the bignum issue is really too much inconsistency for me. It's not just that, but the way the operators are defined on the Int's is also not what an algebraist wants. 1/2 returns 0.5 instead of 1/2. And  div(7//2, 1//3) returns 10. Go figure.

You can't have consistent computer algebra when small integer literals behave differently to bigger ones. It just doesn't work.

This will be an issue when it comes to p-adics for example:

a = 1 + 2*65537 + 3^65537^2 + 11*65537^3 + 13*65537^4 + 12345*65537^5 + O(65537^6) will return who knows what. Certainly not what you entered.

Factorisations of number a = 2*3*5^2*17*71*65537*123456543234565433489*3249861239487619238746032103420132947612384712340813246341 will return who knows what. Certainly not what you entered.

This is just not a tolerable situation for computer algebra. So alas, the input has to be changed to bignum by default.

 

It's also a pain in the neck to take output from one system, such as Pari/GP and cut and paste a long polynomial, having to put ZZ(...) around every bignum. It's just not practical.

There are many good reasons why Sage cannot just be a module for Python. However I remember conversations in this forum years ago when people were musing about how great it would be if, by some miracle, Sage could indeed be a (collection of) module(s): it would be so simple and healthy. (Wasn't the Purple Sage project started with this "simplicity" in mind?) And it would be much easier to explain to people what Sage is. This may be a minor point, but if you've followed the conversation on reddit a few hours ago:

http://www.reddit.com/r/math/comments/2e3qla/william_stein_says_sage_has_overall_failed/

I did.

 

then you will see how many misconceptions people have about Sage, and how scary the technology appears to some.

Again there are very good reasons why Sage is what it is, but I see no reason why Nemo should be compelled to take over Julia like Sage has had to take over Python. You even say that no preparsing is needed! I'm not sure if, for example, you want to insist that integer constants be interpreted automatically as elements of ZZ. It seems to me a minor nuisance to have to write a few extra ZZ's, if in exchange the system can be that much simpler. I'm wondering what other people think.

I wish it were that simple. It might be possible to have Julia load in bare-module mode and still make Nemo available as a Julia Pkg though. I'll talk to the Julia developers about it. It would certainly be great if that were possible.

 

Oh, and I ain't trying Nemo myself until installation is along the lines of "julia --easy_install nemo", sorry... I would love to try it on SageMathCloud if it's available.

Yeah me too.

Bill.

 

best,
Pierre

[Pierre]

Thanks for the examples. However, given that in Julia you can re-define about everything (cf your example with the function extracting elements from a list, replaced by " print 'test' "), surely you could change, say, the behaviour of the function / upon loading the Nemo module? (without mention of a bare-module mode or anything) This in order to make the ordinary machine words behave consistently in comparison with bignums. (There would be the occasional problem when someone has written code speciafically expecting a/b to be a float, but will you ever mix code involving floats with your own ?) This would be the exact reverse of the "from __future__ import division" which you see at the beginning of so many python+numpy scripts.

I don't think the factorization of a = 2*3*5^2*17*71*65537*123456543234565433489*3249861239487619238746032103420132947612384712340813246341 is a good example. You can expect the user to known that this will fail, and wrap things in ZZ. (Do not underestimate the users: there are zillions of casual C programmers out there who realize that "int" is not the same as "double"; understanding the difference between Int64 and ZZ is no harder.) Of course Python or Sage users don't have to worry about that. But isn't Julia about better performance at the cost of specifying types?

As for :

>It's also a pain in the neck to take output from one system, such as Pari/GP and cut and paste a long polynomial, having to put ZZ(...) around every bignum. It's just not practical.

i tend to disagree. A little function parse_PARI(string) or whatever would be easy to write.

I also have an unrelated question, still about Nemo. I think in Julia's type system, when A <: B, then the intuition is that elements of type A are also of type B, or am I wrong? Like elements of type Int64 are also of type Integer (or whatever it's called). If so, your supertype of ZZ should not be called Ring, but RingElement. For an element of ZZ is not a ring. Well, abstract types are not types, they are more like type classes, but still with the machine types the logic is as I've said.


Pierre

[Bill Hart]

On Friday, 22 August 2014 14:10:10 UTC+2, Pierre wrote:

Thanks for the examples. However, given that in Julia you can re-define about everything (cf your example with the function extracting elements from a list, replaced by " print 'test' "), surely you could change, say, the behaviour of the function / upon loading the Nemo module? (without mention of a bare-module mode or anything) This in order to make the ordinary machine words behave consistently in comparison with bignums. (There would be the occasional problem when someone has written code speciafically expecting a/b to be a float, but will you ever mix code involving floats with your own ?) This would be the exact reverse of the "from __future__ import division" which you see at the beginning of so many python+numpy scripts.

You can redefine things that have already been defined, but it sometimes incurs a compiler warning. I will have to look into what can be done.

 

I don't think the factorization of a = 2*3*5^2*17*71*65537*123456543234565433489*3249861239487619238746032103420132947612384712340813246341 is a good example. You can expect the user to known that this will fail, and wrap things in ZZ. (Do not underestimate the users: there are zillions of casual C programmers out there who realize that "int" is not the same as "double"; understanding the difference between Int64 and ZZ is no harder.) Of course Python or Sage users don't have to worry about that. But isn't Julia about better performance at the cost of specifying types?

As for :

>It's also a pain in the neck to take output from one system, such as Pari/GP and cut and paste a long polynomial, having to put ZZ(...) around every bignum. It's just not practical.

i tend to disagree. A little function parse_PARI(string) or whatever would be easy to write.

I also have an unrelated question, still about Nemo. I think in Julia's type system, when A <: B, then the intuition is that elements of type A are also of type B, or am I wrong? Like elements of type Int64 are also of type Integer (or whatever it's called). If so, your supertype of ZZ should not be called Ring, but RingElement. For an element of ZZ is not a ring. Well, abstract types are not types, they are more like type classes, but still with the machine types the logic is as I've said.

A <: B either tells you if a type A is of abstract type B or if abstract type A is of abstract type B.

So, whether or not you use RingElement or Ring depends on whether you consider an element of ZZ to be a ring element or whether you think the integers ZZ form a ring.

With the machine types, we really want a :: T where T <: Integer <: Ring. Again with the baremodule idea, we might be able to do this, but it's not possible if you use the standard provided modules, as far as I can tell.

The logic from my point of view is that a type is a collection of values, and as such Integer <: Ring, not Integer <: RingElement. But the distinction is arbitrary.

In the case of finite fields, I had an argument with myself about whether I should use FFElem or FField. I decided on the latter. I would have liked to use FiniteField, but this is the name I wanted to give to the function that constructs an FField type (what some people call a factory I think).

It's a shame that Julia allows a composite type and constructor to have the same name, but a bitstype and constructor can't have the same name and a factory and composite type cannot have the same name. Similarly a module name and the filename can be the same, but the module name and the main type it defines can't be.

Those are all things I would do differently in Julia. And perhaps some of these will change in a future release, who knows.

They are minor nuisances at best. And there are probably important reasons for those decisions. For example, if a type and factory could have the same name, there could be issues with overloading the dot operator for both. And that's something a lot of people seem to want Julia to do, so that they can write Python in Julia.

Bill.

 

Pierre

[Bill Hart]

One interesting thing I should mention is that loading Julia in bare module mode and redefining all the basic operators and so on does not break all of Julia's libraries. They all continue to operate as normal.

This was quite a surprise to me when I first found out, but all those modules have to import Base themselves, and so as far as they are concerned, everything operates just as it would in normal Julia. 

The operators don't get redefined for everybody when those modules import Base, but they just get locally redefined for the modules that use them.

Maybe that's also the same with Python. I guess it must be. But it costs nothing, at least, in Julia.

[Bill Hart]

I think I know how I will solve the Julia package issue for Nemo.

I will definitely make Nemo a standard Julia package. People who want to just do Pkg.add("Nemo") from within Julia will be able to do so. If you have Julia on your machine that should be all that is necessary. It'll automatically check out flint for you and install it. This sort of thing happens with a lot of Julia packages, as many of them depend on external native libraries.

For people who want algebraic semantics for the operators and possibly bignum parsing as default (if it turns out to be workable), there will be an additional AlgebraicBase module which people will be able to download and install in bare module mode.

I think this should keep everyone happy.

[rjf]

On Wednesday, August 20, 2014 9:54:02 PM UTC-7, Robert Bradshaw wrote:

On Fri, Aug 8, 2014 at 6:57 PM, rjf <fat...@gmail.com> wrote:
>
>
> On Thursday, August 7, 2014 10:55:37 PM UTC-7, Robert Bradshaw wrote:
>>
>> On Thu, Aug 7, 2014 at 9:02 AM, rjf <fat...@gmail.com> wrote:
>> >
>> >
>> > On Wednesday, August 6, 2014 8:11:21 PM UTC-7, Robert Bradshaw wrote:
>> >>
>> >>
>> >>
>> >> The are two representations of the same canonical object.
>> >
>> >
>> > The (computer algebra) use of the term, as in "simplified to a canonical
>> > form"  means
>> > the representation is canonical.  It doesn't make much sense to claim
>> > that
>> > all these
>> > are canonical:   1+1, 2,  2*x^0,  sin(x)^2+cos(x)^2 + exp(0).
>>
>> The point was that there's a canonical domain in which to do the
>> computation.
>
> I have not previously encountered the term "canonical domain".  There is
> a CAS literature which includes the concept of simplification to a canonical
> form.
> There is also a useful concept of a zero-equivalence test, whereby E1-E2
> can be shown to be zero, although there is not necessarily a simplification
> routine that will "canonically simplify"  E1  to E3 and also E2 to E3.

You have to think beyond just the limited domain of a computer *algebra* system.

Actually I am thinking in terms of computer representation, not just a CAS.

You appear to be thinking in some extra-computational way that bits are not bits.

There is a quote from Lewis Carroll's Humpty Dumpty, to the effect that

words mean whatever he says they mean,… who, after all, is the master.

You and I apparently disagree about the term "canonical".

 

If I want to do arithmetic between a \in Z and b \in Z/nZ, I could
either lift b to Z or push a down to Z/nZ. Only one of these maps is
canonical.

I don't know about canonical maps.  The term "canonical representation"

makes sense to me.

 

>> We also have an object called the ring of integers, but really it's
>> the ring of integers that fits into the memory of your computer.
>> Should we not call it a Ring?
>
> The domain of arbitrary-precision integers is an excellent model of the
> ring of integers.  It is true that one can specify a computation that would
> fill up the memory of all the computers in existence. or even all the atoms
> in the (known?) universe.  Presumably a well-constructed support system
> will give an error message on much smaller examples.   I assume
> that your Real Field  operation of   division would give an error if the
> result is inexact.

Such a system would be pedantic to the point of being unuseful.

Quite the contrary. IEEE 754 specifies an "inexact" flag.

 

….snip...

 

> or  log(-1)  when you were first introduced to log?
>>
>> Only being able to store 53 significant bits is completely analogous
>> to only being able to read 3 significant (decimal) figures.
>
>
> Actually this analogy is false.  The 3 digits (sometimes 4) from a slide
> rule are the best that can be read out because of the inherent uncertainty
> in the rulings and construction of the slide rule, the human eye reading
> the lines, etc.   So if I read my slide rule and say 0.25  it is because I
> think
> it is closer to 0.25  than 0.24 or 0.26   There is uncertainty there.
> If a floating point number is computed as 0.25, there is no uncertainty in
> the representation per se.  It is 1/4, exactly a binary fraction, etc.
> Now you could use this representation in various ways, e.g.
> 0.25+-0.01    storing 2 numbers representing a center and a "radius"
> or an interval or ..../   But the floating point number itself is simply
> a computer representation of a particular rational number    aaa x 2^bbb
> Nothing more, nothing less.  And in particular it does NOT mean
> that bits 54,55,56... are uncertain.  Those bits do not exist in the
> representation
> and are  irrelevant for ascertaining the value of the number aaa x 2^bbb.
>
> So the analogy is false.

I would argue that most floating point numbers are either (1)
real-world measurements or (2) intermediate results, both of which are
(again, likely) approximations to the value they're representing.

You could assert this, but what is the point?  You might as well assert

that the computer number system consists of the integers 1,2,3, infinity,  because 

(according to George Gamow) that's what some humans use for counting.

When
they are measured/stored, they are truncated due to the "construction
of the [machine], the [sensor] reading the [values], etc." Thus the
analogy.

Since the computer has no inherent way of recording in a floating-point number anything more than a single exact rational number, that is the starting point for arithmetic.  If you want more information about the possible error, you record TWO numbers.   

> On the other hand, the
>
>>
>> I think
>> the analogy is very suitable for a computer system. It can clearly be
>> made much more rigorous and precise.
>
> What you are referring to is sometimes called significance arithmetic,
> and it has been thoroughly discredited.
> Sadly, Wolfram the physicist put it in Mathematica.

Nope, that's not what I'm referring to.

Can you provide a reference for what you ARE referring to?

 

>> Or are you seriously proposing when adding 3.14159 and 1e-100 it makes
>> more sense, by default, to pad the left hand side with zeros (whether
>> in binary or decimal) and return 3.1415900000...0001 as the result?
>
>
> If you did so, you would preserve the  identity  (a+b)-a   =  b
>
> If you round to some number of bits, say 53,  with a=3.14159  and b=1e-100,
> the left side is 0, and the right side  is 1e-100.  The relative error in
> the answer
> is, um, infinite.
>
> Now if the user specified the kind of arithmetic explicitly, or even
> implicitly
> by saying "use IEEE754 binary floating point arithmetic everywhere" then
> I could go along with that.

You would suggest that this be IEEE754 be requested by the user
(perhaps globally) before using it? Is that how maxima works? (I think
not.)

Numbers that appear with a decimal point are read in as the default float of the underlying lisp system, which is, so far as I know, IEEE754 double, in one form or another in the systems in which Maxima generally runs.   There are options for higher precisions in some lisps.  If a number is written as say 1.3b0   then Mxima's software big floats are used.

 

>> > So it sounds like you actually read the input as  13/10, because only
>> > then
>> > can
>> > you  approximate it to higher precision than 53 bits or whatever.   Why
>> > not
>> > just admit this instead of talking
>> > about 1.3.
>>
>> In this case the user gives us a decimal literal. Yes, this literal is
>> equal to 13/10. We defer interpreting this as a 53-bit binary floating
>> point number long enough for the user to tell us to interpret it
>> differently. This prevents surprises like
>>
>> sage: RealField(100)(float(1.3))
>> 1.3000000000000000444089209850
>>
>> or, more subtly
>>
>> sage: sqrt(RealField(100)(float(1.3)))
>> 1.1401754250991379986106491649
>>
>> instead of
>>
>> sage: sqrt(RealField(100)(1.3))
>> 1.1401754250991379791360490256
>>
>> When you write 1.3, do you really think 5854679515581645 /
>> 4503599627370496, or is your head really thinking "the closest thing
>> to 13/10 that I can get given my choice of floating point
>> representation?" I bet it's the latter, which is why we do what we do.
>
> I suspect it is not what python does.

It is, in the degenerate case that Python only has one native choice
of floating point representation. It's also what (to bring things full
circle), Julia does too.

> It is what Macsyma does if you write 1.3b0   to indicate "bigfloat".

You're still skirting the question of whether *you* mean when you write 1.3.

 It is irrelevant in general what I mean; the questions seem to be what is mathematically appropriate and/or what does the user expect.

Depending on what computer system I am using, I expect different semantics for 1.3 --- FORTRAN/LISP/Mathematica/Maxima/MockMMA.

For example, in MockMMA, a system I wrote, 1.3  is exactly 13/10.

In FORTRAN  1.3d0  and 1.3e0 mean possibly different things.

In Mathematica, 1.3 means different things depending on how many zeros

follow it.

1.3000000000000000000000   vs

1.3000000000000000000

etc

RJF

 

- Robert

[rjf]

On Wednesday, August 20, 2014 9:55:56 PM UTC-7, Robert Bradshaw wrote:

You're (intentionally?) missing the point.

Given the linear system Ax = b, where A and b are given in terms of
floating point numbers, one could

(1) Return x' that is the closest (according to some chosen rounding
scheme) to the x exact solution of Ax = b, treating the entries of A
and b as exact rational numbers.

In other words, find the exact solution and then 'round' it.  Impractical

generally. 

(2) Solve for x approximately, in such a way that x is the exact
solution to a perturbed problem, and one is able to bound the error
wrt the exact solution in terms of condition numbers, etc.

Maybe you could find a quote?  Offhand, I would not phrase it

that way.  More like... Find an exact representable x that is the solution to some perturbed version of the original problem.  

Approach (2) is saying, to quote you, "oh, this data might not be
right so instead of solving the system as best as I can, I will use
some off-hand rule to clip bits off" where the "clipping" here is not
treating the unrepresented lower order bits as all 0. The input is not
taken to be exact (meaning exactly representing the system you really
want to solve).

i thought I was restating YOUR hypothesis that the numbers might not be right so why bother messing with all those bits.    I would prefer to treat all the input numbers as exact rationals.  Unfortunately those numbers may describe a system that does not have an exact representable solution,  so you have this "backward error analysis" approach.

 

All (floating point) linear algebra software that I'm aware of take
approach (2).

Except that your description and mine  of  "sort of"  approach (2) differ.

 

- Robert

[Bill Hart]

He means this. In algebra Z/nZ is actually a ring modulo an ideal. Z is the ring, nZ is the ideal.

The elements of Z/nZ are usually written a + nZ. It means precisely this:

a + nZ = { x in Z such that x = a +nk for some k in Z }

So it is a *set* of numbers (actually, in algebra it is called a coset).

When you write Mod(a, n) you really mean a + nZ if you are an algebraist, i.e. the *set* of numbers congruent to a modulo n.

So Z/nZ consists of (co)sets of the form C = Mod(a, n) = a + nZ. If you like Z/nZ is a set of sets (actually a group of cosets in algebra).

There is a canonical map from Z to Z/nZ: we map b in Z to the (co)set C in Z/nZ that contains b. 

If a = b mod n then C = mod(a, n) is the *only* one of these sets which contains b. So the map is canonical. We only have one choice.

But there is no canonical map going the other way. For any such (co)set C, which element of Z are you going to pick? Any of the elements of C will do.

You can make an arbitrary choice, e.g. pick the unique element of C in the range [0, n). Or you could pick the unique element in the range (-n/2, n/2]. But that is not a canonical choice. You had to make an arbitrary choice. Someone else may have made a different choice. (Indeed some people use the latter choice because some algorithms, such as gcd and Jacobi symbols run faster with this choice.)

So canonical map means there is only one way you could define the map. You don't need to tell anyone what your map is, because they *have* to make the same choice as you, as there are no alternatives. That's the meaning of canonical in mathematics.

It's probably not a terminology they teach in Computer Algebra, but it is taught to undergraduates around the world in pure mathematics.

The whole argument here is that because only one direction gives a canonical map, coercion must only proceed in that direction. Otherwise your computer algebra system is making a choice that someone else's computer algebra system does not.

[Jason Grout]

On 8/21/14, 18:08, Bill Hart wrote:
> In theory, I can do generic programming in Julia that is faster (at
> runtime) than you can do in *any* C compiler. And I don't just mean a
> little bit. I mean a lot. I don't even think the Julia people fully
> realise this yet (I might be wrong about that, I don't know). And it's
> not been realised in practice for many cases. But it is true.

Are you talking about the sort of things illustrated in this julia issue
with the various macros they talk about?

https://github.com/JuliaLang/julia/issues/7033

Thanks,

Jason

P.S. I've been following Julia for a while and I'm very intrigued by it.
I've also started keeping an eye on Nemo recently---very interesting!

[rjf]

On Friday, August 22, 2014 2:04:35 PM UTC-7, Bill Hart wrote:

RJF said...

I don't know about canonical maps.  The term "canonical representation"

makes sense to me.

He means this. In algebra Z/nZ is actually a ring modulo an ideal. Z is the ring, nZ is the ideal.

The elements of Z/nZ are usually written a + nZ. It means precisely this:

a + nZ = { x in Z such that x = a +nk for some k in Z }

So it is a *set* of numbers (actually, in algebra it is called a coset).

When you write Mod(a, n) you really mean a + nZ if you are an algebraist, i.e. the *set* of numbers congruent to a modulo n.

So Z/nZ consists of (co)sets of the form C = Mod(a, n) = a + nZ. If you like Z/nZ is a set of sets (actually a group of cosets in algebra).

There is a canonical map from Z to Z/nZ: we map b in Z to the (co)set C in Z/nZ that contains b. 

If a = b mod n then C = mod(a, n) is the *only* one of these sets which contains b. So the map is canonical. We only have one choice.

But there is no canonical map going the other way. For any such (co)set C, which element of Z are you going to pick? Any of the elements of C will do.

You can make an arbitrary choice, e.g. pick the unique element of C in the range [0, n). Or you could pick the unique element in the range (-n/2, n/2]. But that is not a canonical choice. You had to make an arbitrary choice. Someone else may have made a different choice. (Indeed some people use the latter choice because some algorithms, such as gcd and Jacobi symbols run faster with this choice.)

So canonical map means there is only one way you could define the map. You don't need to tell anyone what your map is, because they *have* to make the same choice as you, as there are no alternatives. That's the meaning of canonical in mathematics.

It's probably not a terminology they teach in Computer Algebra,

There are few courses in Computer Algebra.   But the term canonical is used differently.   In fact, two systems could
have different canonical forms   (e.g. in your example, [0,n-1]    or [-(n-1)/2, (n-1)/2]   for representatives of the numbers mod n.

Just as one system could have  1+x+x^2   and another x^2+x+1.

But for each system with a canonical form for polynomials, all polynomials in its canonical form would be
identical bit patterns.   Some systems e.g. Maple go further and have all canonical and equal forms stored
at the same machine address.

 

but it is taught to undergraduates around the world in pure mathematics.

Um, some English speaking students who take, um, certain electives?

 

The whole argument here is that because only one direction gives a canonical map, coercion must only proceed in that direction. Otherwise your computer algebra system is making a choice that someone else's computer algebra system does not.

But that is perfectly all right.  canoncality is with respect to a particular computer system.  At least
my kind of canonicality.

Otherwise we could not live with BigEndian and LittleEndian computers.  But we do.

  

 

... big snip
RJF
 

[Jason Grout]

On 8/21/14, 11:36, Bill Hart wrote:
> You can define the A.b syntax in Julia if you should so desire. It's
> essentially just another kind of method overload in Julia.
>
> And then, it supports A.<tab> giving a list of all the things that could
> follow the dot, just as Python would.

Are you saying you can overload the dot operator? I thought that was
still an open issue (https://github.com/JuliaLang/julia/issues/1974).

Harold: Julia does have a module system, which Bill illustrated, that
uses the dot to access members of the module. And a dot also accesses
fields of a type instance.

Thanks,

Jason

[Bill Hart]

On Saturday, 23 August 2014 04:26:10 UTC+2, jason wrote:

On 8/21/14, 18:08, Bill Hart wrote:
> In theory, I can do generic programming in Julia that is faster (at
> runtime) than you can do in *any* C compiler. And I don't just mean a
> little bit. I mean a lot. I don't even think the Julia people fully
> realise this yet (I might be wrong about that, I don't know). And it's
> not been realised in practice for many cases. But it is true.

Are you talking about the sort of things illustrated in this julia issue
with the various macros they talk about?

https://github.com/JuliaLang/julia/issues/7033

That's not an example of what I have in mind, but it is an example of the sort of thing that is possible with type inference and macros.

Very nice example. Thanks for pointing it out.

(So the Julia developers really do know what they are onto. :-))

 

Thanks,

Jason

P.S. I've been following Julia for a while and I'm very intrigued by it.
  I've also started keeping an eye on Nemo recently---very interesting!

Cool.

 

[Bill Hart]

On Saturday, 23 August 2014 06:09:00 UTC+2, jason wrote:

On 8/21/14, 11:36, Bill Hart wrote:
> You can define the A.b syntax in Julia if you should so desire. It's
> essentially just another kind of method overload in Julia.
>
> And then, it supports A.<tab> giving a list of all the things that could
> follow the dot, just as Python would.

Are you saying you can overload the dot operator?  I thought that was
still an open issue (https://github.com/JuliaLang/julia/issues/1974).

I don't think it has been implemented just yet. I thought it would make 0.3, but it looks like nothing has happened in that direction yet.

But there's absolutely nothing preventing it, except for the fact that the Julia developers don't seem to want a profusion of dots everywhere.

It's not such a useful thing in a language with multimethods, I think.

[Pierre]

> A <: B either tells you if a type A is of abstract type B or if abstract type A is of abstract type B.

So, whether or not you use RingElement or Ring depends on whether you consider an element of ZZ to be a ring element or whether you think the integers ZZ form a ring.

With the machine types, we really want a :: T where T <: Integer <: Ring. Again with the baremodule idea, we might be able to do this, but it's not possible if you use the standard provided modules, as far as I can tell.

I guess to me the only names that would really make sense would be ZZElement <: RREelement <: RingElement (or else rename the builtin types to something else: it's "Integer" that's not the same of a set!! Rename it Integers and everything makes sense). Anyway, I probably need to learn to think differently about this.

I do wonder about one thing. Since it is not possible in Julia to have instances of abstract types, and you have ZZ <: RR with the possibility of doing ZZ(3) as well as RR(3), does it mean that you've created supplementary concrete types, say ZZ_concrete <: ZZ and RR_concrete <: RR or similar, together with a bunch of constructors to hide everything?

If so, I find this quite complicated. Here good old OOP treats this much more elegantly. But perhaps I am mistaken?

Pierre

[Bill Hart]

On Monday, 25 August 2014 09:48:12 UTC+2, Pierre wrote:

> A <: B either tells you if a type A is of abstract type B or if abstract type A is of abstract type B.

So, whether or not you use RingElement or Ring depends on whether you consider an element of ZZ to be a ring element or whether you think the integers ZZ form a ring.

With the machine types, we really want a :: T where T <: Integer <: Ring. Again with the baremodule idea, we might be able to do this, but it's not possible if you use the standard provided modules, as far as I can tell.

I guess to me the only names that would really make sense would be ZZElement <: RREelement <: RingElement (or else rename the builtin types to something else: it's "Integer" that's not the same of a set!! Rename it Integers and everything makes sense). Anyway, I probably need to learn to think differently about this.

Integer is an abstract type in Julia. For example, Int <: Integer.

Int is a concrete type (actually, it's a typealias for either the concrete type Int32 or Int64, depending on what sort of computer you have). Note it's not IntElement! And Integer is not a concrete type, despite Integer being the name of an element of the ring of Integers. There is no Integers in Julia.

You are probably being confused by the fact that the <: notation in Julia can be used with either concrete types or abstract types on both sides.

This is to prevent needing

Level 1: values

Level 2: types

Level 3: typeclasses

Level 4: typesuperclasses

Level 5: typesupersuperclasses

Level 6: typesupersupersuperclasses

...

etc.

All types, no matter what level in Julia are just DataType's. It's just that some can be instantiated and some are just abstract. The abstract ones serve the purpose of typeclasses, mostly. But they can be arbitrary levels in your hierarchy. It's up to you what interpretation you put on your abstract types.

 

The only thing that can't go on the left side of <: is an actual value, like 1. (There's an argument that it should be possible since types can be parameterised on Int's, but that's another story.)

I do wonder about one thing. Since it is not possible in Julia to have instances of abstract types, and you have ZZ <: RR with the possibility of doing ZZ(3) as well as RR(3), does it mean that you've created supplementary concrete types, say ZZ_concrete <: ZZ and RR_concrete <: RR or similar, together with a bunch of constructors to hide everything?

ZZ is not an abstract type. It's a concrete type. Concrete doesn't mean system defined. It means it has values. A ZZ in Nemo is nothing more and nothing less than a flint fmpz_t integer. So it is very concrete. I could very easily have called it fmpz_t if I wanted. I'd have used Integer to name the type, actually, but that was already taken. As was Int, BigInt and numerous other things. So ZZ was what I decided to use.

I used ZZ because fundamentally, a type is a set of values. ZZ is the name of a set of values. Int is the name of a set of values. It's just a name. I used the first sensible name that popped into my head. Sage uses ZZ for the same purpose roughly, though in Sage there is supposedly a separation between types and mathematical domains. And ZZ is technically a mathematical domain (called a parent in Sage) which has an underlying representation whose values have types.

Aldor is another language that conceptually distinguises mathematical domains from their representations.

On the other hand, you can think of ZZ as being blackboard bold Z, and Z stands for Zahl (number) in German.

 

There is no RR type in Nemo yet. But it can just as well stand for Real if you want, instead of Reals. After all Real is already used for the name of a typeclass in Julia.

If so, I find this quite complicated. Here good old OOP treats this much more elegantly. But perhaps I am mistaken?

How can you even express traits elegantly in good old OOP? It's a template metaprogramming nightmare in C++! You have to model traits yourself in Python. In fact, name one "good old OOP" language that supports traits smoothly. Maybe Scala? Good luck fighting with the type checker though.

Now I will admit that after you first brought this up, I have begun to think that FField should be FFElem. This comes about because we don't have a special name for an element of a finite field, like Polynomial for an element of a Polynomial Ring, or Integer for an element of the ring of Integers or Real Number for an element of the Reals. 

I might change that one. But this has nothing to do with how elegant Julia is or how smooth it's OOP paradigm is. It's just a name I made up. I could have called it anything, including Fred!

And given that I was essentially forced to use ZZ instead of my preferred Integer for a bigint, I think there is no great harm done in using FField instead of FFElem. It's just a name.

I certainly won't be changing Ring to RingElement. That would be wrong. Ring is a typeclass, as is Field. Think of it as an interface or trait in other parlance. Roughly speaking, a typeclass is a collection of types. But the name doesn't have to reflect that. It just has to convey some meaning. And Ring does that perfectly well.

T <: Ring

conveys that T is a type that belongs to the Ring typeclass, i.e. it notionally implements the interface required for doing ring operations.

So in Nemo we will have

Level 1: (values) 1, 1/2, x^2 + x + 1, etc

Level 2: (types) ZZ, QQ, Fraction, Residue, FFElem, Poly, PowerSeries

Level 3: (typeclasses) Ring, Field

So typically function signatures will look like

fnname{T <: Z}(a :: T)

where we have:

a is a value

T is a type

Z is a typeclass

Of course nothing is stopping us having more levels in the hierarchy, e.g. Field <: EuclideanDomain <: PID <: UFD <: IntegralDomain <: CommutativeRing <: Ring

I'd really like to see that done in standard OOP. I don't think it can be even remotely done smoothly. In Julia it is this simple

abstract Ring

abstract CommutativeRing <: Ring

abstract IntegralDomain <: CommutativeRing

abstract UniqueFactorisationDomain <: IntegralDomain

abstract PrincipalIdealDomain <: UniqueFactorisationDomain

abstract EuclideanDomain <: PrincipalIdealDomain

abstract Field <: EuclideanDomain

And then I can even do

abstract EuclideanRing <: CommutativeRing

if I want.

Bill. 

 

Pierre

[Pierre]

Thanks again for the explanations. I realize that there were many things about Julia/Nemo which I had misunderstood. Let me explain a little, so I can (i) embarrass myself a little more, and (ii) perhaps give an indication of common mistakes, for anyone who would write a tutorial about Julia/Nemo :-)

Mistake (1): taking it for OOP. When hearing about the way Julia worked, I jumped when I saw it wasn't OO, and then naturally was looking for ways to emulate OOP with Julia (rather than make the effort of a paradigm shift). And with multi-dispatch, i thought OK, if I want to write code that applies to objects of type A, B and C at once, just declare them A <: X, B <: X and C <: X and define f(x :: X). However this gave me the wrong intuition " A <: X is the same as A inherits from X".

At this point, the fact that abstract types could be instantiated was seen as a limitation compared to OOP. And somehow, it is. Except at the same time, the type system does many new things. Apples and oranges, really.

Mistake (2): taking types for sets. This little discussion about names reveals another mistake of mine: thinking of types as sets. I thought ZZ <: QQ for example. That's not even consistent with the previous mistake. In Sage say, i'm guessing QQ is of type Field or something, which inherits from a type Ring or something -- but not the other way around! So mistake (1) should have led me to believe QQ <: ZZ. Oh well.

Mistake (3): thinking about Sage too much. To my discharge, ZZ in Sage is an object, like 1/2, not a type or anything. This has unconsciously contributed to my confusion greatly. Also the fact that there is a coercion from ZZ to QQ, and somewhere I must have taken A <: B to mean "there is a coercion from A to B".

Pierre

[Pierre]

At this point, the fact that abstract types could be instantiated was seen as a limitation compared to OOP. And somehow, it is. Except at the same time, the type system does many new things.

typo : read "could NOT be instantiated".

[Bill Hart]

Hmm. Some more explaining is necessary.

On Monday, 25 August 2014 12:04:16 UTC+2, Pierre wrote:

Thanks again for the explanations. I realize that there were many things about Julia/Nemo which I had misunderstood. Let me explain a little, so I can (i) embarrass myself a little more, and (ii) perhaps give an indication of common mistakes, for anyone who would write a tutorial about Julia/Nemo :-)

Mistake (1): taking it for OOP.

Julia is OOP!

There are two ways to do OOP. The first is with classes and usually single dispatch and the other is with multi-methods. The latter is the Julia way and it is more flexible.

There is a dual notion for data structures as well, and Julia supports that too, giving even greater flexibility. The combination of the two is more powerful than standard OOP, which can be very restrictive in practice.

 

When hearing about the way Julia worked, I jumped when I saw it wasn't OO, and then naturally was looking for ways to emulate OOP with Julia (rather than make the effort of a paradigm shift). And with multi-dispatch, i thought OK, if I want to write code that applies to objects of type A, B and C at once, just declare them A <: X, B <: X and C <: X and define f(x :: X). However this gave me the wrong intuition " A <: X is the same as A inherits from X".

None of that makes sense.

f(x :: X) is the way you write a function prototype in Julia. It corresponds to the following prototype in C.

void f(X x);

X is just a type, like int or float or double or MyType, etc. In Julia, you just write the type on the right, after a double colon, instead of on the left like C.

A <: X does not mean A is an object of type X. It means A is a subtype of X, where both A and X are types (one or both may be abstract types).

The way to write "A is an object of type X" is A :: X.

The way to write "X is a subtype of Y" is X <: Y.

But remember that subtype here *may* mean X is in the typeclass Y.

At this point, the fact that abstract types could be instantiated was seen as a limitation compared to OOP.

No. There is no similar concept in OOP, except for languages like scala which have traits/interfaces/typeclasses. It's an extension over and above what is found in standard OOP languages.

You can also write all your Julia programs without any <: symbols if you want. And in fact you can leave all types out too and let Julia figure them out.

Types and subtyping simply allows you to restrict what is allowed. Otherwise everything is allowed, just like Python.

 

And somehow, it is.

No it is not. You never need to use the <: notation in Julia if you don't want to. You can just stick to values (objects) and types (classes) if you want. Then you have the same thing as standard OOP. You can even restrict yourself to single dispatch if you want, in which case you have OOP exactly. But there's no reason to restrict yourself in this way, because Julia offers all the other goodies on top of that.

 

Except at the same time, the type system does many new things. Apples and oranges, really.

Mistake (2): taking types for sets.

They sort of are. Formally a type is just a set of values. Well, not in a mathematical sense. But in some intuitive sense.

 

Types of course satisfy some calculus just as sets do and the rules of the calculus are somewhat similar, but not quite the same. Technically dependent types correspond exactly to predicates in the predicate logic, via Type theory (an extension of the Curry-Howard isomorphism). But that's all a bit complicated, so we informally think of types as sets.

One big difference is there are plenty of non-computable functions between sets. But the functions between types are always computable. But again, it's not even vaguely necessary to understand this to get the concept of a type. It's intuitively just a name given to the collection of all values which an object of that type can have.

This little discussion about names reveals another mistake of mine: thinking of types as sets. I thought ZZ <: QQ for example.

Yeah, the <: notation doesn't mean contained in. It's not a mathematical operator. It means "subtype of". There are three things the notation can mean:

1) If A <: B where A and B are types, then it means that if a function can take an argument of type B, then an argument of type A will work just as well in its place. In some sense, the set of values of type A is a restriction of the set of values of type B. This is pretty close to a subset meaning though technically it has more to do with the representation in bits on the machine than anything mathematical.

2) If A <: B where A is a type and B is an abstract type/typeclass/interface/trait, it means that a function which requires a type in class B can take the type A. This is not a subset meaning. It means the type A belongs to the typeclass B.

3) If A <: B where A and B are both abstract types/typeclasses/interfaces/traits, it means that a function which requires a type in class B can take a type in class A. This is a essentially a subset meaning, but the sets here are sets of types, not sets of values as in (1).

There aren't many examples of (1) in Julia that I can think of. Trivial ones such as Int <: Int are the only ones that come to mind. I really can't think of any other examples.

In fact, I think it is a general principle in Julia that concrete types cannot subtype one another. This isn't the case in some other languages, but it is the case in Julia. So (1) above essentially never happens, except for T <: T for concrete types T.

 

This is maybe an oddity, but it's not a restriction. You can still cast between types to your heart's content.

That's not even consistent with the previous mistake. In Sage say, i'm guessing QQ is of type Field or something, which inherits from a type Ring or something -- but not the other way around!

Yeah in Julia you'd have QQ <: Field <: Ring.

The first <: in that is of sort (2) above. The second <: is of sort (3).

So mistake (1) should have led me to believe QQ <: ZZ. Oh well.

We have neither QQ <: ZZ nor ZZ <: QQ. Both are concrete types and may not subtype each other. This is not essential, and in fact is not the case in some other languages. But it is the case in Julia.

 

Mistake (3): thinking about Sage too much. To my discharge, ZZ in Sage is an object,

Correct. But classes are essentially objects in Python because everything is an object, pretty much. This really muddies the water. Better to think of classes and objects as being separate concepts.

Also better to think of ZZ as a class in Sage and elements of ZZ as objects, even though actually ZZ happens to be an object too because Python.

Roughly speaking a Julia value has some similarity to a Python object and a Julia type has some similarity to a Python class. That's a very misleading identification, but there is some truth to it.

Sage is a different beast to Python. As I mentioned, ZZ is not meant to be either an object or a class. It is meant to be a mathematical domain. And mathematical domains are being *modeled* in Python using classes, (which just happen to be objects in Python for the LOL's). That's of course a bastardisation of all that is good and pure. But both identifications occur for somewhat valid reasons.

Classes are objects in Python so that you can do things to classes just as you would any object in the language. The official reason is so that you can rename and import classes. But I rather think there's more to it than that.

And mathematical domains in Sage are modeled as classes because that is the most convenient thing in Python to model them with. Unlike Aldor, Python doesn't have a native domain or category concept. So you have to pick something that works well enough.

It is worth noting that mathematical categories cannot in general be modeled accurately with computer science types. There's a fundamental mismatch between category theory and types. This has led some mathematicians to completely redo the foundations of mathematics in terms of homotopy classes instead of category theory. And this can be modeled accurately in terms of Martin-L\"{o}f Type theory. This new foundation for mathematics is called the Univalent Foundation for Mathematics, and was developed primarily by Vladimir Veovodsky based on the Type theory of Per Martin-L\"{o}f.

I would suggest that a good dependent type system is just as good as a Martin-L\"{o}f Type Theory any day. So for most practical purposes we can model mathematics reasonably well with dependent types. We can't quite push the carpet down at the edges. But we can get most of the wrinkles out.

Sage attempts to model a dependent type system using classes in Python. I won't start a huge digression into the manifold problems with this approach. But it involves a lot of work, for one thing.

like 1/2, not a type or anything.

It's better to think of it as a type or class if you want to think of it as anything computational. It happens to be an object too. But that is an accident forced on you by Smalltalk, err I mean Python.

 

This has unconsciously contributed to my confusion greatly. Also the fact that there is a coercion from ZZ to QQ, and somewhere I must have taken A <: B to mean "there is a coercion from A to B".

Ah, no. Coercion is another thing entirely. That is primarily a mathematical concept. Though there is a related computer science concept called conversion and promotion which has nothing to do with the mathematical concept of coercion, but which is the nearest related concept.

Mathematical coercion is a myth. There is no consistent formal foundation for coercion whatsoever.

The mathematical concept of coercion is nothing more than a bunch of incompatible and inconsistent informal notions that mathematicians carry around in their collective heads to relate the informal, inconsistent and incompatible definitions the carry around in their heads.

Any formal and consistent axiomatisation of computer algebra will deviate arbitrarily far from what any given mathematician thinks is reasonable.

Anyway, best to think of A <: B as meaning A is a subtype of B, where in Julia at least B needs to be an abstract type, except in the trivial A <: A case.

In a sense, you do want to have coercions from types in class A to types in class B where A <: B, usually. But it is much more subtle than that, especially in Julia where types can be parametric. Even if A <: B you don't automatically have T{A} <: T{B} in Julia. This would be called a covariant type system, whereas Julia is invariant.

So even if you had ZZ <: QQ in Julia (which you don't), you wouldn't automatically have Poly{ZZ} <: Poly{QQ}. And even if you did, that would have nothing to do with conversions and promotions because the conversion and promotion system can have arbitrary types in Julia and is completely decoupled from the subtyping system.

Decoupling is good because it gives you more flexibility to do whatever you want to do. You give your own meaning to things, and you only define the promotions and conversions you want.

 

Bill.

Pierre

[Pierre]

Julia is OOP!

If it is, how do you do the following in Julia: define a class (i'm using Python parlance) Foo with some attributes, say if x is an instance of Foo it has an integer x.n ; and then define a class Bar that inherits from Foo, that is, has all the attributes from Foo, and some more, for example if y is an instance of Bar it has y.n because it is an instance of Foo, and say it also has a string y.name. In python

class Foo:
    n= 0

class Bar(Foo):
    name= ""

These classes don't even have methods, so it's not exactly complicated OOP. Yet, I'm not sure how to do this (emulate this?) the Julia way. Presumably, not with the <: relation at all. I'm wondering how multi-dispatch comes in, given there are no methods involved.

(btw I guess i'm only trying to supply another example of a novice's misunderstanding of what he's read here and there. With zero actual Julia practice (who installs a language when reading a tutorial?) May it help you deal with users' questions in the future.)

[Robert Bradshaw]

I think this highlights the difference between your perspective and
mine. When I write "the map Z -> Z/nZ in Sage" you seem to be
interpreting that as "a map from representatives of the integers to
representatives of the integers mod n" and I'm thinking "a
representation of the (canonical) map Z -> Z/nZ." Same for Rings,
Fields, etc.

>> but it is taught to undergraduates around the world in pure mathematics.
>
> Um, some English speaking students who take, um, certain electives?

Many more than take courses in computer algebra. I'd say.

>> The whole argument here is that because only one direction gives a
>> canonical map, coercion must only proceed in that direction. Otherwise your
>> computer algebra system is making a choice that someone else's computer
>> algebra system does not.
>
>
> But that is perfectly all right. canoncality is with respect to a
> particular computer system. At least
> my kind of canonicality.

Clearly, despite some overlap, you're a computer scientist, and I'm a
mathematician:

http://en.wikipedia.org/wiki/Canonical

[Robert Bradshaw]

I, for one, find it very convenient (especially compared to using
reflection and other hacks that are needed to reason about type at
runtime in C++ or Java). Not that some languages aren't even better.

> Also better to think of ZZ as a class in Sage and elements of ZZ as objects,
> even though actually ZZ happens to be an object too because Python.

No. ZZ is an object because you might want to do things with it, like
ask for it's properties (is it infinite?) or use it as a parameter
(e.g. the domain of a Map object). Sometimes the Objects (in the
categorical sense) are more interesting than their elements.

If your Objects are parameterized types, I'm curious how you handle in
Julia the fact that R[x] is a Euclidean domain iff R is a field. (As
you mention, trying to use the type hierarchy in Python to model
categorical relationships in Sage has generally lead to pain which is
why we're moving away from that...)

> --
> You received this message because you are subscribed to the Google Groups
> "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-devel+...@googlegroups.com.
> To post to this group, send email to sage-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-devel.
> For more options, visit https://groups.google.com/d/optout.

[Bill Hart]

Sorry if I have that wrong. But when I type ZZ?? into Sage, it pulls up the definition of a class.

Anyway, it's irrelevant because Python.

In Julia types are also able to be passed around as first class objects. Which means you can define multimethods which take types as parameters. I do this all the time in Nemo.

So I can send a type to a function which gives me some properties of the type, if I so wish. I thought about adding properties like is_infinite to types in Nemo. But I found that I had been thinking the wrong way about things all along.

 

or use it as a parameter
(e.g. the domain of a Map object). Sometimes the Objects (in the
categorical sense) are more interesting than their elements.

If you are a category theorist yes. Most mathematicians think categories are just sets.

 

If your Objects are parameterized types, 

I'm curious how you handle in
Julia the fact that R[x] is a Euclidean domain iff R is a field.

I'm curious, is the category of Euclidean domains a full subcategory of Rings?

 

(As
you mention, trying to use the type hierarchy in Python to model
categorical relationships in Sage has generally lead to pain which is
why we're moving away from that...)

Oh I can't think why.

[Bill Hart]

On Monday, 25 August 2014 17:42:46 UTC+2, Pierre wrote:

Julia is OOP!

If it is,

It is. 

how do you do the following in Julia: define a class (i'm using Python parlance) Foo

 type Foo

    data1::T1

    data2::T2

end

with some attributes,

is_infinite(::Type{Foo}) = yes

is_pierrian(::Type{Foo}) = no

 

say if x is an instance of Foo it has an integer x.n ; and then define a class Bar that inherits from Foo, that is, has all the attributes from Foo, and some more, for example if y is an instance of Bar it has y.n because it is an instance of Foo, and say it also has a string y.name. In python

class Foo:
    n= 0

class Bar(Foo):
    name= ""

Tongue-in-cheek solution:

abstract FooBar

 

type Foo <: FooBar

end

type Bar <: FooBar

end

n{T <: FooBar}(::Type{T}) = 0

name(::Type{Bar}) = ""

Though quite obviously the above is clearly not what you meant. 

These classes don't even have methods, so it's not exactly complicated OOP. Yet, I'm not sure how to do this (emulate this?) the Julia way. Presumably, not with the <: relation at all. I'm wondering how multi-dispatch comes in, given there are no methods involved.

You want to be able to inherit fields, but that highlights the difference between Julia and the traditional OOP perspective. In the traditional perspective the focus is on the data. In the Julia paradigm, the focus is on the behaviour, i.e. the functions.

In a very well defined sense, functions and types are dual. Just as you can inherit structure in traditional OOP, you can inherit behaviour/functionality in the multimethod approach.

The key feature of Julia in this regard is the ability to add new methods to a generic function. So instead of classes and attributes becoming the focus of inheritance, in Julia, Generic functions and their (multi)methods become the focus.

The Julia authors claim that this is more flexible than traditional OOP. A lot of the real world restrictions of OOP come from their servitude to the structure inheritance paradigm.

So in short, if you want to write bad Python or C++ in Julia, you can't. But you also can't write good Julia in Python or C++.

In order to see why the other approach to OOP is not a restriction, you have to go back to what you were trying to model in the first place.

You have in mind something like this:

A real world thing, such as a Rectangle and inheriting from it a ColouredRectangle.

Perhaps the Rectangles are part of a windowing system, and as such any object of type Rectangle (I will avoid using the word class because it may be confused with typeclass, which has nothing to do with the concept of a class) will have a width and length and any object of type ColouredRectangle has the same fields as Rectangle, but also has associated with it a colour.

So that describes the structure of Rectangles and ColouredRectangles.

Now I think about what behaviour I want.

I clearly want to be able to resize any rectangle. I also want to be able to compute the area of a rectangle. I'd like to be able to enquire what the current colour of a coloured rectangle is.

Oh and by the way, some time whilst this project is underway, Circles are invented. We want all the same behaviour for those too, according to management.

So in Julia:

abstract Shape 

abstract RecShape <: Shape

type Rectangle <: RecShape

   width::Int

   height::Int

end

type ColouredRectangle <: RecShape

   width::Int

   height::Int

   colour::Int

end

function resize{T <: RecShape}(a :: T, width :: Int, height :: Int)

   a.width = width

   a.height = height

end

getcolour(a::ColouredRectangle) = a.colour

Oh damn, just got the fax from management. They want Circles too. But Circles don't inherit the structure of rectangles at all. They don't even have sides! What to do?

abstract CircShape <: Shape

type Circle <: CircShape

   radius::Int

end

type ColouredCircle <: CircShape

   radius::Int

   colour::Int

end

resize{T <: CircShape}(a::T, radius::Int)

    a.radius = radius

end

getcolour(a::ColouredCircle) = a.colour

Notice how I inherited the behaviour called resize between the different kinds of rectangles. But later, when a type came along with a different structure, I could still add to the same behaviour resize by adding another method to the generic function resize.

And I also didn't need to rewrite the function getcolour either. I only really needed to go back to the original definition and change it to

getcolour(a::Union(ColouredCircle, ColouredRectangle)) = a.colour

Or if there was multiple inheritance (Julia doesn't have it yet) I could add some more nodes in my abstract type graph so that there is ColouredShape <: Shape. Then just have

getcolour{T <: ColouredShape}(a::T) = a.colour

 

Now, one final observation will demonstrate why this is such a useful paradigm:

Most real world computer programs define *far* fewer classes than methods. The ability to inherit structure is not all that important. The ability to inherit behaviour is much, much more important.

The dual notion of OOP focuses on the latter.

And then comes parameterised types. Which is a story for another time.

[Julien Puydt]

Le 25/08/2014 23:26, Bill Hart a écrit :
> You want to be able to inherit fields, but that highlights the difference
> between Julia and the traditional OOP perspective. In the traditional
> perspective the focus is on the data. In the Julia paradigm, the focus is
> on the behaviour, i.e. the functions.

In the Eiffel language, when you define a class foo with a "member
variable" bar, the api just says that bar has such-type -- only the
class writer knows if it's a constant or a function returning the value
after a computation.

Is that the same principle?

Snark on #sagemath

PS: it's off topic on that list, but I'm still glad following the
discussion.

[Robert Bradshaw]

On Mon, Aug 25, 2014 at 1:24 PM, Bill Hart <goodwi...@googlemail.com> wrote:

>> > Correct. But classes are essentially objects in Python because
>> > everything is
>> > an object, pretty much. This really muddies the water. Better to think
>> > of
>> > classes and objects as being separate concepts.
>>
>> I, for one, find it very convenient (especially compared to using
>> reflection and other hacks that are needed to reason about type at
>> runtime in C++ or Java). Not that some languages aren't even better.
>>
>> > Also better to think of ZZ as a class in Sage and elements of ZZ as
>> > objects,
>> > even though actually ZZ happens to be an object too because Python.
>>
>> No. ZZ is an object because you might want to do things with it, like
>> ask for it's properties (is it infinite?)
>
> Sorry if I have that wrong. But when I type ZZ?? into Sage, it pulls up the
> definition of a class.

If a is an instance of class A, typing a?? will give you the definition of A.

> Anyway, it's irrelevant because Python.
>
> In Julia types are also able to be passed around as first class objects.
> Which means you can define multimethods which take types as parameters. I do
> this all the time in Nemo.
>
> So I can send a type to a function which gives me some properties of the
> type, if I so wish. I thought about adding properties like is_infinite to
> types in Nemo. But I found that I had been thinking the wrong way about
> things all along.

Good, because I was worried there for a minute you were praising Julia
for not treating types as first class objects.

>> or use it as a parameter
>> (e.g. the domain of a Map object). Sometimes the Objects (in the
>> categorical sense) are more interesting than their elements.
>
> If you are a category theorist yes. Most mathematicians think categories are
> just sets.
>
>>
>> If your Objects are parameterized types,
>>
>> I'm curious how you handle in
>> Julia the fact that R[x] is a Euclidean domain iff R is a field.
>
> I'm curious, is the category of Euclidean domains a full subcategory of
> Rings?

I think so; the Euclidean property doesn't place any constraint on the
morphisms.

I am still curious how you would express that
UnivariatePolynomialRing{QQ} <: EuclideanDomans but the same doesn't
hold for UnivariatePolynomialRing{ZZ} in Julia.

>> (As
>> you mention, trying to use the type hierarchy in Python to model
>> categorical relationships in Sage has generally lead to pain which is
>> why we're moving away from that...)
>
> Oh I can't think why.

I buy that the type system is more powerful in Julia (yay), but I'm
still skeptical that it is the best way to express the relationship
between objects and categories, and at the same time elements and
objects, and attach useful information to those categories and
objects.

- Robert

[Bill Hart]

On Tuesday, 26 August 2014 00:40:34 UTC+2, Robert Bradshaw wrote:

On Mon, Aug 25, 2014 at 1:24 PM, Bill Hart <goodwi...@googlemail.com> wrote:

>> > Correct. But classes are essentially objects in Python because
>> > everything is
>> > an object, pretty much. This really muddies the water. Better to think
>> > of
>> > classes and objects as being separate concepts.
>>
>> I, for one, find it very convenient (especially compared to using
>> reflection and other hacks that are needed to reason about type at
>> runtime in C++ or Java). Not that some languages aren't even better.
>>
>> > Also better to think of ZZ as a class in Sage and elements of ZZ as
>> > objects,
>> > even though actually ZZ happens to be an object too because Python.
>>
>> No. ZZ is an object because you might want to do things with it, like
>> ask for it's properties (is it infinite?)
>
> Sorry if I have that wrong. But when I type ZZ?? into Sage, it pulls up the
> definition of a class.

If a is an instance of class A, typing a?? will give you the definition of A.

Ah thanks for pointing that out. I now see that 1 is an object of class Integer, whereas ZZ is an object of class Integer_ring.

So where I spoke of ZZ in Sage, please substitute Integer.

(Except that I was using that as an argument for using ZZ as a name for my bignum integer type in Julia; which argument you have now destroyed. I would have liked to use Integer though. But it was already taken.)

 

> Anyway, it's irrelevant because Python.
>
> In Julia types are also able to be passed around as first class objects.
> Which means you can define multimethods which take types as parameters. I do
> this all the time in Nemo.
>
> So I can send a type to a function which gives me some properties of the
> type, if I so wish. I thought about adding properties like is_infinite to
> types in Nemo. But I found that I had been thinking the wrong way about
> things all along.

Good, because I was worried there for a minute you were praising Julia
for not treating types as first class objects.

Nah. But I don't have anything against not making types first class objects. So long as their constructors are.

It also depends somewhat on your definition of first class.

 

>> or use it as a parameter
>> (e.g. the domain of a Map object). Sometimes the Objects (in the
>> categorical sense) are more interesting than their elements.
>
> If you are a category theorist yes. Most mathematicians think categories are
> just sets.
>
>>
>> If your Objects are parameterized types,
>>
>> I'm curious how you handle in
>> Julia the fact that R[x] is a Euclidean domain iff R is a field.
>
> I'm curious, is the category of Euclidean domains a full subcategory of
> Rings?

I think so; the Euclidean property doesn't place any constraint on the
morphisms.

I should have asked whether Fields were. Some references say yes, some no.

It's also really funny that if you look up Category of Euclidean Domains on the intertubes you get two types of pages in general. The first is Wiki articles that use the words category and subcategory in the common speech sense, not in the mathematical sense. Then there is the Sage documentation, which speaks of the Category of Euclidean Domains.

As a concept, it just isn't that useful to category theorists (apparently; I wouldn't know, I'm not a category theorist).

My point is really that having a typeclass hierarchy that revolves around category theory isn't necessarily useful.

Euclidean rings (and more specifically Euclidean domains) are a difficult concept to model in a type hierarchy, because it is not that easy to tell when a ring is a Euclidean ring, and if so, what the Euclidean function is on that ring.

Also, the Euclidean function may only be used implicitly in any algorithms used to compute a gcd. So the Euclidean function is not that useful. And there may be more than one Euclidean structure on a ring, e.g. Z/nZ for n composite.

The upshot of this is that whether or not a ring has a Euclidean function defined on it is probably not something you want to represent in an "interface". Thus, it is not that useful a concept when defining a typeclass or interface for types to belong to.

You could make your interface depend on whether there is a gcd function. But such a function could just return an ideal, for example in a Dedekind domain (Sage even does this; inconsistently). So whether there is a gcd function is neither here nor there.

So I would argue right from the outset that EuclideanDomain is neither a useful concept category theory-wise, nor computationally.

Note that there is no EuclideanDomain typeclass currently in Nemo (I know I did provide it in my big long list of typeclasses that could be added).

I am still curious how you would express that
UnivariatePolynomialRing{QQ} <: EuclideanDomans but the same doesn't
hold for UnivariatePolynomialRing{ZZ} in Julia.

You wouldn't do this. But suppose you would.

You might have:

abstract Ring

abstract EuclideanDomain <: Ring 

type Poly{T <: Ring} <: Ring

    # blah blah

end

type EuclideanPoly{T <: Field} <: EuclideanDomain

   # more blah blah

end

The type construction function PolynomialRing can then be written such

PolynomialRing{T <: Ring}(::Type{T}, S::string) # this is the generic function which gets called if T is a ring but not a field

   # return a Poly

end

PolynomialRing{T <: Field}(::Type{T}, S::string) # this function specialises the previous one, when T happens to be a field

   # return a EuclideanPoly

end

There are various ways you could handle functions that take either an EuclideanPoly or a Poly. For example

myfun(a :: Union(EuclideanPoly, Poly))

or

typealias AnyPoly Union(EuclideanPoly, Poly)

myfun(a :: AnyPoly)

And there are various ways you could handle functions that take elements of a EuclideanDomain or a general Ring. E.g.

myfun{T <: EuclideanDomain}(a :: T)

myfun{T <: Ring}(a :: T)

But I don't think there is any compelling reason to go to all of this trouble. A much better and cleaner way is the following.

abstract Ring

type Poly{T <: Ring} <: Ring

  # blah blah

end

abstract GCDDomain <: Ring

abstract Field <: GCDDomain # it's questionable whether you really want to do this instead of Field <: Ring, but that's another story

typealias EuclideanDomain{T <: Field} Union(GCDDomain, Poly{T})

function myfun{T <: EuclideanDomain}(a :: T)

   # blah

end

This lot all compiles by the way. And you can check it really works:

First define a type which belongs to Field:

type Bill <: Field

end

julia> myfun(Bill())

julia> myfun(Poly{Bill}())

No complaints.

Now define a type which is not a field, but just a ring:

type Fred <: Ring

end

And then:

julia> myfun(Poly{Fred}())

ERROR: `myfun` has no method matching myfun(::Poly{Fred})

julia> myfun(Fred())

ERROR: `myfun` has no method matching myfun(::Fred)

However, in all of the above, you see that I never actually used GCDDomain at all. So it can really be removed if so desired. It's likely you defined a gcd function for some clean set of rings in the first place, which can all be delineated nicely without having some special typeclass for it.

And then, in reality, EuclideanDomain is just a union of various types, just as the particular rings you can define a gcd function for is an eclectic mix, with no particular rhyme or reason to it.

>> (As
>> you mention, trying to use the type hierarchy in Python to model
>> categorical relationships in Sage has generally lead to pain which is
>> why we're moving away from that...)
>
> Oh I can't think why.

I buy that the type system is more powerful in Julia (yay), but I'm
still skeptical that it is the best way to express the relationship
between objects and categories, and at the same time elements and
objects, and attach useful information to those categories and
objects.

That was kind of the point of my earlier post. There is no such best way. Categories don't match any concept in computer science (well, that is not strictly true, but that's a long diversion into n-categories and homotopy types). Usually we can't push the carpet down around all the corners of the room, but we can get it sufficiently flat locally.

 

Bill.

- Robert

[Bill Hart]

Sorry, I had some errors in this next bit.

 

myfun(a :: Union(EuclideanPoly, Poly))

or

typealias AnyPoly Union(EuclideanPoly, Poly)

myfun(a :: AnyPoly)

It should have been

myfun{T <: Union(EuclideanPoly, Poly)}(a :: T)

or 

typealias AnyPoly Union(EuclideanPoly, Poly)

myfun{T <: AnyPoly}(a :: T)

I didn't compile this lot to check. It's a silly way to do things. Trying too hard to be clever.

[Bill Hart]

I am still curious how you would express that
UnivariatePolynomialRing{QQ} <: EuclideanDomans but the same doesn't
hold for UnivariatePolynomialRing{ZZ} in Julia. 

 

So it doesn't get lost in all that noise. Here is an actual Julia session expressing this:

[Bill Hart]

Stoopid google. 

julia> abstract Ring

julia> type ZZ <: Ring

       end

julia> type Poly{T <: Ring} <: Ring

          #blah

       end

julia> abstract Field <: Ring

julia> type QQ <: Field

       end

julia> typealias EuclideanDomain{T <: Field} Union(Field, Poly{T}, ZZ)

Union(Poly{T<:Field},Field,ZZ)

With all that defined, here is what Julia replies in response to various questions about this:

julia> Field <: EuclideanDomain

true

julia> Ring <: EuclideanDomain

false

julia> ZZ <: EuclideanDomain

true

julia> QQ <: EuclideanDomain

true

julia> Poly{ZZ} <: EuclideanDomain

false

julia> Poly{QQ} <: EuclideanDomain

true

[Robert Bradshaw]

Well, some of the places you wrote ZZ at least, as you interchange the
two concepts quite a bit. That explains at least why you sounded so
confused.

So is this Union type extensible, e.g. if I had a third type that was
also a EuclideanDomain but neither a GCDDomain nor a Poly{T}, could I
add it in later? Or mus they all be enumerated upfront?

I do have to admit I was hoping for something a bit more elegant...
Python has union types too (multiple inheritance with an empty body,
though yes the mix of interface and implementation makes that messy).

A more concrete example is Z/nZ for n composite or prime. What is the
type hierarchy for 3 in Z/5Z? In Z/10Z? Do you have values for Z/nZ,
or is Z/nZ just a parameterized type (the way ZZ is the time of 3 \in
Z)?

Yet you seem to be trying to tie the two together. Or is your point
that Nemo gives you a sufficiently strong type system that the locally
flat place is as big as needed?

- Robert