/*

Copyright (c) 2019, Ben Buhrow

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this

list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,

this list of conditions and the following disclaimer in the documentation

and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND

ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED

WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE

DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR

ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES

(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;

LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND

ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT

(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those

of the authors and should not be interpreted as representing official policies,

either expressed or implied, of the FreeBSD Project.

*/

/* Elliptic Curve Method: toplevel and stage 1 routines.

Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011,

2012, 2016 Paul Zimmermann, Alexander Kruppa, Cyril Bouvier, David Cleaver.

This file is part of the ECM Library.

The ECM Library is free software; you can redistribute it and/or modify

it under the terms of the GNU Lesser General Public License as published by

the Free Software Foundation; either version 3 of the License, or (at your

option) any later version.

The ECM Library 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 Lesser General Public

License for more details.

You should have received a copy of the GNU Lesser General Public License

along with the ECM Library; see the file COPYING.LIB. If not, see

http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,

51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#include "avx_ecm.h"

#include "omp.h"

#include "eratosthenes/soe.h"

//#define DEBUG 1

#define DO_STAGE2_INV

// local functions

void vec_add(monty *mdata, ecm_work *work, ecm_pt *Pin, ecm_pt *Pout);

void vec_duplicate(monty *mdata, ecm_work *work, bignum *insum, bignum *indiff, ecm_pt *P);

void next_pt_vec(monty *mdata, ecm_work *work, ecm_pt *P, uint64_t c);

void euclid(monty *mdata, ecm_work *work, ecm_pt *P, uint64_t c);

void prac(monty *mdata, ecm_work *work, ecm_pt *P, uint64_t c);

int check_factor(mpz_t Z, mpz_t n, mpz_t f);

void build_one_curve(thread_data_t *tdata, mpz_t X, mpz_t Z, mpz_t A, uint64_t sigma);

void ecm_stage1(monty *mdata, ecm_work *work, ecm_pt *P, base_t b1, base_t *primes, int verbose);

void ecm_stage2_init_ref(ecm_pt* P, monty* mdata, ecm_work* work, int verbose);

void ecm_stage2_pair_ref(ecm_pt* P, monty* mdata, ecm_work* work, int verbose);

int ecm_stage2_init(ecm_pt* P, monty* mdata, ecm_work* work, int verbose);

void ecm_stage2_pair(uint32_t pairmap_steps, uint32_t* pm_v, uint32_t* pm_u,

ecm_pt* P, monty* mdata, ecm_work* work, int verbose);

uint32_t pair(uint32_t* pairmap_v, uint32_t* pairmap_u,

ecm_work* work, Queue_t** Q, uint32_t* Qrmap, uint32_t* Qmap,

uint64_t* primes, uint64_t B1, uint64_t B2, int verbose);

#include "threadpool.h"

void ecm_sync(void *vptr);

void ecm_dispatch(void *vptr);

void ecm_build_curve_work_fcn(void *vptr);

void ecm_stage1_work_fcn(void *vptr);

void ecm_stage2_work_fcn(void *vptr);

void ecm_sync(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *udata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

if (udata[tid].ecm_phase == 0)

{

udata[tid].phase_done = 1;

}

else if (udata->ecm_phase == 1)

{

udata[tid].phase_done = 1;

}

else if (udata->ecm_phase == 2)

{

udata[tid].phase_done = 1;

}

else if (udata->ecm_phase == 3)

{

udata[tid].phase_done = 1;

}

return;

}

void ecm_dispatch(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *udata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

if (udata[tid].ecm_phase == 0)

{

if (udata[tid].phase_done == 0)

{

tpdata->work_fcn_id = 0;

}

else

{

tpdata->work_fcn_id = tpdata->num_work_fcn;

}

}

else if (udata[tid].ecm_phase == 1)

{

if (udata[tid].phase_done == 0)

{

tpdata->work_fcn_id = 1;

}

else

{

tpdata->work_fcn_id = tpdata->num_work_fcn;

}

}

else if (udata[tid].ecm_phase == 2)

{

if (udata[tid].phase_done == 0)

{

tpdata->work_fcn_id = 2;

}

else

{

tpdata->work_fcn_id = tpdata->num_work_fcn;

}

}

else if (udata[tid].ecm_phase == 3)

{

if (udata[tid].phase_done == 0)

{

tpdata->work_fcn_id = 3;

}

else

{

tpdata->work_fcn_id = tpdata->num_work_fcn;

}

}

return;

}

void ecm_stage1_work_fcn(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *udata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

ecm_stage1(udata[tid].mdata, udata[tid].work, udata[tid].P, udata[tid].b1, NULL, tid == 0);

return;

}

void ecm_stage2_init_work_fcn(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *udata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

ecm_stage2_init(udata[tid].P, udata[tid].mdata, udata[tid].work, tid == 0);

return;

}

void ecm_stage2_work_fcn(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *udata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

ecm_stage2_pair(udata[tid].pairmap_steps, udata[tid].pairmap_v, udata[tid].pairmap_u,

udata[tid].P, udata[tid].mdata, udata[tid].work, tid == 0);

return;

}

void ecm_build_curve_work_fcn(void *vptr)

{

tpool_t *tpdata = (tpool_t *)vptr;

thread_data_t *tdata = (thread_data_t *)tpdata->user_data;

uint32_t tid = tpdata->tindex;

mpz_t X, Z, A; // , t, n, one, r;

int i;

mpz_init(X);

mpz_init(Z);

mpz_init(A);

for (i = 0; i < VECLEN; i++)

{

build_one_curve(&tdata[tid], X, Z, A, tdata[tid].sigma[i]);

insert_mpz_to_vec(tdata[tid].P->X, X, i);

insert_mpz_to_vec(tdata[tid].P->Z, Z, i);

insert_mpz_to_vec(tdata[tid].work->s, A, i);

tdata[tid].sigma[i] = tdata[tid].work->sigma;

}

//print_vechexbignum(tdata[tid].P->X, "POINT x: ");

//print_vechexbignum(tdata[tid].P->Z, "POINT z: ");

//print_vechexbignum(tdata[tid].work->s, "POINT a: ");

tdata[tid].work->diff1->size = tdata[tid].work->n->size;

tdata[tid].work->diff2->size = tdata[tid].work->n->size;

tdata[tid].work->sum1->size = tdata[tid].work->n->size;

tdata[tid].work->sum2->size = tdata[tid].work->n->size;

tdata[tid].work->pt1.X->size = tdata[tid].work->n->size;

tdata[tid].work->pt1.Z->size = tdata[tid].work->n->size;

tdata[tid].work->pt2.X->size = tdata[tid].work->n->size;

tdata[tid].work->pt2.Z->size = tdata[tid].work->n->size;

tdata[tid].work->tt1->size = tdata[tid].work->n->size;

tdata[tid].work->tt2->size = tdata[tid].work->n->size;

tdata[tid].work->tt3->size = tdata[tid].work->n->size;

tdata[tid].work->tt4->size = tdata[tid].work->n->size;

mpz_clear(X);

mpz_clear(Z);

mpz_clear(A);

return;

}

void ecm_work_init(ecm_work *work)

{

int i, j, m;

uint32_t U = work->U;

uint32_t L = work->L;

uint32_t D = work->D;

uint32_t R = work->R;

work->ptadds = 0;

work->ptdups = 0;

work->numinv = 0;

work->numprimes = 0;

work->paired = 0;

work->diff1 = vecInit();

work->diff2 = vecInit();

work->sum1 = vecInit();

work->sum2 = vecInit();

work->tt1 = vecInit();

work->tt2 = vecInit();

work->tt3 = vecInit();

work->tt4 = vecInit();

work->tt5 = vecInit();

work->s = vecInit();

work->n = vecInit();

work->Paprod = (bignum **)malloc((2 * L) * sizeof(bignum *));

work->Pa_inv = (bignum * *)malloc((2 * L) * sizeof(bignum*));

work->Pa = (ecm_pt *)malloc((2 * L) * sizeof(ecm_pt));

for (j = 0; j < (2 * L); j++)

{

work->Paprod[j] = vecInit();

work->Pa_inv[j] = vecInit();

ecm_pt_init(&work->Pa[j]);

}

work->Pad = (ecm_pt *)malloc(sizeof(ecm_pt));

ecm_pt_init(work->Pad);

work->Pdnorm = (ecm_pt*)malloc(sizeof(ecm_pt));

ecm_pt_init(work->Pdnorm);

// build an array to hold values of f(b).

// will need to be of size U*R, allowed values will be

// up to a multiple of U of the residues mod D

work->Pb = (ecm_pt *)malloc(U * (R + 1) * sizeof(ecm_pt));

work->Pbprod = (bignum **)malloc(U * (R + 1) * sizeof(bignum *));

for (j = 0; j < U * (R + 1); j++)

{

ecm_pt_init(&work->Pb[j]);

work->Pbprod[j] = vecInit();

}

work->map = (uint32_t *)calloc(U * (D + 1) + 3, sizeof(uint32_t));

work->map[0] = 0;

work->map[1] = 1;

work->map[2] = 2;

m = 3;

for (i = 0; i < U; i++)

{

if (i == 0)

j = 3;

else

j = 1;

for (; j < D; j++)

{

if (spGCD(j, D) == 1)

{

work->map[i*D + j] = m++;

}

else

{

work->map[i*D + j] = 0;

}

}

if (i == 0)

work->map[i*D + j] = m++;

}

ecm_pt_init(&work->pt1);

ecm_pt_init(&work->pt2);

ecm_pt_init(&work->pt3);

ecm_pt_init(&work->pt4);

ecm_pt_init(&work->pt5);

work->stg2acc = vecInit();

return;

}

void ecm_work_free(ecm_work *work)

{

int i, j;

uint32_t U = work->U;

uint32_t L = work->L;

uint32_t D = work->D;

uint32_t R = work->R;

vecFree(work->diff1);

vecFree(work->diff2);

vecFree(work->sum1);

vecFree(work->sum2);

vecFree(work->tt1);

vecFree(work->tt2);

vecFree(work->tt3);

vecFree(work->tt4);

vecFree(work->tt5);

vecFree(work->s);

vecFree(work->n);

ecm_pt_free(&work->pt1);

ecm_pt_free(&work->pt2);

ecm_pt_free(&work->pt3);

ecm_pt_free(&work->pt4);

ecm_pt_free(&work->pt5);

for (i = 0; i < (2 * L); i++)

{

ecm_pt_free(&work->Pa[i]);

vecFree(work->Paprod[i]);

vecFree(work->Pa_inv[i]);

}

ecm_pt_free(work->Pad);

ecm_pt_free(work->Pdnorm);

free(work->Pa);

free(work->Pad);

free(work->Pdnorm);

free(work->map);

for (i = 0; i < U * (R + 1); i++)

{

ecm_pt_free(&work->Pb[i]);

vecFree(work->Pbprod[i]);

}

free(work->Pbprod);

free(work->Paprod);

free(work->Pb);

vecFree(work->stg2acc);

return;

}

void ecm_pt_init(ecm_pt *pt)

{

pt->X = vecInit();

pt->Z = vecInit();

}

void ecm_pt_free(ecm_pt *pt)

{

vecFree(pt->X);

vecFree(pt->Z);

}

void vec_add(monty *mdata, ecm_work *work, ecm_pt *Pin, ecm_pt *Pout)

{

// compute:

//x+ = z- * [(x1-z1)(x2+z2) + (x1+z1)(x2-z2)]^2

//z+ = x- * [(x1-z1)(x2+z2) - (x1+z1)(x2-z2)]^2

// where:

//x- = original x

//z- = original z

// given the sums and differences of the original points (stored in work structure).

vecmulmod_ptr(work->diff1, work->sum2, work->tt1, work->n, work->tt4, mdata); //U

vecmulmod_ptr(work->sum1, work->diff2, work->tt2, work->n, work->tt4, mdata); //V

vecaddsubmod_ptr(work->tt1, work->tt2, work->tt3, work->tt4, mdata);

vecsqrmod_ptr(work->tt3, work->tt1, work->n, work->tt5, mdata); //(U + V)^2

vecsqrmod_ptr(work->tt4, work->tt2, work->n, work->tt5, mdata); //(U - V)^2

// choosing the initial point Pz0 = 1 means that z_p-q = 1 and this mul isn't necessary...

// but that involves a different way to initialize curves, so for now

// we can't assume Z=1

if (Pin->X == Pout->X)

{

base_t *swap;

vecmulmod_ptr(work->tt1, Pin->Z, Pout->Z, work->n, work->tt4, mdata); //Z * (U + V)^2

vecmulmod_ptr(work->tt2, Pin->X, Pout->X, work->n, work->tt4, mdata); //x * (U - V)^2

swap = Pout->Z->data;

Pout->Z->data = Pout->X->data;

Pout->X->data = swap;

}

else

{

vecmulmod_ptr(work->tt1, Pin->Z, Pout->X, work->n, work->tt4, mdata); //Z * (U + V)^2

vecmulmod_ptr(work->tt2, Pin->X, Pout->Z, work->n, work->tt4, mdata); //x * (U - V)^2

}

work->ptadds++;

return;

}

void vec_duplicate(monty *mdata, ecm_work *work, bignum *insum, bignum *indiff, ecm_pt *P)

{

vecsqrmod_ptr(indiff, work->tt1, work->n, work->tt4, mdata); // V=(x1 - z1)^2

vecsqrmod_ptr(insum, work->tt2, work->n, work->tt4, mdata); // U=(x1 + z1)^2

vecmulmod_ptr(work->tt1, work->tt2, P->X, work->n, work->tt4, mdata); // x=U*V

vecsubmod_ptr(work->tt2, work->tt1, work->tt3, mdata); // w = U-V

vecmulmod_ptr(work->tt3, work->s, work->tt2, work->n, work->tt4, mdata); // t = (A+2)/4 * w

vecaddmod_ptr(work->tt2, work->tt1, work->tt2, mdata); // t = t + V

vecmulmod_ptr(work->tt2, work->tt3, P->Z, work->n, work->tt4, mdata); // Z = t*w

work->ptdups++;

return;

}

#define ADD 5.5 // counting squares as 0.75 of a mul

#define DUP 4.5 // counting squares as 0.75 of a mul

// PRAC with all of the conditions specified. Prime95 found

// that removing some of them slightly improves performance for

// very large inputs. Testing also revealed a very slight

// improvement for the smaller numbers this routine can deal with.

// To use the original version, define this.

//#define ORIG_PRAC

#define NV 10

/* 1/val[0] = the golden ratio (1+sqrt(5))/2, and 1/val[i] for i>0

is the real number whose continued fraction expansion is all 1s

except for a 2 in i+1-st place */

static double val[NV] =

{ 0.61803398874989485, 0.72360679774997897, 0.58017872829546410,

0.63283980608870629, 0.61242994950949500, 0.62018198080741576,

0.61721461653440386, 0.61834711965622806, 0.61791440652881789,

0.61807966846989581};

static double

lucas_cost(uint64_t n, double v)

{

uint64_t d, e, r;

double c; /* cost */

d = n;

r = (uint64_t)((double)d * v + 0.5);

if (r >= n)

return (ADD * (double)n);

d = n - r;

e = 2 * r - n;

c = DUP + ADD; /* initial duplicate and final addition */

while (d != e)

{

if (d < e)

{

r = d;

d = e;

e = r;

}

#ifdef ORIG_PRAC

if (d - e <= e / 4 && ((d + e) % 3) == 0)

{ /* condition 1 */

d = (2 * d - e) / 3;

e = (e - d) / 2;

c += 3.0 * ADD; /* 3 additions */

}

else if (d - e <= e / 4 && (d - e) % 6 == 0)

{ /* condition 2 */

d = (d - e) / 2;

c += ADD + DUP; /* one addition, one duplicate */

}

else if ((d + 3) / 4 <= e)

#else

if ((d + 3) / 4 <= e)

#endif

{ /* condition 3 */

d -= e;

c += ADD; /* one addition */

}

else if ((d + e) % 2 == 0)

{ /* condition 4 */

d = (d - e) / 2;

c += ADD + DUP; /* one addition, one duplicate */

}

/* now d+e is odd */

else if (d % 2 == 0)

{ /* condition 5 */

d /= 2;

c += ADD + DUP; /* one addition, one duplicate */

}

#ifdef ORIG_PRAC

/* now d is odd and e is even */

else if (d % 3 == 0)

{ /* condition 6 */

d = d / 3 - e;

c += 3.0 * ADD + DUP; /* three additions, one duplicate */

}

else if ((d + e) % 3 == 0)

{ /* condition 7 */

d = (d - 2 * e) / 3;

c += 3.0 * ADD + DUP; /* three additions, one duplicate */

}

else if ((d - e) % 3 == 0)

{ /* condition 8 */

d = (d - e) / 3;

c += 3.0 * ADD + DUP; /* three additions, one duplicate */

}

#endif

else /* necessarily e is even: catches all cases */

{ /* condition 9 */

e /= 2;

c += ADD + DUP; /* one addition, one duplicate */

}

}

if (d != 1)

{

return 999999999.;

}

return c;

}

void prac(monty *mdata, ecm_work *work, ecm_pt *P, uint64_t c)

{

uint64_t d, e, r;

double cmin, cost;

int i;

bignum *s1, *s2, *d1, *d2;

base_t *sw_x, *sw_z;

/* chooses the best value of v */

for (i = d = 0, cmin = ADD * (double)c; d < NV; d++)

{

cost = lucas_cost(c, val[d]);

if (cost < cmin)

{

cmin = cost;

i = d;

}

}

d = c;

r = (uint64_t)((double)d * val[i] + 0.5);

s1 = work->sum1;

s2 = work->sum2;

d1 = work->diff1;

d2 = work->diff2;

/* first iteration always begins by Condition 3, then a swap */

d = c - r;

e = 2 * r - c;

// mpres_set(xB, xA, n);

// mpres_set(zB, zA, n); /* B=A */

// mpres_set(xC, xA, n);

// mpres_set(zC, zA, n); /* C=A */

// duplicate(xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */

// the first one is always a doubling

// point1 is [1]P

vecCopy(P->X, work->pt1.X);

vecCopy(P->Z, work->pt1.Z);

vecCopy(P->X, work->pt2.X);

vecCopy(P->Z, work->pt2.Z);

vecCopy(P->X, work->pt3.X);

vecCopy(P->Z, work->pt3.Z);

vecsubmod_ptr(work->pt1.X, work->pt1.Z, d1, mdata);

vecaddmod_ptr(work->pt1.X, work->pt1.Z, s1, mdata);

// point2 is [2]P

vec_duplicate(mdata, work, s1, d1, &work->pt1);

while (d != e)

{

if (d < e)

{

r = d;

d = e;

e = r;

//mpres_swap(xA, xB, n);

//mpres_swap(zA, zB, n);

sw_x = work->pt1.X->data;

sw_z = work->pt1.Z->data;

work->pt1.X->data = work->pt2.X->data;

work->pt1.Z->data = work->pt2.Z->data;

work->pt2.X->data = sw_x;

work->pt2.Z->data = sw_z;

}

/* do the first line of Table 4 whose condition qualifies */

#ifdef ORIG_PRAC

if (d - e <= e / 4 && ((d + e) % 3) == 0)

{ /* condition 1 */

d = (2 * d - e) / 3;

e = (e - d) / 2;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt4); // T = A + B (C)

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt2, &work->pt5); // T2 = T + A (B)

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt1, &work->pt2); // B = B + T (A)

