1 /* $NetBSD: bn_mp_is_square.c,v 1.1.1.2 2014/04/24 12:45:31 pettai Exp $ */
2
3 #include <tommath.h>
4 #ifdef BN_MP_IS_SQUARE_C
5 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6 *
7 * LibTomMath is a library that provides multiple-precision
8 * integer arithmetic as well as number theoretic functionality.
9 *
10 * The library was designed directly after the MPI library by
11 * Michael Fromberger but has been written from scratch with
12 * additional optimizations in place.
13 *
14 * The library is free for all purposes without any express
15 * guarantee it works.
16 *
17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
18 */
19
20 /* Check if remainders are possible squares - fast exclude non-squares */
21 static const char rem_128[128] = {
22 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
23 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
24 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
25 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
26 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
27 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
28 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
29 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
30 };
31
32 static const char rem_105[105] = {
33 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
34 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
35 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
36 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
37 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
38 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
39 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
40 };
41
42 /* Store non-zero to ret if arg is square, and zero if not */
mp_is_square(mp_int * arg,int * ret)43 int mp_is_square(mp_int *arg,int *ret)
44 {
45 int res;
46 mp_digit c;
47 mp_int t;
48 unsigned long r;
49
50 /* Default to Non-square :) */
51 *ret = MP_NO;
52
53 if (arg->sign == MP_NEG) {
54 return MP_VAL;
55 }
56
57 /* digits used? (TSD) */
58 if (arg->used == 0) {
59 return MP_OKAY;
60 }
61
62 /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
63 if (rem_128[127 & DIGIT(arg,0)] == 1) {
64 return MP_OKAY;
65 }
66
67 /* Next check mod 105 (3*5*7) */
68 if ((res = mp_mod_d(arg,105,&c)) != MP_OKAY) {
69 return res;
70 }
71 if (rem_105[c] == 1) {
72 return MP_OKAY;
73 }
74
75
76 if ((res = mp_init_set_int(&t,11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
77 return res;
78 }
79 if ((res = mp_mod(arg,&t,&t)) != MP_OKAY) {
80 goto ERR;
81 }
82 r = mp_get_int(&t);
83 /* Check for other prime modules, note it's not an ERROR but we must
84 * free "t" so the easiest way is to goto ERR. We know that res
85 * is already equal to MP_OKAY from the mp_mod call
86 */
87 if ( (1L<<(r%11)) & 0x5C4L ) goto ERR;
88 if ( (1L<<(r%13)) & 0x9E4L ) goto ERR;
89 if ( (1L<<(r%17)) & 0x5CE8L ) goto ERR;
90 if ( (1L<<(r%19)) & 0x4F50CL ) goto ERR;
91 if ( (1L<<(r%23)) & 0x7ACCA0L ) goto ERR;
92 if ( (1L<<(r%29)) & 0xC2EDD0CL ) goto ERR;
93 if ( (1L<<(r%31)) & 0x6DE2B848L ) goto ERR;
94
95 /* Final check - is sqr(sqrt(arg)) == arg ? */
96 if ((res = mp_sqrt(arg,&t)) != MP_OKAY) {
97 goto ERR;
98 }
99 if ((res = mp_sqr(&t,&t)) != MP_OKAY) {
100 goto ERR;
101 }
102
103 *ret = (mp_cmp_mag(&t,arg) == MP_EQ) ? MP_YES : MP_NO;
104 ERR:mp_clear(&t);
105 return res;
106 }
107 #endif
108
109 /* Source: /cvs/libtom/libtommath/bn_mp_is_square.c,v */
110 /* Revision: 1.4 */
111 /* Date: 2006/12/28 01:25:13 */
112