Quaternions, how to speed up computation

[Peter Luschny]

How can I speed up this computation?

    H.<i,j,k> = QuaternionAlgebra(SR, -1, -1)

    def Q(a, b, c, d): return H(a + b*i + c*j + d*k)

    def P(n): return Q(x+1,1,1,1)*P(n-1) if n > 0 else Q(1,0,0,0)

    def p(n): return P(n)[0].list()

    for n in (0..20): print [n], p(n)

    [0] [1]

    [1] [1, 1]

    [2] [-2, 2, 1]

    [3] [-8, -6, 3, 1]

    [4] [-8, -32, -12, 4, 1]

    [5] [16, -40, -80, -20, 5, 1]

    [6] [64, 96, -120, -160, -30, 6, 1]

    ...

With Mathematica this takes 6 sec, with Sage it takes 

hours, (in fact I interrupted after n=15). 

Thanks, Peter

[William Stein]

On Fri, Sep 14, 2018 at 9:00 AM, Peter Luschny <peter....@gmail.com> wrote:
> How can I speed up this computation?
>
> H.<i,j,k> = QuaternionAlgebra(SR, -1, -1)

Do NOT use SR. Instead use QQ(x). Then it is 100x faster than Mathematica...

R.<x> = QQ[]
K = R.fraction_field()
H.<i,j,k> = QuaternionAlgebra(K, -1, -1)

def Q(a, b, c, d): return H(a + b*i + c*j + d*k)
def P(n): return Q(x+1,1,1,1)*P(n-1) if n > 0 else Q(1,0,0,0)

def p(n): return P(n)[0].numerator().list()

for n in (0..20): print [n], p(n)

This takes 0.05s walltime. See

https://share.cocalc.com/share/4a5f0542-5873-4eed-a85c-a18c706e8bcd/support/2018-09-14-091151-quat.sagews?viewer=share

Also, if you want to actually push things further, note that your
function P depends on all smaller values of
P, so you could make this massively faster by making P use
@cached_function, like this:

R.<x> = QQ[]
K = R.fraction_field()
H.<i,j,k> = QuaternionAlgebra(K, -1, -1)

def Q(a, b, c, d): return H(a + b*i + c*j + d*k)

@cached_function

def P(n):
return Q(x+1,1,1,1)*P(n-1) if n > 0 else Q(1,0,0,0)

def p(n): return P(n)[0].numerator().list()

for n in (0..20): print [n], p(n)

This takes 0.01s walltime, and will work for 20 replaced by a much
larger number.
Read about this cached_function decorator:

http://doc.sagemath.org/html/en/reference/misc/sage/misc/cachefunc.html

> def Q(a, b, c, d): return H(a + b*i + c*j + d*k)
> def P(n): return Q(x+1,1,1,1)*P(n-1) if n > 0 else Q(1,0,0,0)
> def p(n): return P(n)[0].list()
> for n in (0..20): print [n], p(n)
>
> [0] [1]
> [1] [1, 1]
> [2] [-2, 2, 1]
> [3] [-8, -6, 3, 1]
> [4] [-8, -32, -12, 4, 1]
> [5] [16, -40, -80, -20, 5, 1]
> [6] [64, 96, -120, -160, -30, 6, 1]
> ...
>
> With Mathematica this takes 6 sec, with Sage it takes
> hours, (in fact I interrupted after n=15).
>
> Thanks, Peter
>
>

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



