Repository: sympy/sympy

Issue:
pi_hex_digits doesn't print leading 0 and errors
```python
>>> from sympy.ntheory.bbp_pi import pi_hex_digits
>>> pi_hex_digits(13)
'8d313198a2e037'  <-- should be '08d313198a2e03'
```

In addition the following starting pts lead to errors in the last of the 14 digits (and only the differeing digits are shown). There are more errors, but these are the first instances where a progressively larger trailing error is found:

```
pi_hex_digts(14) -> ...6 should be ...7
pi_hex_digts(381) -> ...af should be ...b0
pi_hex_digts(722) -> ...1ff should be ...200
```

Discussion:
smichr: Confirmed with the reference added:

```python
>>> from time import time
>>> t=time();pi_hex_digits(10**2-10 + 1, 10), time()-t
('e90c6cc0ac', 0.0)
>>> t=time();pi_hex_digits(10**4-10 + 1, 10), time()-t
('26aab49ec6', 0.17100000381469727)
>>> t=time();pi_hex_digits(10**5-10 + 1, 10), time()-t
('a22673c1a5', 4.7109999656677246)
>>> t=time();pi_hex_digits(10**6-10 + 1, 10), time()-t
('9ffd342362', 59.985999822616577)
>>> t=time();pi_hex_digits(10**7-10 + 1, 10), time()-t
('c1a42e06a1', 689.51800012588501)
```

Patch:
diff --git a/sympy/ntheory/bbp_pi.py b/sympy/ntheory/bbp_pi.py
index 993360675a..38832e0dbb 100644
--- a/sympy/ntheory/bbp_pi.py
+++ b/sympy/ntheory/bbp_pi.py
@@ -23,59 +23,99 @@
 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 Modifications:
