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:  sage7.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 )
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/sagedevel/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 followup: ↓ 2 Changed 5 years ago by
comment:2 in reply to: ↑ 1 Changed 5 years ago by
 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
 Description modified (diff)
 Type changed from enhancement to defect
comment:4 Changed 5 years ago by
 Branch set to u/rws/implement_harmonic_number_function_h_n_m_
comment:5 followup: ↓ 6 Changed 5 years ago by
 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/sagedevel/SM_FLL33t2g as to argument order of the generalized version.
New commits:
5310ca1  16671: implement harmonic number function H(n,m) (modulo argument order)

comment:6 in reply to: ↑ 5 Changed 5 years ago by
Replying to rws:
It is already implemented, I'm just waiting for https://groups.google.com/forum/?hl=en#!topic/sagedevel/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 argumentswapping 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
 Commit changed from 5310ca10a440c2536da3f043deba80f5c35d041a to ebc93ecacf652b5b83077ba20c22b37f8178ca78
comment:8 Changed 5 years ago by
 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 followup: ↓ 13 Changed 5 years ago by
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 languagethey 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
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
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
 Commit changed from ebc93ecacf652b5b83077ba20c22b37f8178ca78 to 0c93bfe4da1bd2b3a77865ae0a6878007726fff8
comment:13 in reply to: ↑ 9 ; followup: ↓ 14 Changed 5 years ago by
Replying to nbruin:
Regular languages are not sufficient to describe expression languages (the nested parentheses can't be expressed in a regular languagethey 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
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
 Commit changed from 0c93bfe4da1bd2b3a77865ae0a6878007726fff8 to fd0ccd6ea2e7f6a0abf92ee85207c011250b9496
comment:16 Changed 5 years ago by
Thanks again. I completely missed this trick, although I already created a dummy function there.
comment:17 Changed 5 years ago by
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
 Commit changed from fd0ccd6ea2e7f6a0abf92ee85207c011250b9496 to 4fa8dd871063de9d658a11135861d3564d215b4c
comment:19 Changed 5 years ago by
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
 Commit changed from 4fa8dd871063de9d658a11135861d3564d215b4c to e263ac41caa14a9ba4b537f0f5091cf012d64c3b
Branch pushed to git repo; I updated commit sha1. New commits:
e263ac4  Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:21 followup: ↓ 24 Changed 5 years ago by
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
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
See #16785 for differential operator translation fix.
comment:24 in reply to: ↑ 21 ; followup: ↓ 25 Changed 5 years ago by
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
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
 Milestone changed from sage6.3 to sage6.4
comment:27 Changed 5 years ago by
 Commit changed from e263ac41caa14a9ba4b537f0f5091cf012d64c3b to 860820fd5d454bfd4e2e8f40af82e724537f45dc
Branch pushed to git repo; I updated commit sha1. New commits:
9c0f66d  trac #16785: improve differential operator translation to maxima

ef6d620  Merge commit '9c0f66dca5d59aef571c69fe7561b6446ec23d71' of trac.sagemath.org:sage into t/16671/implement_harmonic_number_function_h_n_m_

860820f  Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:28 Changed 5 years ago by
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
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
 Commit changed from 860820fd5d454bfd4e2e8f40af82e724537f45dc to 26652819bf557f1daed5e01d892d7ad982e09d29
Branch pushed to git repo; I updated commit sha1. New commits:
2665281  16671: doctest fix (appeared due to 16785 merge)

comment:31 Changed 5 years ago by
Ah thanks, the $X
should have rung a bell much earlier.
comment:32 Changed 5 years ago by
 Commit changed from 26652819bf557f1daed5e01d892d7ad982e09d29 to fd8c09c4be9334f4fd9af3225aa3a2898dd8366d
Branch pushed to git repo; I updated commit sha1. New commits:
fd8c09c  Merge branch 'develop' into t/16671/implement_harmonic_number_function_h_n_m_

comment:33 Changed 5 years ago by
 Commit changed from fd8c09c4be9334f4fd9af3225aa3a2898dd8366d to 87597bde4c041766a26674a3a8be1d4b076180d5
comment:34 Changed 5 years ago by
 Cc eviatarbach added
comment:35 Changed 5 years ago by
 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
 Branch changed from u/rws/implement_harmonic_number_function_h_n_m_ to u/rws/16671
comment:37 Changed 5 years ago by
 Commit changed from 87597bde4c041766a26674a3a8be1d4b076180d5 to c26cd0b047aef2446737587bd784e38851277ed4
 Status changed from needs_work to needs_review
Done. Squashed it all into one commit.
New commits:
c26cd0b  16671: implement harmonic number function H(n,m)

comment:38 Changed 5 years ago by
 Status changed from needs_review to needs_work
In flint/arith.pyx
, you're missing actual examples.
comment:39 Changed 5 years ago by
 Commit changed from c26cd0b047aef2446737587bd784e38851277ed4 to 3ec7fbba6b2c095ef2aa38cec5ffc47c9ba0cae6
Branch pushed to git repo; I updated commit sha1. New commits:
3ec7fbb  16671: add missing doctest

comment:41 Changed 5 years ago by
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
 Commit changed from 3ec7fbba6b2c095ef2aa38cec5ffc47c9ba0cae6 to fabf86397adf8c07ce9b63b64c64e5f724702d21
Branch pushed to git repo; I updated commit sha1. New commits:
fabf863  16671: code tightened

comment:43 Changed 5 years ago by
What's with the condition elif z < 2**20:
??
comment:44 Changed 5 years ago by
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
 Status changed from needs_review to needs_work
comment:46 Changed 5 years ago by
Why doesn't this work?
sage: harmonic_number(x).diff(x) D[0]harmonic_number(x)
comment:47 followups: ↓ 54 ↓ 65 Changed 5 years ago by
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
Also please fix \ No newline at end of file
comment:49 Changed 5 years ago by
I don't know what SR('x')
does, but I think it's better to use SR.var('x')
instead.
comment:50 followup: ↓ 51 Changed 5 years ago by
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 ; followup: ↓ 52 Changed 5 years ago by
Replying to jdemeyer:
Same comment as above: are you sure you want
elif z in QQ:
and notelif 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
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
The general rule for functions should be: whenever the answer is nonsymbolic, the answer should have a parent which "matches" the parent of the input.
comment:54 in reply to: ↑ 47 ; followup: ↓ 55 Changed 5 years ago by
Replying to jdemeyer:
Artificial example (this could be a doctest):
sage: harmonic_number._eval_(10, Qp(5)(1)) 7381/2520but 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 "<ipythoninput903dbaa59749c>", line 1, in <module> harmonic_number(Integer(5),Qp(Integer(5))(Integer(8))) File "/home/ralf/sage/local/lib/python2.7/sitepackages/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 5adic Field with capped relative precision 20 to Symbolic Ring
comment:55 in reply to: ↑ 54 ; followup: ↓ 57 Changed 5 years ago by
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
 Commit changed from fabf86397adf8c07ce9b63b64c64e5f724702d21 to 309660a11c0929f223087f5f208a2ea8f14e7625
Branch pushed to git repo; I updated commit sha1. New commits:
309660a  16671: generalize evaluation of H(n,m); add its derivative; several fixes; document FLINT computation limit

comment:57 in reply to: ↑ 55 ; followup: ↓ 61 Changed 5 years ago by
 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 followup: ↓ 60 Changed 5 years ago by
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 highprecision floating point computation using `.n()`
comment:59 Changed 5 years ago by
 Commit changed from 309660a11c0929f223087f5f208a2ea8f14e7625 to ca696281c6d2425b72d72b6408e4a787b32d637a
Branch pushed to git repo; I updated commit sha1. New commits:
ca69628  16671: remove restriction

comment:60 in reply to: ↑ 58 Changed 5 years ago by
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
 Dependencies set to #17790
comment:62 Changed 4 years ago by
 Status changed from needs_review to needs_work
 Work issues set to merge fix for #17790
comment:63 Changed 4 years ago by
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 2015Mar07) H(k,n) = binomial(n+k1,k1)*(harmonic_number(n+k1)harmonic_number(k1)) (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.
comment:64 Changed 4 years ago by
 Commit changed from ca696281c6d2425b72d72b6408e4a787b32d637a to d60311eda7cce822b74480550498da09b34b7a63
Branch pushed to git repo; I updated commit sha1. New commits:
5eb502d  Merge branch 'develop' into t/16671/16671

a808a0e  Merge branch 'develop' into t/16671/16671

f6e2113  17790: add padics to the conversion workaround for symbolic function arguments

05ad72b  Merge branch 'u/rws/builtinfunction_doesn_t_pass_non_sr_coercible_arguments_to_function_code' of trac.sagemath.org:sage into t/16671/16671

ca1a03d  Merge branch 'develop' into t/16671/16671

d60311e  16671: handle more argument types; fix num. H(n,m); cosmetics

comment:65 in reply to: ↑ 47 Changed 4 years ago by
 Milestone changed from sage6.4 to sage6.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/2520but 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 4 years ago by
 Milestone changed from sage6.6 to sagepending
comment:67 Changed 4 years ago by
 Cc dkrenn added
comment:68 Changed 4 years ago by
 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/sitepackages/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/sitepackages/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
 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 3 years ago by
 Status changed from needs_review to needs_work
comment:71 Changed 3 years ago by
 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 sagepending to sage7.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 padics 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
 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
 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
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
I would ignore both but let's wait for the release manager.
comment:76 Changed 3 years ago by
 Branch changed from u/arminstraub/ticket/16671/harmonic_numbers to 9042104196f2d3b9822b3827aa2bbd70097f7b95
 Resolution set to fixed
 Status changed from positive_review to closed
For exact values, FLINT has arith_harmonic_number()