[BACK]Return to ddM.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / asm

Annotation of OpenXM_contrib2/asir2000/asm/ddM.c, Revision 1.9

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.9     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/asm/ddM.c,v 1.8 2009/03/02 19:01:43 ohara Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "base.h"
                     52: #include "inline.h"
                     53:
                     54: /*
                     55:  * mod is declared as 'int', because several xxxum functions contains signed
                     56:  * integer addition/subtraction. So mod should be less than 2^31.
                     57:  */
                     58:
1.5       noro       59: void mulum(int mod,UM p1,UM p2,UM pr)
1.1       noro       60: {
1.9     ! noro       61:   int *pc1,*pcr;
        !            62:   int *c1,*c2,*cr;
        !            63:   unsigned int mul;
        !            64:   int i,j,d1,d2;
        !            65:
        !            66:   if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
        !            67:     DEG(pr) = -1;
        !            68:     return;
        !            69:   }
        !            70:   c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
        !            71:   bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
        !            72:   for ( i = 0; i <= d2; i++, cr++ )
        !            73:     if ( mul = *c2++ )
        !            74:       for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ ) {
        !            75:         DMAR(*pc1,mul,*pcr,mod,*pcr)
        !            76:       }
        !            77:   DEG(pr) = d1 + d2;
1.1       noro       78: }
                     79:
1.5       noro       80: void mulsum(int mod,UM p,int n,UM pr)
1.1       noro       81: {
1.9     ! noro       82:   int *sp,*dp;
        !            83:   int i;
1.1       noro       84:
1.9     ! noro       85:   for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
        !            86:       i >= 0; i--, dp--, sp-- ) {
        !            87:     DMAR(*sp,n,0,mod,*dp)
        !            88:   }
1.1       noro       89: }
1.9     ! noro       90:
1.5       noro       91: int divum(int mod,UM p1,UM p2,UM pq)
1.1       noro       92: {
1.9     ! noro       93:   int *pc1,*pct;
        !            94:   int *c1,*c2,*ct;
        !            95:   unsigned int inv,hd,tmp;
        !            96:   int i,j, d1,d2,dd;
        !            97:
        !            98:   if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
        !            99:     DEG(pq) = -1;
        !           100:     return d1;
        !           101:   }
        !           102:   c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
        !           103:   if ( ( hd = c2[d2] ) != 1 ) {
        !           104:     inv = invm(hd,mod);
        !           105:     for ( pc1 = c2 + d2; pc1 >= c2; pc1-- ) {
        !           106:       DMAR(*pc1,inv,0,mod,*pc1)
        !           107:     }
        !           108:   } else
        !           109:     inv = 1;
        !           110:   for ( i = dd, ct = c1+d1; i >= 0; i-- )
        !           111:     if ( tmp = *ct-- ) {
        !           112:       tmp = mod - tmp;
        !           113:       for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- ) {
        !           114:         DMAR(*pc1,tmp,*pct,mod,*pct)
        !           115:       }
        !           116:     }
        !           117:   if ( inv != 1 ) {
        !           118:     for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ ) {
        !           119:       DMAR(*pc1,inv,0,mod,*pc1)
        !           120:     }
        !           121:     for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ ) {
        !           122:       DMAR(*pc1,inv,0,mod,*pc1)
        !           123:     }
        !           124:   }
        !           125:   for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
        !           126:   for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
        !           127:     *pct-- = *pc1--;
        !           128:   return i;
1.1       noro      129: }
                    130:
1.5       noro      131: void diffum(int mod,UM f,UM fd)
1.1       noro      132: {
1.9     ! noro      133:   int *dp,*sp;
        !           134:   int i;
1.1       noro      135:
1.9     ! noro      136:   for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
        !           137:     i >= 1; i--, dp--, sp-- ) {
        !           138:     DMAR(*sp,i,0,mod,*dp)
        !           139:   }
        !           140:   degum(fd,DEG(f) - 1);
1.1       noro      141: }
                    142:
1.5       noro      143: unsigned int pwrm(int mod,int a,int n)
1.1       noro      144: {
1.9     ! noro      145:   unsigned int s,t;
1.1       noro      146:
1.9     ! noro      147:   if ( !n )
        !           148:     return 1;
        !           149:   else if ( n == 1 )
        !           150:     return a;
        !           151:   else {
        !           152:     t = pwrm(mod,a,n/2);
        !           153:     DMAR(t,t,0,mod,s)
        !           154:     if ( n % 2 ) {
        !           155:       DMAR(s,a,0,mod,t)
        !           156:       return t;
        !           157:     } else
        !           158:       return s;
        !           159:   }
1.1       noro      160: }
                    161:
