1 /* mpn_toom2_sqr -- Square {ap,an}.
2
3 Contributed to the GNU project by Torbjorn Granlund.
4
5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
6 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
7 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
8
9 Copyright 2006-2010, 2012, 2014, 2018 Free Software Foundation, Inc.
10
11 This file is part of the GNU MP Library.
12
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of either:
15
16 * the GNU Lesser General Public License as published by the Free
17 Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
19
20 or
21
22 * the GNU General Public License as published by the Free Software
23 Foundation; either version 2 of the License, or (at your option) any
24 later version.
25
26 or both in parallel, as here.
27
28 The GNU MP Library is distributed in the hope that it will be useful, but
29 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
30 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
31 for more details.
32
33 You should have received copies of the GNU General Public License and the
34 GNU Lesser General Public License along with the GNU MP Library. If not,
35 see https://www.gnu.org/licenses/. */
36
37
38 #include "gmp-impl.h"
39
40 /* Evaluate in: -1, 0, +inf
41
42 <-s--><--n-->
43 ____ ______
44 |_a1_|___a0_|
45
46 v0 = a0 ^2 # A(0)^2
47 vm1 = (a0- a1)^2 # A(-1)^2
48 vinf= a1 ^2 # A(inf)^2
49 */
50
51 #if TUNE_PROGRAM_BUILD || WANT_FAT_BINARY
52 #define MAYBE_sqr_toom2 1
53 #else
54 #define MAYBE_sqr_toom2 \
55 (SQR_TOOM3_THRESHOLD >= 2 * SQR_TOOM2_THRESHOLD)
56 #endif
57
58 #define TOOM2_SQR_REC(p, a, n, ws) \
59 do { \
60 if (! MAYBE_sqr_toom2 \
61 || BELOW_THRESHOLD (n, SQR_TOOM2_THRESHOLD)) \
62 mpn_sqr_basecase (p, a, n); \
63 else \
64 mpn_toom2_sqr (p, a, n, ws); \
65 } while (0)
66
67 void
mpn_toom2_sqr(mp_ptr pp,mp_srcptr ap,mp_size_t an,mp_ptr scratch)68 mpn_toom2_sqr (mp_ptr pp,
69 mp_srcptr ap, mp_size_t an,
70 mp_ptr scratch)
71 {
72 const int __gmpn_cpuvec_initialized = 1;
73 mp_size_t n, s;
74 mp_limb_t cy, cy2;
75 mp_ptr asm1;
76
77 #define a0 ap
78 #define a1 (ap + n)
79
80 s = an >> 1;
81 n = an - s;
82
83 ASSERT (0 < s && s <= n && s >= n - 1);
84
85 asm1 = pp;
86
87 /* Compute asm1. */
88 if (s == n)
89 {
90 if (mpn_cmp (a0, a1, n) < 0)
91 {
92 mpn_sub_n (asm1, a1, a0, n);
93 }
94 else
95 {
96 mpn_sub_n (asm1, a0, a1, n);
97 }
98 }
99 else /* n - s == 1 */
100 {
101 if (a0[s] == 0 && mpn_cmp (a0, a1, s) < 0)
102 {
103 mpn_sub_n (asm1, a1, a0, s);
104 asm1[s] = 0;
105 }
106 else
107 {
108 asm1[s] = a0[s] - mpn_sub_n (asm1, a0, a1, s);
109 }
110 }
111
112 #define v0 pp /* 2n */
113 #define vinf (pp + 2 * n) /* s+s */
114 #define vm1 scratch /* 2n */
115 #define scratch_out scratch + 2 * n
116
117 /* vm1, 2n limbs */
118 TOOM2_SQR_REC (vm1, asm1, n, scratch_out);
119
120 /* vinf, s+s limbs */
121 TOOM2_SQR_REC (vinf, a1, s, scratch_out);
122
123 /* v0, 2n limbs */
124 TOOM2_SQR_REC (v0, ap, n, scratch_out);
125
126 /* H(v0) + L(vinf) */
127 cy = mpn_add_n (pp + 2 * n, v0 + n, vinf, n);
128
129 /* L(v0) + H(v0) */
130 cy2 = cy + mpn_add_n (pp + n, pp + 2 * n, v0, n);
131
132 /* L(vinf) + H(vinf) */
133 cy += mpn_add (pp + 2 * n, pp + 2 * n, n, vinf + n, s + s - n);
134
135 cy -= mpn_sub_n (pp + n, pp + n, vm1, 2 * n);
136
137 ASSERT (cy + 1 <= 3);
138 ASSERT (cy2 <= 2);
139
140 if (LIKELY (cy <= 2)) {
141 MPN_INCR_U (pp + 2 * n, s + s, cy2);
142 MPN_INCR_U (pp + 3 * n, s + s - n, cy);
143 } else { /* cy is negative */
144 /* The total contribution of v0+vinf-vm1 can not be negative. */
145 #if WANT_ASSERT
146 /* The borrow in cy stops the propagation of the carry cy2, */
147 ASSERT (cy2 == 1);
148 cy += mpn_add_1 (pp + 2 * n, pp + 2 * n, n, cy2);
149 ASSERT (cy == 0);
150 #else
151 /* we simply fill the area with zeros. */
152 MPN_FILL (pp + 2 * n, n, 0);
153 #endif
154 }
155 }
156