Opened 5 years ago

Closed 3 years ago

#16671 closed defect (fixed)

implement harmonic number function H(n,m)

Reported by: rws Owned by:
Priority: major Milestone: sage-7.4
Component: symbolics Keywords: special, log
Cc: eviatarbach, dkrenn Merged in:
Authors: Ralf Stephan, Armin Straub Reviewers: Ralf Stephan, Armin Straub
Report Upstream: N/A Work issues:
Branch: 9042104 (Commits) Commit: 9042104196f2d3b9822b3827aa2bbd70097f7b95
Dependencies: Stopgaps:

Description (last modified by rws)

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.

Change History (76)

comment:1 follow-up: Changed 5 years ago by fredrik.johansson

For exact values, FLINT has arith_harmonic_number()

comment:2 in reply to: ↑ 1 Changed 5 years ago by rws

  • Description modified (diff)
  • Summary changed from implement harmonic number function H_n to implement harmonic number function H(n,m)

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.

comment:3 Changed 5 years ago by rws

  • Description modified (diff)
  • Type changed from enhancement to defect

comment:4 Changed 5 years ago by rws

  • Branch set to u/rws/implement_harmonic_number_function_h_n_m_

comment:5 follow-up: Changed 5 years ago by rws

  • Authors set to Ralf Stephan
  • Commit set to 5310ca10a440c2536da3f043deba80f5c35d041a
  • Status changed from new to needs_info

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:

5310ca116671: implement harmonic number function H(n,m) (modulo argument order)

comment:6 in reply to: ↑ 5 Changed 5 years ago by nbruin

Replying to rws:

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:

5310ca116671: 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.

comment:7 Changed 5 years ago by git

  • Commit changed from 5310ca10a440c2536da3f043deba80f5c35d041a to ebc93ecacf652b5b83077ba20c22b37f8178ca78

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

61572c916671: swap arguments wrt maxima(_lib); add more documentation
ebc93ec16671: fix doc typo

comment:8 Changed 5 years ago by rws

  • Keywords special log added
  • Status changed from needs_info to needs_review

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

Please review.

comment:9 follow-up: Changed 5 years ago by nbruin

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.

comment:10 Changed 5 years ago by nbruin

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.

comment:11 Changed 5 years ago by rws

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.

comment:12 Changed 5 years ago by git

  • Commit changed from ebc93ecacf652b5b83077ba20c22b37f8178ca78 to 0c93bfe4da1bd2b3a77865ae0a6878007726fff8

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

0b15dd1Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_
0c93bfe16671: fix maxima_lib interface

comment:13 in reply to: ↑ 9 ; follow-up: Changed 5 years ago by rws

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.

comment:14 in reply to: ↑ 13 Changed 5 years ago by nbruin

Replying to rws:

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.

comment:15 Changed 5 years ago by git

  • Commit changed from 0c93bfe4da1bd2b3a77865ae0a6878007726fff8 to fd0ccd6ea2e7f6a0abf92ee85207c011250b9496

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

e90d77b16671: replace maxima interface arg swap with pynac conversion
fd0ccd616671: improve previous patch; handle maple swap too

comment:16 Changed 5 years ago by rws

Thanks again. I completely missed this trick, although I already created a dummy function there.

comment:17 Changed 5 years ago by nbruin

Sorry, this isn't really your fault, but as it turns out the whole _maxima_init_evaled_ mechanism gets used in a way that's broken: parameters do not get processed recursively. See #16732. As an illustration:

sage:sage: maxima_calculus(polylog(3,polylog(3,x)))
li[3](polylog(3,x))

Note that the second argument doesn't get appropriately processed. Your implementation suffers from the same problem.

