Return to N.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Changed a tab to two space charaters.

/*
 * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
 * All rights reserved.
 * 
 * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
 * non-exclusive and royalty-free license to use, copy, modify and
 * redistribute, solely for non-commercial and non-profit purposes, the
 * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
 * conditions of this Agreement. For the avoidance of doubt, you acquire
 * only a limited right to use the SOFTWARE hereunder, and FLL or any
 * third party developer retains all rights, including but not limited to
 * copyrights, in and to the SOFTWARE.
 * 
 * (1) FLL does not grant you a license in any way for commercial
 * purposes. You may use the SOFTWARE only for non-commercial and
 * non-profit purposes only, such as academic, research and internal
 * business use.
 * (2) The SOFTWARE is protected by the Copyright Law of Japan and
 * international copyright treaties. If you make copies of the SOFTWARE,
 * with or without modification, as permitted hereunder, you shall affix
 * to all such copies of the SOFTWARE the above copyright notice.
 * (3) An explicit reference to this SOFTWARE and its copyright owner
 * shall be made on your publication or presentation in any form of the
 * results obtained by use of the SOFTWARE.
 * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
 * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
 * for such modification or the source code of the modified part of the
 * SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
 * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
 * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
 * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
 * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
 * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
 * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
 * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
 * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
 * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
 * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
 * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
 * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
 * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
 * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
 * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
 * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
 *
 * $OpenXM: OpenXM_contrib2/asir2000/engine/N.c,v 1.12 2018/03/29 01:32:51 noro Exp $ 
*/
#include "ca.h"
#include "base.h"

#if (defined(_M_IX86) || defined(i386)) && !defined(__MINGW32__)
void addn(N n1,N n2,N *nr)
{
  unsigned int *m1,*m2,*mr;
  unsigned int c;
  N r;
  int i,d1,d2;
  unsigned int tmp;

  if ( !n1 )
    COPY(n2,*nr);
  else if ( !n2 )
    COPY(n1,*nr);
  else {
    if ( PL(n1) > PL(n2) ) {  
      d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
    } else {
      d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
    }
    *nr = r = NALLOC(d1 + 1); INITRC(r); mr = BD(r);

#if defined(_M_IX86)
    __asm {
    push  esi
    push  edi
    mov esi,m1
    mov edi,m2
    mov ebx,mr
    mov ecx,d2
    xor  eax,eax
    Lstart_addn:
    mov eax,DWORD PTR [esi]
    mov edx,DWORD PTR [edi]
    adc eax,edx
    mov DWORD PTR [ebx],eax
    lea esi,DWORD PTR [esi+4]
    lea edi,DWORD PTR [edi+4]
    lea ebx,DWORD PTR [ebx+4]
    dec ecx
    jnz Lstart_addn
    pop  edi
    pop  esi
    mov eax,0
    adc eax,eax
    mov c,eax
    }
#else
    asm volatile("\
    pushl  %%ebx;\
    movl  %1,%%esi;\
    movl  %2,%%edi;\
    movl  %3,%%ebx;\
    movl  %4,%%ecx;\
    testl  %%eax,%%eax;\
    Lstart_addn:;\
    movl  (%%esi),%%eax;\
    movl  (%%edi),%%edx;\
    adcl  %%edx,%%eax;\
    movl  %%eax,(%%ebx);\
    leal  4(%%esi),%%esi;\
    leal  4(%%edi),%%edi;\
    leal  4(%%ebx),%%ebx;\
    decl  %%ecx;\
    jnz Lstart_addn;\
    movl  $0,%%eax;\
    adcl  %%eax,%%eax;\
    movl  %%eax,%0;\
    popl  %%ebx"\
    :"=m"(c)\
    :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
    :"eax","ecx","edx","esi","edi");
#endif
    for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
      tmp = *m1++ + c;
      c = tmp < c ? 1 : 0;
      *mr++ = tmp;
    }
    for ( ; i < d1; i++ ) 
      *mr++ = *m1++;
    *mr = c;
    PL(r) = (c?d1+1:d1);
  }
}