--
William (http://wstein.org)

[Kolen Cheung]

Hi, I’m trying to translate this Sage syntax to Python syntax (i.e. using sage as a Python library.) But I got stuck even on the first command.

In Sage,

>>> R.<x> = QQ[]
>>> type(R)

<class 'sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category'>

Then I thought I can import it in Python like this:

import sage.rings

# OK

sage.rings.polynomial.polynomial_ring.PolynomialRing_field

# AttributeError

sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category

How would you write the same program in Python using sage as a library? And in general I see many unfamiliar syntax (from Python’s point of view) like R.<x>, QQ[], (0..20), etc. Do you think it is realistic to use sage as a Python library and completely not using sage (as an interpreter) itself?

c.f. Documentation on using sage as a library?

Thanks.

​

[Dima Pasechnik]

Sage uses a preparer to translate things like R.<X>=QQ[] into usual Python 

R=PolynomialRing(...).

In the library code the preparer is not used.

--

[Simon King]

Hi

On 2018-11-19, Kolen Cheung <christi...@gmail.com> wrote:
> Then I thought I can import it in Python like this:
>
> import sage.rings
> # OK
>
> sage.rings.polynomial.polynomial_ring.PolynomialRing_field
> # AttributeError

Admittedly the following is not an ideal solution, but you can do
>>> from sage import all
>>> from sage.rings.polynomial.polynomial_ring import PolynomialRing_field

> sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category

The "..._with_category" classes cannot be imported: These are
dynamically created classes.

> How would you write the same program in Python using sage as a library? And
> in general I see many unfamiliar syntax (from Python’s point of view) like
> R.<x>, QQ[], (0..20), etc.

That's a feature. Syntax such as
R.<x> = QQ[]
was inspired from Magma (if I recall correctly). This and other
syntactical features are very useful in interactive sessions.

It works because of a preparser. If you want to know how to translate an
interactive session into module code, the "preparse" function is your
friend:
sage: preparse("R.<x> = QQ[]")
"R = QQ['x']; (x,) = R._first_ngens(1)"

So, the output of preparse(`interactive sage command`) is `valid python
command`.

Concerning imports of constants, there is the "import_statements"
function:
sage: R.<x> = QQ[]
sage: import_statements(R)
from sage.rings.qqbar import QQx
sage: R.<x,y> = QQ[]
sage: import_statements(R)
from sage.rings.qqbar import QQxy
sage: R.<x,y,z> = QQ[]
sage: import_statements(R)
Traceback (most recent call last):
...
ValueError: no import statement found for 'Multivariate Polynomial
Ring in x, y, z over Rational Field'.

So, "import_statements" is for stuff that is defined in some Sage
module, not for interactively created objects

If I recall correctly, there is a function that for *many* (not all)
interactively created objects in Sage tells how they can be constructed,
but I don't recall the name of that function.

> Do you think it is realistic to use sage as a
> Python library and completely not using sage (as an interpreter) itself?

Sure. If you write code (in the sense of "python or cython module") for
Sage, then you are in fact using Sage as a Python library. To execute it
in Python, I think "from sage import all" should be enough to make it
work (but a less thorough import statement might work too).

Best regards,
Simon

[slelievre]

Mon 2018-11-19 09:41:03 UTC+1, Simon King:

>
> If I recall correctly, there is a function that for *many* (not all)
> interactively created objects in Sage tells how they can be constructed,
> but I don't recall the name of that function.

That is sage_input, which can be used as follows:

    sage: R.<x> = QQ[]

    sage: sage_input(R)

    QQ['x']

    sage: sage_input(x)

    R.<x> = QQ[]

    x

Note that it does not always work:

    sage: K = R.fraction_field()

    sage: sage_input(K)

    Traceback (most recent call last)

    ...

    ValueError: Can't convert Fraction Field of Univariate Polynomial Ring in x over Rational Field to sage_input form

[Kolen Cheung]

But if you're reading carefully this isn't what I'm asking.

[Kolen Cheung]

Thank you both for the answers. However, I’m still stuck:

Focusing on just translating the first line: R.<x> = QQ[]

In sage,

>>> preparse("R.<x> = QQ[]")
"R = QQ['x']; (x,) = R._first_ngens(1)"

>>> import_statements(QQ)
# ** Warning **: several names for that object: Q, QQ
from sage.rings.rational_field import Q

>>> import_statements(R)
from sage.rings.qqbar import QQx

Immediately I have a question: how come the import_statements of QQ is ... import Q, and the import statement of R is ... import QQx? In either case the namespace in question is not imported. Does it mean ... import Q as QQ and ... import QQx as R respectively?

Then I’ve a problem: In sage, I can import them using these import statements. But when entering these 2 import statements in Python,

>>> from sage.rings.rational_field import Q
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-0fb4c40d1a13> in <module>()
----> 1 from sage.rings.rational_field import Q

/usr/lib/python2.7/site-packages/sage/rings/rational_field.py in <module>()
     62     _long_type = int
     63 
---> 64 from .rational import Rational
     65 from .integer import Integer
     66 

/usr/lib/python2.7/site-packages/sage/rings/rational.pyx in init sage.rings.rational (build/cythonized/sage/rings/rational.c:40976)()
     94 
     95 
---> 96 import sage.rings.real_mpfr
     97 import sage.rings.real_double
     98 from libc.stdint cimport uint64_t

/usr/lib/python2.7/site-packages/sage/rings/real_mpfr.pyx in init sage.rings.real_mpfr (build/cythonized/sage/rings/real_mpfr.c:44298)()
----> 1 r"""
      2 Arbitrary Precision Real Numbers
      3 
      4 AUTHORS:
      5 

/usr/lib/python2.7/site-packages/sage/rings/complex_number.pxd in init sage.libs.mpmath.utils (build/cythonized/sage/libs/mpmath/utils.c:8831)()
      4 from .real_mpfr cimport RealNumber
      5 
----> 6 cdef class ComplexNumber(sage.structure.element.FieldElement):
      7     cdef mpfr_t __re
      8     cdef mpfr_t __im

/usr/lib/python2.7/site-packages/sage/rings/complex_double.pxd in init sage.rings.complex_number (build/cythonized/sage/rings/complex_number.c:24212)()
      8 
      9 
---> 10 cdef class ComplexDoubleField_class(sage.rings.ring.Field):
     11     pass
     12 

/usr/lib/python2.7/site-packages/sage/rings/complex_double.pyx in init sage.rings.complex_double (build/cythonized/sage/rings/complex_double.c:24230)()
     96 from cypari2.convert cimport new_gen_from_double, new_t_COMPLEX_from_double
     97 
---> 98 from . import complex_number
     99 
    100 from .complex_field import ComplexField

ImportError: cannot import name complex_number
>>> from sage.rings.qqbar import QQx
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-4809457c5588> in <module>()
----> 1 from sage.rings.qqbar import QQx

/usr/lib/python2.7/site-packages/sage/rings/qqbar.py in <module>()
    512                                     rich_to_bool, richcmp_not_equal,
    513                                     op_EQ, op_NE)
--> 514 from sage.rings.real_mpfr import RR
    515 from sage.rings.real_mpfi import RealIntervalField, RIF, is_RealIntervalFieldElement, RealIntervalField_class
    516 from sage.rings.complex_field import ComplexField

/usr/lib/python2.7/site-packages/sage/rings/real_mpfr.pyx in init sage.rings.real_mpfr (build/cythonized/sage/rings/real_mpfr.c:44298)()
----> 1 r"""
      2 Arbitrary Precision Real Numbers
      3 
      4 AUTHORS:
      5 

/usr/lib/python2.7/site-packages/sage/rings/complex_number.pxd in init sage.libs.mpmath.utils (build/cythonized/sage/libs/mpmath/utils.c:8831)()
      4 from .real_mpfr cimport RealNumber
      5 
----> 6 cdef class ComplexNumber(sage.structure.element.FieldElement):
      7     cdef mpfr_t __re
      8     cdef mpfr_t __im

/usr/lib/python2.7/site-packages/sage/rings/complex_double.pxd in init sage.rings.complex_number (build/cythonized/sage/rings/complex_number.c:24212)()
      8 
      9 
---> 10 cdef class ComplexDoubleField_class(sage.rings.ring.Field):
     11     pass
     12 

/usr/lib/python2.7/site-packages/sage/rings/complex_double.pyx in init sage.rings.complex_double (build/cythonized/sage/rings/complex_double.c:24230)()
     96 from cypari2.convert cimport new_gen_from_double, new_t_COMPLEX_from_double
     97 
---> 98 from . import complex_number
     99 
    100 from .complex_field import ComplexField

ImportError: cannot import name complex_number

Thanks.

​

[John Cremona]

I recommend importing anything you need from sage.all since the details of where everything is might change in time.  This works perfectly well:

$ sage -python  # so we use Sage's python not my system-wide one

Python 2.7.15 (default, Nov  2 2018, 14:32:42) 

[GCC 4.8.5 20150623 (Red Hat 4.8.5-28)] on linux2

Type "help", "copyright", "credits" or "license" for more information.

>>> from sage.all import PolynomialRing, QQ

>>> R = PolynomialRing(QQ,'x')

>>> x = R.gen()

>>> f = x^3+1 # no good as no preparser

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "sage/structure/element.pyx", line 944, in sage.structure.element.Element.__xor__ (build/cythonized/sage/structure/element.c:9006)

RuntimeError: Use ** for exponentiation, not '^', which means xor

in Python, and has the wrong precedence.

# Now we can work with this:

>>> f = x**3+1  # that's better

>>> f.factor()

(x + 1) * (x^2 - x + 1)

Another thing I do is to put my commands into a .py file and run pyflakes on it, as that reveals what things need to be imported.  Then import them all from sage.all.

John Cremona

--

[Dima Pasechnik]

On Mon, Nov 19, 2018 at 11:58 AM John Cremona <john.c...@gmail.com> wrote:
>
> I recommend importing anything you need from sage.all since the details of where everything is might change in time. This works perfectly well:
>
> $ sage -python # so we use Sage's python not my system-wide one
> Python 2.7.15 (default, Nov 2 2018, 14:32:42)
> [GCC 4.8.5 20150623 (Red Hat 4.8.5-28)] on linux2
> Type "help", "copyright", "credits" or "license" for more information.
> >>> from sage.all import PolynomialRing, QQ
> >>> R = PolynomialRing(QQ,'x')
> >>> x = R.gen()
> >>> f = x^3+1 # no good as no preparser
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> File "sage/structure/element.pyx", line 944, in sage.structure.element.Element.__xor__ (build/cythonized/sage/structure/element.c:9006)
> RuntimeError: Use ** for exponentiation, not '^', which means xor
> in Python, and has the wrong precedence.
>
> # Now we can work with this:
>
> >>> f = x**3+1 # that's better
> >>> f.factor()
> (x + 1) * (x^2 - x + 1)
>
> Another thing I do is to put my commands into a .py file and run pyflakes on it, as that reveals what things need to be imported. Then import them all from sage.all.

you can also find out exactly what you need to import:

sage: import_statements("QQ")
from sage.rings.rational_field import QQ

[Kolen Cheung]

I wonder why for me the result is

>>> import_statements(QQ)

# ** Warning **: several names for that object: Q, QQ from sage.rings.rational_field import Q

[Kolen Cheung]

Thanks. This answer my question and I put the tip in https://groups.google.com/d/msg/sage-support/NFtI5XqjQWg/sz5WPcFMAgAJ

[Peter Luschny]

> Hi, I’m trying to translate this Sage syntax to Python syntax (i.e. using sage as a Python library.) But I got stuck even on the first command.

Why do you hijack this thread with a completely different topic?

I'm sure that your question and the answers of the experts are of
interest to many, but they will not expect them under this subject.

Please open in the future new threads for your questions!  Thanks!

[Kolen Cheung]

Because I use exactly the code shown here as an example. And I did open a thread posting a direct question on how to use sage as a library but the question is too general. So I figure I should use an example and so happen I saw this thread with an interesting example.

So this is not hijack. And not a completely different topic. The piece of code is exactly the same.

[Kolen Cheung]

And if anyone is dictating this, that one is surely not you, even if you're the OP. Shut up.

On Monday, November 19, 2018 at 4:43:44 AM UTC-8, Peter Luschny wrote:

[Peter Luschny]

Kolen Cheung:

And if anyone is dictating this, that one is surely not you, even if you're the OP. Shut up.

You reconfirm the obvious, namely that you cannot behave. 

[Dima Pasechnik]

On Mon, Nov 19, 2018 at 12:53 PM Kolen Cheung <christi...@gmail.com> wrote:
>
> And if anyone is dictating this, that one is surely not you, even if you're the OP. Shut up.

Language!!! You might get banned for this...

>
> On Monday, November 19, 2018 at 4:43:44 AM UTC-8, Peter Luschny wrote:
>>
>> > Hi, I’m trying to translate this Sage syntax to Python syntax (i.e. using sage as a Python library.) But I got stuck even on the first command.
>>
>> Why do you hijack this thread with a completely different topic?
>>
>> I'm sure that your question and the answers of the experts are of
>> interest to many, but they will not expect them under this subject.
>>
>> Please open in the future new threads for your questions! Thanks!
>>

[Kolen Cheung]

"shut up" is language issue? You dont know whats bad language.

[Kolen Cheung]

How's your behavior so good? you are a selfish guy that only want to consumea others time and ask question and cannot stand people asking following up question on the same piece of code. You should really have shutted up and may be mark it as spam but not speak up.

[Kolen Cheung]

Do ban me if you think "shut up" is language issue. It must be a joke.

[Kolen Cheung]

And do point out where's the guideline saying I can ask follow up questions like this.

[Kolen Cheung]

Snd if someone say something trying to dictate what can or cannot be asked without backing up with an authority to dictate (such as a guideline), he should really shut up literally and I wonder why this is a 'language issue". Did you guys not literally means shut up when you say shut up? Do you guys use shut up as foul language like you don't mean it? (I wonder why guys will say f**k up to guys if they are not gay, you know, like they don't mean it.)

[Kolen Cheung]

If one don't have an authority to dictate, why can't he told to be shut up? Instead ask him to continue to dictate? What kind of community is this?

[Kolen Cheung]

*can't

[Kolen Cheung]

And if anyone has language issue it surely is OP. hijack vs shut up, which one is more serious? Accusing someone hijacking something is a very serious accusation.

[Dima Pasechnik]

On Monday, November 19, 2018 at 1:14:38 PM UTC, Kolen Cheung wrote:

"shut up" is language issue? You dont know whats bad language.

I bet I can swear in more languages than you do: English, Russian, Dutch, German, Ukrainian, Polish, Yiddish, French, Italian... 

And when I worked in Singapore I even picked up some Hokkein and Tamil bad words, and other local words such as "ang moh" :-)

So, yes, I do consider "shut up" rude. Please don't use it here.

[John Cremona]

Seconded (by a polite Englishman).

[John H Palmieri]

On Monday, November 19, 2018 at 5:32:39 AM UTC-8, Kolen Cheung wrote:

And if anyone has language issue it surely is OP. hijack vs shut up, which one is more serious? Accusing someone hijacking something is a very serious accusation.

Since you asked what kind of community this is, from https://wiki.sagemath.org/Community?action=show&redirect=CodeOfConduct (abridged):

The Sage community is comprised of an international mixture of mathematicians, computer scientists, engineers, researchers, teachers, amateurs, and others with varied backgrounds. This diversity is one of our strengths, but it can also lead to communication problems and unhappiness. People who love working on Sage can more effectively collaborate with others if they follow this code.

1) Be friendly and patient.

