/*
 *    MIRACL  C++ Implementation file zzn2.cpp
 *
 *    AUTHOR  : M. Scott
 *  
 *    PURPOSE : Implementation of class ZZn2  (Arithmetic over n^2)
 *
 * WARNING: This class has been cobbled together for a specific use with
 * the MIRACL library. It is not complete, and may not work in other 
 * applications
 *
 *    Note: This code assumes that -1 is a Quadratic Non-Residue,
 *          In other words the modulus p is a prime = 3 mod 4
 *          OR -2 is a QNR if p=5 mod 8
 *
 *    Copyright (c) 2001-2004 Shamus Software Ltd.
 */

#include "zzn2.h"

using namespace std;

void ZZn2::get(Big& x,Big& y)  
{x=Big(a); y=Big(b);} 

void ZZn2::get(Big& x) 
{x=Big(a); }


ZZn2& ZZn2::operator*=(const ZZn2& x)
{ // optimized to reduce constructor/destructor calls

 int qnr=get_mip()->qnr;

 if (&x==this)
 { // a = (a+b)(a-b), b=2ab  (2 muls)
   ZZn t=a;  t+=b;
   ZZn t2=a; t2-=b; if (qnr==-2) t2-=b;
   t*=t2;
   b*=a;
   if (qnr==-2) t+=b;
   b+=b;
   a=t;
 }
 else
 { // a=a*x.a-b*x.b), b=(a+b)(x.a+x.b) - a*x.a - b*x.b  (3 muls)
   ZZn t=a; t*=x.a;
   ZZn t2=b; t2*=x.b;
   ZZn t3=x.a; t3+=x.b;
   b+=a; b*=t3; b-=t; b-=t2;
   if (qnr==-2) t2+=t2;
   t-=t2; a=t;  
 }

 return *this;
}

ZZn2& ZZn2::operator/=(const ZZn2& x)
{
    *this*=inverse(x);
     return *this;
}

ZZn2& ZZn2::operator/=(const ZZn& x)
{
    ZZn t=(ZZn)1/x;
    a*=t; b*=t;
    return *this;
}

ZZn2& ZZn2::operator/=(int i)
{
    ZZn t=(ZZn)1/i;
    a*=t; b*=t;
    return *this;
}

ZZn2 inverse(const ZZn2& w)
{
    ZZn2 y=conj(w);
    ZZn u=w.a;
    ZZn v=w.b;
    ZZn i=(ZZn)1;
    int qnr=get_mip()->qnr;
    u*=u;
    v*=v;
    if (qnr== -2) v+=v;
    u+=v;
    i/=u;
    y*=i;
    return y;
}

ZZn2 operator+(const ZZn2& x,const ZZn2& y) 
{ZZn2 w=x; w.a+=y.a; w.b+=y.b; return w; } 

ZZn2 operator+(const ZZn2& x,const ZZn& y) 
{ZZn2 w=x; w.a+=y; return w; } 

ZZn2 operator-(const ZZn2& x,const ZZn2& y) 
{ZZn2 w=x; w.a-=y.a; w.b-=y.b; return w; } 

ZZn2 operator-(const ZZn2& x,const ZZn& y) 
{ZZn2 w=x; w.a-=y; return w; } 

ZZn2 operator-(const ZZn2& x) 
{ZZn2 w; w.a=-x.a; w.b=-x.b; return w; } 

ZZn2 operator*(const ZZn2& x,const ZZn2& y)
{
 ZZn2 w=x;  
 if (&x==&y) w*=w;
 else        w*=y;  
 return w;
}

ZZn2 operator*(const ZZn2& x,const ZZn& y)
{ZZn2 w=x; w*=y; return w;}

ZZn2 operator*(const ZZn& y,const ZZn2& x)
{ZZn2 w=x; w*=y; return w;}

ZZn2 operator*(const ZZn2& x,int y)
{ZZn2 w=x; w*=y; return w;}

ZZn2 operator*(int y,const ZZn2& x)
{ZZn2 w=x; w*=y;return w;}

ZZn2 operator/(const ZZn2& x,const ZZn2& y)
{ZZn2 w=x; w/=y; return w;}

ZZn2 operator/(const ZZn2& x,const ZZn& y)
{ZZn2 w=x; w/=y; return w;}

ZZn2 operator/(const ZZn2& x,int i)
{ZZn2 w=x; w/=i; return w;}

ZZn2 randn2(void)
{ZZn2 w; w.a=randn(); w.b=randn(); return w;}

BOOL is_on_curve(const ZZn2& x)
{
    ZZn2 w;
    int qnr=get_mip()->qnr;
    BOOL twist=get_mip()->TWIST;

    if (twist)
    {
        w=x*x*x+qnr*getA()*x+qnr*tx((ZZn2)getB());
    }
    else
    {
        w=x*x*x+getA()*x+getB();
    }

    if (qr(w)) return TRUE;
    return FALSE;
}

BOOL qr(const ZZn2& x)
{
    ZZn s;
    int qnr=get_mip()->qnr;
    if (x.iszero()) return TRUE;
    if (x.b.iszero()) return TRUE;
    if (qnr==-1)
    {
        if (x.a.iszero()) return TRUE;
    }

    s=x.b;
    s*=s;
    if (qnr== -2) s+=s;
    if (!qr(x.a*x.a+s)) return FALSE;
    s=sqrt(x.a*x.a+s);
    if (qr((x.a+s)/2) || qr((x.a-s)/2)) return TRUE;
  
    return FALSE;
}