1.5       noro      162: unsigned int invm(unsigned int s,int mod)
1.1       noro      163: {
1.9     ! noro      164:   unsigned int r,a2,q;
        !           165:   unsigned int f1,f2,a1;
1.1       noro      166:
1.9     ! noro      167:   for ( f1 = s, f2 = mod, a1 = 1, a2 = 0; ; ) {
        !           168:     q = f1/f2; r = f1 - f2*q; f1 = f2;
        !           169:     if ( !(f2 = r) )
        !           170:       break;
        !           171:     DMAR(a2,q,0,mod,r)
        !           172: /*    r = ( a1 - r + mod ) % mod; */
        !           173:     if ( a1 >= r )
        !           174:       r = a1 - r;
        !           175:     else {
        !           176:       r = (mod - r) + a1;
        !           177:     }
        !           178:     a1 = a2; a2 = r;
        !           179:   }
        !           180: /*  return( ( a2 >= 0 ? a2 : a2 + mod ) ); */
        !           181:   return a2;
1.1       noro      182: }
                    183:
1.5       noro      184: unsigned int rem(N n,int m)
1.1       noro      185: {
1.9     ! noro      186:   unsigned int *x;
        !           187:   unsigned int t,r;
        !           188:   int i;
        !           189:
        !           190:   if ( !n )
        !           191:     return 0;
        !           192:   for ( i = PL(n)-1, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
1.7       ohara     193: #if defined(sparc) && !defined(__sparcv9)
1.9     ! noro      194:     r = dsar(m,r,*x);
1.1       noro      195: #else
1.9     ! noro      196:     DSAB(m,r,*x,t,r)
1.1       noro      197: #endif
1.9     ! noro      198:   }
        !           199:   return r;
1.1       noro      200: }
                    201:
1.6       noro      202: #if !defined(sparc) || defined(__sparcv9)
1.5       noro      203: void addpadic(int mod,int n,unsigned int *n1,unsigned int *n2)
1.1       noro      204: {
1.9     ! noro      205:   unsigned int carry,tmp;
        !           206:   int i;
1.1       noro      207:
1.9     ! noro      208:   for ( i = 0, carry = 0; i < n; i++ ) {
        !           209:     tmp = *n1++ + *n2 + carry;
        !           210:     DQR(tmp,mod,carry,*n2++)
1.1       noro      211: /*
1.9     ! noro      212:     carry = tmp / mod;
        !           213:     *n2++ = tmp - ( carry * mod );
1.1       noro      214: */
1.9     ! noro      215:   }
1.1       noro      216: }
                    217: #endif
                    218:
1.5       noro      219: void mulpadic(int mod,int n,unsigned int *n1,unsigned int *n2,unsigned int *nr)
1.1       noro      220: {
1.9     ! noro      221:   unsigned int *pn1,*pnr;
        !           222:   unsigned int carry,mul;
        !           223:   int i,j;
        !           224:
        !           225:   bzero((char *)nr,(int)(n*sizeof(int)));
        !           226:   for ( j = 0; j < n; j++, n2++, nr++ )
        !           227:     if ( mul = *n2 )
        !           228:       for ( i = j, carry = 0, pn1 = n1, pnr = nr;
        !           229:         i < n; i++, pn1++, pnr++ ) {
        !           230:         carry += *pnr;
        !           231:         DMAB(mod,*pn1,mul,carry,carry,*pnr)
        !           232:       }
1.1       noro      233: }
                    234:
                    235: extern up_kara_mag;
                    236:
1.5       noro      237: void kmulum(int mod,UM n1,UM n2,UM nr)
1.1       noro      238: {
1.9     ! noro      239:   UM n,t,s,m,carry;
        !           240:   int d,d1,d2,len,i,l;
        !           241:   unsigned int *r,*r0;
        !           242:
        !           243:   if ( !n1 || !n2 ) {
        !           244:     nr->d = -1; return;
        !           245:   }
        !           246:   d1 = DEG(n1)+1; d2 = DEG(n2)+1;
        !           247:   if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {
        !           248:     mulum(mod,n1,n2,nr); return;
        !           249:   }
        !           250:   if ( d1 < d2 ) {
        !           251:     n = n1; n1 = n2; n2 = n;
        !           252:     d = d1; d1 = d2; d2 = d;
        !           253:   }
        !           254:   if ( d2 > (d1+1)/2 ) {
        !           255:     kmulummain(mod,n1,n2,nr); return;
        !           256:   }
        !           257:   d = (d1/d2)+((d1%d2)!=0);
        !           258:   len = (d+1)*d2;
        !           259:   r0 = (unsigned int *)ALLOCA(len*sizeof(int));
        !           260:   bzero((char *)r0,len*sizeof(int));
        !           261:   m = W_UMALLOC(d2+1);
        !           262:   carry = W_UMALLOC(d2+1);
        !           263:   t = W_UMALLOC(d1+d2+1);
        !           264:   s = W_UMALLOC(d1+d2+1);
        !           265:   for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {
        !           266:     extractum(n1,i*d2,d2,m);
        !           267:     if ( m ) {
        !           268:       kmulum(mod,m,n2,t);
        !           269:       addum(mod,t,carry,s);
        !           270:       c_copyum(s,d2,r);
        !           271:       extractum(s,d2,d2,carry);
        !           272:     } else {
        !           273:       c_copyum(carry,d2,r);
        !           274:       carry = 0;
        !           275:     }
        !           276:   }
        !           277:   c_copyum(carry,d2,r);
        !           278:   for ( l = len - 1; !r0[l]; l-- );
        !           279:   l++;
        !           280:   DEG(nr) = l-1;
        !           281:   bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
1.1       noro      282: }
                    283:
