implement harmonic number function H(n,m) #16671

rwst · 2014-07-17T16:03:35Z

H(n,m) belongs to the holonomic repertoir and should be available in Sage, both numerically and symbolically. For some info, see https://groups.google.com/forum/#!topic/sage-devel/uaA0yU7dZHU and https://en.wikipedia.org/wiki/Harmonic_number

This is actually a defect because we leave a Maxima result undefined:

sage: sum(1/k^3,k,1,n)
gen_harmonic_number(3, n)
sage: gen_harmonic_number?
Object `gen_harmonic_number` not found.

CC: @eviatarbach @dkrenn

Component: symbolics

Keywords: special, log

Author: Ralf Stephan, Armin Straub

Branch/Commit: 9042104

Reviewer: Ralf Stephan, Armin Straub

Issue created by migration from https://trac.sagemath.org/ticket/16671

The text was updated successfully, but these errors were encountered:

fredrik-johansson · 2014-07-17T16:37:59Z

comment:1

For exact values, FLINT has arith_harmonic_number()

rwst · 2014-07-26T13:54:54Z

comment:2

Replying to @fredrik-johansson:

For exact values, FLINT has arith_harmonic_number()

Unfortunately, neither FLINT nor mpmath implements the generalized H(m,n) returned by Maxima in the code example.

rwst · 2014-07-27T16:04:02Z

Branch: u/rws/implement_harmonic_number_function_h_n_m_

rwst · 2014-07-27T16:06:23Z

Commit: 5310ca1

rwst · 2014-07-27T16:06:23Z

Author: Ralf Stephan

rwst · 2014-07-27T16:06:23Z

comment:5

It is already implemented, I'm just waiting for https://groups.google.com/forum/?hl=en#!topic/sage-devel/SM_FLL33t2g as to argument order of the generalized version.

New commits:

`5310ca1`	`16671: implement harmonic number function H(n,m) (modulo argument order)`

nbruin · 2014-07-27T18:32:11Z

comment:6

Replying to @rwst:

It is already implemented, I'm just waiting for https://groups.google.com/forum/?hl=en#!topic/sage-devel/SM_FLL33t2g as to argument order of the generalized version.

New commits:

`5310ca1`	`16671: implement harmonic number function H(n,m) (modulo argument order)`

For maxima_lib you definitely have to implement a special conversion rule (in both special_max_to_sage and special_sage_to_max). It'll be a very simple conversion rule, though, so it should be quick to write. There are plenty of examples of entries there already that have to do more complicated manipulations, so just swapping arguments will be easy to figure out.

For conversions to other maxima interfaces a good start may the special examples in maxima_lib. It seems that implementing _maxima_init_evaled_ might do the trick for conversion to maxima.

The workhorse for the other direction is sage.calculus.calculus.symbolic_expression_from_maxima_string which defines regexes for polylog etc., so it seems that any special conversion rules would need to be added there.

Your conversion is not that strange though. You just need to swap arguments, which means no special conversion at all as long as you can get it converted to a sage function that expects the arguments in the required order and produces the desired symbolic expression.

It seems that as long as you can get an entry into sage.symbolic.pynac.symbol_table.get('maxima', {}) that has a argument-swapping function associated with it, you might be good to go. sage.symbolic.pynac.register_symbol might do the trick for that.

sagetrac-git · 2014-07-28T09:05:56Z

Changed commit from 5310ca1 to ebc93ec

sagetrac-git · 2014-07-28T09:05:56Z

Branch pushed to git repo; I updated commit sha1. New commits:

`61572c9`	`16671: swap arguments wrt maxima(_lib); add more documentation`
`ebc93ec`	`16671: fix doc typo`

rwst · 2014-07-28T09:08:58Z

comment:8

Thanks, that was absolutely helpful. I'll put a link on the symbolics wiki page.

Please review.

rwst · 2014-07-28T09:08:58Z

Changed keywords from none to special, log

nbruin · 2014-07-28T16:19:08Z

comment:9

A few remarks.

in the maxima_harmonic regular expression you seem to be trying to find the argument by matching subexpressions. I'm afraid that won't work reliably. If you put in arguments that use parentheses and commas itself and you'll find the regex gets it wrong. Indeed, the example you've taken it from is equally broken:

sage: maxima_calculus(psi(psi(x,x),x))._sage_()
TypeError: unable to make sense of Maxima expression 'psi[psi(x,x)](x)' in Sage

(never mind that this is a ValueError and not a TypeError). The maxima_lib specific translator does not have trouble with this:

sage: max_to_sr(maxima_calculus(psi(psi(x,x),x)).ecl())
psi(psi(x, x), x)