2) Be welcoming.

3) Be considerate.

4) Be respectful and polite.

Telling someone to "shut up" and calling someone selfish do not adhere to these guidelines. Maybe accusing someone of hijacking also does not, but in any case, your response was not "friendly and patient" nor "welcoming", etc.

[HG]

sage-support is one of the best list. English is not my mother language and sometimes I don't explain myself well but always people here are very kind and trying to help... It's very important this can long !

Le vendredi 14 septembre 2018 18:00:55 UTC+2, Peter Luschny a écrit :

How can I speed up this computation?

    H.<i,j,k> = QuaternionAlgebra(SR, -1, -1)

    def Q(a, b, c, d): return H(a + b*i + c*j + d*k)

    def P(n): return Q(x+1,1,1,1)*P(n-1) if n > 0 else Q(1,0,0,0)

    def p(n): return P(n)[0].list()

    for n in (0..20): print [n], p(n)

    [0] [1]

    [1] [1, 1]

    [2] [-2, 2, 1]

    [3] [-8, -6, 3, 1]

    [4] [-8, -32, -12, 4, 1]

    [5] [16, -40, -80, -20, 5, 1]

    [6] [64, 96, -120, -160, -30, 6, 1]

    ...

With Mathematica this takes 6 sec, with Sage it takes 

hours, (in fact I interrupted after n=15). 

Thanks, Peter

[Peter Luschny]

You say: "hijack vs shut up, which one is more serious?  Accusing 

someone hijacking something is a very serious accusation."

I certainly didn't mean it as a 'very serious accusation' like

hijacking people or planes. I used the term "hijacked" as a 

technical term like it is defined in the urban dictionary: 

https://www.urbandictionary.com/define.php?term=Thread%20Hijacking

I'm sorry if it reached you differently, in a stronger sense. 

You say: "And not a completely different topic." 

I deny this. It didn't add anything to the topic. My question 

was excellently answered, and in all your posts and those of 

the experts not the slightest new information has been added 

to the original question. Why? Because you have a completely

different topic. In your own words:

"Hi, I’m trying to translate this Sage syntax to Python syntax 

Do you think it is realistic to use sage as a Python library 

and completely not using sage (as an interpreter) itself?"

Does this has anything to do with the question of how to 

calculate efficiently with the Quaternion algebra?

In your own words, a recommendation you gave above to someone 

who wanted to help you: