Tcl Source Code

Check-in [3f35b52355]
Login
Bounty program for improvements to Tcl and certain Tcl packages.
Tcl 2019 Conference, Houston/TX, US, Nov 4-8
Send your abstracts to [email protected]
or submit via the online form by Sep 9.

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:Only use special mp_sqrt() code when double format/tommath format are exactly what's expected. Otherwise, use original always-working tommath code. Simplify overflow check in bignum expononent code, not using bignums where it's not necessary. Don't overallocate bignums when using wideint's only.
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | core-8-6-branch
Files: files | file ages | folders
SHA3-256: 3f35b523558b49db296fa34756dad9bde0df6734eb67006112cdc017408fec64
User & Date: jan.nijtmans 2019-04-11 20:09:27
Context
2019-04-17
14:28
Isolate tests of [info frame] results from testing environment. check-in: 572f113bbb user: dgp tags: core-8-6-branch
2019-04-11
20:37
Merge 8.6 check-in: 2a6c012bff user: jan.nijtmans tags: core-8-branch
20:09
Only use special mp_sqrt() code when double format/tommath format are exactly what's expected. Other... check-in: 3f35b52355 user: jan.nijtmans tags: core-8-6-branch
2019-04-09
19:39
merge 8.5 check-in: 4cb9044dfa user: sebres tags: core-8-6-branch
Changes
Hide Diffs Unified Diffs Ignore Whitespace Patch

Changes to generic/tclExecute.c.

6028
6029
6030
6031
6032
6033
6034
6035
6036
6037
6038
6039
6040
6041
6042
6043
6044
6045
6046
6047
6048
6049
6050
6051
6052
6053
6054
6055
6056
6057
6058
6059
....
8980
8981
8982
8983
8984
8985
8986

8987
8988
8989
8990
8991
8992
8993
8994
8995
8996
8997
8998
8999
9000
9001
9002
9003
9004
9005
9006
9007
9008
9009
	if (GetNumberFromObj(NULL, OBJ_AT_TOS, &ptr1, &type1) != TCL_OK) {
	    type1 = 0;
	} else if (type1 == TCL_NUMBER_LONG) {
	    /* value is between LONG_MIN and LONG_MAX */
	    /* [string is integer] is -UINT_MAX to UINT_MAX range */
	    int i;

	    if (Tcl_GetIntFromObj(NULL, OBJ_AT_TOS, &i) != TCL_OK) {
		type1 = TCL_NUMBER_WIDE;
	    }
#ifndef TCL_WIDE_INT_IS_LONG
	} else if (type1 == TCL_NUMBER_WIDE) {
	    /* value is between WIDE_MIN and WIDE_MAX */
	    /* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
	    int i;
	    if (Tcl_GetIntFromObj(NULL, OBJ_AT_TOS, &i) == TCL_OK) {
		type1 = TCL_NUMBER_LONG;
	    }
#endif
	} else if (type1 == TCL_NUMBER_BIG) {
	    /* value is an integer outside the WIDE_MIN to WIDE_MAX range */
	    /* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
	    Tcl_WideInt w;

	    if (Tcl_GetWideIntFromObj(NULL, OBJ_AT_TOS, &w) == TCL_OK) {
		type1 = TCL_NUMBER_WIDE;
	    }
	}
	TclNewIntObj(objResultPtr, type1);
	TRACE(("\"%.20s\" => %d\n", O2S(OBJ_AT_TOS), type1));
	NEXT_INST_F(1, 1, 1);

................................................................................
		wResult = oddExponent ? -Exp64Value[base] : Exp64Value[base];
		WIDE_RESULT(wResult);
	    }
	}
#endif

    overflowExpon:

	Tcl_TakeBignumFromObj(NULL, value2Ptr, &big2);
	if ((big2.used > 1)
#if DIGIT_BIT > 28
		    || ((big2.used == 1) && (big2.dp[0] >= (1<<28)))
#endif
	) {
	    mp_clear(&big2);
	    Tcl_SetObjResult(interp, Tcl_NewStringObj(
		    "exponent too large", -1));
	    return GENERAL_ARITHMETIC_ERROR;
	}
	Tcl_TakeBignumFromObj(NULL, valuePtr, &big1);
	mp_init(&bigResult);
	mp_expt_d(&big1, big2.dp[0], &bigResult);
	mp_clear(&big1);
	mp_clear(&big2);
	BIG_RESULT(&bigResult);
    }

    case INST_ADD:
    case INST_SUB:
    case INST_MULT:
    case INST_DIV:






|







|








|







 







>
|
<
|
|
<
<
<






|

<







6028
6029
6030
6031
6032
6033
6034
6035
6036
6037
6038
6039
6040
6041
6042
6043
6044
6045
6046
6047
6048
6049
6050
6051
6052
6053
6054
6055
6056
6057
6058
6059
....
8980
8981
8982
8983
8984
8985
8986
8987
8988

8989
8990



8991
8992
8993
8994
8995
8996
8997
8998

8999
9000
9001
9002
9003
9004
9005
	if (GetNumberFromObj(NULL, OBJ_AT_TOS, &ptr1, &type1) != TCL_OK) {
	    type1 = 0;
	} else if (type1 == TCL_NUMBER_LONG) {
	    /* value is between LONG_MIN and LONG_MAX */
	    /* [string is integer] is -UINT_MAX to UINT_MAX range */
	    int i;

	    if (TclGetIntFromObj(NULL, OBJ_AT_TOS, &i) != TCL_OK) {
		type1 = TCL_NUMBER_WIDE;
	    }
#ifndef TCL_WIDE_INT_IS_LONG
	} else if (type1 == TCL_NUMBER_WIDE) {
	    /* value is between WIDE_MIN and WIDE_MAX */
	    /* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
	    int i;
	    if (TclGetIntFromObj(NULL, OBJ_AT_TOS, &i) == TCL_OK) {
		type1 = TCL_NUMBER_LONG;
	    }
#endif
	} else if (type1 == TCL_NUMBER_BIG) {
	    /* value is an integer outside the WIDE_MIN to WIDE_MAX range */
	    /* [string is wideinteger] is -UWIDE_MAX to UWIDE_MAX range */
	    Tcl_WideInt w;

	    if (TclGetWideIntFromObj(NULL, OBJ_AT_TOS, &w) == TCL_OK) {
		type1 = TCL_NUMBER_WIDE;
	    }
	}
	TclNewIntObj(objResultPtr, type1);
	TRACE(("\"%.20s\" => %d\n", O2S(OBJ_AT_TOS), type1));
	NEXT_INST_F(1, 1, 1);

