forked from mdickinson/bigfloat
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmpfr.pyx
2952 lines (2348 loc) · 97.6 KB
/
mpfr.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: utf-8 -*-
# cython: embedsignature = True
# Copyright 2009--2015 Mark Dickinson.
#
# This file is part of the bigfloat package.
#
# The bigfloat package is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# The bigfloat package is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
# for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with the bigfloat package. If not, see <http://www.gnu.org/licenses/>.
cimport cmpfr
import sys
cdef extern from "limits.h":
cdef long LONG_MAX
cdef long LONG_MIN
###############################################################################
# Various constants exported to Python
###############################################################################
# Make LONG_MAX and LONG_MIN available to Python. These are the limits of
# values accepted by functions like e.g., mpfr_set_si.
_LONG_MAX = LONG_MAX
_LONG_MIN = LONG_MIN
# Make precision limits available to Python
MPFR_PREC_MIN = cmpfr.MPFR_PREC_MIN
MPFR_PREC_MAX = cmpfr.MPFR_PREC_MAX
# Make rounding mode values available to Python
MPFR_RNDN = cmpfr.MPFR_RNDN
MPFR_RNDZ = cmpfr.MPFR_RNDZ
MPFR_RNDU = cmpfr.MPFR_RNDU
MPFR_RNDD = cmpfr.MPFR_RNDD
MPFR_RNDA = cmpfr.MPFR_RNDA
# Default values for Emax and Emin
MPFR_EMAX_DEFAULT = cmpfr.MPFR_EMAX_DEFAULT
MPFR_EMIN_DEFAULT = cmpfr.MPFR_EMIN_DEFAULT
###############################################################################
# Helper functions, not exposed to Python
###############################################################################
# Forward declaration
cdef class Mpfr_t
# Checks for valid parameter ranges
cdef int check_rounding_mode(cmpfr.mpfr_rnd_t rnd) except -1:
"""
Check that the given rounding mode is valid. Raise ValueError if not.
"""
if MPFR_RNDN <= rnd <= MPFR_RNDA:
return 0
else:
raise ValueError("invalid rounding mode {}".format(rnd))
cdef int check_base(int b, int allow_zero) except -1:
"""
Check that the given base (for string conversion) is valid.
Raise ValueError if not.
"""
if allow_zero:
if 2 <= b <= 62 or b == 0:
return 0
else:
raise ValueError(
"base should be zero or in the range 2 to 62 (inclusive)"
)
else:
if 2 <= b <= 62:
return 0
else:
raise ValueError("base should be in the range 2 to 62 (inclusive)")
cdef int check_get_str_n(int b, size_t n) except -1:
"""
Check that the given number of requested digits is valid
for the given base.
Raise ValueError if not.
"""
if n == 0 or 2 <= n:
return 0
# Undocumented MPFR feature: mpfr_get_str works for n = 1
# in non-power-of-two bases. See comments in vasprintf.c and
# get_str.c.
if n == 1 and (b & (b-1) != 0):
return 0
raise ValueError("n should be either 0 or at least 2")
cdef int check_precision(cmpfr.mpfr_prec_t precision) except -1:
"""
Check that the given precision is valid. Raise ValueError if not.
"""
if MPFR_PREC_MIN <= precision <= MPFR_PREC_MAX:
return 0
else:
raise ValueError(
"precision should be between {} and {}".format(
MPFR_PREC_MIN, MPFR_PREC_MAX
)
)
cdef int check_initialized(Mpfr_t x) except -1:
"""
Check that the given Mpfr_t x instance has been initialized.
Raise ValueError if not.
"""
if not cmpfr_initialized_p(&x._value):
raise ValueError(
"Mpfr_t instance {} should be initialized before use".format(x)
)
cdef int check_not_initialized(Mpfr_t x) except -1:
"""
Check that the given Mpfr_t x instance has *not* been initialized.
Raise ValueError if it has. This function is used by the mpfr_init
and mpfr_init2 functions.
"""
if cmpfr_initialized_p(&x._value):
raise ValueError(
"Mpfr_t instance {} is already initialized.".format(x)
)
cdef decode_ternary_pair(int ternary_pair):
"""
Decode an encoded pair of ternary values.
Some MPFR functions with two outputs (mpfr_sin_cos, mpfr_sinh_cosh,
mpfr_modf) also return a pair of ternary values encoded into a single int.
This function decodes that ternary pair, returning a Python pair of
the corresponding ternary values.
"""
cdef int first_ternary, second_ternary
first_ternary = ternary_pair & 3
if first_ternary == 2:
first_ternary = -1
second_ternary = ternary_pair >> 2
if second_ternary == 2:
second_ternary = -1
return first_ternary, second_ternary
cdef int cmpfr_initialized_p(cmpfr.mpfr_ptr op):
"""
Return non-zero if op is initialized. Return zero otherwise.
"""
return op._mpfr_d != NULL
###############################################################################
# The main Python extension type, based on mpfr_t.
###############################################################################
cdef class Mpfr_t:
"""
Mutable arbitrary-precision binary floating-point numbers.
Mpfr_t() -> new, uninitialized Mpfr object
Mpfr_t() creates a new, uninitialized Mpfr_t object. This object must be
initialized before use, for example by using the mpfr_init2 function.
However, unlike the underlying MPFR library, it's not necessary to clear
the object when it's no longer used.
"""
cdef cmpfr.__mpfr_struct _value
def __dealloc__(self):
if cmpfr_initialized_p(&self._value):
cmpfr.mpfr_clear(&self._value)
##############################################################################
# Additional functions provided by this extension module
###############################################################################
def mpfr_initialized_p(Mpfr_t op not None):
"""
Return True if op has been initialized. Return False otherwise.
"""
return bool(cmpfr_initialized_p(&op._value))
##############################################################################
# 5.1 Initialization Functions
###############################################################################
def mpfr_init2(Mpfr_t x not None, cmpfr.mpfr_prec_t prec):
"""
Initialize x, set its precision to be prec bits and its value to NaN.
Normally, a variable should be initialized once only or at least be
cleared, using mpfr_clear, between initializations. To change the precision
of a variable which has already been initialized, use mpfr_set_prec. The
precision prec must be an integer between MPFR_PREC_MIN and MPFR_PREC_MAX
(otherwise the behavior is undefined).
"""
check_not_initialized(x)
check_precision(prec)
cmpfr.mpfr_init2(&x._value, prec)
def mpfr_inits2(cmpfr.mpfr_prec_t prec, *args):
"""
Initialize each of the variables in args: set its precision to prec bits
and its value to NaN.
"""
check_precision(prec)
for arg in args:
mpfr_init2(arg, prec)
def mpfr_clear(Mpfr_t x not None):
"""
Free the space occupied by the significand of x.
It's not usually necessary to call this function directly: it will be
called automatically when x is garbage-collected.
"""
check_initialized(x)
cmpfr.mpfr_clear(&x._value)
def mpfr_clears(*args):
"""
Free the space occupied by each of the variables in args. See mpfr_clear
for more details.
"""
for arg in args:
mpfr_clear(arg)
def mpfr_init(Mpfr_t x not None):
"""
Initialize x, set its precision to the default precision, and set its value
to NaN. The default precision can be changed by a call to
mpfr_set_default_prec.
Warning! In a given program, some other libraries might change the default
precision and not restore it. Thus it is safer to use mpfr_init2.
"""
check_not_initialized(x)
cmpfr.mpfr_init(&x._value)
def mpfr_inits(*args):
"""
Initialize each of the variables in args: set its precision to the default
precision, and set its value to NaN. The default precision can be changed
by a call to mpfr_set_default_prec.
Warning! In a given program, some other libraries might change the default
precision and not restore it. Thus it is safer to use mpfr_inits2.
"""
for arg in args:
mpfr_init(arg)
def mpfr_set_default_prec(cmpfr.mpfr_prec_t prec):
"""
Set the default precision to be exactly prec bits.
prec can be any integer between MPFR_PREC_MIN and MPFR_PREC_MAX. The
precision of a variable means the number of bits used to store its
significand. All subsequent calls to mpfr_init or mpfr_inits will use this
precision, but previously initialized variables are unaffected. The default
precision is set to 53 bits initially.
"""
check_precision(prec)
cmpfr.mpfr_set_default_prec(prec)
def mpfr_get_default_prec():
"""
Return the current default MPFR precision in bits.
"""
return cmpfr.mpfr_get_default_prec()
def mpfr_set_prec(Mpfr_t x not None, cmpfr.mpfr_prec_t prec):
"""
Reset precision of x.
Reset the precision of x to be exactly prec bits, and set its value to
NaN. The previous value stored in x is lost. It is equivalent to a call to
mpfr_clear(x) followed by a call to mpfr_init2(x, prec), but more efficient
as no allocation is done in case the current allocated space for the
significand of x is enough. The precision prec can be any integer between
MPFR_PREC_MIN and MPFR_PREC_MAX. In case you want to keep the previous
value stored in x, use mpfr_prec_round instead.
"""
check_initialized(x)
check_precision(prec)
cmpfr.mpfr_set_prec(&x._value, prec)
def mpfr_get_prec(Mpfr_t x not None):
"""
Return the precision of x
Returns the number of bits used to store the significand of x.
"""
check_initialized(x)
return cmpfr.mpfr_get_prec(&x._value)
###############################################################################
# 5.2 Assignment Functions
###############################################################################
def mpfr_set(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop from op, rounded in the direction rnd.
Set the value of rop from the value of the Mpfr_t object op, rounded toward
the given direction rnd.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_set(&rop._value, &op._value, rnd)
def mpfr_set_si(Mpfr_t rop not None, long int op, cmpfr.mpfr_rnd_t rnd):
"""
Set the value of rop from a Python int, rounded in the direction rnd.
Set the value of rop from op, rounded toward the given direction rnd. Note
that the input 0 is converted to +0, regardless of the rounding mode.
"""
check_initialized(rop)
check_rounding_mode(rnd)
return cmpfr.mpfr_set_si(&rop._value, op, rnd)
def mpfr_set_d(Mpfr_t rop not None, double op, cmpfr.mpfr_rnd_t rnd):
"""
Set the value of rop from a Python float op, rounded in the direction rnd.
Set the value of rop from op, rounded toward the given direction rnd. If
the system does not support the IEEE 754 standard, mpfr_set_d might not
preserve the signed zeros.
Note: If you want to store a floating-point constant to an Mpfr_t object, you
should use mpfr_set_str (or one of the MPFR constant functions, such as
mpfr_const_pi for Pi) instead of mpfr_set_d. Otherwise the floating-point
constant will be first converted into a reduced-precision (e.g., 53-bit)
binary number before MPFR can work with it.
"""
check_initialized(rop)
check_rounding_mode(rnd)
return cmpfr.mpfr_set_d(&rop._value, op, rnd)
def mpfr_set_si_2exp(Mpfr_t rop not None, long int op,
cmpfr.mpfr_exp_t e, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op multiplied by a power of 2.
Set the value of rop from op multiplied by two to the power e, rounded
toward the given direction rnd. Note that the input 0 is converted to +0.
"""
check_initialized(rop)
check_rounding_mode(rnd)
return cmpfr.mpfr_set_si_2exp(&rop._value, op, e, rnd)
def mpfr_set_str(Mpfr_t rop not None, object s, int base, cmpfr.mpfr_rnd_t rnd):
"""
Set rop from a string s.
Set rop to the value of the string s in base base, rounded in the direction
rnd. See the documentation of mpfr_strtofr for a detailed description of
the valid string formats. Contrary to mpfr_strtofr, mpfr_set_str requires
the whole string to represent a valid floating-point number. This function
returns 0 if the entire string up to the final null character is a valid
number in base base; otherwise it returns −1, and rop may have
changed. Note: it is preferable to use mpfr_strtofr if one wants to
distinguish between an infinite rop value coming from an infinite s or from
an overflow.
"""
cdef bytes bytes_s
bytes_s = s.encode('ascii')
check_initialized(rop)
check_base(base, False)
check_rounding_mode(rnd)
return cmpfr.mpfr_set_str(&rop._value, bytes_s, base, rnd)
def mpfr_strtofr(Mpfr_t rop not None, object s, int base, cmpfr.mpfr_rnd_t rnd):
"""
Read a floating-point number from a string.
Read a floating-point number from a string s in base base, rounded in
the direction rnd; base must be either 0 (to detect the base, as described
below) or a number from 2 to 62 (otherwise the behavior is undefined).
If s starts with valid data, the result is stored in rop and the function
returns a pair (ternary, endindex) where ternary is the usual ternary
return value and endindex gives the index of the character just after the
valid data. Otherwise rop is set to zero (for consistency with strtod) and
endindex is 0.
Parsing follows the standard C strtod function with some extensions. After
optional leading whitespace, one has a subject sequence consisting of an
optional sign (+ or -), and either numeric data or special data. The
subject sequence is defined as the longest initial subsequence of the input
string, starting with the first non-whitespace character, that is of the
expected form.
The form of numeric data is a non-empty sequence of significand digits with
an optional decimal point, and an optional exponent consisting of an
exponent prefix followed by an optional sign and a non-empty sequence of
decimal digits. A significand digit is either a decimal digit or a Latin
letter (62 possible characters), with A = 10, B = 11, ..., Z = 35; case is
ignored in bases less or equal to 36, in bases larger than 36, a = 36, b =
37, ..., z = 61. The value of a significand digit must be strictly less
than the base. The decimal point can be either the one defined by the
current locale or the period (the first one is accepted for consistency
with the C standard and the practice, the second one is accepted to allow
the programmer to provide MPFR numbers from strings in a way that does not
depend on the current locale). The exponent prefix can be e or E for bases
up to 10, or @ in any base; it indicates a multiplication by a power of the
base. In bases 2 and 16, the exponent prefix can also be p or P, in which
case the exponent, called binary exponent, indicates a multiplication by a
power of 2 instead of the base (there is a difference only for base 16); in
base 16 for example 1p2 represents 4 whereas 1@2 represents 256. The value
of an exponent is always written in base 10.
If the argument base is 0, then the base is automatically detected as
follows. If the significand starts with 0b or 0B, base 2 is assumed. If the
significand starts with 0x or 0X, base 16 is assumed. Otherwise base 10 is
assumed.
Note: The exponent (if present) must contain at least a digit. Otherwise
the possible exponent prefix and sign are not part of the number (which
ends with the significand). Similarly, if 0b, 0B, 0x or 0X is not followed
by a binary/hexadecimal digit, then the subject sequence stops at the
character 0, thus 0 is read.
Special data (for infinities and NaN) can be @inf@ or
@nan@(n-char-sequence-opt), and if base <= 16, it can also be infinity,
inf, nan or nan(n-char-sequence-opt), all case insensitive. A
n-char-sequence-opt is a possibly empty string containing only digits,
Latin letters and the underscore (0, 1, 2, ..., 9, a, b, ..., z, A, B, ...,
Z, _). Note: one has an optional sign for all data, even NaN. For example,
-@nAn@(This_Is_Not_17) is a valid representation for NaN in base 17.
"""
cdef char* endptr
cdef char* startptr
cdef bytes bytes_s
bytes_s = s.encode('ascii')
startptr = bytes_s
check_initialized(rop)
check_base(base, True)
check_rounding_mode(rnd)
ternary = cmpfr.mpfr_strtofr(
&rop._value,
bytes_s,
&endptr,
base,
rnd,
)
endindex = endptr - startptr
return ternary, endindex
def mpfr_set_nan(Mpfr_t op not None):
""" Set x to a NaN.
Set the variable x to NaN (Not-a-Number). The sign bit of the result is
unspecified.
"""
check_initialized(op)
cmpfr.mpfr_set_nan(&op._value)
def mpfr_set_inf(Mpfr_t op not None, int sign):
""" Set x to an infinity.
Set the variable x to infinity. x is set to positive infinity if the sign
is nonnegative, and negative infinity otherwise.
Note the unusual sign convention here: most MPFR functions that deal with
signs use a nonzero value (or True) to indicate a negative number, and zero
to indiciate a positive number.
"""
check_initialized(op)
cmpfr.mpfr_set_inf(&op._value, sign)
def mpfr_set_zero(Mpfr_t op not None, int sign):
""" Set x to a zero.
Set the variable x to zero. x is set to positive zero if the sign is
nonnegative, and negative zero otherwise.
Note the unusual sign convention here: most MPFR functions that deal with
signs use a nonzero value (or True) to indicate a negative number, and zero
to indiciate a positive number.
"""
check_initialized(op)
cmpfr.mpfr_set_zero(&op._value, sign)
def mpfr_swap(Mpfr_t x not None, Mpfr_t y not None):
"""
Swap the values of x and y efficiently.
Warning: the precisions are exchanged too; in case the precisions are
different, mpfr_swap is thus not equivalent to three mpfr_set calls using a
third auxiliary variable.
"""
check_initialized(x)
check_initialized(y)
cmpfr.mpfr_swap(&x._value, &y._value)
###############################################################################
# 5.4 Conversion Functions
###############################################################################
def mpfr_get_d(Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Convert op to a Python float.
Convert op to a Python float using the rounding mode rnd. If op is NaN,
some fixed NaN (either quiet or signaling) or the result of 0.0/0.0 is
returned. If op is ±Inf, an infinity of the same sign or the result of
±1.0/0.0 is returned. If op is zero, this function returns a zero, trying
to preserve its sign, if possible.
"""
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_get_d(&op._value, rnd)
def mpfr_get_si(Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Convert op to a Python int.
Convert op to a Python int after rounding it with respect to rnd. If op is
NaN, 0 is returned and the erange flag is set. If op is too big for a
Python int, the function returns the maximum or the minimum representable
int, depending on the direction of the overflow; the erange flag is set
too.
"""
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_get_si(&op._value, rnd)
def mpfr_get_d_2exp(Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Convert op to a Python float and an exponent.
Return a pair (d, exp) consisting of a Python float d and an exponent exp
such that 0.5<=abs(d)<1 and d times 2 raised to exp equals op rounded to
double (resp. long double) precision, using the given rounding mode. If op
is zero, then a zero of the same sign (or an unsigned zero, if the
implementation does not have signed zeros) is returned, and exp is set to
0. If op is NaN or an infinity, then the corresponding double precision
(resp. long-double precision) value is returned, and exp is undefined.
"""
cdef long int exp
cdef double d
check_initialized(op)
check_rounding_mode(rnd)
d = cmpfr.mpfr_get_d_2exp(&exp, &op._value, rnd)
return d, exp
def mpfr_get_str(int b, size_t n, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Compute a base 'b' string representation for 'op'.
Convert op to a string of digits in base b, with rounding in the direction
rnd, where n is either zero (see below) or the number of significant digits
output in the string; in the latter case, n must be greater or equal to
2. The base may vary from 2 to 62. Returns a pair (digits, exp) where
digits gives the base-b digits of op, and for an ordinary number, exp is
the exponent (for input 0, the current minimal exponent is written).
The generated string is a fraction, with an implicit radix point
immediately to the left of the first digit. For example, the number −3.1416
would be returned as ("−31416", 1). If rnd is to nearest, and op is exactly
in the middle of two consecutive possible outputs, the one with an even
significand is chosen, where both significands are considered with the
exponent of op. Note that for an odd base, this may not correspond to an
even last digit: for example with 2 digits in base 7, (14) and a half is
rounded to (15) which is 12 in decimal, (16) and a half is rounded to (20)
which is 14 in decimal, and (26) and a half is rounded to (26) which is 20
in decimal.
If n is zero, the number of digits of the significand is chosen large
enough so that re-reading the printed value with the same precision,
assuming both output and input use rounding to nearest, will recover the
original value of op. More precisely, in most cases, the chosen precision
of str is the minimal precision m depending only on p = PREC(op) and b that
satisfies the above property, i.e., m = 1 + ceil(p*log(2)/log(b)), with p
replaced by p−1 if b is a power of 2, but in some very rare cases, it might
be m+1 (the smallest case for bases up to 62 is when p equals 186564318007
for bases 7 and 49).
Space for the digit string is automatically allocated, and freed by Python
when no longer needed. There's no requirement to free this space manually.
RuntimeError is raised on error.
"""
cdef cmpfr.mpfr_exp_t exp
cdef bytes digits
cdef char *c_digits
check_base(b, False)
check_get_str_n(b, n)
check_initialized(op)
check_rounding_mode(rnd)
c_digits = cmpfr.mpfr_get_str(NULL, &exp, b, n, &op._value, rnd)
if c_digits == NULL:
raise RuntimeError("Error during string conversion.")
# It's possible for the conversion from c_digits to digits to raise, so use
# a try-finally block to ensure that c_digits always gets freed.
try:
digits = bytes(c_digits)
finally:
cmpfr.mpfr_free_str(c_digits)
if sys.version_info < (3,):
return digits, exp
else:
return digits.decode('ascii'), exp
def mpfr_fits_slong_p(Mpfr_t x not None, cmpfr.mpfr_rnd_t rnd):
"""
Return True if op would fit into a Python int.
Return True if op would fit into a Python int when rounded to an integer
in the direction rnd.
"""
check_initialized(x)
check_rounding_mode(rnd)
return bool(cmpfr.mpfr_fits_slong_p(&x._value, rnd))
###############################################################################
# 5.5 Basic Arithmetic Functions
###############################################################################
def mpfr_add(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op1 + op2 rounded in the direction rnd.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_add(&rop._value, &op1._value, &op2._value, rnd)
def mpfr_sub(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op1 - op2, rounded in the direction rnd.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_sub(&rop._value, &op1._value, &op2._value, rnd)
def mpfr_mul(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op1 times op2, rounded in the direction rnd.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_mul(&rop._value, &op1._value, &op2._value, rnd)
def mpfr_sqr(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the square of op, rounded in the direction rnd.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_sqr(&rop._value, &op._value, rnd)
def mpfr_div(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op1 divided by op2, rounded in the direction rnd.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_div(&rop._value, &op1._value, &op2._value, rnd)
def mpfr_sqrt(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the square root of op, rounded in the direction rnd.
Set rop to −0 if op is −0, to be consistent with the IEEE 754 standard. Set
rop to NaN if op is negative.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_sqrt(&rop._value, &op._value, rnd)
def mpfr_rec_sqrt(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the reciprocal square root of op, rounded in the direction rnd.
Set rop to +Inf if op is ±0, +0 if op is +Inf, and NaN if op is negative.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_rec_sqrt(&rop._value, &op._value, rnd)
def mpfr_cbrt(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the cube root of op rounded in the direction rnd.
For op negative, set rop to a negative number. The cube root of -0 is
defined to be -0.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_cbrt(&rop._value, &op._value, rnd)
def mpfr_root(Mpfr_t rop not None, Mpfr_t op not None,
unsigned long int k, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the kth root of op, rounding in the direction rnd.
For k odd (resp. even) and op negative (including −Inf), set rop to a
negative number (resp. NaN). The kth root of −0 is defined to be −0,
whatever the parity of k.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_root(&rop._value, &op._value, k, rnd)
def mpfr_pow(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to op1 raised to the power op2, rounded in the direction rnd.
Special values are handled as described in the ISO C99 and IEEE 754-2008
standards for the pow function.
* pow(±0, y) returns plus or minus infinity for y a negative odd integer.
* pow(±0, y) returns plus infinity for y negative and not an odd integer.
* pow(±0, y) returns plus or minus zero for y a positive odd integer.
* pow(±0, y) returns plus zero for y positive and not an odd integer.
* pow(-1, ±Inf) returns 1.
* pow(+1, y) returns 1 for any y, even a NaN.
* pow(x, ±0) returns 1 for any x, even a NaN.
* pow(x, y) returns NaN for finite negative x and finite non-integer y.
* pow(x, -Inf) returns plus infinity for 0 < abs(x) < 1, and plus zero
for abs(x) > 1.
* pow(x, +Inf) returns plus zero for 0 < abs(x) < 1, and plus infinity
for abs(x) > 1.
* pow(-Inf, y) returns minus zero for y a negative odd integer.
* pow(-Inf, y) returns plus zero for y negative and not an odd integer.
* pow(-Inf, y) returns minus infinity for y a positive odd integer.
* pow(-Inf, y) returns plus infinity for y positive and not an odd
integer.
* pow(+Inf, y) returns plus zero for y negative, and plus infinity for y
positive.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_pow(&rop._value, &op1._value, &op2._value, rnd)
def mpfr_neg(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to -op, rounded in the direction rnd.
This function just changes or adjusts the sign if rop and op are the same
variable, otherwise a rounding might occur if the precision of rop is less
than that of op.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_neg(&rop._value, &op._value, rnd)
def mpfr_abs(Mpfr_t rop not None, Mpfr_t op not None, cmpfr.mpfr_rnd_t rnd):
"""
Set rop to the absolute value of op, rounded in the direction rnd.
This function just changes or adjusts the sign if rop and op are the same
variable, otherwise a rounding might occur if the precision of rop is less
than that of op.
"""
check_initialized(rop)
check_initialized(op)
check_rounding_mode(rnd)
return cmpfr.mpfr_abs(&rop._value, &op._value, rnd)
def mpfr_dim(Mpfr_t rop not None, Mpfr_t op1 not None, Mpfr_t op2 not None,
cmpfr.mpfr_rnd_t rnd):
"""
Set rop to max(op1 - op2, 0), rounded in the direction rnd.
Set rop to op1 - op2 rounded in the direction rnd if op1 > op2, +0 if op1
<= op2, and NaN if op1 or op2 is NaN.
"""
check_initialized(rop)
check_initialized(op1)
check_initialized(op2)
check_rounding_mode(rnd)
return cmpfr.mpfr_dim(&rop._value, &op1._value, &op2._value, rnd)
###############################################################################
# 5.6 Comparison Functions
###############################################################################
def mpfr_cmp(Mpfr_t op1 not None, Mpfr_t op2 not None):
"""
Perform a three-way comparison of op1 and op2.
Return a positive value if op1 > op2, zero if op1 = op2, and a negative
value if op1 < op2. Both op1 and op2 are considered to their full own
precision, which may differ. If one of the operands is NaN, set the erange
flag and return zero.
Note: This function may be useful to distinguish the three possible
cases. If you need to distinguish two cases only, it is recommended to use
the predicate functions (e.g., mpfr_equal_p for the equality) described
below; they behave like the IEEE 754 comparisons, in particular when one or
both arguments are NaN.
"""
check_initialized(op1)
check_initialized(op2)
return cmpfr.mpfr_cmp(&op1._value, &op2._value)
def mpfr_cmpabs(Mpfr_t op1 not None, Mpfr_t op2 not None):
"""
Compare the absolute values of op1 and op2.
Compare |op1| and |op2|. Return a positive value if |op1| > |op2|, zero if
|op1| == |op2|, and a negative value if |op1| < |op2|. If one of the
operands is NaN, set the erange flag and return zero.
"""
check_initialized(op1)
check_initialized(op2)
return cmpfr.mpfr_cmpabs(&op1._value, &op2._value)
def mpfr_nan_p(Mpfr_t op not None):
"""
Return True if op is a NaN. Return False otherwise.
"""
check_initialized(op)
return bool(cmpfr.mpfr_nan_p(&op._value))
def mpfr_inf_p(Mpfr_t op not None):
"""
Return True if op is an infinity. Return False otherwise.
"""
check_initialized(op)
return bool(cmpfr.mpfr_inf_p(&op._value))
def mpfr_number_p(Mpfr_t op not None):
"""
Return True if op is an ordinary number. Return False otherwise.
An ordinary number is a number which is neither a NaN nor an infinity.
"""
check_initialized(op)
return bool(cmpfr.mpfr_number_p(&op._value))
def mpfr_zero_p(Mpfr_t op not None):
"""
Return True if op is zero. Return False otherwise.
"""
check_initialized(op)
return bool(cmpfr.mpfr_zero_p(&op._value))
def mpfr_regular_p(Mpfr_t op not None):
"""
Return True if op is a regular number. Return False otherwise.
A regular number is a number which is neither a NaN, nor an infinity, nor a
zero.
"""
check_initialized(op)
return bool(cmpfr.mpfr_regular_p(&op._value))
def mpfr_sgn(Mpfr_t op not None):
"""
Return the sign of op.
Return a positive value if op > 0, zero if op = 0, and a negative value if
op < 0. If the operand is NaN, set the erange flag and return zero. This is
equivalent to mpfr_cmp_ui (op, 0), but more efficient.
"""
check_initialized(op)
return cmpfr.mpfr_sgn(&op._value)
def mpfr_greater_p(Mpfr_t op1 not None, Mpfr_t op2 not None):
"""
Return True if op1 > op2 and False otherwise.
This function returns False whenever op1 and/or op2 is a NaN.
"""
check_initialized(op1)
check_initialized(op2)
return bool(cmpfr.mpfr_greater_p(&op1._value, &op2._value))