Regular languages are not sufficient to describe expression languages (the nested parentheses can't be expressed in a regular language--they need a stack or something similar). Can't you get that conversion by giving an appropriate sage.symbolic.pynac.register_symbol? then you get to use the parentheses parser that is already in place.

for this line, you took a way too complicated example.

sage.functions.log.harmonic_number : lambda N,X : [[mqapply],[[max_harmo, max_array],X],N],

from

sage: maxima_calculus("gen_harmonic_number(a,b)").ecl()
<ECL: (($GEN_HARMONIC_NUMBER SIMP) $A $B)>

you see that this should simply be

sage.functions.log.harmonic_number : lambda N,X : [[max_harmo],X,N],

With your version one gets:

sage: maxima_calculus(EclObject([["mqapply"],[["$gen_harmonic_number","array"],10],20]))
gen_harmonic_number[10](20)

which is not the desired form, and you'll find this form can't be converted back either.

nbruin · 2014-07-28T19:50:05Z

comment:10

By the way, what's your most convincing argument to use an argument order opposite to that of maxima? I think your ordering agrees with what maple and mathematica have chosen, so it may well be that that is indeed the ordering that has most consensus. The wikipedia page is ambiguous since it lists H_{n,m}=H_n^m=H_m(n). Referencing explicitly which convention we follow will make for a better motivated and more stable interface.

rwst · 2014-07-28T20:01:41Z

comment:11

The paper I looked at had H_n^(m), and Wikipedia has mostly H_n,m and H_n^(m), but the paper would not be relevant. WP is, and it's listed as reference.

sagetrac-git · 2014-07-29T14:45:59Z

Branch pushed to git repo; I updated commit sha1. New commits:

`0b15dd1`	`Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_`
`0c93bfe`	`16671: fix maxima_lib interface`

sagetrac-git · 2014-07-29T14:45:59Z

Changed commit from ebc93ec to 0c93bfe

rwst · 2014-07-29T14:46:40Z

comment:13

Replying to @nbruin:

Regular languages are not sufficient to describe expression languages (the nested parentheses can't be expressed in a regular language--they need a stack or something similar). Can't you get that conversion by giving an appropriate sage.symbolic.pynac.register_symbol? then you get to use the parentheses parser that is already in place.

Is there ready mechanics for this? I can't find one. Do you expect me to implement it?

you see that this should simply be
sage.functions.log.harmonic_number : lambda N,X : [[max_harmo],X,N],

Thanks, included.

nbruin · 2014-07-29T15:33:34Z

comment:14

Replying to @rwst:

Replying to @nbruin:

Can't you get that conversion by giving an appropriate sage.symbolic.pynac.register_symbol?

Is there ready mechanics for this? I can't find one. Do you expect me to implement it?

No, just use it. If you do search_src("register_symbol") you get a few hits on its use. The following seems to basically work without your patch so I hope you can adjust it for your purposes:

sage: var('x,t,s')
(x, t, s)
sage: function('g')
g
sage: def swap(a,b): return g(b,a)
sage: sage.symbolic.pynac.register_symbol(swap,{'maxima':'gen_harmonic_number'})
sage: sum(1/t^2,t,1,s)
g(s, 2)

That way, the parsing won't be more broken than for other functions (and it seems we're doing fine for most expressions nowadays). The regex for psi isn't really wrong, because it specifies there should be only digits between the square brackets. It is perhaps a little overly restrictive, though:

sage: psi(x,s)
psi(x, s)
sage: maxima_calculus(psi(x,s))
psi[x](s)
sage: maxima_calculus(psi(x,s))._sage_()
TypeError: unable to make sense of Maxima expression 'psi[x](s)' in Sage
sage: from sage.interfaces.maxima_lib import max_to_sr
sage: max_to_sr(maxima_calculus(psi(x,s)).ecl()) #this doesn't rely on parsing or regexes
psi(x, s)

I think "parsing" psi via regex was just a quick fix around the fact that the expression parser doesn't know what to do with square brackets (and since they're rare, you can get away with not really parsing them properly). In your case, round brackets are so common that you'll quickly get caught if you don't deal with them properly. Luckily, the parser knows how to deal with them, you just have to swap the arguments.

dkrenn · 2015-11-29T11:29:07Z

comment:68

sage -t src/sage/symbolic/function.pyx
**********************************************************************
File "src/sage/symbolic/function.pyx", line 394, in sage.symbolic.function.Function.__call__
Failed example:
    binomial(Qp(2)(9),5)