//add3(xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */

//add3(xT2, zT2, xT, zT, xA, zA, xB, zB, n, u, v, w); /* T2 = f(T,A,B) */

//add3(xB, zB, xB, zB, xT, zT, xA, zA, n, u, v, w); /* B = f(B,T,A) */

//mpres_swap(xA, xT2, n);

//mpres_swap(zA, zT2, n); /* swap A and T2 */

sw_x = work->pt1.X->data;

sw_z = work->pt1.Z->data;

work->pt1.X->data = work->pt5.X->data;

work->pt1.Z->data = work->pt5.Z->data;

work->pt5.X->data = sw_x;

work->pt5.Z->data = sw_z;

}

else if (d - e <= e / 4 && (d - e) % 6 == 0)

{ /* condition 2 */

d = (d - e) / 2;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt2); // B = A + B (C)

vec_duplicate(mdata, work, s1, d1, &work->pt1); // A = 2A

//add3(xB, zB, xA, zA, xB, zB, xC, zC, n, u, v, w); /* B = f(A,B,C) */

//duplicate(xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */

}

else if ((d + 3) / 4 <= e)

#else

if ((d + 3) / 4 <= e)

#endif

{ /* condition 3 */

d -= e;

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt4); // T = B + A (C)

//add3(xT, zT, xB, zB, xA, zA, xC, zC, n, u, v, w); /* T = f(B,A,C) */