int subn(N n1,N n2,N *nr)
{
  N r;
  unsigned int *m1,*m2,*mr,br;
  unsigned int tmp,t;
  int d1,d2,sgn,i;

  if ( !n1 ) {
    if ( n2 ) {
      COPY(n2,*nr);
      return -1;
    } else {
      *nr = 0;
      return 0;
    }
  } else if ( !n2 ) {
    COPY(n1,*nr);
    return 1;
  } else {
    d1 = PL(n1); d2 = PL(n2);
    m1 = BD(n1); m2 = BD(n2);
    if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
      sgn = 1;
    else if ( d1 < d2 ) {
      d1 = PL(n2); d2 = PL(n1);
      m1 = BD(n2); m2 = BD(n1);
      sgn = -1;
    } else {
      for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
      if ( i < 0 ) {
        *nr = 0;
        return 0;
      }
      d1 = d2 = i+1;
      if ( m1[i] > m2[i] )
        sgn = 1;
      else {
        m1 = BD(n2); m2 = BD(n1);
        sgn = -1;
      }
    }
    *nr = r = NALLOC(d1); INITRC(r); mr = BD(r);

#if defined(_M_IX86)
    __asm {
    push  esi
    push  edi
    mov esi,m1
    mov edi,m2
    mov ebx,mr
    mov ecx,d2
    xor  eax,eax
    Lstart_subn:
    mov eax,DWORD PTR [esi]
    mov edx,DWORD PTR [edi]
    sbb eax,edx
    mov DWORD PTR [ebx],eax
    lea esi,DWORD PTR [esi+4]
    lea edi,DWORD PTR [edi+4]
    lea ebx,DWORD PTR [ebx+4]
    dec ecx
    jnz Lstart_subn
    pop  edi
    pop  esi
    mov eax,0
    adc eax,eax
    mov br,eax
    }
#else
    asm volatile("\
    pushl  %%ebx;\
    movl  %1,%%esi;\
    movl  %2,%%edi;\
    movl  %3,%%ebx;\
    movl  %4,%%ecx;\
    testl  %%eax,%%eax;\
    Lstart_subn:;\
    movl  (%%esi),%%eax;\
    movl  (%%edi),%%edx;\
    sbbl  %%edx,%%eax;\
    movl  %%eax,(%%ebx);\
    leal  4(%%esi),%%esi;\
    leal  4(%%edi),%%edi;\
    leal  4(%%ebx),%%ebx;\
    decl  %%ecx;\
    jnz Lstart_subn;\
    movl  $0,%%eax;\
    adcl  %%eax,%%eax;\
    movl  %%eax,%0;\
    popl  %%ebx"\
    :"=m"(br)\
    :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
    :"eax","ecx","edx","esi","edi");
#endif
    for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
      t = *m1++;
      tmp = t - br;
      br = tmp > t ? 1 : 0;
      *mr++ = tmp;
    }
    for ( ; i < d1; i++ )
      *mr++ = *m1++;
    for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
    PL(r) = i + 1;
    return sgn;
  }
}

void _addn(N n1,N n2,N nr)
{
  unsigned int *m1,*m2,*mr;
  unsigned int c;
  int i,d1,d2;
  unsigned int tmp;

  if ( !n1 || !PL(n1) )
    dupn(n2,nr);
  else if ( !n2 || !PL(n2) )
    dupn(n1,nr);
  else {
    if ( PL(n1) > PL(n2) ) {  
      d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
    } else {
      d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
    }
    mr = BD(nr);

#if defined(_M_IX86)
    __asm {
    push  esi
    push  edi
    mov esi,m1
    mov edi,m2
    mov ebx,mr
    mov ecx,d2
    xor  eax,eax
    Lstart__addn:
    mov eax,DWORD PTR [esi]
    mov edx,DWORD PTR [edi]
    adc eax,edx
    mov DWORD PTR [ebx],eax
    lea esi,DWORD PTR [esi+4]
    lea edi,DWORD PTR [edi+4]
    lea ebx,DWORD PTR [ebx+4]
    dec ecx
    jnz Lstart__addn
    pop  edi
    pop  esi
    mov eax,0
    adc eax,eax
    mov c,eax
    }
#else
    asm volatile("\
    pushl  %%ebx;\
    movl  %1,%%esi;\
    movl  %2,%%edi;\
    movl  %3,%%ebx;\
    movl  %4,%%ecx;\
    testl  %%eax,%%eax;\
    Lstart__addn:;\
    movl  (%%esi),%%eax;\
    movl  (%%edi),%%edx;\
    adcl  %%edx,%%eax;\
    movl  %%eax,(%%ebx);\
    leal  4(%%esi),%%esi;\
    leal  4(%%edi),%%edi;\
    leal  4(%%ebx),%%ebx;\
    decl  %%ecx;\
    jnz Lstart__addn;\
    movl  $0,%%eax;\
    adcl  %%eax,%%eax;\
    movl  %%eax,%0;\
    popl  %%ebx"\
    :"=m"(c)\
    :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
    :"eax","ecx","edx","esi","edi");
#endif
    for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
      tmp = *m1++ + c;
      c = tmp < c ? 1 : 0;
      *mr++ = tmp;
    }
    for ( ; i < d1; i++ ) 
      *mr++ = *m1++;
    *mr = c;
    PL(nr) = (c?d1+1:d1);
  }
}

int _subn(N n1,N n2,N nr)
{
  unsigned int *m1,*m2,*mr,br;
  unsigned int tmp,t;
  int d1,d2,sgn,i;

  if ( !n1 || !PL(n1) ) {
    if ( n2 && PL(n2) ) {
      dupn(n2,nr);
      return -1;
    } else {
      PL(nr) = 0;
      return 0;
    }
  } else if ( !n2 || !PL(n2) ) {
    dupn(n1,nr);
    return 1;
  } else {
    d1 = PL(n1); d2 = PL(n2);
    m1 = BD(n1); m2 = BD(n2);
    if ( (d1 = PL(n1)) > (d2 = PL(n2)) )
      sgn = 1;
    else if ( d1 < d2 ) {
      d1 = PL(n2); d2 = PL(n1);
      m1 = BD(n2); m2 = BD(n1);
      sgn = -1;
    } else {
      for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
      if ( i < 0 ) {
        PL(nr) = 0;
        return 0;
      }
      d1 = d2 = i+1;
      if ( m1[i] > m2[i] )
        sgn = 1;
      else {
        m1 = BD(n2); m2 = BD(n1);
        sgn = -1;
      }
    }
    mr = BD(nr);

#if defined(_M_IX86)
    __asm {
    push  esi
    push  edi
    mov esi,m1
    mov edi,m2
    mov ebx,mr
    mov ecx,d2
    xor  eax,eax
    Lstart__subn:
    mov eax,DWORD PTR [esi]
    mov edx,DWORD PTR [edi]
    sbb eax,edx
    mov DWORD PTR [ebx],eax
    lea esi,DWORD PTR [esi+4]
    lea edi,DWORD PTR [edi+4]
    lea ebx,DWORD PTR [ebx+4]
    dec ecx
    jnz Lstart__subn
    pop  edi
    pop  esi
    mov eax,0
    adc eax,eax
    mov br,eax
    }
#else
    asm volatile("\
    pushl  %%ebx;\
    movl  %1,%%esi;\
    movl  %2,%%edi;\
    movl  %3,%%ebx;\
    movl  %4,%%ecx;\
    testl  %%eax,%%eax;\
    Lstart__subn:;\
    movl  (%%esi),%%eax;\
    movl  (%%edi),%%edx;\
    sbbl  %%edx,%%eax;\
    movl  %%eax,(%%ebx);\
    leal  4(%%esi),%%esi;\
    leal  4(%%edi),%%edi;\
    leal  4(%%ebx),%%ebx;\
    decl  %%ecx;\
    jnz Lstart__subn;\
    movl  $0,%%eax;\
    adcl  %%eax,%%eax;\
    movl  %%eax,%0;\
    popl  %%ebx"\
    :"=m"(br)\
    :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
    :"eax","ecx","edx","esi","edi");
#endif
    for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
      t = *m1++;
      tmp = t - br;
      br = tmp > t ? 1 : 0;
      *mr++ = tmp;
    }
    for ( ; i < d1; i++ )
      *mr++ = *m1++;
    for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
    PL(nr) = i + 1;
    return sgn;
  }
}
#else

void addn(N n1,N n2,N *nr)
{
  unsigned int *m1,*m2,*mr,i,c;
  N r;
  int d1,d2;
  unsigned int tmp;

  if ( !n1 )
    COPY(n2,*nr);
  else if ( !n2 )
    COPY(n1,*nr);
  else {
    if ( PL(n1) > PL(n2) ) {  
      d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
    } else {
      d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
    }
    *nr = r = NALLOC(d1 + 1); INITRC(r);
    for ( i = 0, c = 0, mr = BD(r); i < d2; i++, m1++, m2++, mr++ ) {
      tmp = *m1 + *m2;
      if ( tmp < *m1 ) {
        tmp += c;
        c = 1;
      } else {
        tmp += c;
        c = tmp < c ? 1 : 0;
      }
      *mr = tmp;
    }
    for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
      tmp = *m1 + c;
      c = tmp < c ? 1 : 0;
      *mr = tmp;
    }
    for ( ; i < d1; i++ ) 
      *mr++ = *m1++;
    *mr = c;
    PL(r) = (c?d1+1:d1);
  }
}

