amcl/
fp12.rs

1/*
2Licensed to the Apache Software Foundation (ASF) under one
3or more contributor license agreements.  See the NOTICE file
4distributed with this work for additional information
5regarding copyright ownership.  The ASF licenses this file
6to you under the Apache License, Version 2.0 (the
7"License"); you may not use this file except in compliance
8with the License.  You may obtain a copy of the License at
9
10  http://www.apache.org/licenses/LICENSE-2.0
11
12Unless required by applicable law or agreed to in writing,
13software distributed under the License is distributed on an
14"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
15KIND, either express or implied.  See the License for the
16specific language governing permissions and limitations
17under the License.
18*/
19
20use super::big;
21use super::fp2::FP2;
22use super::fp4::FP4;
23use super::big::BIG;
24use super::rom;
25use types::SexticTwist;
26use std::str::SplitWhitespace;
27
28#[derive(Copy, Clone)]
29pub struct FP12 {
30    a: FP4,
31    b: FP4,
32    c: FP4,
33}
34
35impl PartialEq for FP12 {
36	fn eq(&self, other: &FP12) -> bool {
37		self.equals(other)
38	}
39}
40
41impl FP12 {
42    pub fn new() -> FP12 {
43        FP12 {
44            a: FP4::new(),
45            b: FP4::new(),
46            c: FP4::new(),
47        }
48    }
49
50    pub fn new_int(a: isize) -> FP12 {
51        let mut f = FP12::new();
52        f.a.copy(&FP4::new_int(a));
53        f.b.zero();
54        f.c.zero();
55        return f;
56    }
57
58    pub fn new_copy(x: &FP12) -> FP12 {
59        let mut f = FP12::new();
60        f.a.copy(&x.a);
61        f.b.copy(&x.b);
62        f.c.copy(&x.c);
63        return f;
64    }
65
66    pub fn new_fp4s(d: &FP4, e: &FP4, f: &FP4) -> FP12 {
67        let mut g = FP12::new();
68        g.a.copy(d);
69        g.b.copy(e);
70        g.c.copy(f);
71        return g;
72    }
73
74    pub fn new_fp4(d: &FP4) -> FP12 {
75        let mut g = FP12::new();
76        g.a.copy(d);
77        g.b.zero();
78        g.c.zero();
79        return g;
80    }
81
82    /* reduce components mod Modulus */
83    pub fn reduce(&mut self) {
84        self.a.reduce();
85        self.b.reduce();
86        self.c.reduce();
87    }
88
89    /* normalise components of w */
90    pub fn norm(&mut self) {
91        self.a.norm();
92        self.b.norm();
93        self.c.norm();
94    }
95
96    /* test self=0 ? */
97    pub fn iszilch(&self) -> bool {
98        //self.reduce();
99        return self.a.iszilch() && self.b.iszilch() && self.c.iszilch();
100    }
101
102    /* Conditional move of g to self dependant on d */
103    pub fn cmove(&mut self, g: &FP12, d: isize) {
104        self.a.cmove(&g.a, d);
105        self.b.cmove(&g.b, d);
106        self.c.cmove(&g.c, d);
107    }
108
109    /* return 1 if b==c, no branching */
110    fn teq(b: i32, c: i32) -> isize {
111        let mut x = b ^ c;
112        x -= 1; // if x=0, x now -1
113        return ((x >> 31) & 1) as isize;
114    }
115
116    /* Constant time select from pre-computed table */
117    pub fn selector(&mut self, g: &[FP12], b: i32) {
118        let m = b >> 31;
119        let mut babs = (b ^ m) - m;
120
121        babs = (babs - 1) / 2;
122
123        self.cmove(&g[0], FP12::teq(babs, 0)); // conditional move
124        self.cmove(&g[1], FP12::teq(babs, 1));
125        self.cmove(&g[2], FP12::teq(babs, 2));
126        self.cmove(&g[3], FP12::teq(babs, 3));
127        self.cmove(&g[4], FP12::teq(babs, 4));
128        self.cmove(&g[5], FP12::teq(babs, 5));
129        self.cmove(&g[6], FP12::teq(babs, 6));
130        self.cmove(&g[7], FP12::teq(babs, 7));
131
132        let mut invf = FP12::new_copy(self);
133        invf.conj();
134        self.cmove(&invf, (m & 1) as isize);
135    }
136
137    /* test self=1 ? */
138    pub fn isunity(&self) -> bool {
139        let one = FP4::new_int(1);
140        return self.a.equals(&one) && self.b.iszilch() && self.c.iszilch();
141    }
142
143    /* test self=x */
144    pub fn equals(&self, x: &FP12) -> bool {
145        return self.a.equals(&x.a) && self.b.equals(&x.b) && self.c.equals(&x.c);
146    }
147
148    pub fn geta(&mut self) -> FP4 {
149        let f = FP4::new_copy(&self.a);
150        return f;
151    }
152
153    pub fn getb(&mut self) -> FP4 {
154        let f = FP4::new_copy(&self.b);
155        return f;
156    }
157
158    pub fn getc(&mut self) -> FP4 {
159        let f = FP4::new_copy(&self.c);
160        return f;
161    }
162
163    /* copy self=x */
164    pub fn copy(&mut self, x: &FP12) {
165        self.a.copy(&x.a);
166        self.b.copy(&x.b);
167        self.c.copy(&x.c);
168    }
169
170    /* set self=1 */
171    pub fn one(&mut self) {
172        self.a.one();
173        self.b.zero();
174        self.c.zero();
175    }
176
177    /* this=conj(this) */
178    pub fn conj(&mut self) {
179        self.a.conj();
180        self.b.nconj();
181        self.c.conj();
182    }
183
184    /* Granger-Scott Unitary Squaring */
185    pub fn usqr(&mut self) {
186        let mut a = FP4::new_copy(&self.a);
187        let mut b = FP4::new_copy(&self.c);
188        let mut c = FP4::new_copy(&self.b);
189        let mut d = FP4::new();
190
191        self.a.sqr();
192        d.copy(&self.a);
193        d.add(&self.a);
194        self.a.add(&d);
195
196        self.a.norm();
197        a.nconj();
198
199        a.dbl();
200        self.a.add(&a);
201        b.sqr();
202        b.times_i();
203
204        d.copy(&b);
205        d.add(&b);
206        b.add(&d);
207        b.norm();
208
209        c.sqr();
210        d.copy(&c);
211        d.add(&c);
212        c.add(&d);
213        c.norm();
214
215        self.b.conj();
216        self.b.dbl();
217        self.c.nconj();
218
219        self.c.dbl();
220        self.b.add(&b);
221        self.c.add(&c);
222        self.reduce();
223    }
224
225    /* Chung-Hasan SQR2 method from http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf */
226    pub fn sqr(&mut self) {
227        let mut a = FP4::new_copy(&self.a);
228        let mut b = FP4::new_copy(&self.b);
229        let mut c = FP4::new_copy(&self.c);
230        let mut d = FP4::new_copy(&self.a);
231
232        a.sqr();
233        b.mul(&self.c);
234        b.dbl();
235        b.norm();
236        c.sqr();
237        d.mul(&self.b);
238        d.dbl();
239
240        self.c.add(&self.a);
241        self.c.add(&self.b);
242        self.c.norm();
243        self.c.sqr();
244
245        self.a.copy(&a);
246        a.add(&b);
247        a.norm();
248        a.add(&c);
249        a.add(&d);
250        a.norm();
251
252        a.neg();
253        b.times_i();
254        c.times_i();
255
256        self.a.add(&b);
257
258        self.b.copy(&c);
259        self.b.add(&d);
260        self.c.add(&a);
261        self.norm();
262    }
263
264    /* FP12 full multiplication self=self*y */
265    pub fn mul(&mut self, y: &FP12) {
266        let mut z0 = FP4::new_copy(&self.a);
267        let mut z1 = FP4::new();
268        let mut z2 = FP4::new_copy(&mut self.b);
269        let mut z3 = FP4::new();
270        let mut t0 = FP4::new_copy(&self.a);
271        let mut t1 = FP4::new_copy(&y.a);
272
273        z0.mul(&y.a);
274        z2.mul(&y.b);
275
276        t0.add(&self.b);
277        t1.add(&y.b);
278
279        t0.norm();
280        t1.norm();
281
282        z1.copy(&t0);
283        z1.mul(&t1);
284        t0.copy(&self.b);
285        t0.add(&self.c);
286        t1.copy(&y.b);
287        t1.add(&y.c);
288
289        t0.norm();
290        t1.norm();
291
292        z3.copy(&t0);
293        z3.mul(&t1);
294
295        t0.copy(&z0);
296        t0.neg();
297        t1.copy(&z2);
298        t1.neg();
299
300        z1.add(&t0);
301        self.b.copy(&z1);
302        self.b.add(&t1);
303
304        z3.add(&t1);
305        z2.add(&t0);
306
307        t0.copy(&self.a);
308        t0.add(&self.c);
309        t0.norm();
310        t1.copy(&y.a);
311        t1.add(&y.c);
312        t1.norm();
313        t0.mul(&t1);
314        z2.add(&t0);
315
316        t0.copy(&self.c);
317        t0.mul(&y.c);
318        t1.copy(&t0);
319        t1.neg();
320
321        self.c.copy(&z2);
322        self.c.add(&t1);
323        z3.add(&t1);
324        t0.times_i();
325        self.b.add(&t0);
326        z3.norm();
327
328        z3.times_i();
329        self.a.copy(&z0);
330        self.a.add(&z3);
331        self.norm();
332    }
333
334    /* Special case of multiplication arises from special form of ATE pairing line function */
335    pub fn smul(&mut self, y: &FP12, twist: usize) {
336        if twist == SexticTwist::D_TYPE.into() {
337            let mut z0 = FP4::new_copy(&self.a);
338            let mut z2 = FP4::new_copy(&self.b);
339            let mut z3 = FP4::new_copy(&self.b);
340            let mut t0 = FP4::new();
341            let mut t1 = FP4::new_copy(&y.a);
342
343            z0.mul(&y.a);
344            z2.pmul(&y.b.real());
345            self.b.add(&self.a);
346            t1.padd(&y.b.real());
347
348            self.b.norm();
349            t1.norm();
350
351            self.b.mul(&t1);
352            z3.add(&self.c);
353            z3.norm();
354            z3.pmul(&y.b.real());
355
356            t0.copy(&z0);
357            t0.neg();
358            t1.copy(&z2);
359            t1.neg();
360
361            self.b.add(&t0);
362
363            self.b.add(&t1);
364            z3.add(&t1);
365            z2.add(&t0);
366
367            t0.copy(&self.a);
368            t0.add(&self.c);
369            t0.norm();
370            z3.norm();
371
372            t0.mul(&y.a);
373            self.c.copy(&z2);
374            self.c.add(&t0);
375
376            z3.times_i();
377            self.a.copy(&z0);
378            self.a.add(&z3);
379        }
380        if twist == SexticTwist::M_TYPE.into() {
381            let mut z0 = FP4::new_copy(&self.a);
382            let mut z1 = FP4::new();
383            let mut z2 = FP4::new();
384            let mut z3 = FP4::new();
385            let mut t0 = FP4::new_copy(&self.a);
386            let mut t1 = FP4::new();
387
388            z0.mul(&y.a);
389            t0.add(&self.b);
390            t0.norm();
391
392            z1.copy(&t0);
393            z1.mul(&y.a);
394            t0.copy(&self.b);
395            t0.add(&self.c);
396            t0.norm();
397
398            z3.copy(&t0); //z3.mul(y.c);
399            z3.pmul(&y.c.getb());
400            z3.times_i();
401
402            t0.copy(&z0);
403            t0.neg();
404
405            z1.add(&t0);
406            self.b.copy(&z1);
407            z2.copy(&t0);
408
409            t0.copy(&self.a);
410            t0.add(&self.c);
411            t1.copy(&y.a);
412            t1.add(&y.c);
413
414            t0.norm();
415            t1.norm();
416
417            t0.mul(&t1);
418            z2.add(&t0);
419
420            t0.copy(&self.c);
421
422            t0.pmul(&y.c.getb());
423            t0.times_i();
424
425            t1.copy(&t0);
426            t1.neg();
427
428            self.c.copy(&z2);
429            self.c.add(&t1);
430            z3.add(&t1);
431            t0.times_i();
432            self.b.add(&t0);
433            z3.norm();
434            z3.times_i();
435            self.a.copy(&z0);
436            self.a.add(&z3);
437        }
438        self.norm();
439    }
440
441    /* self=1/self */
442    pub fn inverse(&mut self) {
443        let mut f0 = FP4::new_copy(&self.a);
444        let mut f1 = FP4::new_copy(&self.b);
445        let mut f2 = FP4::new_copy(&self.a);
446        let mut f3 = FP4::new();
447
448        self.norm();
449        f0.sqr();
450        f1.mul(&self.c);
451        f1.times_i();
452        f0.sub(&f1);
453        f0.norm();
454
455        f1.copy(&self.c);
456        f1.sqr();
457        f1.times_i();
458        f2.mul(&self.b);
459        f1.sub(&f2);
460        f1.norm();
461
462        f2.copy(&self.b);
463        f2.sqr();
464        f3.copy(&self.a);
465        f3.mul(&self.c);
466        f2.sub(&f3);
467        f2.norm();
468
469        f3.copy(&self.b);
470        f3.mul(&f2);
471        f3.times_i();
472        self.a.mul(&f0);
473        f3.add(&self.a);
474        self.c.mul(&f1);
475        self.c.times_i();
476
477        f3.add(&self.c);
478        f3.norm();
479        f3.inverse();
480        self.a.copy(&f0);
481        self.a.mul(&f3);
482        self.b.copy(&f1);
483        self.b.mul(&f3);
484        self.c.copy(&f2);
485        self.c.mul(&f3);
486    }
487
488    /* self=self^p using Frobenius */
489    pub fn frob(&mut self, f: &FP2) {
490        let mut f2 = FP2::new_copy(f);
491        let mut f3 = FP2::new_copy(f);
492
493        f2.sqr();
494        f3.mul(&f2);
495
496        self.a.frob(&f3);
497        self.b.frob(&f3);
498        self.c.frob(&f3);
499
500        self.b.pmul(f);
501        self.c.pmul(&f2);
502    }
503
504    /* trace function */
505    pub fn trace(&mut self) -> FP4 {
506        let mut t = FP4::new();
507        t.copy(&self.a);
508        t.imul(3);
509        t.reduce();
510        return t;
511    }
512
513    /* convert from byte array to FP12 */
514    pub fn frombytes(w: &[u8]) -> FP12 {
515        let mut t: [u8; big::MODBYTES as usize] = [0; big::MODBYTES as usize];
516        let mb = big::MODBYTES as usize;
517
518        for i in 0..mb {
519            t[i] = w[i]
520        }
521        let mut a = BIG::frombytes(&t);
522        for i in 0..mb {
523            t[i] = w[i + mb]
524        }
525        let mut b = BIG::frombytes(&t);
526        let mut c = FP2::new_bigs(&a, &b);
527
528        for i in 0..mb {
529            t[i] = w[i + 2 * mb]
530        }
531        a = BIG::frombytes(&t);
532        for i in 0..mb {
533            t[i] = w[i + 3 * mb]
534        }
535        b = BIG::frombytes(&t);
536        let mut d = FP2::new_bigs(&a, &b);
537
538        let e = FP4::new_fp2s(&c, &d);
539
540        for i in 0..mb {
541            t[i] = w[i + 4 * mb]
542        }
543        a = BIG::frombytes(&t);
544        for i in 0..mb {
545            t[i] = w[i + 5 * mb]
546        }
547        b = BIG::frombytes(&t);
548        c = FP2::new_bigs(&a, &b);
549
550        for i in 0..mb {
551            t[i] = w[i + 6 * mb]
552        }
553        a = BIG::frombytes(&t);
554        for i in 0..mb {
555            t[i] = w[i + 7 * mb]
556        }
557        b = BIG::frombytes(&t);
558        d = FP2::new_bigs(&a, &b);
559
560        let f = FP4::new_fp2s(&c, &d);
561
562        for i in 0..mb {
563            t[i] = w[i + 8 * mb]
564        }
565        a = BIG::frombytes(&t);
566        for i in 0..mb {
567            t[i] = w[i + 9 * mb]
568        }
569        b = BIG::frombytes(&t);
570
571        c = FP2::new_bigs(&a, &b);
572
573        for i in 0..mb {
574            t[i] = w[i + 10 * mb]
575        }
576        a = BIG::frombytes(&t);
577        for i in 0..mb {
578            t[i] = w[i + 11 * mb]
579        }
580        b = BIG::frombytes(&t);
581        d = FP2::new_bigs(&a, &b);
582
583        let g = FP4::new_fp2s(&c, &d);
584
585        return FP12::new_fp4s(&e, &f, &g);
586    }
587
588    /* convert this to byte array */
589    pub fn tobytes(&mut self, w: &mut [u8]) {
590        let mut t: [u8; big::MODBYTES as usize] = [0; big::MODBYTES as usize];
591        let mb = big::MODBYTES as usize;
592
593        self.a.geta().geta().tobytes(&mut t);
594        for i in 0..mb {
595            w[i] = t[i]
596        }
597        self.a.geta().getb().tobytes(&mut t);
598        for i in 0..mb {
599            w[i + mb] = t[i]
600        }
601        self.a.getb().geta().tobytes(&mut t);
602        for i in 0..mb {
603            w[i + 2 * mb] = t[i]
604        }
605        self.a.getb().getb().tobytes(&mut t);
606        for i in 0..mb {
607            w[i + 3 * mb] = t[i]
608        }
609
610        self.b.geta().geta().tobytes(&mut t);
611        for i in 0..mb {
612            w[i + 4 * mb] = t[i]
613        }
614        self.b.geta().getb().tobytes(&mut t);
615        for i in 0..mb {
616            w[i + 5 * mb] = t[i]
617        }
618        self.b.getb().geta().tobytes(&mut t);
619        for i in 0..mb {
620            w[i + 6 * mb] = t[i]
621        }
622        self.b.getb().getb().tobytes(&mut t);
623        for i in 0..mb {
624            w[i + 7 * mb] = t[i]
625        }
626
627        self.c.geta().geta().tobytes(&mut t);
628        for i in 0..mb {
629            w[i + 8 * mb] = t[i]
630        }
631        self.c.geta().getb().tobytes(&mut t);
632        for i in 0..mb {
633            w[i + 9 * mb] = t[i]
634        }
635        self.c.getb().geta().tobytes(&mut t);
636        for i in 0..mb {
637            w[i + 10 * mb] = t[i]
638        }
639        self.c.getb().getb().tobytes(&mut t);
640        for i in 0..mb {
641            w[i + 11 * mb] = t[i]
642        }
643    }
644
645    /* output to hex string */
646    pub fn tostring(&mut self) -> String {
647        return format!(
648            "[{},{},{}]",
649            self.a.tostring(),
650            self.b.tostring(),
651            self.c.tostring()
652        );
653    }
654
655    pub fn to_hex(&self) -> String {
656        format!("{} {} {}", self.a.to_hex(), self.b.to_hex(), self.c.to_hex())
657    }
658
659    pub fn from_hex_iter(iter: &mut SplitWhitespace) -> FP12 {
660        FP12 {
661            a: FP4::from_hex_iter(iter),
662            b: FP4::from_hex_iter(iter),
663            c: FP4::from_hex_iter(iter)
664        }
665    }
666
667    pub fn from_hex(val: String) -> FP12 {
668        let mut iter = val.split_whitespace();
669        return FP12::from_hex_iter(&mut iter);
670    }
671
672    /* self=self^e */
673    pub fn pow(&self, e: &BIG) -> FP12 {
674        let mut r = FP12::new_copy(self);
675        r.norm();
676        let mut e1 = BIG::new_copy(e);
677        e1.norm();
678        let mut e3 = BIG::new_copy(&e1);
679        e3.pmul(3);
680        e3.norm();
681        let mut w = FP12::new_copy(&r);
682
683        let nb = e3.nbits();
684        for i in (1..nb - 1).rev() {
685            w.usqr();
686            let bt = e3.bit(i) - e1.bit(i);
687            if bt == 1 {
688                w.mul(&r);
689            }
690            if bt == -1 {
691                r.conj();
692                w.mul(&r);
693                r.conj();
694            }
695        }
696
697        w.reduce();
698        return w;
699    }
700
701    /* constant time powering by small integer of max length bts */
702    pub fn pinpow(&mut self, e: i32, bts: i32) {
703        let mut r: [FP12; 2] = [FP12::new_int(1), FP12::new_copy(self)];
704        let mut t = FP12::new();
705
706        for i in (0..bts).rev() {
707            let b: usize = ((e >> i) & 1) as usize;
708            t.copy(&r[b]);
709            r[1 - b].mul(&t);
710            r[b].usqr();
711        }
712        self.copy(&r[0]);
713    }
714
715    pub fn compow(&mut self, e: &BIG, r: &BIG) -> FP4 {
716        let f = FP2::new_bigs(&BIG::new_ints(&rom::FRA), &BIG::new_ints(&rom::FRB));
717        let q = BIG::new_ints(&rom::MODULUS);
718
719        let mut g1 = FP12::new_copy(self);
720        let mut g2 = FP12::new_copy(self);
721
722        let mut m = BIG::new_copy(&q);
723        m.rmod(&r);
724
725        let mut a = BIG::new_copy(&e);
726        a.rmod(&mut m);
727
728        let mut b = BIG::new_copy(&e);
729        b.div(&mut m);
730
731        let mut c = g1.trace();
732
733        if b.iszilch() {
734            c = c.xtr_pow(&mut a);
735            return c;
736        }
737
738        g2.frob(&f);
739        let cp = g2.trace();
740        g1.conj();
741        g2.mul(&g1);
742        let cpm1 = g2.trace();
743        g2.mul(&g1);
744        let cpm2 = g2.trace();
745
746        c = c.xtr_pow2(&cp, &cpm1, &cpm2, &mut a, &mut b);
747
748        return c;
749    }
750
751    /* p=q0^u0.q1^u1.q2^u2.q3^u3 */
752    // Bos & Costello https://eprint.iacr.org/2013/458.pdf
753    // Faz-Hernandez & Longa & Sanchez  https://eprint.iacr.org/2013/158.pdf
754    // Side channel attack secure
755    pub fn pow4(q: &[FP12], u: &[BIG]) -> FP12 {
756        let mut g: [FP12; 8] = [
757            FP12::new(),
758            FP12::new(),
759            FP12::new(),
760            FP12::new(),
761            FP12::new(),
762            FP12::new(),
763            FP12::new(),
764            FP12::new(),
765        ];
766
767        let mut r = FP12::new();
768        let mut p = FP12::new();
769        const CT: usize = 1 + big::NLEN * (big::BASEBITS as usize);
770        let mut w: [i8; CT] = [0; CT];
771        let mut s: [i8; CT] = [0; CT];
772
773        let mut mt = BIG::new();
774        let mut t: [BIG; 4] = [
775            BIG::new_copy(&u[0]),
776            BIG::new_copy(&u[1]),
777            BIG::new_copy(&u[2]),
778            BIG::new_copy(&u[3]),
779        ];
780
781        for i in 0..4 {
782            t[i].norm();
783        }
784
785        // precomputation
786        g[0].copy(&q[0]);
787        r.copy(&g[0]);
788        g[1].copy(&r);
789        g[1].mul(&q[1]); // q[0].q[1]
790        g[2].copy(&r);
791        g[2].mul(&q[2]);
792        r.copy(&g[1]); // q[0].q[2]
793        g[3].copy(&r);
794        g[3].mul(&q[2]);
795        r.copy(&g[0]); // q[0].q[1].q[2]
796        g[4].copy(&r);
797        g[4].mul(&q[3]);
798        r.copy(&g[1]); // q[0].q[3]
799        g[5].copy(&r);
800        g[5].mul(&q[3]);
801        r.copy(&g[2]); // q[0].q[1].q[3]
802        g[6].copy(&r);
803        g[6].mul(&q[3]);
804        r.copy(&g[3]); // q[0].q[2].q[3]
805        g[7].copy(&r);
806        g[7].mul(&q[3]); // q[0].q[1].q[2].q[3]
807
808        // Make it odd
809        let pb = 1 - t[0].parity();
810        t[0].inc(pb);
811        t[0].norm();
812
813        // Number of bits
814        mt.zero();
815        for i in 0..4 {
816            mt.or(&t[i]);
817        }
818
819        let nb = 1 + mt.nbits();
820
821        // Sign pivot
822        s[nb - 1] = 1;
823        for i in 0..nb - 1 {
824            t[0].fshr(1);
825            s[i] = (2 * t[0].parity() - 1) as i8;
826            //println!("s={}",s[i]);
827        }
828
829        // Recoded exponent
830        for i in 0..nb {
831            w[i] = 0;
832            let mut k = 1;
833            for j in 1..4 {
834                let bt = s[i] * (t[j].parity() as i8);
835                t[j].fshr(1);
836                t[j].dec((bt >> 1) as isize);
837                t[j].norm();
838                w[i] += bt * (k as i8);
839                k = 2 * k;
840            }
841        }
842
843        // Main loop
844        p.selector(&g, (2 * w[nb - 1] + 1) as i32);
845        for i in (0..nb - 1).rev() {
846            p.usqr();
847            r.selector(&g, (2 * w[i] + s[i]) as i32);
848            p.mul(&r);
849        }
850
851        // apply correction
852        r.copy(&q[0]);
853        r.conj();
854        r.mul(&p);
855        p.cmove(&r, pb);
856        p.reduce();
857        return p;
858    }
859}