ZZn2 sqrt(const ZZn2& x)
{
// sqrt(a+ib) = sqrt(a+sqrt(a*a+n*b*b)/2)+ib/(2*sqrt(a+sqrt(a*a+n*b*b)/2))
// where i*i=n

    ZZn2 w;
    ZZn a,s,t;
    int qnr=get_mip()->qnr;

    if (x.iszero()) return w; 
    if (x.b.iszero())
    {
        if (qr(x.a))
        {
            s=sqrt(x.a);
            w.a=s; w.b=0;
        }
        else
        {
            s=sqrt(x.a/qnr);
            w.a=0; w.b=s;
        }
        return w;
    }

    if (qnr==-1)
    {
        if (x.a.iszero())
        {
            if (qr(x.b/2))
            {
                s=sqrt(x.b/2);
                w.a=w.b=s;
            }
            else
            {
                s=sqrt(-x.b/2);
                w.a=-s; w.b=s;
            }
            return w;
        }
    }
    s=x.b;
    s*=s;
    if (qnr== -2) s+=s;
    s=sqrt(x.a*x.a+s);

    if (s.iszero()) return w;

    if (qr((x.a+s)/2))
    {
        a=sqrt((x.a+s)/2);
    }    
    else
    {
        a=sqrt((x.a-s)/2);
        if (a.iszero()) return w;
    } 
    
    w.a=a;
    w.b=x.b/(2*a);

    return w;
}

ZZn2 conj(const ZZn2& x)
{
    ZZn2 u=x;
    u.conj();
    return u;
}

// for use by ZZn4...

ZZn2 txx(const ZZn2& w)
{ // multiply w by t^2 where x^2-t is irreducible polynomial for ZZn4
  //
  // for p=5 mod 8 t=sqrt(sqrt(-2)), qnr=-2
  // for p=3 mod 8 t=sqrt(1+sqrt(-1)), qnr=-1
  // for p=7 mod 8 and p=1 mod 3 t=sqrt(1+sqrt(-2)), qnr=-2

    switch (get_mip()->pmod8)
    {
    case 5:
        return tx(w);
    case 3:
    case 7:
        return w+tx(w);
    default: return w;
    }
}

ZZn2 txd(const ZZn2& w)
{ // divide w by t^2 where x^2-t is irreducible polynomial for ZZn4
  //
  // for p=5 mod 8 t=sqrt(sqrt(-2)), qnr=-2
  // for p=3 mod 8 t=sqrt(1+sqrt(-1)), qnr=-1
  // for p=7 mod 8 and p=1 mod 3 t=sqrt(1+sqrt(-2)), qnr=-2
    ZZn2 u;
    switch (get_mip()->pmod8)
    {
    case 5:
        u.set(w.b,-w.a/2);
        return u;
    case 3:
        u.set(w.a+w.b,w.b-w.a);
        return (u/2);
    case 7:
        u.set(w.a+w.b+w.b,w.b-w.a);
        return (u/3);
        break;
        
    default: return u;
    }
}

ZZn2 tx(const ZZn2& w)
{ // multiply w by i, the square root of the qnr
    ZZn t=-w.b;
    int qnr=get_mip()->qnr;
    if (qnr==-2) t+=t;
    ZZn2 u(t,w.a);
    return u;
}

// regular ZZn2 powering - but see powl function in zzn.h

ZZn2 pow(const ZZn2& x,const Big& k)
{
    int i,j,nb,n,nbw,nzs;
    ZZn2 u,u2,t[16];

    if (x.iszero()) return (ZZn2)0;
    if (k==0) return (ZZn2)1;
    u=x;
    if (k==1) return u;
//
// Prepare table for windowing
//
    u2=(u*u);
    t[0]=u;
   
    for (i=1;i<16;i++)
        t[i]=u2*t[i-1];

// Left to right method - with windows

    nb=bits(k);
    if (nb>1) for (i=nb-2;i>=0;)
    {
        n=window(k,i,&nbw,&nzs);
        for (j=0;j<nbw;j++) u*=u;
        if (n>0) u*=t[n/2];
        i-=nbw;
        if (nzs)
        {
            for (j=0;j<nzs;j++) u*=u;
            i-=nzs;
        }
    }
    return u;
}

ZZn2 powl(const ZZn2& x,const Big& k)
{
    ZZn2 w8,w9,two,y;
    int i,nb;

    two=(ZZn)2;
    y=two*x;
    if (k==0)  return (ZZn2)two;
    if (k==1)  return y;

    w8=two;
    w9=y;
    nb=bits(k);
    for (i=nb-1;i>=0;i--)
    {
        if (bit(k,i))
        {
            w8*=w9; w8-=y; w9*=w9; w9-=two;
        }
        else
        {
            w9*=w8; w9-=y; w8*=w8; w8-=two;
        }
    }
    return (w8/2);
}

#ifndef MR_NO_STANDARD_IO

ostream& operator<<(ostream& s,const ZZn2& xx)
{
    ZZn2 b=xx;
    Big x,y;
    b.get(x,y);
    s << "[" << x << "," << y << "]";
    return s;    
}

#endif

