123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604 |
- /*
- * File: dclxvi-20130329/fp2e.c
- * Author: Ruben Niederhagen, Peter Schwabe
- * Public Domain
- */
- #include <stdio.h>
- #include <assert.h>
- #include <math.h>
- #include "cpucycles.h"
- #ifdef NEW_PARAMETERS
- #include "scalar_512.h"
- #else
- #include "scalar.h"
- #endif
- #include "mul.h"
- extern "C" {
- #include "fpe.h"
- #include "fp2e.h"
- }
- extern const double bn_v;
- extern const double bn_v6;
- extern const scalar_t bn_pminus2;
-
- #ifdef N_OPS
- unsigned long long mulfp2ctr;
- unsigned long long mulfp2fpctr;
- unsigned long long sqfp2ctr;
- unsigned long long invfp2ctr;
- unsigned long long double2fp2ctr;
- unsigned long long doublefp2ctr;
- unsigned long long triple2fp2ctr;
- unsigned long long triplefp2ctr;
- unsigned long long mul_scalarfp2ctr;
- unsigned long long add2fp2ctr;
- unsigned long long addfp2ctr;
- unsigned long long sub2fp2ctr;
- unsigned long long subfp2ctr;
- unsigned long long neg2fp2ctr;
- unsigned long long negfp2ctr;
- unsigned long long mulxifp2ctr;
- unsigned long long conjugatefp2ctr;
- unsigned long long short_coeffredfp2ctr;
- #endif
- #ifndef QHASM
- void fp2e_short_coeffred_c(fp2e_t rop)
- {
- #ifdef N_OPS
- short_coeffredfp2ctr++;
- #endif
- mydouble carry11 = round(rop->v[22]/bn_v);
- mydouble carry11b = round(rop->v[23]/bn_v);
- rop->v[22] = remround(rop->v[22],bn_v);
- rop->v[23] = remround(rop->v[23],bn_v);
- rop->v[0] = rop->v[0] - carry11;
- rop->v[1] = rop->v[1] - carry11b;
- rop->v[6] = rop->v[6] - carry11;
- rop->v[7] = rop->v[7] - carry11b;
- rop->v[12] = rop->v[12] - 4*carry11;
- rop->v[13] = rop->v[13] - 4*carry11b;
- rop->v[18] = rop->v[18] - carry11;
- rop->v[19] = rop->v[19] - carry11b;
- mydouble carry1 = round(rop->v[2]/bn_v);
- mydouble carry1b = round(rop->v[3]/bn_v);
- rop->v[2] = remround(rop->v[2],bn_v);
- rop->v[3] = remround(rop->v[3],bn_v);
- rop->v[4] += carry1;
- rop->v[5] += carry1b;
-
- mydouble carry3 = round(rop->v[6]/bn_v);
- mydouble carry3b = round(rop->v[7]/bn_v);
- rop->v[6] = remround(rop->v[6],bn_v);
- rop->v[7] = remround(rop->v[7],bn_v);
- rop->v[8] += carry3;
- rop->v[9] += carry3b;
- mydouble carry5 = round(rop->v[10]/bn_v);
- mydouble carry5b = round(rop->v[11]/bn_v);
- rop->v[10] = remround(rop->v[10],bn_v);
- rop->v[11] = remround(rop->v[11],bn_v);
- rop->v[12] += carry5;
- rop->v[13] += carry5b;
- mydouble carry7 = round(rop->v[14]/bn_v);
- mydouble carry7b = round(rop->v[15]/bn_v);
- rop->v[14] = remround(rop->v[14],bn_v);
- rop->v[15] = remround(rop->v[15],bn_v);
- rop->v[16] += carry7;
- rop->v[17] += carry7b;
- mydouble carry9 = round(rop->v[18]/bn_v);
- mydouble carry9b = round(rop->v[19]/bn_v);
- rop->v[18] = remround(rop->v[18],bn_v);
- rop->v[19] = remround(rop->v[19],bn_v);
- rop->v[20] += carry9;
- rop->v[21] += carry9b;
- mydouble carry0 = round(rop->v[0]/bn_v6);
- mydouble carry0b = round(rop->v[1]/bn_v6);
- rop->v[0] = remround(rop->v[0],bn_v6);
- rop->v[1] = remround(rop->v[1],bn_v6);
- rop->v[2] += carry0;
- rop->v[3] += carry0b;
- mydouble carry2 = round(rop->v[4]/bn_v);
- mydouble carry2b = round(rop->v[5]/bn_v);
- rop->v[4] = remround(rop->v[4],bn_v);
- rop->v[5] = remround(rop->v[5],bn_v);
- rop->v[6] += carry2;
- rop->v[7] += carry2b;
- mydouble carry4 = round(rop->v[8]/bn_v);
- mydouble carry4b = round(rop->v[9]/bn_v);
- rop->v[8] = remround(rop->v[8],bn_v);
- rop->v[9] = remround(rop->v[9],bn_v);
- rop->v[10] += carry4;
- rop->v[11] += carry4b;
- mydouble carry6 = round(rop->v[12]/bn_v6);
- mydouble carry6b = round(rop->v[13]/bn_v6);
- rop->v[12] = remround(rop->v[12],bn_v6);
- rop->v[13] = remround(rop->v[13],bn_v6);
- rop->v[14] += carry6;
- rop->v[15] += carry6b;
- mydouble carry8 = round(rop->v[16]/bn_v);
- mydouble carry8b = round(rop->v[17]/bn_v);
- rop->v[16] = remround(rop->v[16],bn_v);
- rop->v[17] = remround(rop->v[17],bn_v);
- rop->v[18] += carry8;
- rop->v[19] += carry8b;
- mydouble carry10 = round(rop->v[20]/bn_v);
- mydouble carry10b = round(rop->v[21]/bn_v);
- rop->v[20] = remround(rop->v[20],bn_v);
- rop->v[21] = remround(rop->v[21],bn_v);
- rop->v[22] += carry10;
- rop->v[23] += carry10b;
- }
- #endif
- void fp2e_to_2fpe(fpe_t ropa, fpe_t ropb, const fp2e_t op)//ai+b
- {
- int i;
- for(i=0;i<12;i++)
- {
- ropb->v[i] = op->v[2*i];
- ropa->v[i] = op->v[2*i+1];
- }
- }
- void _2fpe_to_fp2e(fp2e_t rop, const fpe_t opa, const fpe_t opb)//ai+b
- {
- int i;
- for(i=0;i<12;i++)
- {
- rop->v[2*i] = opb->v[i];
- rop->v[2*i+1] = opa->v[i];
- }
- }
- // Set fp2e_t rop to given value:
- void fp2e_set(fp2e_t rop, const fp2e_t op)
- {
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = op->v[i];
- }
- /* Communicate the fact that the fp2e is reduced (and that we don't know anything more about it) */
- void fp2e_isreduced(fp2e_t rop)
- {
- setmax(rop->v[0],(long)bn_v6/2);
- setmax(rop->v[1],(long)bn_v6/2);
- setmax(rop->v[12],(long)bn_v6/2);
- setmax(rop->v[13],(long)bn_v6/2);
- setmax(rop->v[2],(long)bn_v/2);
- setmax(rop->v[3],(long)bn_v/2);
- setmax(rop->v[6],(long)bn_v/2);
- setmax(rop->v[7],(long)bn_v/2);
- setmax(rop->v[8],(long)bn_v/2);
- setmax(rop->v[9],(long)bn_v/2);
- setmax(rop->v[14],(long)bn_v/2);
- setmax(rop->v[15],(long)bn_v/2);
- setmax(rop->v[18],(long)bn_v/2);
- setmax(rop->v[19],(long)bn_v/2);
- setmax(rop->v[20],(long)bn_v/2);
- setmax(rop->v[21],(long)bn_v/2);
- //XXX: Change additive constant
- setmax(rop->v[4],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[5],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[10],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[11],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[16],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[17],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[22],(long)bn_v/2+2331); /* TODO: Think about value */
- setmax(rop->v[23],(long)bn_v/2+2331); /* TODO: Think about value */
- }
- // Set fp2e_t rop to given value:
- void fp2e_set_fpe(fp2e_t rop, const fpe_t op)
- {
- int i;
- for(i=0;i<12;i++)
- {
- rop->v[2*i] = op->v[i];
- rop->v[2*i+1] = 0;
- }
- }
- // Set rop to one
- void fp2e_setone(fp2e_t rop)
- {
- int i;
- for(i=1;i<24;i++)
- rop->v[i] = 0;
- rop->v[0] = 1.;
- }
- // Set rop to zero
- void fp2e_setzero(fp2e_t rop)
- {
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = 0;
- }
- // Compare for equality:
- int fp2e_iseq(const fp2e_t op1, const fp2e_t op2)
- {
- fpe_t a1,b1,a2,b2;
- fp2e_to_2fpe(a1,b1,op1);
- fp2e_to_2fpe(a2,b2,op2);
- return fpe_iseq(a1,a2) && fpe_iseq(b1,b2);
- }
- int fp2e_isone(const fp2e_t op)
- {
- fpe_t ta, tb;
- fp2e_to_2fpe(ta, tb, op);
- int ret = fpe_iszero(ta);
- ret = ret && fpe_isone(tb);
- return ret;
- }
- int fp2e_iszero(const fp2e_t op)
- {
- fpe_t ta, tb;
- fp2e_to_2fpe(ta, tb, op);
- int ret = fpe_iszero(ta);
- ret = ret && fpe_iszero(tb);
- return ret;
- }
- void fp2e_cmov(fp2e_t rop, const fp2e_t op, int c)
- {
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = (1-c)*rop->v[i] + c*op->v[i];
- }
- #ifndef QHASM
- // Double an fp2e:
- void fp2e_double2_c(fp2e_t rop)
- {
- #ifdef N_OPS
- double2fp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = 2*rop->v[i];
- }
- #endif
- #ifndef QHASM
- // Double an fp2e:
- void fp2e_double_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- doublefp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = 2*op->v[i];
- }
- #endif
- #ifndef QHASM
- // Triple an fp2e:
- void fp2e_triple2_c(fp2e_t rop)
- {
- #ifdef N_OPS
- triple2fp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = 3*rop->v[i];
- }
- #endif
- #ifndef QHASM
- // Triple an fp2e:
- void fp2e_triple_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- triplefp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = 3*op->v[i];
- }
- #endif
- void fp2e_mul_scalar(fp2e_t rop, const fp2e_t op, const int s)
- {
- #ifdef N_OPS
- mul_scalarfp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = s*op->v[i];
- }
- // Add two fp2e, store result in op1:
- #ifndef QHASM
- void fp2e_add2_c(fp2e_t op1, const fp2e_t op2)
- {
- #ifdef N_OPS
- add2fp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- op1->v[i] += op2->v[i];
- }
- #endif
- #ifndef QHASM
- // Add two fp2e, store result in rop:
- void fp2e_add_c(fp2e_t rop, const fp2e_t op1, const fp2e_t op2)
- {
- #ifdef N_OPS
- addfp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = op1->v[i] + op2->v[i];
- }
- #endif
- #ifndef QHASM
- // Subtract op2 from op1, store result in op1:
- void fp2e_sub2_c(fp2e_t op1, const fp2e_t op2)
- {
- #ifdef N_OPS
- sub2fp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- op1->v[i] -= op2->v[i];
- }
- #endif
- #ifndef QHASM
- // Subtract op2 from op1, store result in rop:
- void fp2e_sub_c(fp2e_t rop, const fp2e_t op1, const fp2e_t op2)
- {
- #ifdef N_OPS
- subfp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = op1->v[i] - op2->v[i];
- }
- #endif
- #ifndef QHASM
- // Negate op
- void fp2e_neg2_c(fp2e_t op)
- {
- #ifdef N_OPS
- neg2fp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- op->v[i] = -op->v[i];
- }
- #endif
- #ifndef QHASM
- // Negate op
- void fp2e_neg_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- negfp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i++)
- rop->v[i] = -op->v[i];
- }
- #endif
- #ifndef QHASM
- // Conjugates: aX+b to -aX+b
- void fp2e_conjugate_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- conjugatefp2ctr++;
- #endif
- int i;
- for(i=0;i<24;i+=2)
- {
- rop->v[i] = op->v[i];
- rop->v[i+1] = op->v[i+1] * (-1);
- }
- }
- #endif
- #ifndef QHASM
- // Multiply two fp2e, store result in rop:
- void fp2e_mul_c(fp2e_t rop, const fp2e_t op1, const fp2e_t op2)
- {
- #ifdef N_OPS
- mulfp2ctr += 1;
- #endif
- fpe_t a1, b1, a2, b2, r1, r2;
- mydouble a3[24], b3[24];
- int i;
- mydouble t0[24], t1[24], t2[24], t3[24];
- fp2e_to_2fpe(a1, b1, op1);
- fp2e_to_2fpe(a2, b2, op2);
- polymul(t1, a1->v, b2->v); // t1 = a1*b2
- polymul(t2, b1->v, a2->v); // t2 = b1*a2
- for(i=0; i<12; i++) // t3 = 1*a1
- {
- t3[i] = 1*a1->v[i];
- }
- polymul(t3, t3, a2->v); // t3 = 1*a1*a2
- polymul(t0, b1->v, b2->v); // t0 = b1*b2
- for(i=0; i<23; i++)
- {
- a3[i] = t1[i] + t2[i]; // a3 = a1*b2 + b1*a2
- b3[i] = t0[i] - t3[i]; // b3 = b1*b2 - 1*a1*a2
- }
- degred(a3);
- degred(b3);
- coeffred_round_par(a3);
- coeffred_round_par(b3);
- fpe_set_doublearray(r1, a3);
- fpe_set_doublearray(r2, b3);
- _2fpe_to_fp2e(rop, r1, r2);
- }
- #endif
- #ifndef QHASM
- // Square an fp2e, store result in rop:
- void fp2e_square_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- sqfp2ctr += 1;
- #endif
- fpe_t a1, b1, r1, r2;
- mydouble ropa[24], ropb[24];
- fp2e_to_2fpe(a1, b1, op);
- int i;
- /* CheckDoubles are not smart enough to recognize
- * binomial formula to compute b^2-a^2 */
- #ifdef CHECK
- mydouble d1[24];
- polymul(d1, a1->v, a1->v);
- polymul(ropb, b1->v, b1->v);
- polymul(ropa, b1->v, a1->v);
- for(i=0;i<23;i++)
- {
- ropb[i] -= d1[i];
- ropa[i] *= 2;
- }
- #else
- fpe_t t1, t2, t3;
- for(i=0;i<12;i++)
- {
- t1->v[i] = a1->v[i] + b1->v[i];
- t2->v[i] = b1->v[i] - a1->v[i];
- t3->v[i] = 2*b1->v[i];
- }
- polymul(ropa, a1->v, t3->v);
- polymul(ropb, t1->v, t2->v);
- #endif
- degred(ropa);
- degred(ropb);
- coeffred_round_par(ropa);
- coeffred_round_par(ropb);
- fpe_set_doublearray(r1, ropa);
- fpe_set_doublearray(r2, ropb);
- _2fpe_to_fp2e(rop, r1, r2);
- }
- #endif
- #ifndef NEW_PARAMETERS
- #ifndef QHASM
- // Multiply by xi=i+3 which is used to construct F_p^6
- // (a*i + b)*(i + 3) = (3*b - 1*a) + (3*a + b)*i
- void fp2e_mulxi_c(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- mulxifp2ctr++;
- #endif
- fpe_t a, b, t1, t2, t3, t4, t5;
- fp2e_to_2fpe(a, b, op);
- int i;
- for(i=0; i<12; i++)
- {
- t1->v[i] = 3*a->v[i]; // t1 = 3*a
- t2->v[i] = 3*b->v[i]; // t2 = 3*b
- t3->v[i] = 1*a->v[i]; // t3 = 1*a
- }
- fpe_add(t4, t1, b); // t4 = 3*a + b
- fpe_sub(t5, t2, t3); // t5 = 3*b - 1*a
- _2fpe_to_fp2e(rop, t4, t5);
- }
- #endif
- #endif
- // Scalar multiple of an fp2e, store result in rop:
- #ifndef QHASM
- void fp2e_mul_fpe_c(fp2e_t rop, const fp2e_t op1, const fpe_t op2)
- {
- #ifdef N_OPS
- mulfp2fpctr += 1;
- #endif
- fpe_t a1,b1;
- fp2e_to_2fpe(a1,b1,op1);
- fpe_mul(a1,a1,op2);
- fpe_mul(b1,b1,op2);
- _2fpe_to_fp2e(rop,a1,b1);
- }
- #endif
- #ifndef QHASM
- /* computes (op1->m_a*op2->m_a, op1->m_b*op2->m_b) */
- void fp2e_parallel_coeffmul_c(fp2e_t rop, const fp2e_t op1, const fp2e_t op2)
- {
- fpe_t a1, b1, a2, b2; // Needed for intermediary results
- fp2e_to_2fpe(a1,b1,op1);
- fp2e_to_2fpe(a2,b2,op2);
- fpe_mul(a1, a1, a2);
- fpe_mul(b1, b1, b2);
- _2fpe_to_fp2e(rop, a1, b1);
- }
- #endif
- // Inverse multiple of an fp2e, store result in rop:
- void fp2e_invert(fp2e_t rop, const fp2e_t op)
- {
- #ifdef N_OPS
- invfp2ctr += 1;
- #endif
- /* New version */
- fp2e_t d1, d2;
- int i;
- fp2e_parallel_coeffmul(d1, op, op);
- for(i=0;i<24;i+=2)
- d1->v[i] = d1->v[i+1] = d1->v[i] + d1->v[i+1];
- fp2e_short_coeffred(d1);
- for(i=0;i<24;i+=2)
- {
- d2->v[i] = op->v[i];
- d2->v[i+1] = -op->v[i+1];
- }
- fp2e_set(rop,d1);
- for(i = 254; i >= 0; i--)
- {
- fp2e_parallel_coeffmul(rop,rop,rop);
- if(scalar_getbit(bn_pminus2, i))
- fp2e_parallel_coeffmul(rop,rop,d1);
- }
- fp2e_parallel_coeffmul(rop,rop,d2);
- }
- // Print the fp2e:
- void fp2e_print(FILE *outfile, const fp2e_t op)
- {
- fpe_t a,b;
- fp2e_to_2fpe(a,b,op);
- fpe_print(outfile, a);
- fprintf(outfile," * X + ");
- fpe_print(outfile, b);
- }
|