/* circular permutation (B,T,C) */

//tmp = xB;

//xB = xT;

//xT = xC;

//xC = tmp;

//tmp = zB;

//zB = zT;

//zT = zC;

//zC = tmp;

sw_x = work->pt2.X->data;

sw_z = work->pt2.Z->data;

work->pt2.X->data = work->pt4.X->data;

work->pt2.Z->data = work->pt4.Z->data;

work->pt4.X->data = work->pt3.X->data;

work->pt4.Z->data = work->pt3.Z->data;

work->pt3.X->data = sw_x;

work->pt3.Z->data = sw_z;

}

else if ((d + e) % 2 == 0)

{ /* condition 4 */

d = (d - e) / 2;

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt2); // B = B + A (C)

vec_duplicate(mdata, work, s2, d2, &work->pt1); // A = 2A

//add3(xB, zB, xB, zB, xA, zA, xC, zC, n, u, v, w); /* B = f(B,A,C) */

//duplicate(xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */

}

/* now d+e is odd */

else if (d % 2 == 0)

{ /* condition 5 */

d /= 2;

vecaddsubmod_ptr(work->pt3.X, work->pt3.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt2, &work->pt3); // C = C + A (B)

vec_duplicate(mdata, work, s2, d2, &work->pt1); // A = 2A

//add3(xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(C,A,B) */

//duplicate(xA, zA, xA, zA, n, b, u, v, w); /* A = 2*A */

}

