File: gf2x.c

package info (click to toggle)
gf2x 1.1-2
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd, stretch
  • size: 2,612 kB
  • ctags: 725
  • sloc: sh: 11,112; ansic: 7,549; cpp: 1,369; makefile: 490; perl: 151
file content (217 lines) | stat: -rw-r--r-- 6,431 bytes parent folder | download
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
/* This file is part of the gf2x library.

   Copyright 2007, 2008, 2009
   Richard Brent, Pierrick Gaudry, Emmanuel Thome', Paul Zimmermann

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2 of the License, or (at your
   option) any later version.

   This program 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 General Public License for
   more details.

   You should have received a copy of the GNU General Public License along
   with this program; see the file COPYING.  If not, write to the Free
   Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
   02111-1307, USA.
*/

#include <string.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#include "gf2x.h"
#include "gf2x/gf2x-impl.h"

/* We _do_ need the inlines here. */
#include "gf2x/gf2x-small.h"

void gf2x_mul_basecase(unsigned long * c, const unsigned long * a,
			 long na, const unsigned long * b, long nb)
{
    assert(c != a);
    assert(c != b);

    if (na == nb) {
	switch (na) {
	case 0:		// This can occur in call from KarMul return;
	case 1: gf2x_mul1(c, a[0], b[0]); return;
	case 2: gf2x_mul2(c, a, b); return;
	case 3: gf2x_mul3(c, a, b); return;
	case 4: gf2x_mul4(c, a, b); return;
	case 5: gf2x_mul5(c, a, b); return;
	case 6: gf2x_mul6(c, a, b); return;
	case 7: gf2x_mul7(c, a, b); return;
	case 8: gf2x_mul8(c, a, b); return;
	case 9: gf2x_mul9(c, a, b); return;
	default:
	    fprintf(stderr, "basecase.c: ran off end of switch\n"
		    "na=nb=%ld ; decrease GF2X_MUL_KARA_THRESHOLD\n", na);
	    exit(1);
	}
    } else if (na < nb) {
        /* FIXME -- this does not seem efficient */
        long i;
        /* No need to clear c[], it's filled up progressively */
        c[nb] = gf2x_mul_1_n(c, b, nb, a[0]);
        for (i = 1; i < na; i++) {
            c[nb + i] = gf2x_addmul_1_n(c + i, c + i, b, nb, a[i]);
        }
    } else {
        gf2x_mul_basecase(c, b, nb, a, na);
    }
}

#ifdef GF2X_MUL_FFT_TABLE
int64_t T_FFT_TAB[][2] = GF2X_MUL_FFT_TABLE;
#endif

/* This is the toplevel multiplication routine. It handles the temporary
 * storage if necessary.
 */

/* Having this static ensures initialization to zero */
static gf2x_mul_pool_t global_pool;

void gf2x_mul_pool_init(gf2x_mul_pool_t p)
{
    memset(p, 0, sizeof(gf2x_mul_pool_t));
}

void gf2x_mul_pool_clear(gf2x_mul_pool_t p)
{
    free(p->stk);
    p->stk_size = 0;
}

void gf2x_mul(unsigned long * c,
        const unsigned long * a, unsigned int sa,
        const unsigned long * b, unsigned int sb)
{
    gf2x_mul_r(c, a, sa, b, sb, global_pool);
}

void gf2x_mul_r(unsigned long * c,
        const unsigned long * a, unsigned int sa,
        const unsigned long * b, unsigned int sb, gf2x_mul_pool_t pool)
{
    unsigned int sc = sa + sb;
    /* As a starting guess, assume that c may alias a or b */
    unsigned long * dst = c;

    if (sa > sb) {
        gf2x_mul_r(c, b, sb, a, sa, pool);
        return;
    }
    // now sa <= sb (note: sa and sb are interchanged in Toom3uMul etc)

    if (sa < GF2X_MUL_KARA_THRESHOLD) {
        /* This calls the hand-crafted code if sa == sb */
        if (dst == a || dst == b) {
            dst = malloc(sc * sizeof(unsigned long));
        }
        gf2x_mul_basecase(dst, a, sa, b, sb);
        goto end_of_gf2x_mul_r;
    }

    /* This ugly cpp block entirely disables the FFT if it has
     * not yet been tuned */
#ifdef GF2X_MUL_FFT_TABLE
    long ix, K, sab = sc / 2;
    long max_ix = sizeof(T_FFT_TAB)/sizeof(T_FFT_TAB[0]);
    for (ix = 0; T_FFT_TAB[ix + 1][0] <= sab && ix + 1 < max_ix; ix++);
    /* now T_FFT_TAB[ix][0] <= sab < T_FFT_TAB[ix+1][0] */
    K = T_FFT_TAB[ix][1];

    /* gf2x_mul_fft can handle unbalanced operands if not too
     * small: return the result in {dst, sa+sb} */

    /* Note that dst aliasing a or b works with fft. */

    if (sc >= GF2X_MUL_FFT_THRESHOLD && K && K != 1) {
        gf2x_mul_fft(dst, a, sa, b, sb, K);
        goto end_of_gf2x_mul_r;
    }
#endif

    unsigned int sp1, sp2, sp3, sp;
    sp1 = gf2x_toomspace(sa);	// Space for balanced TC routines
    sp2 = gf2x_toomuspace(2 * sa);	// Space for unbalanced TC routines
    sp3 = 2 * sa + gf2x_toomspace(sa);	// Space for unbalanced TC routines w/ lazy cut

    sp = sp1;
    if (sp < sp2)
        sp = sp2;
    if (sp < sp3)
        sp = sp3;		// Worst-case space required

    if (pool->stk_size < sp) {
        pool->stk = realloc(pool->stk, sp * sizeof(unsigned long));
        pool->stk_size = sp;
    }

    /* None of the alternatives below supports c aliasing a or b */
    if (dst == a || dst == b)
        dst = malloc(sc * sizeof(unsigned long));

    if (sa == sb) {
        // Avoid copy in common case
        gf2x_mul_toom(dst, a, b, sa, pool->stk);
    } else if ((sa == (sb + 1) / 2) && gf2x_best_utoom(sb)) {
        // Another common case
        // due to GCD algorithm
        gf2x_mul_tc3u(dst, b, sb, a, pool->stk);
    } else {
        unsigned long *v = pool->stk + gf2x_toomspace(sa);

        unsigned int i, j;

        memset(dst, 0, sc * sizeof(unsigned long));

        unsigned long * ptr = dst;

        for (;;) {
            if (sa == 0)
                break;

            if (sa == 1) {
                ptr[sb] ^= gf2x_addmul_1_n(ptr, ptr, b, sb, a[0]);
                break;
            }

            // finally: the general case
            for (i = 0; i + sa <= sb; i += sa) {
                // Generic (balanced) Toom-Cook mult.
                gf2x_mul_toom(v, a, b + i, sa, pool->stk);
                for (j = 0; j < 2 * sa; j++)
                    ptr[i + j] ^= v[j];
            }

            {
                const unsigned long *t;
                unsigned int st;

                /* Swap a and b, and go for the next spin */
                t = a;
                st = sa;
                a = b + i;
                sa = sb - i;
                b = t;
                sb = st;
            }
            ptr = ptr + i;
        }
    }
end_of_gf2x_mul_r:
    if (dst && dst != c) {
        /* Then we have allocated a temp buffer */
        memcpy(c, dst, sc * sizeof(unsigned long));
        free(dst);
    }
}

/* vim: set sw=4 sta et: */