private BigInteger oddModPow(BigInteger y,BigInteger z){
  if (y.equals(ONE))   return this;
  if (signum == 0)   return ZERO;
  int[] base=mag.clone();
  int[] exp=y.mag;
  int[] mod=z.mag;
  int modLen=mod.length;
  if ((modLen & 1) != 0) {
    int[] x=new int[modLen + 1];
    System.arraycopy(mod,0,x,1,modLen);
    mod=x;
    modLen++;
  }
  int wbits=0;
  int ebits=bitLength(exp,exp.length);
  if ((ebits != 17) || (exp[0] != 65537)) {
    while (ebits > bnExpModThreshTable[wbits]) {
      wbits++;
    }
  }
  int tblmask=1 << wbits;
  int[][] table=new int[tblmask][];
  for (int i=0; i < tblmask; i++)   table[i]=new int[modLen];
  long n0=(mod[modLen - 1] & LONG_MASK) + ((mod[modLen - 2] & LONG_MASK) << 32);
  long inv=-MutableBigInteger.inverseMod64(n0);
  int[] a=leftShift(base,base.length,modLen << 5);
  MutableBigInteger q=new MutableBigInteger(), a2=new MutableBigInteger(a), b2=new MutableBigInteger(mod);
  b2.normalize();
  MutableBigInteger r=a2.divide(b2,q);
  table[0]=r.toIntArray();
  if (table[0].length < modLen) {
    int offset=modLen - table[0].length;
    int[] t2=new int[modLen];
    System.arraycopy(table[0],0,t2,offset,table[0].length);
    table[0]=t2;
  }
  int[] b=montgomerySquare(table[0],mod,modLen,inv,null);
  int[] t=Arrays.copyOf(b,modLen);
  for (int i=1; i < tblmask; i++) {
    table[i]=montgomeryMultiply(t,table[i - 1],mod,modLen,inv,null);
  }
  int bitpos=1 << ((ebits - 1) & (32 - 1));
  int buf=0;
  int elen=exp.length;
  int eIndex=0;
  for (int i=0; i <= wbits; i++) {
    buf=(buf << 1) | (((exp[eIndex] & bitpos) != 0) ? 1 : 0);
    bitpos>>>=1;
    if (bitpos == 0) {
      eIndex++;
      bitpos=1 << (32 - 1);
      elen--;
    }
  }
  int multpos=ebits;
  ebits--;
  boolean isone=true;
  multpos=ebits - wbits;
  while ((buf & 1) == 0) {
    buf>>>=1;
    multpos++;
  }
  int[] mult=table[buf >>> 1];
  buf=0;
  if (multpos == ebits)   isone=false;
  while (true) {
    ebits--;
    buf<<=1;
    if (elen != 0) {
      buf|=((exp[eIndex] & bitpos) != 0) ? 1 : 0;
      bitpos>>>=1;
      if (bitpos == 0) {
        eIndex++;
        bitpos=1 << (32 - 1);
        elen--;
      }
    }
    if ((buf & tblmask) != 0) {
      multpos=ebits - wbits;
      while ((buf & 1) == 0) {
        buf>>>=1;
        multpos++;
      }
      mult=table[buf >>> 1];
      buf=0;
    }
    if (ebits == multpos) {
      if (isone) {
        b=mult.clone();
        isone=false;
      }
 else {
        t=b;
        a=montgomeryMultiply(t,mult,mod,modLen,inv,a);
        t=a;
        a=b;
        b=t;
      }
    }
    if (ebits == 0)     break;
    if (!isone) {
      t=b;
      a=montgomerySquare(t,mod,modLen,inv,a);
      t=a;
      a=b;
      b=t;
    }
  }
  int[] t2=new int[2 * modLen];
  System.arraycopy(b,0,t2,modLen,modLen);
  b=montReduce(t2,mod,modLen,(int)inv);
  t2=Arrays.copyOf(b,modLen);
  return new BigInteger(1,t2);
}