int subn(N n1,N n2,N *nr)
{
  N r;
  unsigned int *m1,*m2,*mr,i,br;
  L tmp;
  int d1,d2,nz,sgn;

  if ( !n1 ) {
    if ( n2 ) {
      COPY(n2,*nr);
      return -1;
    } else {
      *nr = 0;
      return 0;
    }
  } else if ( !n2 ) {
    COPY(n1,*nr);
    return 1;
  } else {
    switch ( cmpn(n1,n2) ) {
      case 1:
        d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
        sgn = 1; break;
      case -1:
        d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
        sgn = -1; break;
      case 0:
      default:
        *nr = 0; return ( 0 ); break;
    }
    *nr = r = NALLOC(d1); INITRC(r);
    for ( i = 0, br = 0, nz = -1, mr = BD(r); 
      i < d2; i++ ) {
      if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) ) 
        nz = i;
      if ( tmp < 0 ) {
        br = 1; *mr++ = (unsigned int)(tmp + LBASE);
      } else {
        br = 0; *mr++ = (unsigned int)tmp;
      }
    }
    for ( ; (i < d1) && br; i++ ) {
      if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) ) 
        nz = i;
      if ( tmp < 0 ) {
        br = 1; *mr++ = (unsigned int)(tmp + LBASE);
      } else {
        br = 0; *mr++ = (unsigned int)tmp;
      }
    }
    for ( ; i < d1; i++ )
      if ( *mr++ = *m1++ )
        nz = i;
    PL(r) = nz + 1;
    return sgn;
  }
}

void _addn(N n1,N n2,N nr)
{
  unsigned int *m1,*m2,*mr,i,c;
  int d1,d2;
  unsigned int tmp;

  if ( !n1 || !PL(n1) )
    dupn(n2,nr);
  else if ( !n2 || !PL(n2) )
    dupn(n1,nr);
  else {
    if ( PL(n1) > PL(n2) ) {  
      d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
    } else {
      d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
    }
    for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {
      tmp = *m1 + *m2;
      if ( tmp < *m1 ) {
        tmp += c;
        c = 1;
      } else {
        tmp += c;
        c = tmp < c ? 1 : 0;
      }
      *mr = tmp;
    }
    for ( ; (i < d1) && c ; i++, m1++, mr++ ) {
      tmp = *m1 + c;
      c = tmp < c ? 1 : 0;
      *mr = tmp;
    }
    for ( ; i < d1; i++ ) 
      *mr++ = *m1++;
    *mr = c;
    PL(nr) = (c?d1+1:d1);
  }
}

int _subn(N n1,N n2,N nr)
{
  N r;
  unsigned int *m1,*m2,*mr,i,br;
  L tmp;
  int d1,d2,nz,sgn;

  if ( !n1 || !PL(n1) ) {
    if ( n2 && PL(n2) ) {
      dupn(n2,nr);
      return -1;
    } else {
      PL(nr) = 0;
      return 0;
    }
  } else if ( !n2 || !PL(n2) ) {
    dupn(n1,nr);
    return 1;
  } else {
    switch ( cmpn(n1,n2) ) {
      case 1:
        d1 = PL(n1); d2 = PL(n2); m1 = BD(n1); m2 = BD(n2);
        sgn = 1; break;
      case -1:
        d1 = PL(n2); d2 = PL(n1); m1 = BD(n2); m2 = BD(n1);
        sgn = -1; break;
      case 0:
      default:
        PL(nr) = 0; return ( 0 ); break;
    }
    for ( i = 0, br = 0, nz = -1, mr = BD(nr); 
      i < d2; i++ ) {
      if ( (tmp = (L)*m1++ - (L)*m2++ - (L)br) && ( tmp > -LBASE ) ) 
        nz = i;
      if ( tmp < 0 ) {
        br = 1; *mr++ = (unsigned int)(tmp + LBASE);
      } else {
        br = 0; *mr++ = (unsigned int)tmp;
      }
    }
    for ( ; (i < d1) && br; i++ ) {
      if ( (tmp = (L)*m1++ - (L)br) && ( tmp > -LBASE ) ) 
        nz = i;
      if ( tmp < 0 ) {
        br = 1; *mr++ = (unsigned int)(tmp + LBASE);
      } else {
        br = 0; *mr++ = (unsigned int)tmp;
      }
    }
    for ( ; i < d1; i++ )
      if ( *mr++ = *m1++ )
        nz = i;
    PL(nr) = nz + 1;
    return sgn;
  }
}
#endif