................................................................................
		wResult = oddExponent ? -Exp64Value[base] : Exp64Value[base];
		WIDE_RESULT(wResult);
	    }
	}
#endif

    overflowExpon:

	if ((TclGetWideIntFromObj(NULL, value2Ptr, &w2) != TCL_OK)

		|| (value2Ptr->typePtr != &tclIntType)
		|| (Tcl_WideUInt)w2 >= (1<<28)) {



	    Tcl_SetObjResult(interp, Tcl_NewStringObj(
		    "exponent too large", -1));
	    return GENERAL_ARITHMETIC_ERROR;
	}
	Tcl_TakeBignumFromObj(NULL, valuePtr, &big1);
	mp_init(&bigResult);
	mp_expt_d(&big1, w2, &bigResult);
	mp_clear(&big1);

	BIG_RESULT(&bigResult);
    }

    case INST_ADD:
    case INST_SUB:
    case INST_MULT:
    case INST_DIV:

Changes to generic/tclTomMathInterface.c.

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
...
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
    unsigned long v;
    mp_digit *p;

    /*
     * Allocate enough memory to hold the largest possible long
     */

    status = mp_init_size(a,
	    (CHAR_BIT * sizeof(long) + DIGIT_BIT - 1) / DIGIT_BIT);
    if (status != MP_OKAY) {
	Tcl_Panic("initialization failure in TclBNInitBignumFromLong");
    }

    /*
     * Convert arg to sign and magnitude.
     */
................................................................................
    int status;
    mp_digit *p;

    /*
     * Allocate enough memory to hold the largest possible Tcl_WideUInt.
     */

    status = mp_init_size(a,
	    (CHAR_BIT * sizeof(Tcl_WideUInt) + DIGIT_BIT - 1) / DIGIT_BIT);
    if (status != MP_OKAY) {
	Tcl_Panic("initialization failure in TclBNInitBignumFromWideUInt");
    }

    a->sign = 0;

    /*






|
<







 







|
<







115
116
117
118
119
120
121
122

123
124
125
126
127
128
129
...
201
202
203
204
205
206
207
208

209
210
211
212
213
214
215
    unsigned long v;
    mp_digit *p;

    /*
     * Allocate enough memory to hold the largest possible long
     */

    status = mp_init(a);

    if (status != MP_OKAY) {
	Tcl_Panic("initialization failure in TclBNInitBignumFromLong");
    }

    /*
     * Convert arg to sign and magnitude.
     */
................................................................................
    int status;
    mp_digit *p;

    /*
     * Allocate enough memory to hold the largest possible Tcl_WideUInt.
     */

    status = mp_init(a);

    if (status != MP_OKAY) {
	Tcl_Panic("initialization failure in TclBNInitBignumFromWideUInt");
    }

    a->sign = 0;

    /*

Changes to libtommath/bn_mp_set_double.c.

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 *
 * SPDX-License-Identifier: Unlicense
 */