/* now d is odd, e is even */

#ifdef ORIG_PRAC

else if (d % 3 == 0)

{ /* condition 6 */

d = d / 3 - e;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vec_duplicate(mdata, work, s1, d1, &work->pt4); // T = 2A

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt5); // T2 = A + B (C)

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt1, &work->pt1); // A = T + A (A)

vecaddsubmod_ptr(work->pt5.X, work->pt5.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt4); // T = T + T2 (C)

//duplicate(xT, zT, xA, zA, n, b, u, v, w); /* T = 2*A */

//add3(xT2, zT2, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T2 = f(A,B,C) */

//add3(xA, zA, xT, zT, xA, zA, xA, zA, n, u, v, w); /* A = f(T,A,A) */

//add3(xT, zT, xT, zT, xT2, zT2, xC, zC, n, u, v, w); /* T = f(T,T2,C) */

/* circular permutation (C,B,T) */

//tmp = xC;

//xC = xB;

//xB = xT;

//xT = tmp;

//tmp = zC;

//zC = zB;

//zB = zT;

//zT = tmp;

sw_x = work->pt3.X->data;

sw_z = work->pt3.Z->data;

work->pt3.X->data = work->pt2.X->data;

work->pt3.Z->data = work->pt2.Z->data;

work->pt2.X->data = work->pt4.X->data;

work->pt2.Z->data = work->pt4.Z->data;

work->pt4.X->data = sw_x;

work->pt4.Z->data = sw_z;

}

else if ((d + e) % 3 == 0)

{ /* condition 7 */

d = (d - 2 * e) / 3;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt4); // T = A + B (C)

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt2, &work->pt2); // B = T + A (B)

vec_duplicate(mdata, work, s2, d2, &work->pt4); // T = 2A

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt1, &work->pt1); // A = A + T (A) = 3A

//add3(xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */

//add3(xB, zB, xT, zT, xA, zA, xB, zB, n, u, v, w); /* B = f(T,A,B) */

//duplicate(xT, zT, xA, zA, n, b, u, v, w);

//add3(xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */

}

else if ((d - e) % 3 == 0)

{ /* condition 8 */

d = (d - e) / 3;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt3, &work->pt4); // T = A + B (C)

vecaddsubmod_ptr(work->pt3.X, work->pt3.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt2, &work->pt3); // C = C + A (B)

//add3(xT, zT, xA, zA, xB, zB, xC, zC, n, u, v, w); /* T = f(A,B,C) */

//add3(xC, zC, xC, zC, xA, zA, xB, zB, n, u, v, w); /* C = f(A,C,B) */

//mpres_swap(xB, xT, n);

//mpres_swap(zB, zT, n); /* swap B and T */

sw_x = work->pt2.X->data;

sw_z = work->pt2.Z->data;

work->pt2.X->data = work->pt4.X->data;

work->pt2.Z->data = work->pt4.Z->data;

work->pt4.X->data = sw_x;

work->pt4.Z->data = sw_z;

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s2, d2, mdata);

vec_duplicate(mdata, work, s2, d2, &work->pt4); // T = 2A

vecaddsubmod_ptr(work->pt1.X, work->pt1.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt4.X, work->pt4.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt1, &work->pt1); // A = A + T (A) = 3A

//duplicate(xT, zT, xA, zA, n, b, u, v, w);

//add3(xA, zA, xA, zA, xT, zT, xA, zA, n, u, v, w); /* A = 3*A */

}

#endif

else /* necessarily e is even here */