/* a2 += a1; n2 >= n1 */

void addarray_to(unsigned int *a1,int n1,unsigned int *a2,int n2)
{
  int i;
  unsigned int c,tmp;

  for ( i = 0, c = 0; i < n2; i++, a1++, a2++ ) {
    tmp = *a1 + *a2;
    if ( tmp < *a1 ) {
      tmp += c;
      c = 1;
    } else {
      tmp += c;
      c = tmp < c ? 1 : 0;
    }
    *a2 = tmp;
  }
  for ( ; (i < n2) && c ; i++, a2++ ) {
    tmp = *a2 + c;
    c = tmp < c ? 1 : 0;
    *a2 = tmp;
  }
  if ( i == n2 && c )
    *a2 = c;
}

void pwrn(N n,int e,N *nr)
{
  N nw,nw1;

  if ( e == 1 ) {
    COPY(n,*nr);
    return;
  }
  pwrn(n,e / 2,&nw);
  if ( e % 2 == 0 ) 
    kmuln(nw,nw,nr);
  else {
    kmuln(nw,nw,&nw1); kmuln(nw1,n,nr); FREEN(nw1);
  }
  FREEN(nw);
}



extern int igcd_algorithm;

void gcdEuclidn(), gcdn_HMEXT();

void lcmn(N n1,N n2,N *nr)
{
  N g,t;

  gcdn(n1,n2,&g);
  divsn(n1,g,&t);
  muln(t,n2,nr);
}

void gcdn(N n1,N n2,N *nr)
{
  if ( !igcd_algorithm )
    gcdEuclidn(n1,n2,nr);
  else {
    gcdn_HMEXT(n1,n2,nr);
  }
}

#include "Ngcd.c"

void gcdEuclidn(N n1,N n2,N *nr)
{
  N m1,m2,q,r;
  unsigned int i1,i2,ir;

  if ( !n1 )
    COPY(n2,*nr);
  else if ( !n2 )
    COPY(n1,*nr);
  else {
    if ( PL(n1) > PL(n2) ) {
      COPY(n1,m1); COPY(n2,m2);
    } else {
      COPY(n1,m2); COPY(n2,m1);
    }
    while ( PL(m1) > 1 ) {
      divn(m1,m2,&q,&r); FREEN(m1); FREEN(q);
      if ( !r ) {
        *nr = m2;
        return;
      } else {
        m1 = m2; m2 = r;
      }
    }
    for ( i1 = BD(m1)[0], i2 = BD(m2)[0]; ir = i1 % i2; i2 = ir )
      i1 = i2;
    if ( i2 == 1 )
      COPY(ONEN,*nr);
    else {
      *nr = r = NALLOC(1); INITRC(r); PL(r) = 1; BD(r)[0] = i2;
    }
  }
}

int cmpn(N n1,N n2)
{
  int i;
  unsigned int *m1,*m2;

  if ( !n1 ) 
    if ( !n2 ) 
      return 0;
    else 
      return -1;
  else if ( !n2 ) 
    return 1;
  else if ( PL(n1) > PL(n2) )
    return 1;
  else if ( PL(n1) < PL(n2) ) 
    return -1;
  else {  
    for ( i = PL(n1)-1, m1 = BD(n1)+i, m2 = BD(n2)+i;
      i >= 0; i--, m1--, m2-- ) 
      if ( *m1 > *m2 ) 
        return 1;
      else if ( *m1 < *m2 )
        return -1;
    return 0;
  }
}