private BigInteger oddModPow ( BigInteger y , BigInteger z ) { if ( y . equals ( ONE ) ) return this ; if ( signum == 0 ) return ZERO ; int [ ] base = mag . clone ( ) ; int [ ] exp = y . mag ; int [ ] mod = z . mag ; int modLen = mod . length ; if ( ( modLen & 1 ) != 0 ) { int [ ] x = new int [ modLen + 1 ] ; System . arraycopy ( mod , 0 , x , 1 , modLen ) ; mod = x ; modLen ++ ; } int wbits = 0 ; int ebits = bitLength ( exp , exp . length ) ; if ( ( ebits != 17 ) || ( exp [ 0 ] != 65537 ) ) { while ( ebits > bnExpModThreshTable [ wbits ] ) { wbits ++ ; } } int tblmask = 1 << wbits ; int [ ] [ ] table = new int [ tblmask ] [ ] ; for ( int i = 0 ; i < tblmask ; i ++ ) table [ i ] = new int [ modLen ] ; long n0 = ( mod [ modLen - 1 ] & LONG_MASK ) + ( ( mod [ modLen - 2 ] & LONG_MASK ) << 32 ) ; long inv = - MutableBigInteger . inverseMod64 ( n0 ) ; int [ ] a = leftShift ( base , base . length , modLen << 5 ) ; MutableBigInteger q = new MutableBigInteger ( ) , a2 = new MutableBigInteger ( a ) , b2 = new MutableBigInteger ( mod ) ; b2 . normalize ( ) ; MutableBigInteger r = a2 . divide ( b2 , q ) ; table [ 0 ] = r . toIntArray ( ) ; if ( table [ 0 ] . length < modLen ) { int offset = modLen - table [ 0 ] . length ; int [ ] t2 = new int [ modLen ] ; System . arraycopy ( table [ 0 ] , 0 , t2 , offset , table [ 0 ] . length ) ; table [ 0 ] = t2 ; } int [ ] b = montgomerySquare ( table [ 0 ] , mod , modLen , inv , null ) ; int [ ] t = Arrays . copyOf ( b , modLen ) ; for ( int i = 1 ; i < tblmask ; i ++ ) { table [ i ] = montgomeryMultiply ( t , table [ i - 1 ] , mod , modLen , inv , null ) ; } int bitpos = 1 << ( ( ebits - 1 ) & ( 32 - 1 ) ) ; int buf = 0 ; int elen = exp . length ; int eIndex = 0 ; for ( int i = 0 ; i <= wbits ; i ++ ) { buf = ( buf << 1 ) | ( ( ( exp [ eIndex ] & bitpos ) != 0 ) ? 1 : 0 ) ; bitpos >>>= 1 ; if ( bitpos == 0 ) { eIndex ++ ; bitpos = 1 << ( 32 - 1 ) ; elen -- ; } } int multpos = ebits ; ebits -- ; boolean isone = true ; multpos = ebits - wbits ; while ( ( buf & 1 ) == 0 ) { buf >>>= 1 ; multpos ++ ; } int [ ] mult = table [ buf > > > 1 ] ; buf = 0 ; if ( multpos == ebits ) isone = false ; while ( true ) { ebits -- ; buf <<= 1 ; if ( elen != 0 ) { buf |= ( ( exp [ eIndex ] & bitpos ) != 0 ) ? 1 : 0 ; bitpos >>>= 1 ; if ( bitpos == 0 ) { eIndex ++ ; bitpos = 1 << ( 32 - 1 ) ; elen -- ; } } if ( ( buf & tblmask ) != 0 ) { multpos = ebits - wbits ; while ( ( buf & 1 ) == 0 ) { buf >>>= 1 ; multpos ++ ; } mult = table [ buf > > > 1 ] ; buf = 0 ; } if ( ebits == multpos ) { if ( isone ) { b = mult . clone ( ) ; isone = false ; } else { t = b ; a = montgomeryMultiply ( t , mult , mod , modLen , inv , a ) ; t = a ; a = b ; b = t ; } } if ( ebits == 0 ) break ; if ( ! isone ) { t = b ; a = montgomerySquare ( t , mod , modLen , inv , a ) ; t = a ; a = b ; b = t ; } } int [ ] t2 = new int [ 2 * modLen ] ; System . arraycopy ( b , 0 , t2 , modLen , modLen ) ; b = montReduce ( t2 , mod , modLen , ( int ) inv ) ; t2 = Arrays . copyOf ( b , modLen ) ; return new BigInteger ( 1 , t2 ) ; }
returns a biginteger whose value is x to the power of y mod z . assumes : z is odd & & x < z .