-1.Once the nth digit is selected the number of digits of working
-precision is calculated to ensure that the 14 Hexadecimal representation
-of that region is accurate. This was found empirically to be
-int((math.log10(n//1000))+18). This was found by searching for a value
-of working precision for the n = 0 and n = 1 then n was increased until
-the result was less precise, therefore increased again this was repeated
-for increasing n and an effective fit was found between n and the
-working precision value.
-
-2. The while loop to evaluate whether the series has converged has be
-replaced with a fixed for loop, that option was selected because in a
-very large number of cases the loop converged to a point where no
-difference can be detected in less than 15 iterations. (done for more
-accurate memory and time banking).
-
-3. output hex string constrained to 14 characters (accuracy assured to be
-n = 10**7)
+
+1.Once the nth digit and desired number of digits is selected, the
+number of digits of working precision is calculated to ensure that
+the hexadecimal digits returned are accurate. This is calculated as
+
+    int(math.log(start + prec)/math.log(16) + prec + 3)
+    ---------------------------------------   --------
+                      /                          /
+    number of hex digits           additional digits
+
+This was checked by the following code which completed without
+errors (and dig are the digits included in the test_bbp.py file):
+
+    for i in range(0,1000):
+     for j in range(1,1000):
+      a, b = pi_hex_digits(i, j), dig[i:i+j]
+      if a != b:
+        print('%s\n%s'%(a,b))
+
+Deceasing the additional digits by 1 generated errors, so '3' is
+the smallest additional precision needed to calculate the above
+loop without errors. The following trailing 10 digits were also
+checked to be accurate (and the times were slightly faster with
+some of the constant modifications that were made):
+
+    >> from time import time
+    >> t=time();pi_hex_digits(10**2-10 + 1, 10), time()-t
+    ('e90c6cc0ac', 0.0)
+    >> t=time();pi_hex_digits(10**4-10 + 1, 10), time()-t
+    ('26aab49ec6', 0.17100000381469727)
+    >> t=time();pi_hex_digits(10**5-10 + 1, 10), time()-t
+    ('a22673c1a5', 4.7109999656677246)
+    >> t=time();pi_hex_digits(10**6-10 + 1, 10), time()-t
+    ('9ffd342362', 59.985999822616577)
+    >> t=time();pi_hex_digits(10**7-10 + 1, 10), time()-t
+    ('c1a42e06a1', 689.51800012588501)
+
+2. The while loop to evaluate whether the series has converged quits
+when the addition amount `dt` has dropped to zero.
+
+3. the formatting string to convert the decimal to hexidecimal is
+calculated for the given precision.
 
 4. pi_hex_digits(n) changed to have coefficient to the formula in an
 array (perhaps just a matter of preference).
 
 '''
+
 from __future__ import print_function, division
 
 import math
-from sympy.core.compatibility import range
+from sympy.core.compatibility import range, as_int
 
 
-def _series(j, n):
+def _series(j, n, prec=14):
 
     # Left sum from the bbp algorithm
     s = 0
-    D = _dn(n)
+    D = _dn(n, prec)
+    D4 = 4 * D
+    k = 0
+    d = 8 * k + j
     for k in range(n + 1):
-        r = 8*k + j
-        s += (pow(16, n - k, r) << 4 * D) // r
+        s += (pow(16, n - k, d) << D4) // d
+        d += 8
 
-    # Right sum. should iterate to infinty, but now just iterates to the point where
-    # one iterations change is beyond the resolution of the data type used
+    # Right sum iterates to infinity for full precision, but we
+    # stop at the point where one iteration is beyond the precision
+    # specified.
 
     t = 0
-    for k in range(n + 1, n + 15):
-        xp = int(16**(n - k) * 16**D)
-        t += xp // (8 * k + j)
+    k = n + 1
+    e = 4*(D + n - k)
+    d = 8 * k + j
+    while True:
+        dt = (1 << e) // d
+        if not dt:
+            break
+        t += dt
+        # k += 1
+        e -= 4
+        d += 8
     total = s + t
 
     return total
 
 
-def pi_hex_digits(n):
-    """Returns a string containing 14 digits after the nth value of pi in hex
-       The decimal has been taken out of the number, so
-       n = 0[0] = 3 # First digit of pi in hex, 3
+def pi_hex_digits(n, prec=14):
+    """Returns a string containing ``prec`` (default 14) digits
+    starting at the nth digit of pi in hex. Counting of digits
+    starts at 0 and the decimal is not counted, so for n = 0 the
+    returned value starts with 3; n = 1 corresponds to the first
+    digit past the decimal point (which in hex is 2).
 
     Examples
     ========
@@ -83,9 +123,19 @@ def pi_hex_digits(n):
     >>> from sympy.ntheory.bbp_pi import pi_hex_digits
     >>> pi_hex_digits(0)
     '3243f6a8885a30'
-    >>> pi_hex_digits(10)
-    '5a308d313198a2'
+    >>> pi_hex_digits(0, 3)
+    '324'
+
+    References
+    ==========
+
+    .. [1] http://www.numberworld.org/digits/Pi/
     """
+    n, prec = as_int(n), as_int(prec)
+    if n < 0:
+        raise ValueError('n cannot be negative')
+    if prec == 0:
+        return ''
 
     # main of implementation arrays holding formulae coefficients
     n -= 1
@@ -93,18 +143,19 @@ def pi_hex_digits(n):
     j = [1, 4, 5, 6]
 
     #formulae
-    x = + (a[0]*_series(j[0], n)
-         - a[1]*_series(j[1], n)
-         - a[2]*_series(j[2], n)
-         - a[3]*_series(j[3], n)) & (16**(_dn(n)) - 1)
+    D = _dn(n, prec)
+    x = + (a[0]*_series(j[0], n, prec)
+         - a[1]*_series(j[1], n, prec)
+         - a[2]*_series(j[2], n, prec)
+         - a[3]*_series(j[3], n, prec)) & (16**D - 1)
 
-    s = ("%014x" % x)
-    #s is constrained between 0 and 14
-    return s[:14]
+    s = ("%0" + "%ix" % prec) % (x // 16**(D - prec))
+    return s
 
 
-def _dn(n):
+def _dn(n, prec):
     # controller for n dependence on precision
-    if (n < 1000):
-        return 16
-    return int(math.log10(n//1000) + 18)
+    # n = starting digit index
+    # prec = the number of total digits to compute
+    n += 1  # because we subtract 1 for _series
+    return int(math.log(n + prec)/math.log(16) + prec + 3)


Test patch:
diff --git a/sympy/ntheory/tests/test_bbp_pi.py b/sympy/ntheory/tests/test_bbp_pi.py
index d063079ed4..2b0918e48a 100644
--- a/sympy/ntheory/tests/test_bbp_pi.py
+++ b/sympy/ntheory/tests/test_bbp_pi.py
@@ -1,7 +1,133 @@
+from random import randint
+
 from sympy.ntheory.bbp_pi import pi_hex_digits
+from sympy.utilities.pytest import raises
+
+
+# http://www.herongyang.com/Cryptography/Blowfish-First-8366-Hex-Digits-of-PI.html
+# There are actually 8336 listed there; with the preppended 3 there are 8337
+# below
+dig=''.join('''
+3243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89452821e638d013
+77be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b54709179216d5d98979fb1bd1310ba698dfb5
+ac2ffd72dbd01adfb7b8e1afed6a267e96ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc
+16636920d871574e69a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59
+b59c30d5392af26013c5d1b023286085f0ca417918b8db38ef8e79dcb0603a180e6c9e0e8bb01e8a
+3ed71577c1bd314b2778af2fda55605c60e65525f3aa55ab945748986263e8144055ca396a2aab10
+b6b4cc5c341141e8cea15486af7c72e993b3ee1411636fbc2a2ba9c55d741831f6ce5c3e169b8793
+1eafd6ba336c24cf5c7a325381289586773b8f48986b4bb9afc4bfe81b6628219361d809ccfb21a9
+91487cac605dec8032ef845d5de98575b1dc262302eb651b8823893e81d396acc50f6d6ff383f442
+392e0b4482a484200469c8f04a9e1f9b5e21c66842f6e96c9a670c9c61abd388f06a51a0d2d8542f
+68960fa728ab5133a36eef0b6c137a3be4ba3bf0507efb2a98a1f1651d39af017666ca593e82430e
+888cee8619456f9fb47d84a5c33b8b5ebee06f75d885c12073401a449f56c16aa64ed3aa62363f77
+061bfedf72429b023d37d0d724d00a1248db0fead349f1c09b075372c980991b7b25d479d8f6e8de
+f7e3fe501ab6794c3b976ce0bd04c006bac1a94fb6409f60c45e5c9ec2196a246368fb6faf3e6c53
+b51339b2eb3b52ec6f6dfc511f9b30952ccc814544af5ebd09bee3d004de334afd660f2807192e4b
+b3c0cba85745c8740fd20b5f39b9d3fbdb5579c0bd1a60320ad6a100c6402c7279679f25fefb1fa3
+cc8ea5e9f8db3222f83c7516dffd616b152f501ec8ad0552ab323db5fafd23876053317b483e00df
+829e5c57bbca6f8ca01a87562edf1769dbd542a8f6287effc3ac6732c68c4f5573695b27b0bbca58
+c8e1ffa35db8f011a010fa3d98fd2183b84afcb56c2dd1d35b9a53e479b6f84565d28e49bc4bfb97
+90e1ddf2daa4cb7e3362fb1341cee4c6e8ef20cada36774c01d07e9efe2bf11fb495dbda4dae9091
+98eaad8e716b93d5a0d08ed1d0afc725e08e3c5b2f8e7594b78ff6e2fbf2122b648888b812900df0
+1c4fad5ea0688fc31cd1cff191b3a8c1ad2f2f2218be0e1777ea752dfe8b021fa1e5a0cc0fb56f74
+e818acf3d6ce89e299b4a84fe0fd13e0b77cc43b81d2ada8d9165fa2668095770593cc7314211a14
+77e6ad206577b5fa86c75442f5fb9d35cfebcdaf0c7b3e89a0d6411bd3ae1e7e4900250e2d2071b3
+5e226800bb57b8e0af2464369bf009b91e5563911d59dfa6aa78c14389d95a537f207d5ba202e5b9
+c5832603766295cfa911c819684e734a41b3472dca7b14a94a1b5100529a532915d60f573fbc9bc6
+e42b60a47681e6740008ba6fb5571be91ff296ec6b2a0dd915b6636521e7b9f9b6ff34052ec58556
+6453b02d5da99f8fa108ba47996e85076a4b7a70e9b5b32944db75092ec4192623ad6ea6b049a7df
+7d9cee60b88fedb266ecaa8c71699a17ff5664526cc2b19ee1193602a575094c29a0591340e4183a
+3e3f54989a5b429d656b8fe4d699f73fd6a1d29c07efe830f54d2d38e6f0255dc14cdd20868470eb
+266382e9c6021ecc5e09686b3f3ebaefc93c9718146b6a70a1687f358452a0e286b79c5305aa5007
+373e07841c7fdeae5c8e7d44ec5716f2b8b03ada37f0500c0df01c1f040200b3ffae0cf51a3cb574
+b225837a58dc0921bdd19113f97ca92ff69432477322f547013ae5e58137c2dadcc8b576349af3dd
+a7a94461460fd0030eecc8c73ea4751e41e238cd993bea0e2f3280bba1183eb3314e548b384f6db9
+086f420d03f60a04bf2cb8129024977c795679b072bcaf89afde9a771fd9930810b38bae12dccf3f
+2e5512721f2e6b7124501adde69f84cd877a5847187408da17bc9f9abce94b7d8cec7aec3adb851d
+fa63094366c464c3d2ef1c18473215d908dd433b3724c2ba1612a14d432a65c45150940002133ae4
+dd71dff89e10314e5581ac77d65f11199b043556f1d7a3c76b3c11183b5924a509f28fe6ed97f1fb
+fa9ebabf2c1e153c6e86e34570eae96fb1860e5e0a5a3e2ab3771fe71c4e3d06fa2965dcb999e71d
+0f803e89d65266c8252e4cc9789c10b36ac6150eba94e2ea78a5fc3c531e0a2df4f2f74ea7361d2b
+3d1939260f19c279605223a708f71312b6ebadfe6eeac31f66e3bc4595a67bc883b17f37d1018cff
+28c332ddefbe6c5aa56558218568ab9802eecea50fdb2f953b2aef7dad5b6e2f841521b628290761
+70ecdd4775619f151013cca830eb61bd960334fe1eaa0363cfb5735c904c70a239d59e9e0bcbaade
+14eecc86bc60622ca79cab5cabb2f3846e648b1eaf19bdf0caa02369b9655abb5040685a323c2ab4
+b3319ee9d5c021b8f79b540b19875fa09995f7997e623d7da8f837889a97e32d7711ed935f166812
+810e358829c7e61fd696dedfa17858ba9957f584a51b2272639b83c3ff1ac24696cdb30aeb532e30
+548fd948e46dbc312858ebf2ef34c6ffeafe28ed61ee7c3c735d4a14d9e864b7e342105d14203e13
+e045eee2b6a3aaabeadb6c4f15facb4fd0c742f442ef6abbb5654f3b1d41cd2105d81e799e86854d
+c7e44b476a3d816250cf62a1f25b8d2646fc8883a0c1c7b6a37f1524c369cb749247848a0b5692b2
+85095bbf00ad19489d1462b17423820e0058428d2a0c55f5ea1dadf43e233f70613372f0928d937e
+41d65fecf16c223bdb7cde3759cbee74604085f2a7ce77326ea607808419f8509ee8efd85561d997
+35a969a7aac50c06c25a04abfc800bcadc9e447a2ec3453484fdd567050e1e9ec9db73dbd3105588
+cd675fda79e3674340c5c43465713e38d83d28f89ef16dff20153e21e78fb03d4ae6e39f2bdb83ad
+f7e93d5a68948140f7f64c261c94692934411520f77602d4f7bcf46b2ed4a20068d40824713320f4
+6a43b7d4b7500061af1e39f62e9724454614214f74bf8b88404d95fc1d96b591af70f4ddd366a02f
+45bfbc09ec03bd97857fac6dd031cb850496eb27b355fd3941da2547e6abca0a9a28507825530429
+f40a2c86dae9b66dfb68dc1462d7486900680ec0a427a18dee4f3ffea2e887ad8cb58ce0067af4d6
+b6aace1e7cd3375fecce78a399406b2a4220fe9e35d9f385b9ee39d7ab3b124e8b1dc9faf74b6d18
+5626a36631eae397b23a6efa74dd5b43326841e7f7ca7820fbfb0af54ed8feb397454056acba4895
+2755533a3a20838d87fe6ba9b7d096954b55a867bca1159a58cca9296399e1db33a62a4a563f3125
+f95ef47e1c9029317cfdf8e80204272f7080bb155c05282ce395c11548e4c66d2248c1133fc70f86
+dc07f9c9ee41041f0f404779a45d886e17325f51ebd59bc0d1f2bcc18f41113564257b7834602a9c
+60dff8e8a31f636c1b0e12b4c202e1329eaf664fd1cad181156b2395e0333e92e13b240b62eebeb9
+2285b2a20ee6ba0d99de720c8c2da2f728d012784595b794fd647d0862e7ccf5f05449a36f877d48
+fac39dfd27f33e8d1e0a476341992eff743a6f6eabf4f8fd37a812dc60a1ebddf8991be14cdb6e6b
+0dc67b55106d672c372765d43bdcd0e804f1290dc7cc00ffa3b5390f92690fed0b667b9ffbcedb7d
+9ca091cf0bd9155ea3bb132f88515bad247b9479bf763bd6eb37392eb3cc1159798026e297f42e31
+2d6842ada7c66a2b3b12754ccc782ef11c6a124237b79251e706a1bbe64bfb63501a6b101811caed
+fa3d25bdd8e2e1c3c9444216590a121386d90cec6ed5abea2a64af674eda86a85fbebfe98864e4c3
+fe9dbc8057f0f7c08660787bf86003604dd1fd8346f6381fb07745ae04d736fccc83426b33f01eab
+71b08041873c005e5f77a057bebde8ae2455464299bf582e614e58f48ff2ddfda2f474ef388789bd
+c25366f9c3c8b38e74b475f25546fcd9b97aeb26618b1ddf84846a0e79915f95e2466e598e20b457
+708cd55591c902de4cb90bace1bb8205d011a862487574a99eb77f19b6e0a9dc09662d09a1c43246
+33e85a1f0209f0be8c4a99a0251d6efe101ab93d1d0ba5a4dfa186f20f2868f169dcb7da83573906
+fea1e2ce9b4fcd7f5250115e01a70683faa002b5c40de6d0279af88c27773f8641c3604c0661a806
+b5f0177a28c0f586e0006058aa30dc7d6211e69ed72338ea6353c2dd94c2c21634bbcbee5690bcb6
+deebfc7da1ce591d766f05e4094b7c018839720a3d7c927c2486e3725f724d9db91ac15bb4d39eb8
+fced54557808fca5b5d83d7cd34dad0fc41e50ef5eb161e6f8a28514d96c51133c6fd5c7e756e14e
+c4362abfceddc6c837d79a323492638212670efa8e406000e03a39ce37d3faf5cfabc277375ac52d
+1b5cb0679e4fa33742d382274099bc9bbed5118e9dbf0f7315d62d1c7ec700c47bb78c1b6b21a190
+45b26eb1be6a366eb45748ab2fbc946e79c6a376d26549c2c8530ff8ee468dde7dd5730a1d4cd04d
+c62939bbdba9ba4650ac9526e8be5ee304a1fad5f06a2d519a63ef8ce29a86ee22c089c2b843242e
+f6a51e03aa9cf2d0a483c061ba9be96a4d8fe51550ba645bd62826a2f9a73a3ae14ba99586ef5562
+e9c72fefd3f752f7da3f046f6977fa0a5980e4a91587b086019b09e6ad3b3ee593e990fd5a9e34d7
+972cf0b7d9022b8b5196d5ac3a017da67dd1cf3ed67c7d2d281f9f25cfadf2b89b5ad6b4725a88f5
+4ce029ac71e019a5e647b0acfded93fa9be8d3c48d283b57ccf8d5662979132e28785f0191ed7560
+55f7960e44e3d35e8c15056dd488f46dba03a161250564f0bdc3eb9e153c9057a297271aeca93a07
+2a1b3f6d9b1e6321f5f59c66fb26dcf3197533d928b155fdf5035634828aba3cbb28517711c20ad9
+f8abcc5167ccad925f4de817513830dc8e379d58629320f991ea7a90c2fb3e7bce5121ce64774fbe
+32a8b6e37ec3293d4648de53696413e680a2ae0810dd6db22469852dfd09072166b39a460a6445c0
+dd586cdecf1c20c8ae5bbef7dd1b588d40ccd2017f6bb4e3bbdda26a7e3a59ff453e350a44bcb4cd
+d572eacea8fa6484bb8d6612aebf3c6f47d29be463542f5d9eaec2771bf64e6370740e0d8de75b13
+57f8721671af537d5d4040cb084eb4e2cc34d2466a0115af84e1b0042895983a1d06b89fb4ce6ea0
+486f3f3b823520ab82011a1d4b277227f8611560b1e7933fdcbb3a792b344525bda08839e151ce79
+4b2f32c9b7a01fbac9e01cc87ebcc7d1f6cf0111c3a1e8aac71a908749d44fbd9ad0dadecbd50ada
+380339c32ac69136678df9317ce0b12b4ff79e59b743f5bb3af2d519ff27d9459cbf97222c15e6fc
+2a0f91fc719b941525fae59361ceb69cebc2a8645912baa8d1b6c1075ee3056a0c10d25065cb03a4
+42e0ec6e0e1698db3b4c98a0be3278e9649f1f9532e0d392dfd3a0342b8971f21e1b0a74414ba334
+8cc5be7120c37632d8df359f8d9b992f2ee60b6f470fe3f11de54cda541edad891ce6279cfcd3e7e
+6f1618b166fd2c1d05848fd2c5f6fb2299f523f357a632762393a8353156cccd02acf081625a75eb
+b56e16369788d273ccde96629281b949d04c50901b71c65614e6c6c7bd327a140a45e1d006c3f27b
+9ac9aa53fd62a80f00bb25bfe235bdd2f671126905b2040222b6cbcf7ccd769c2b53113ec01640e3
+d338abbd602547adf0ba38209cf746ce7677afa1c52075606085cbfe4e8ae88dd87aaaf9b04cf9aa
+7e1948c25c02fb8a8c01c36ae4d6ebe1f990d4f869a65cdea03f09252dc208e69fb74e6132ce77e2
+5b578fdfe33ac372e6'''.split())
 
 
 def test_hex_pi_nth_digits():
     assert pi_hex_digits(0) == '3243f6a8885a30'
-    assert pi_hex_digits(1) == '243f6a8885a308'
+    assert pi_hex_digits(1) ==  '243f6a8885a308'
     assert pi_hex_digits(10000) == '68ac8fcfb8016c'
+    assert pi_hex_digits(13) == '08d313198a2e03'
+    assert pi_hex_digits(0, 3) == '324'
+    assert pi_hex_digits(0, 0) == ''
+    raises(ValueError, lambda: pi_hex_digits(-1))
+    raises(ValueError, lambda: pi_hex_digits(3.14))
+
+    # this will pick a random segment to compute every time
+    # it is run. If it ever fails, there is an error in the
+    # computation.
+    n = randint(0, len(dig))
+    prec = randint(0, len(dig) - n)
+    assert pi_hex_digits(n, prec) == dig[n: n + prec]