{ /* condition 9 */

e /= 2;

vecaddsubmod_ptr(work->pt3.X, work->pt3.Z, s1, d1, mdata);

vecaddsubmod_ptr(work->pt2.X, work->pt2.Z, s2, d2, mdata);

vec_add(mdata, work, &work->pt1, &work->pt3); // C = C + B (A)

vec_duplicate(mdata, work, s2, d2, &work->pt2); // B = 2B

//add3(xC, zC, xC, zC, xB, zB, xA, zA, n, u, v, w); /* C = f(C,B,A) */

//duplicate(xB, zB, xB, zB, n, b, u, v, w); /* B = 2*B */

}

}

vecsubmod_ptr(work->pt1.X, work->pt1.Z, d1, mdata);

vecaddmod_ptr(work->pt1.X, work->pt1.Z, s1, mdata);

vecsubmod_ptr(work->pt2.X, work->pt2.Z, d2, mdata);

vecaddmod_ptr(work->pt2.X, work->pt2.Z, s2, mdata);

vec_add(mdata, work, &work->pt3, P); // A = A + B (C)

//add3(xA, zA, xA, zA, xB, zB, xC, zC, n, u, v, w);

if (d != 1)

{

printf("problem: d != 1\n");

}

return;

}

void next_pt_vec(monty *mdata, ecm_work *work, ecm_pt *P, uint64_t c)

{

uint64_t mask;

bignum *x1, *z1, *x2, *z2, *s1, *s2, *d1, *d2;

uint64_t e, d;

x1 = work->pt1.X;

z1 = work->pt1.Z;

x2 = work->pt2.X;

z2 = work->pt2.Z;

s1 = work->sum1;

s2 = work->sum2;

d1 = work->diff1;

d2 = work->diff2;

if (c == 1)

{

return;

}

vecCopy(P->X, x1);

vecCopy(P->Z, z1);

vecsubmod_ptr(P->X, P->Z, d1, mdata);

vecaddmod_ptr(P->X, P->Z, s1, mdata);

vec_duplicate(mdata, work, s1, d1, &work->pt2);

//mulcnt[tid] += 3;

//sqrcnt[tid] += 2;

if (c == 2)

{

vecCopy(x2, P->X);

vecCopy(z2, P->Z);

return;

}

// find the first '1' bit then skip it

#ifdef __INTEL_COMPILER

mask = 1ULL << (64 - _lzcnt_u64((uint64_t)c) - 2);

#else

mask = 1ULL << (64 - __builtin_clzll((uint64_t)c) - 2);

#endif

//goal is to compute x_c, z_c using montgomery's addition

//and multiplication formula's for x and z

//the procedure will be similar to p+1 lucas chain formula's

//but this time we are simultaneously incrementing both x and z

//rather than just Vn. In each bit of the binary expansion of

//c, we need just one addition and one duplication for both x and z.

//compute loop

//for each bit of M to the right of the most significant bit

d = 1;

e = 2;

while (mask > 0)

{

vecaddsubmod_ptr(x2, z2, s2, d2, mdata);

vecaddsubmod_ptr(x1, z1, s1, d1, mdata);

//if the bit is 1

if (c & mask)

{

//add x1,z1, duplicate x2,z2

vec_add(mdata, work, P, &work->pt1);

vec_duplicate(mdata, work, s2, d2, &work->pt2);

d = d + e;

e *= 2;

}

else

{

//add x2,z2, duplicate x1,z1

vec_add(mdata, work, P, &work->pt2);

vec_duplicate(mdata, work, s1, d1, &work->pt1);

e = e + d;

d *= 2;

}

//mulcnt[tid] += 7;

//sqrcnt[tid] += 4;

mask >>= 1;

}

if (d != c)

{

printf("expected %lu, ladder returned %lu\n", c, d);

exit(1);

}

vecCopy(x1, P->X);

vecCopy(z1, P->Z);

return;

}

void array_mul(uint64_t *primes, int num_p, mpz_t piprimes)

{

mpz_t *p;

int i;

int alloc;

uint64_t stg1 = (uint64_t)STAGE1_MAX;

if (num_p & 1)

{

p = (mpz_t *)malloc((num_p + 1) * sizeof(mpz_t));

alloc = num_p + 1;

}

else

{

p = (mpz_t *)malloc((num_p + 0) * sizeof(mpz_t));

alloc = num_p;

}

for (i = 0; i < num_p; i++)

{

uint64_t c = 1;

uint64_t q;