1.5       noro      284: void ksquareum(int mod,UM n1,UM nr)
1.1       noro      285: {
1.9     ! noro      286:   int d1;
1.1       noro      287:
1.9     ! noro      288:   if ( !n1 ) {
        !           289:     nr->d = -1; return;
        !           290:   }
        !           291:   d1 = DEG(n1)+1;
        !           292:   if ( (d1 < up_kara_mag) ) {
        !           293:     pwrum(mod,n1,2,nr); return;
        !           294:   }
        !           295:   ksquareummain(mod,n1,nr);
1.1       noro      296: }
                    297:
1.5       noro      298: void extractum(UM n,int index,int len,UM nr)
1.1       noro      299: {
1.9     ! noro      300:   int *m;
        !           301:   int l;
1.1       noro      302:
1.9     ! noro      303:   if ( !n ) {
        !           304:     nr->d = -1; return;
        !           305:   }
        !           306:   m = COEF(n)+index;
        !           307:   if ( (l = (DEG(n)+1)-index) >= len ) {
        !           308:     for ( l = len - 1; (l >= 0) && !m[l]; l-- );
        !           309:     l++;
        !           310:   }
        !           311:   if ( l <= 0 ) {
        !           312:     nr->d = -1; return;
        !           313:   } else {
        !           314:     DEG(nr) = l-1;
        !           315:     bcopy((char *)m,(char *)COEF(nr),l*sizeof(Q));
        !           316:   }
1.1       noro      317: }
                    318:
1.5       noro      319: void copyum(UM n1,UM n2)
1.1       noro      320: {
1.9     ! noro      321:   n2->d = n1->d;
        !           322:   bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(int));
1.1       noro      323: }
                    324:
1.5       noro      325: void c_copyum(UM n,int len,int *p)
1.1       noro      326: {
1.9     ! noro      327:   if ( n )
        !           328:     bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(int));
1.1       noro      329: }
                    330:
1.5       noro      331: void kmulummain(int mod,UM n1,UM n2,UM nr)
1.1       noro      332: {
1.9     ! noro      333:   int d1,d2,h,len;
        !           334:   UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
1.1       noro      335:
1.9     ! noro      336:   d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
        !           337:   n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
        !           338:   n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);
        !           339:   lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);
        !           340:   mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);
        !           341:   mid = W_UMALLOC(d1+d2+1);
        !           342:   s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);
        !           343:   extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
        !           344:   extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
        !           345:   kmulum(mod,n1hi,n2hi,hi); kmulum(mod,n1lo,n2lo,lo);
        !           346:   len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
        !           347:   bzero((char *)COEF(t1),len*sizeof(int));
        !           348:   if ( lo )
        !           349:     bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
        !           350:   if ( hi )
        !           351:     bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
        !           352:
        !           353:   addum(mod,hi,lo,mid1);
        !           354:   subum(mod,n1hi,n1lo,s1); subum(mod,n2lo,n2hi,s2); kmulum(mod,s1,s2,mid2);
        !           355:   addum(mod,mid1,mid2,mid);
        !           356:   if ( mid ) {
        !           357:     len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
        !           358:     bzero((char *)COEF(t2),len*sizeof(int));
        !           359:     bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
        !           360:     addum(mod,t1,t2,nr);
        !           361:   } else
        !           362:     copyum(t1,nr);
1.1       noro      363: }
                    364:
1.5       noro      365: void ksquareummain(int mod,UM n1,UM nr)
1.1       noro      366: {
1.9     ! noro      367:   int d1,h,len;
        !           368:   UM n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
1.1       noro      369:
1.9     ! noro      370:   d1 = DEG(n1)+1; h = (d1+1)/2;
        !           371:   n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
        !           372:   lo = W_UMALLOC(2*d1+1); hi = W_UMALLOC(2*d1+1);
        !           373:   mid1 = W_UMALLOC(2*d1+1); mid2 = W_UMALLOC(2*d1+1);
        !           374:   mid = W_UMALLOC(2*d1+1);
        !           375:   s1 = W_UMALLOC(2*d1+1);
        !           376:   extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
        !           377:   ksquareum(mod,n1hi,hi); ksquareum(mod,n1lo,lo);
        !           378:   len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
        !           379:   bzero((char *)COEF(t1),len*sizeof(int));
        !           380:   if ( lo )
        !           381:     bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
        !           382:   if ( hi )
        !           383:     bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
        !           384:
        !           385:   addum(mod,hi,lo,mid1);
        !           386:   subum(mod,n1hi,n1lo,s1); ksquareum(mod,s1,mid2);
        !           387:   subum(mod,mid1,mid2,mid);
        !           388:   if ( mid ) {
        !           389:     len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
        !           390:     bzero((char *)COEF(t2),len*sizeof(int));
        !           391:     bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
        !           392:     addum(mod,t1,t2,nr);
        !           393:   } else
        !           394:     copyum(t1,nr);
1.1       noro      395: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>