Exception raised:
    Traceback (most recent call last):
      File "/local/dakrenn/sage/6.9/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/local/dakrenn/sage/6.9/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.symbolic.function.Function.__call__[20]>", line 1, in <module>
        binomial(Qp(Integer(2))(Integer(9)),Integer(5))
      File "sage/symbolic/function.pyx", line 851, in sage.symbolic.function.GinacFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:9858)
        res = super(GinacFunction, self).__call__(*args, **kwds)
      File "sage/symbolic/function.pyx", line 998, in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:11421)
        res = super(BuiltinFunction, self).__call__(
      File "sage/symbolic/function.pyx", line 511, in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:7288)
        res = g_function_eval2(self._serial, (<Expression>args[0])._gobj,
      File "sage/symbolic/pynac.pyx", line 187, in sage.symbolic.pynac.pyExpression_to_ex (build/cythonized/sage/symbolic/pynac.cpp:4108)
        raise TypeError("function did not return a symbolic expression or an element that can be coerced into a symbolic expression")
    TypeError: function did not return a symbolic expression or an element that can be coerced into a symbolic expression

rwst · 2015-11-29T13:37:55Z

comment:69

Yes, this is because of #17790, and for that reason this ticket is pending (but needs review for everything else).

arminstraub · 2016-08-17T22:33:37Z

Changed branch from u/rws/16671 to u/arminstraub/ticket/16671/harmonic_numbers

arminstraub · 2016-08-17T22:33:37Z

comment:71

I have merged the code with the current codebase and simplified Function_harmonic_number_generalized._eval_ and Function_harmonic_number._eval_ (please check that I didn't oversimplify).

When merging the code into the current codebase, there was the following conflict in symbolic/function.pyx:

<<<<<<< HEAD
=======
                # There is no natural coercion from QQbar to the symbolic ring
                # in order to support
                #     sage: QQbar(sqrt(2)) + sqrt(3)
                #     3.146264369941973?
                # to work around this limitation, we manually convert
                # elements of QQbar to symbolic expressions here
                from sage.rings.qqbar import QQbar, AA
                from sage.rings.padics.padic_generic_element import pAdicGenericElement
                nargs = [None]*len(args)
                for i in range(len(args)):
                    carg = args[i]
                    if (isinstance(carg, Element) and
                            ((<Element>carg)._parent is QQbar or
                            (<Element>carg)._parent is AA or
                            isinstance(carg, pAdicGenericElement))):
                        nargs[i] = SR(carg)
                    else:
                        try:
                            nargs[i] = SR.coerce(carg)
                        except Exception:
                            raise TypeError, "cannot coerce arguments: %s"%(err)
                args = nargs
>>>>>>> FETCH_HEAD

It was unclear to me to which degree the underlying code had changed, which is why this code for working with QQbar over the symbolic ring is not merged into the branch I posted.

In any case, it seems to me that any such change might be better suited for a separate ticket because it is not specific to harmonic numbers. Consequently, I also removed from symbolic/function.pyx the doctest

            sage: binomial(Qp(2)(9),5)
            126

as well as from functions/log.py the doctests:

            sage: harmonic_number(Qp(5)(10),1)
            4*5^-1 + 2 + 2*5 + 4*5^2 + 3*5^3 + 5^4 + ...
            sage: harmonic_number(Qp(5)(10),2)
            4*5^-1 + 3 + 5 + 3*5^2 + 2*5^3 + 3*5^5 + ...
...
            sage: harmonic_number(Qp(5)(3))
            1 + 5 + 4*5^2 + 4*5^4 + 4*5^6 + ...

I think that a separate ticket should be used for uniformly fixing the coercion/conversion of p-adics to the symbolic ring. Since the present code doesn't touch these issues, I have removed the dependency on #17790.

By the way, I liked the symbolic sum you were using for computing harmonic_number(n) when n is a rational number. While testing this, I noticed that sum(1/x/(1/5+x),x,1,infinity) results in an error. On the other hand, using the present code, we have:

sage: harmonic_number(1/5)
euler_gamma + psi(6/5)
sage: _.full_simplify()
-1/4*(sqrt(5) + 1)*log(1/2*sqrt(5) + 5/2) + 1/4*(sqrt(5) - 1)*log(-1/2*sqrt(5) + 5/2) - 1/10*pi*sqrt(10*sqrt(5) + 25) - log(5) + 5

The only other changes I made were small touches to the doc strings. A quick sage -t src/sage/functions did not reveal any troubles.

Thank you for the nice work! It would be great to finally have harmonic numbers in Sage.

arminstraub · 2016-08-17T22:33:37Z

Changed commit from d60311e to 9042104

arminstraub · 2016-08-17T22:33:37Z

Changed dependencies from #17790 to none

rwst · 2016-08-18T06:49:06Z

Reviewer: Ralf Stephan

rwst · 2016-08-18T06:49:06Z

comment:72

I agree with your changes, especially the removal of the dependency on #17790, this specialty is stuff for another ticket. I'm assuming you reviewed my code, so I set positive now, after more extensive tests pass here.

Please add your name in the Author and Reviewer fields. Thanks.

arminstraub · 2016-08-18T10:35:16Z

Changed author from Ralf Stephan to Ralf Stephan, Armin Straub

arminstraub · 2016-08-18T10:35:16Z

comment:73

Thanks for checking!

Yes, I did also review your code. I don't know enough about maxima and the intricacies of interfacing it, but it is my understanding that that part of the code was checked by you and Nils, and I didn't run into any additional trouble.

arminstraub · 2016-08-18T10:35:16Z

Changed reviewer from Ralf Stephan to Ralf Stephan, Armin Straub

arminstraub · 2016-08-19T21:25:11Z

comment:74

The patchbot complains about _swap_harmonic not having a doctest. Do we need to add one despite the underscore in the name?

Also, the bot indicates an issue with startup_modules, but I couldn't make sense of it...

rwst · 2016-08-20T06:18:50Z

comment:75

I would ignore both but let's wait for the release manager.

vbraun · 2016-08-21T13:14:00Z

Changed branch from u/arminstraub/ticket/16671/harmonic_numbers to 9042104

rwst added this to the sage-6.3 milestone Jul 17, 2014

rwst added c: symbolics labels Jul 17, 2014

This comment has been minimized.

Sign in to view

rwst changed the title ~~implement harmonic number function H_n~~ implement harmonic number function H(n,m) Jul 26, 2014

This comment has been minimized.

Sign in to view

rwst added t: bug and removed t: enhancement labels Jul 26, 2014

rwst added the s: needs info label Jul 27, 2014

rwst added s: needs review and removed s: needs info labels Jul 28, 2014

rwst added the pending label Jun 11, 2015

dkrenn added s: needs work and removed s: needs review labels Nov 29, 2015

rwst added s: needs review and removed s: needs work labels Nov 29, 2015

rwst added s: needs work and removed s: needs review labels Jun 14, 2016

arminstraub added s: needs review and removed s: needs work labels Aug 17, 2016

arminstraub added this to the sage-7.4 milestone Aug 17, 2016

arminstraub removed the pending label Aug 17, 2016

rwst added s: positive review and removed s: needs review labels Aug 18, 2016

vbraun removed the s: positive review label Aug 21, 2016

vbraun closed this as completed in 699518a Aug 21, 2016

rwst mentioned this issue Oct 27, 2014

replace gap.eval with libgap calls in combinat/combinat.py #16719

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

implement harmonic number function H(n,m) #16671

implement harmonic number function H(n,m) #16671

rwst commented Jul 17, 2014

fredrik-johansson commented Jul 17, 2014

rwst commented Jul 26, 2014

This comment has been minimized.

This comment has been minimized.

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

nbruin commented Jul 27, 2014

sagetrac-git mannequin commented Jul 28, 2014

sagetrac-git mannequin commented Jul 28, 2014

rwst commented Jul 28, 2014

rwst commented Jul 28, 2014

nbruin commented Jul 28, 2014

nbruin commented Jul 28, 2014

rwst commented Jul 28, 2014

sagetrac-git mannequin commented Jul 29, 2014

sagetrac-git mannequin commented Jul 29, 2014

rwst commented Jul 29, 2014

nbruin commented Jul 29, 2014

dkrenn commented Nov 29, 2015

rwst commented Nov 29, 2015

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

rwst commented Aug 18, 2016

rwst commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 19, 2016

rwst commented Aug 20, 2016

vbraun commented Aug 21, 2016

implement harmonic number function H(n,m) #16671

implement harmonic number function H(n,m) #16671

Comments

rwst commented Jul 17, 2014

fredrik-johansson commented Jul 17, 2014

rwst commented Jul 26, 2014

This comment has been minimized.

This comment has been minimized.

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

rwst commented Jul 27, 2014

nbruin commented Jul 27, 2014

sagetrac-git mannequin commented Jul 28, 2014

sagetrac-git mannequin commented Jul 28, 2014

rwst commented Jul 28, 2014

rwst commented Jul 28, 2014

nbruin commented Jul 28, 2014

nbruin commented Jul 28, 2014

rwst commented Jul 28, 2014

sagetrac-git mannequin commented Jul 29, 2014

sagetrac-git mannequin commented Jul 29, 2014

rwst commented Jul 29, 2014

nbruin commented Jul 29, 2014

dkrenn commented Nov 29, 2015

rwst commented Nov 29, 2015

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

arminstraub commented Aug 17, 2016

rwst commented Aug 18, 2016

rwst commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 18, 2016

arminstraub commented Aug 19, 2016

rwst commented Aug 20, 2016

vbraun commented Aug 21, 2016