Either we change the way that _maxima_init_evaled_ gets called (i.e., we pass it arguments that are processed by the maxima translator itself) or we must change all the _maxima_init_evaled_ implementations to do that properly (at least a bunch presently don't). An argument against giving _maxima_init_evaled_ preprocessed arguments is that in principle it's conceivable that a function has an argument that needs special treatment. So letting _maxima_init_evaled_ do the recursive argument processing under its own control is the more flexible (but more onerous) scheme.

comment:18 Changed 5 years ago by git

  • Commit changed from fd0ccd6ea2e7f6a0abf92ee85207c011250b9496 to 4fa8dd871063de9d658a11135861d3564d215b4c

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

019257d16671: _maxima_init_evaled_ fix and doctest
4fa8dd816671: fix call with m an expression: allow complex m

comment:19 Changed 5 years ago by rws

Done. I'm not sure if Maxima does the right thing here by expanding the sum. This can get messy for big integers:

sage: maxima_calculus(harmonic_number(3,harmonic_number(x,3),hold=True)).sage()
1/3^harmonic_number(x, 3) + 1/2^harmonic_number(x, 3) + 1

I also removed the restriction on the second argument.

comment:20 Changed 5 years ago by git

  • Commit changed from 4fa8dd871063de9d658a11135861d3564d215b4c to e263ac41caa14a9ba4b537f0f5091cf012d64c3b

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

e263ac4Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:21 follow-up: Changed 5 years ago by nbruin

Ouch, given what we've just run into on #16697 with derivatives, our _maxima_init_evaled_ hack for swapping arguments would fail in the presence of derivatives:

maxima_calculus(harmonic_number(x,y).diff(x))

Should we disallow differentiating harmonic numbers? Otherwise I think we need to rewrite how differential operators get translated. Currently, given:

D[0,1,2](f)(u,v,w)

gets translated basically to

"diff("+"f"+"(u,v,w)"+",u,1,v,1,w,1)"

Instead, we could make the *expression* f(u,v,w) and convert that to a maxima string *as a whole* and compose that with the rest. Then we WOULD be going through _maxima_init_evaled_ etc. There are situations where this is a little more complicated: when "u,v,w" are not pure, distinct, variables, then sage needs to substitute temporary place holder variables in order to express the D[0,1,2] operator in Leibnitz notation. (and then, after differentiation, substitute the values). Example:

sage: function('f')
sage: diff( f(x+1),x)._maxima_init_()
"at(diff(''f(t0), t0, 1), [t0 = x + 1])"
sage: maxima_calculus(diff( f(x+1),x))
?%at('diff(f(t0),t0,1),[t0=x+1])

Presently, those temporary variables t0 do not live in Sage at all. With this change, we would have to instantiate them as symbolic variables in sage to get them converted. (I have checked before, and the at substitutions in maxima happen simultaneously, so symbol clashes shouldn't be an issue, other than the possibility of assumptions interfering)

comment:22 Changed 5 years ago by nbruin

Ouch, this actually shows a missed conversion:

sage: maxima_calculus(f(x).diff(x))
'diff(f(_SAGE_VAR_x),_SAGE_VAR_x,1)
sage: maxima_calculus(f(x+1).diff(x)) #this is bad!
?%at('diff(f(t0),t0,1),[t0=x+1])

and to compare:

sage: from sage.interfaces.maxima_lib import sr_to_max
sage: maxima_calculus(sr_to_max(f(x+1).diff(x)))
?%at('diff(f(_SAGE_VAR_t0),_SAGE_VAR_t0,1),[_SAGE_VAR_t0=_SAGE_VAR_x+1])
sage: maxima_calculus(sr_to_max(f(x).diff(x)))
'diff(f(_SAGE_VAR_x),_SAGE_VAR_x,1)

How the temporary variable t0 gets represented doesn't matter, but x should absolutely be converted to _SAGE_VAR_x, in all cases.

The problem is in line 505. Changing it like this:

-            subs = ["%s = %s"%(t,a) for t,a in zip(temp_args,args)]
+            subs = ["%s = %s"%(t._maxima_init_(),a._maxima_init_()) for t,a in zip(temp_args,args)]

and some more changes of the same nature should do the trick.

This code also shows that the temporary variables are instantiated in sage anyway, so changing the conversion wouldn't be so bad.

comment:23 Changed 5 years ago by nbruin

See #16785 for differential operator translation fix.

comment:24 in reply to: ↑ 21 ; follow-up: Changed 5 years ago by rws

Replying to nbruin:

Ouch, given what we've just run into on #16697 with derivatives, our _maxima_init_evaled_ hack for swapping arguments would fail in the presence of derivatives:

maxima_calculus(harmonic_number(x,y).diff(x))

Should we disallow differentiating harmonic numbers?

I think so, there's really no need for it.

comment:25 in reply to: ↑ 24 Changed 5 years ago by nbruin

Replying to rws:

Should we disallow differentiating harmonic numbers?

I think so, there's really no need for it.

With #16785 (which needs to be merged anyway) the urgency is reduced. Given the complex continuations that exist for the harmonic number function, perhaps we shouldn't go out of our way prohibiting it, even if we don't have any functionality for it (yet).

comment:26 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:27 Changed 5 years ago by git

  • Commit changed from e263ac41caa14a9ba4b537f0f5091cf012d64c3b to 860820fd5d454bfd4e2e8f40af82e724537f45dc

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

9c0f66dtrac #16785: improve differential operator translation to maxima
ef6d620Merge commit '9c0f66dca5d59aef571c69fe7561b6446ec23d71' of trac.sagemath.org:sage into t/16671/implement_harmonic_number_function_h_n_m_
860820fMerge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:28 Changed 5 years ago by rws

This looks good with #16785:

sage: maxima_calculus(harmonic_number(x,y).diff(x))
'diff(gen_harmonic_number(_SAGE_VAR_y,_SAGE_VAR_x),_SAGE_VAR_x,1)

So, thanks for your review so far. Is there anything missing?

comment:29 Changed 5 years ago by nbruin

Patchbot complains about a failing doctest

File "src/sage/interfaces/maxima_lib.py", line 1458, in sage.interfaces.maxima_lib.max_harmonic_to_sage
Failed example:
    c.ecl()
Expected:
    <ECL: (($GEN_HARMONIC_NUMBER SIMP) 2 $X)>
Got:
    <ECL: (($GEN_HARMONIC_NUMBER SIMP) 2 |$_SAGE_VAR_x|)>

comment:30 Changed 5 years ago by git

  • Commit changed from 860820fd5d454bfd4e2e8f40af82e724537f45dc to 26652819bf557f1daed5e01d892d7ad982e09d29

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

266528116671: doctest fix (appeared due to 16785 merge)

comment:31 Changed 5 years ago by rws

Ah thanks, the $X should have rung a bell much earlier.

comment:32 Changed 5 years ago by git

  • Commit changed from 26652819bf557f1daed5e01d892d7ad982e09d29 to fd8c09c4be9334f4fd9af3225aa3a2898dd8366d

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

fd8c09cMerge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:33 Changed 5 years ago by git

  • Commit changed from fd8c09c4be9334f4fd9af3225aa3a2898dd8366d to 87597bde4c041766a26674a3a8be1d4b076180d5

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

fa73ddeMerge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_
87597bd16671: adaptations due to improved BuiltinFunction

comment:34 Changed 5 years ago by eviatarbach

  • Cc eviatarbach added

comment:35 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

You don't need this:

diff --git a/src/sage/libs/flint/arith.pyx b/src/sage/libs/flint/arith.pyx
index 7280536..91caf32 100644
--- a/src/sage/libs/flint/arith.pyx
+++ b/src/sage/libs/flint/arith.pyx
@@ -10,10 +10,23 @@ FLINT Arithmetic Functions
include "sage/ext/interrupt.pxi"
+cdef extern from "flint/fmpq.h":
+ ctypedef void * fmpq_t
+ void fmpq_init(fmpq_t)
+ void fmpq_clear(fmpq_t)
+ void fmpq_get_mpq(mpq_t, fmpq_t)
+ void fmpq_set_mpq(fmpq_t, mpq_t)
+
+cdef extern from "flint/arith.h":
+ void arith_number_of_partitions(fmpz_t x, unsigned long n)
+ void arith_dedekind_sum(fmpq_t, fmpz_t, fmpz_t)
+ void arith_harmonic_number(fmpq_t, unsigned long n)
+
from fmpz cimport *
from fmpq cimport *
from arith cimport *
+
from sage.rings.integer cimport Integer
from sage.rings.rational cimport Rational

Add the functions you need (only arith_harmonic_number in this case) to src/sage/libs/flint/arith.pxd instead.

comment:36 Changed 5 years ago by rws

  • Branch changed from u/rws/implement_harmonic_number_function_h_n_m_ to u/rws/16671

comment:37 Changed 5 years ago by rws

  • Commit changed from 87597bde4c041766a26674a3a8be1d4b076180d5 to c26cd0b047aef2446737587bd784e38851277ed4
  • Status changed from needs_work to needs_review

Done. Squashed it all into one commit.


New commits:

c26cd0b16671: implement harmonic number function H(n,m)

comment:38 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

In flint/arith.pyx, you're missing actual examples.

comment:39 Changed 5 years ago by git

  • Commit changed from c26cd0b047aef2446737587bd784e38851277ed4 to 3ec7fbba6b2c095ef2aa38cec5ffc47c9ba0cae6

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

3ec7fbb16671: add missing doctest

comment:40 Changed 5 years ago by rws

  • Status changed from needs_work to needs_review

Indeed.

comment:41 Changed 5 years ago by jdemeyer

In __call__, why mess with len(args)? Why not simply something like

def __call__(self, z, m=1, **kwds):
    return BuiltinFunction.__call__(self, z, m, **kwds)

comment:42 Changed 5 years ago by git

  • Commit changed from 3ec7fbba6b2c095ef2aa38cec5ffc47c9ba0cae6 to fabf86397adf8c07ce9b63b64c64e5f724702d21

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

fabf86316671: code tightened

comment:43 Changed 5 years ago by jdemeyer

What's with the condition elif z < 2**20:??

comment:44 Changed 5 years ago by jdemeyer

In flint/arith.pyx, you should add

+sig_on()
 arith_harmonic_number(s_fmpq, n)

 fmpq_get_mpq((<Rational>s).value, s_fmpq)
+sig_off()

comment:45 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:46 Changed 5 years ago by jdemeyer

Why doesn't this work?

sage: harmonic_number(x).diff(x)
D[0]harmonic_number(x)

comment:47 follow-ups: Changed 5 years ago by jdemeyer

I don't agree with sum(QQ(1)/(k**m) for k in range(1,z+1)) since the result is always in QQ. I would use srange to respect the input type (foo in ZZ does not imply that the type of foo is ZZ):

sum(1/(k**m) for k in srange(1,z+1))

Artificial example (this could be a doctest):

sage: harmonic_number._eval_(10, Qp(5)(1))
7381/2520

but the answer should be

sage: Qp(5)(harmonic_number(10))
4*5^-1 + 2 + 2*5 + 4*5^2 + 3*5^3 + 5^4 + 4*5^5 + 4*5^6 + 5^7 + 4*5^8 + 3*5^9 + 5^10 + 4*5^11 + 4*5^12 + 5^13 + 4*5^14 + 3*5^15 + 5^16 + 4*5^17 + 4*5^18 + O(5^19)

comment:48 Changed 5 years ago by jdemeyer

Also please fix \ No newline at end of file

comment:49 Changed 5 years ago by jdemeyer

I don't know what SR('x') does, but I think it's better to use SR.var('x') instead.

comment:50 follow-up: Changed 5 years ago by jdemeyer

Same comment as above: are you sure you want elif z in QQ: and not elif isinstance(z, Rational)?

comment:51 in reply to: ↑ 50 ; follow-up: Changed 5 years ago by rws

Replying to jdemeyer:

Same comment as above: are you sure you want elif z in QQ: and not elif isinstance(z, Rational)?

BTW I would think this if var in ZZ/QQ appears more often throughout the Sage code. Is there a manual page that deals with this?

comment:52 in reply to: ↑ 51 Changed 5 years ago by jdemeyer

Replying to rws:

BTW I would think this if var in ZZ/QQ appears more often throughout the Sage code.

Perhaps, but it's not always wrong. It just depends on what you want to do with it.

comment:53 Changed 5 years ago by jdemeyer

The general rule for functions should be: whenever the answer is non-symbolic, the answer should have a parent which "matches" the parent of the input.

comment:54 in reply to: ↑ 47 ; follow-up: Changed 5 years ago by rws

Replying to jdemeyer:

Artificial example (this could be a doctest):

sage: harmonic_number._eval_(10, Qp(5)(1))
7381/2520

but the answer should be

sage: Qp(5)(harmonic_number(10))
4*5^-1 + 2 + 2*5 + 4*5^2 + 3*5^3 + 5^4 + 4*5^5 + 4*5^6 + 5^7 + 4*5^8 + 3*5^9 + 5^10 + 4*5^11 + 4*5^12 + 5^13 + 4*5^14 + 3*5^15 + 5^16 + 4*5^17 + 4*5^18 + O(5^19)

I think this is the only example where H(n,m) would have a useful meaning, and it can be computed, as you show, by staying in the rationals and converting the result. Trying to support this case however will fail in BuiltinFunction.__call__:

sage: harmonic_number(5,Qp(5)(8))
Traceback (most recent call last):
  File "<ipython-input-9-03dbaa59749c>", line 1, in <module>
    harmonic_number(Integer(5),Qp(Integer(5))(Integer(8)))
  File "/home/ralf/sage/local/lib/python2.7/site-packages/sage/functions/log.py", line 886, in __call__
    return BuiltinFunction.__call__(self, z, m, **kwds)
  File "sage/symbolic/function.pyx", line 993, in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:10572)
    res = super(BuiltinFunction, self).__call__(
  File "sage/symbolic/function.pyx", line 487, in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6301)
    raise TypeError, "cannot coerce arguments: %s"%(err)
TypeError: cannot coerce arguments: no canonical coercion from 5-adic Field with capped relative precision 20 to Symbolic Ring

comment:55 in reply to: ↑ 54 ; follow-up: Changed 5 years ago by jdemeyer

Replying to rws:

it can be computed, as you show, by staying in the rationals and converting the result.

But that wouldn't be the optimal way to compute it. I think it's easy to make the formula work for arbitrary types.

Trying to support this case however will fail in BuiltinFunction.__call__:

I know, but that doesn't mean we shouldn't adapt _eval_.

comment:56 Changed 5 years ago by git

  • Commit changed from fabf86397adf8c07ce9b63b64c64e5f724702d21 to 309660a11c0929f223087f5f208a2ea8f14e7625

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

309660a16671: generalize evaluation of H(n,m); add its derivative; several fixes; document FLINT computation limit

comment:57 in reply to: ↑ 55 ; follow-up: Changed 5 years ago by rws

  • Status changed from needs_work to needs_review

Trying to support this case however will fail in BuiltinFunction.__call__:

I know, but that doesn't mean we shouldn't adapt _eval_.

That means it will be fixed elsewhere?

comment:58 follow-up: Changed 5 years ago by jdemeyer

I still don't understand the point of this:

If the argument is an integer greater than `2^20` no rational
evaluation is done, in order to allow for high-precision floating
point computation using `.n()`

comment:59 Changed 5 years ago by git

  • Commit changed from 309660a11c0929f223087f5f208a2ea8f14e7625 to ca696281c6d2425b72d72b6408e4a787b32d637a

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

ca6962816671: remove restriction

comment:60 in reply to: ↑ 58 Changed 5 years ago by rws

Replying to jdemeyer:

I still don't understand the point of this:

OK, me neither. Presumably I had a memory error at the time ...

comment:61 in reply to: ↑ 57 Changed 5 years ago by rws

  • Dependencies set to #17790

Replying to rws:

Trying to support this case however will fail in BuiltinFunction.__call__:

I know, but that doesn't mean we shouldn't adapt _eval_.

That means it will be fixed elsewhere?

See #17790

comment:62 Changed 5 years ago by rws

  • Status changed from needs_review to needs_work
  • Work issues set to merge fix for #17790

comment:63 Changed 5 years ago by rws

Let's document here that both forms for the generalized harmonic numbers

H(q/p,m) = zeta(m) - p^m * sum(1/(q+p*k)^m, k, 1, oo)   (Wikipedia 2015-Mar-07)

H(k,n) = binomial(n+k-1,k-1)*(harmonic_number(n+k-1)-harmonic_number(k-1))
         (http://mathworld.wolfram.com/HarmonicNumber.html eq. 55)

are both wrong, with the counterexamples H(5/2,3)=1/27*sqrt(3) + 1/8*sqrt(2) + 1 and H(5,3)=256103/216000.

EDIT. no, the first is right, my implementation was wrong.

Last edited 5 years ago by rws (previous) (diff)

comment:64 Changed 5 years ago by git

  • Commit changed from ca696281c6d2425b72d72b6408e4a787b32d637a to d60311eda7cce822b74480550498da09b34b7a63

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

5eb502dMerge branch 'develop' into t/16671/16671
a808a0eMerge branch 'develop' into t/16671/16671
f6e211317790: add padics to the conversion workaround for symbolic function arguments
05ad72bMerge branch 'u/rws/builtinfunction_doesn_t_pass_non_sr_coercible_arguments_to_function_code' of trac.sagemath.org:sage into t/16671/16671
ca1a03dMerge branch 'develop' into t/16671/16671
d60311e16671: handle more argument types; fix num. H(n,m); cosmetics

comment:65 in reply to: ↑ 47 Changed 5 years ago by rws

  • Milestone changed from sage-6.4 to sage-6.6
  • Status changed from needs_work to needs_review
  • Work issues merge fix for #17790 deleted

Replying to jdemeyer:

sage: harmonic_number._eval_(10, Qp(5)(1))
7381/2520

but the answer should be

sage: Qp(5)(harmonic_number(10))

Disagree. That should be the result of harmonic_number._eval_(Qp(5)(10), 1). Anything else than a rational as second argument makes no sense.

I have now adapted the code so it gives back the same type as the input, but wrapped as Expression. Also, you can expect that only if it can be converted back, which is not the case for H(p/q) or H(anything, p/q) as the result is not rational.

comment:66 Changed 5 years ago by rws

  • Milestone changed from sage-6.6 to sage-pending

comment:67 Changed 4 years ago by dkrenn

  • Cc dkrenn added

comment:68 Changed 4 years ago by dkrenn

  • Status changed from needs_review to needs_work
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

comment:69 Changed 4 years ago by rws

  • Status changed from needs_work to needs_review

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

comment:70 Changed 4 years ago by rws

  • Status changed from needs_review to needs_work

comment:71 Changed 3 years ago by arminstraub

  • Branch changed from u/rws/16671 to u/arminstraub/ticket/16671/harmonic_numbers
  • Commit changed from d60311eda7cce822b74480550498da09b34b7a63 to 9042104196f2d3b9822b3827aa2bbd70097f7b95
  • Dependencies #17790 deleted
  • Milestone changed from sage-pending to sage-7.4
  • Status changed from needs_work to needs_review

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.

comment:72 Changed 3 years ago by rws

  • Reviewers set to Ralf Stephan
  • Status changed from needs_review to positive_review

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.

comment:73 Changed 3 years ago by arminstraub

  • Authors changed from Ralf Stephan to Ralf Stephan, Armin Straub
  • Reviewers changed from Ralf Stephan to Ralf Stephan, Armin Straub

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.

comment:74 Changed 3 years ago by arminstraub

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...

comment:75 Changed 3 years ago by rws

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

comment:76 Changed 3 years ago by vbraun

  • Branch changed from u/arminstraub/ticket/16671/harmonic_numbers to 9042104196f2d3b9822b3827aa2bbd70097f7b95
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.