void bshiftn(N n,int b,N *r)
{
  int w,l,nl,i,j;
  N z;
  unsigned int msw;
  unsigned int *p,*pz;

  if ( b == 0 ) {
    COPY(n,*r); return;
  }
  if ( b > 0 ) { /* >> */
    w = b / BSH; l = PL(n)-w;
    if ( l <= 0 ) {
      *r = 0; return;
    }
    b %= BSH; p = BD(n)+w;
    if ( !b ) {
      *r = z = NALLOC(l); INITRC(z); PL(z) = l;
      bcopy(p,BD(z),l*sizeof(unsigned int));
      return;
    }
    msw = p[l-1];
    for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
    if ( b >= i ) {
      l--;
      if ( !l ) {
        *r = 0; return;
      }
      *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
      for ( j = 0; j < l; j++, p++ )
        *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
    } else {
      *r = z = NALLOC(l); INITRC(z); PL(z) = l; pz = BD(z);
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
      *pz = *p>>b;
    }
  } else { /* << */
    b = -b;
    w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
    if ( !b ) {
      nl = l+w; *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      bcopy(p,BD(z)+w,l*sizeof(unsigned int));
      return;
    }
    msw = p[l-1];
    for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
    if ( b + i > BSH ) {
      nl = l+w+1;
      *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      *pz++ = *p++<<b;
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
      *pz = *(p-1)>>(BSH-b);
    } else {
      nl = l+w;
      *r = z = NALLOC(nl); INITRC(z); PL(z) = nl; pz = BD(z)+w;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      *pz++ = *p++<<b;
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
    }
  }
}

#if 0
void _bshiftn(N n,int b,N z)
{
  int w,l,nl,i,j;
  unsigned int msw;
  unsigned int *p,*pz;

  if ( b == 0 ) {
    copyn(n,PL(n),BD(z)); PL(z) = PL(n); return;
  }
  if ( b > 0 ) { /* >> */
    w = b / BSH; l = PL(n)-w;
    if ( l <= 0 ) {
      PL(z) = 0; return;
    }
    b %= BSH; p = BD(n)+w;
    if ( !b ) {
      PL(z) = l;
      bcopy(p,BD(z),l*sizeof(unsigned int));
      return;
    }
    msw = p[l-1];
    for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
    if ( b >= i ) {
      l--;
      if ( !l ) {
        PL(z) = 0; return;
      }
      PL(z) = l; pz = BD(z);
      for ( j = 0; j < l; j++, p++ )
        *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
    } else {
      PL(z) = l; pz = BD(z);
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*(p+1)<<(BSH-b))|(*p>>b);
      *pz = *p>>b;
    }
  } else { /* << */
    b = -b;
    w = b / BSH; b %= BSH; l = PL(n); p = BD(n);
    if ( !b ) {
      nl = l+w; PL(z) = nl;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      bcopy(p,BD(z)+w,l*sizeof(unsigned int));
      return;
    }
    msw = p[l-1];
    for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
    if ( b + i > BSH ) {
      nl = l+w+1;
      PL(z) = nl; pz = BD(z)+w;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      *pz++ = *p++<<b;
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
      *pz = *(p-1)>>(BSH-b);
    } else {
      nl = l+w;
      PL(z) = nl; pz = BD(z)+w;
      bzero((char *)BD(z),w*sizeof(unsigned int));
      *pz++ = *p++<<b;
      for ( j = 1; j < l; j++, p++ )
        *pz++ = (*p<<b)|(*(p-1)>>(BSH-b));
    }
  }
}
#endif

void shiftn(N n,int w,N *r)
{
  int l,nl;
  N z;

  if ( w == 0 )
    COPY(n,*r);
  else if ( w > 0 ) { /* >> */
    l = PL(n)-w;
    if ( l <= 0 )
      *r = 0;
    else {
      *r = z = NALLOC(l); INITRC(z); PL(z) = l;
      bcopy(BD(n)+w,BD(z),l*sizeof(unsigned int));
    }
  } else { /* << */
    w = -w;
    l = PL(n); nl = l+w;
    *r = z = NALLOC(nl); INITRC(z); PL(z) = nl;
    bzero((char *)BD(z),w*sizeof(unsigned int));
    bcopy(BD(n),BD(z)+w,l*sizeof(unsigned int));
  }
}

void randomn(int bits,N *r)
{
  int l,i;
  unsigned int *tb;
  N t;

  l = (bits+31)>>5; /* word length */
  *r = t = NALLOC(l);
  tb = BD(t);
  for ( i = 0; i < l; i++ )
    tb[i] = mt_genrand();
  if ( bits&31 )
    tb[l-1] &= (1<<(bits&31))-1;
  for ( i = l-1; i >= 0 && !tb[i]; i-- );
  if ( i < 0 )
    *r = 0;
  else
    PL(t) = i+1;
}

void freen(N n)
{
  if ( n && (n != ONEN) )
    free(n);
}

/* accepts Z */
int n_bits(N n)
{
  unsigned int i,t;
  int l;

  if ( !n )
    return 0;
  l = PL(n);
  if ( l < 0 ) l = -l;
  t = BD(n)[l-1];
  for ( i = 0; t; t>>=1, i++);
  return i + (l-1)*BSH;
}