#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
int mp_set_double(mp_int *a, double b)
{
   uint64_t frac;
   int exp, res;
   union {
      double   dbl;
      uint64_t bits;
   } cast;
   cast.dbl = b;

   exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFU);
   frac = (cast.bits & ((1ULL << 52) - 1ULL)) | (1ULL << 52);

   if (exp == 0x7FF) { /* +-inf, NaN */






|



|







11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 *
 * SPDX-License-Identifier: Unlicense
 */

#if defined(__STDC_IEC_559__) || defined(__GCC_IEC_559)
int mp_set_double(mp_int *a, double b)
{
   unsigned long long frac;
   int exp, res;
   union {
      double   dbl;
      unsigned long long bits;
   } cast;
   cast.dbl = b;

   exp = (int)((unsigned)(cast.bits >> 52) & 0x7FFU);
   frac = (cast.bits & ((1ULL << 52) - 1ULL)) | (1ULL << 52);

   if (exp == 0x7FF) { /* +-inf, NaN */

Changes to libtommath/bn_mp_sqrt.c.

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
..
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
..
92
93
94
95
96
97
98
99

100
101
102

103





104
105
106
107
108
109
110
 * additional optimizations in place.
 *
 * SPDX-License-Identifier: Unlicense
 */

#ifndef NO_FLOATING_POINT
#include <math.h>



#endif

/* this function is less generic than mp_n_root, simpler and faster */
int mp_sqrt(const mp_int *arg, mp_int *ret)
{
   int res;
   mp_int t1, t2;
   int i, j, k;
#ifndef NO_FLOATING_POINT

   volatile double d;
   mp_digit dig;
#endif

   /* must be positive */
   if (arg->sign == MP_NEG) {
      return MP_VAL;
................................................................................
   }

   /* easy out */
   if (mp_iszero(arg) == MP_YES) {
      mp_zero(ret);
      return MP_OKAY;
   }



   i = (arg->used / 2) - 1;
   j = 2 * i;
   if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) {
      return res;
   }

................................................................................
      goto E2;
   }

   for (k = 0; k < i; ++k) {
      t1.dp[k] = (mp_digit) 0;
   }

#ifndef NO_FLOATING_POINT

   /* Estimate the square root using the hardware floating point unit. */

   d = 0.0;
   for (k = arg->used-1; k >= j; --k) {
      d = ldexp(d, DIGIT_BIT) + (double)(arg->dp[k]);
   }

................................................................................
   } else {
      t1.used = i+1;
      t1.dp[i] = ((mp_digit) d) - 1;
   }

#else

   /* Estimate the square root as having 1 in the most significant place. */


   t1.used = i + 2;
   t1.dp[i+1] = (mp_digit) 1;

   t1.dp[i] = (mp_digit) 0;






#endif

   /* t1 > 0  */
   if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
      goto E1;
   }






>
>
>







<

>







 







>
>







 







<
<







 







|
>
|
<
<
>
|
>
>
>
>
>







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
..
53
54
55
56
57
58
59


60
61
62
63
64
65
66
..
95
96
97
98
99
100
101
102
103
104


105
106
107
108
109
110
111
112
113
114
115
116
117
118
 * additional optimizations in place.
 *
 * SPDX-License-Identifier: Unlicense
 */

#ifndef NO_FLOATING_POINT
#include <math.h>
#if (DIGIT_BIT != 28) || (FLT_RADIX != 2) || (DBL_MANT_DIG != 53) || (DBL_MAX_EXP != 1024)
#define NO_FLOATING_POINT
#endif
#endif

/* this function is less generic than mp_n_root, simpler and faster */
int mp_sqrt(const mp_int *arg, mp_int *ret)
{
   int res;
   mp_int t1, t2;

#ifndef NO_FLOATING_POINT
   int i, j, k;
   volatile double d;
   mp_digit dig;
#endif

   /* must be positive */
   if (arg->sign == MP_NEG) {
      return MP_VAL;
................................................................................
   }

   /* easy out */
   if (mp_iszero(arg) == MP_YES) {
      mp_zero(ret);
      return MP_OKAY;
   }

#ifndef NO_FLOATING_POINT

   i = (arg->used / 2) - 1;
   j = 2 * i;
   if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) {
      return res;
   }

................................................................................
      goto E2;
   }

   for (k = 0; k < i; ++k) {
      t1.dp[k] = (mp_digit) 0;
   }



   /* Estimate the square root using the hardware floating point unit. */

   d = 0.0;
   for (k = arg->used-1; k >= j; --k) {
      d = ldexp(d, DIGIT_BIT) + (double)(arg->dp[k]);
   }

................................................................................
   } else {
      t1.used = i+1;
      t1.dp[i] = ((mp_digit) d) - 1;
   }

#else

   if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
      return res;
   }



   if ((res = mp_init(&t2)) != MP_OKAY) {
      goto E2;
   }

   /* First approx. (not very bad for large arg) */
   mp_rshd(&t1, t1.used/2);

#endif

   /* t1 > 0  */
   if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
      goto E1;
   }