1/*
2 * Copyright (c) 2008-2020 Stefan Krah. All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright
9 *    notice, this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright
12 *    notice, this list of conditions and the following disclaimer in the
13 *    documentation and/or other materials provided with the distribution.
14 *
15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25 * SUCH DAMAGE.
26 */
27
28
29#include "mpdecimal.h"
30
31#include <assert.h>
32#include <stdlib.h>
33
34#include "bits.h"
35#include "numbertheory.h"
36#include "umodarith.h"
37
38
39/* Bignum: Initialize the Number Theoretic Transform. */
40
41
42/*
43 * Return the nth root of unity in F(p). This corresponds to e**((2*pi*i)/n)
44 * in the Fourier transform. We have w**n == 1 (mod p).
45 *    n := transform length.
46 *    sign := -1 for forward transform, 1 for backward transform.
47 *    modnum := one of {P1, P2, P3}.
48 */
49mpd_uint_t
50_mpd_getkernel(mpd_uint_t n, int sign, int modnum)
51{
52    mpd_uint_t umod, p, r, xi;
53#ifdef PPRO
54    double dmod;
55    uint32_t dinvmod[3];
56#endif
57
58    SETMODULUS(modnum);
59    r = mpd_roots[modnum]; /* primitive root of F(p) */
60    p = umod;
61    xi = (p-1) / n;
62
63    if (sign == -1)
64        return POWMOD(r, (p-1-xi));
65    else
66        return POWMOD(r, xi);
67}
68
69/*
70 * Initialize and return transform parameters.
71 *    n := transform length.
72 *    sign := -1 for forward transform, 1 for backward transform.
73 *    modnum := one of {P1, P2, P3}.
74 */
75struct fnt_params *
76_mpd_init_fnt_params(mpd_size_t n, int sign, int modnum)
77{
78    struct fnt_params *tparams;
79    mpd_uint_t umod;
80#ifdef PPRO
81    double dmod;
82    uint32_t dinvmod[3];
83#endif
84    mpd_uint_t kernel, w;
85    mpd_uint_t i;
86    mpd_size_t nhalf;
87
88    assert(ispower2(n));
89    assert(sign == -1 || sign == 1);
90    assert(P1 <= modnum && modnum <= P3);
91
92    nhalf = n/2;
93    tparams = mpd_sh_alloc(sizeof *tparams, nhalf, sizeof (mpd_uint_t));
94    if (tparams == NULL) {
95        return NULL;
96    }
97
98    SETMODULUS(modnum);
99    kernel = _mpd_getkernel(n, sign, modnum);
100
101    tparams->modnum = modnum;
102    tparams->modulus = umod;
103    tparams->kernel = kernel;
104
105    /* wtable[] := w**0, w**1, ..., w**(nhalf-1) */
106    w = 1;
107    for (i = 0; i < nhalf; i++) {
108        tparams->wtable[i] = w;
109        w = MULMOD(w, kernel);
110    }
111
112    return tparams;
113}
114
115/* Initialize wtable of size three. */
116void
117_mpd_init_w3table(mpd_uint_t w3table[3], int sign, int modnum)
118{
119    mpd_uint_t umod;
120#ifdef PPRO
121    double dmod;
122    uint32_t dinvmod[3];
123#endif
124    mpd_uint_t kernel;
125
126    SETMODULUS(modnum);
127    kernel = _mpd_getkernel(3, sign, modnum);
128
129    w3table[0] = 1;
130    w3table[1] = kernel;
131    w3table[2] = POWMOD(kernel, 2);
132}
133