[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.4

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.4     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/asm/ddM.c,v 1.3 2000/08/22 05:03:55 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "base.h"
                     52: #include "inline.h"
                     53:
                     54: void ksquareummain(int,UM,UM);
                     55: void kmulummain(int,UM,UM,UM);
                     56: void c_copyum(UM,int,int *);
                     57: void copyum(UM,UM);
                     58: void extractum(UM,int,int,UM);
                     59: void ksquareum(int,UM,UM);
                     60: void kmulum(int,UM,UM,UM);
                     61:
                     62: /*
                     63:  * mod is declared as 'int', because several xxxum functions contains signed
                     64:  * integer addition/subtraction. So mod should be less than 2^31.
                     65:  */
                     66:
                     67: void mulum(mod,p1,p2,pr)
                     68: int mod;
                     69: UM p1,p2,pr;
                     70: {
                     71:        int *pc1,*pcr;
                     72:        int *c1,*c2,*cr;
                     73:        unsigned int mul;
                     74:        int i,j,d1,d2;
                     75:
                     76:        if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
                     77:                DEG(pr) = -1;
                     78:                return;
                     79:        }
                     80:        c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
                     81:        bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
                     82:        for ( i = 0; i <= d2; i++, cr++ )
                     83:                if ( mul = *c2++ )
                     84:                        for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ ) {
                     85:                                DMAR(*pc1,mul,*pcr,mod,*pcr)
                     86:                        }
                     87:        DEG(pr) = d1 + d2;
                     88: }
                     89:
                     90: void mulsum(mod,p,n,pr)
                     91: int mod,n;
                     92: UM p,pr;
                     93: {
                     94:        int *sp,*dp;
                     95:        int i;
                     96:
                     97:        for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
                     98:                  i >= 0; i--, dp--, sp-- ) {
                     99:                DMAR(*sp,n,0,mod,*dp)
                    100:        }
                    101: }
                    102:
                    103: int divum(mod,p1,p2,pq)
                    104: int mod;
                    105: UM p1,p2,pq;
                    106: {
                    107:        int *pc1,*pct;
                    108:        int *c1,*c2,*ct;
                    109:        unsigned int inv,hd,tmp;
                    110:        int i,j, d1,d2,dd;
                    111:
                    112:        if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
                    113:                DEG(pq) = -1;
                    114:                return d1;
                    115:        }
                    116:        c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
                    117:        if ( ( hd = c2[d2] ) != 1 ) {
                    118:                inv = invm(hd,mod);
                    119:                for ( pc1 = c2 + d2; pc1 >= c2; pc1-- ) {
                    120:                        DMAR(*pc1,inv,0,mod,*pc1)
                    121:                }
                    122:        } else
                    123:                inv = 1;
                    124:        for ( i = dd, ct = c1+d1; i >= 0; i-- )
                    125:                if ( tmp = *ct-- ) {
                    126:                        tmp = mod - tmp;
                    127:                        for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- ) {
                    128:                                DMAR(*pc1,tmp,*pct,mod,*pct)
                    129:                        }
                    130:                }
                    131:        if ( inv != 1 ) {
                    132:                for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ ) {
                    133:                        DMAR(*pc1,inv,0,mod,*pc1)
                    134:                }
                    135:                for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ ) {
                    136:                        DMAR(*pc1,inv,0,mod,*pc1)
                    137:                }
                    138:        }
                    139:        for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
                    140:        for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
                    141:                *pct-- = *pc1--;
                    142:        return i;
                    143: }
                    144:
                    145: void diffum(mod,f,fd)
                    146: int mod;
                    147: UM f,fd;
                    148: {
                    149:        int *dp,*sp;
                    150:        int i;
                    151:        UL ltmp;
                    152:
                    153:        for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
                    154:                i >= 1; i--, dp--, sp-- ) {
                    155:                DMAR(*sp,i,0,mod,*dp)
                    156:        }
                    157:        degum(fd,DEG(f) - 1);
                    158: }
                    159:
                    160: unsigned int pwrm(mod,a,n)
                    161: int mod,a;
                    162: int n;
                    163: {
                    164:        unsigned int s,t;
                    165:
                    166:        if ( !n )
                    167:                return 1;
                    168:        else if ( n == 1 )
                    169:                return a;
                    170:        else {
                    171:                t = pwrm(mod,a,n/2);
                    172:                DMAR(t,t,0,mod,s)
                    173:                if ( n % 2 ) {
                    174:                        DMAR(s,a,0,mod,t)
                    175:                        return t;
                    176:                } else
                    177:                        return s;
                    178:        }
                    179: }
                    180:
                    181: unsigned int invm(s,mod)
                    182: unsigned int s;
                    183: int mod;
                    184: {
                    185:        unsigned int r,a2,q;
                    186:        unsigned int f1,f2,a1;
                    187:
                    188:        for ( f1 = s, f2 = mod, a1 = 1, a2 = 0; ; ) {
                    189:                q = f1/f2; r = f1 - f2*q; f1 = f2;
                    190:                if ( !(f2 = r) )
                    191:                        break;
                    192:                DMAR(a2,q,0,mod,r)
                    193: /*             r = ( a1 - r + mod ) % mod; */
                    194:                if ( a1 >= r )
                    195:                        r = a1 - r;
                    196:                else {
                    197:                        r = (mod - r) + a1;
                    198:                }
                    199:                a1 = a2; a2 = r;
                    200:        }
                    201: /*     return( ( a2 >= 0 ? a2 : a2 + mod ) ); */
                    202:        return a2;
                    203: }
                    204:
                    205: unsigned int rem(n,m)
                    206: N n;
                    207: unsigned int m;
                    208: {
                    209:        unsigned int *x;
                    210:        unsigned int t,r;
                    211:        int i;
                    212:
                    213:        if ( !n )
                    214:                return 0;
                    215:        for ( i = PL(n)-1, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
                    216: #if defined(sparc)
                    217:                r = dsar(m,r,*x);
                    218: #else
                    219:                DSAB(m,r,*x,t,r)
                    220: #endif
                    221:        }
                    222:        return r;
                    223: }
                    224:
                    225: #ifndef sparc
                    226: void addpadic(mod,n,n1,n2)
                    227: int mod;
                    228: int n;
                    229: unsigned int *n1,*n2;
                    230: {
                    231:        unsigned int carry,tmp;
                    232:        int i;
                    233:
                    234:        for ( i = 0, carry = 0; i < n; i++ ) {
                    235:                tmp = *n1++ + *n2 + carry;
                    236:                DQR(tmp,mod,carry,*n2++)
                    237: /*
                    238:                carry = tmp / mod;
                    239:                *n2++ = tmp - ( carry * mod );
                    240: */
                    241:        }
                    242: }
                    243: #endif
                    244:
                    245: void mulpadic(mod,n,n1,n2,nr)
                    246: int mod;
                    247: int n;
                    248: unsigned int *n1;
                    249: unsigned int *n2,*nr;
                    250: {
                    251:        unsigned int *pn1,*pnr;
                    252:        unsigned int carry,mul;
                    253:        int i,j;
                    254:
                    255:        bzero((char *)nr,(int)(n*sizeof(int)));
                    256:        for ( j = 0; j < n; j++, n2++, nr++ )
                    257:                if ( mul = *n2 )
                    258:                        for ( i = j, carry = 0, pn1 = n1, pnr = nr;
                    259:                                i < n; i++, pn1++, pnr++ ) {
                    260:                                carry += *pnr;
                    261:                                DMAB(mod,*pn1,mul,carry,carry,*pnr)
                    262:                        }
                    263: }
                    264:
                    265: extern up_kara_mag;
                    266:
                    267: void kmulum(mod,n1,n2,nr)
                    268: UM n1,n2,nr;
                    269: {
                    270:        UM n,t,s,m,carry;
                    271:        int d,d1,d2,len,i,l;
                    272:        unsigned int *r,*r0;
                    273:
                    274:        if ( !n1 || !n2 ) {
                    275:                nr->d = -1; return;
                    276:        }
                    277:        d1 = DEG(n1)+1; d2 = DEG(n2)+1;
                    278:        if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {
                    279:                mulum(mod,n1,n2,nr); return;
                    280:        }
                    281:        if ( d1 < d2 ) {
                    282:                n = n1; n1 = n2; n2 = n;
                    283:                d = d1; d1 = d2; d2 = d;
                    284:        }
                    285:        if ( d2 > (d1+1)/2 ) {
                    286:                kmulummain(mod,n1,n2,nr); return;
                    287:        }
                    288:        d = (d1/d2)+((d1%d2)!=0);
                    289:        len = (d+1)*d2;
                    290:        r0 = (unsigned int *)ALLOCA(len*sizeof(int));
                    291:        bzero((char *)r0,len*sizeof(int));
                    292:        m = W_UMALLOC(d2+1);
                    293:        carry = W_UMALLOC(d2+1);
                    294:        t = W_UMALLOC(d1+d2+1);
                    295:        s = W_UMALLOC(d1+d2+1);
1.4     ! noro      296:        for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {
1.1       noro      297:                extractum(n1,i*d2,d2,m);
                    298:                if ( m ) {
                    299:                        kmulum(mod,m,n2,t);
                    300:                        addum(mod,t,carry,s);
                    301:                        c_copyum(s,d2,r);
                    302:                        extractum(s,d2,d2,carry);
                    303:                } else {
                    304:                        c_copyum(carry,d2,r);
                    305:                        carry = 0;
                    306:                }
                    307:        }
                    308:        c_copyum(carry,d2,r);
                    309:        for ( l = len - 1; !r0[l]; l-- );
                    310:        l++;
                    311:        DEG(nr) = l-1;
                    312:        bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
                    313: }
                    314:
                    315: void ksquareum(mod,n1,nr)
                    316: int mod;
                    317: UM n1,nr;
                    318: {
                    319:        int d1;
                    320:
                    321:        if ( !n1 ) {
                    322:                nr->d = -1; return;
                    323:        }
                    324:        d1 = DEG(n1)+1;
                    325:        if ( (d1 < up_kara_mag) ) {
                    326:                pwrum(mod,n1,2,nr); return;
                    327:        }
                    328:        ksquareummain(mod,n1,nr);
                    329: }
                    330:
                    331: void extractum(n,index,len,nr)
                    332: UM n;
                    333: int index,len;
                    334: UM nr;
                    335: {
                    336:        int *m;
                    337:        int l;
                    338:
                    339:        if ( !n ) {
                    340:                nr->d = -1; return;
                    341:        }
                    342:        m = COEF(n)+index;
                    343:        if ( (l = (DEG(n)+1)-index) >= len ) {
                    344:                for ( l = len - 1; (l >= 0) && !m[l]; l-- );
                    345:                l++;
                    346:        }
                    347:        if ( l <= 0 ) {
                    348:                nr->d = -1; return;
                    349:        } else {
                    350:                DEG(nr) = l-1;
                    351:                bcopy((char *)m,(char *)COEF(nr),l*sizeof(Q));
                    352:        }
                    353: }
                    354:
                    355: void copyum(n1,n2)
                    356: UM n1,n2;
                    357: {
                    358:        n2->d = n1->d;
                    359:        bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(int));
                    360: }
                    361:
                    362: void c_copyum(n,len,p)
                    363: UM n;
                    364: int len;
                    365: int *p;
                    366: {
                    367:        if ( n )
                    368:                bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(int));
                    369: }
                    370:
                    371: void kmulummain(mod,n1,n2,nr)
                    372: int mod;
                    373: UM n1,n2,nr;
                    374: {
                    375:        int d1,d2,h,len;
                    376:        UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
                    377:
                    378:        d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
                    379:        n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
                    380:        n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);
                    381:        lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);
                    382:        mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);
                    383:        mid = W_UMALLOC(d1+d2+1);
                    384:        s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);
                    385:        extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
                    386:        extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
                    387:        kmulum(mod,n1hi,n2hi,hi); kmulum(mod,n1lo,n2lo,lo);
                    388:        len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
                    389:        bzero((char *)COEF(t1),len*sizeof(int));
                    390:        if ( lo )
                    391:                bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
                    392:        if ( hi )
                    393:                bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
                    394:
                    395:        addum(mod,hi,lo,mid1);
                    396:        subum(mod,n1hi,n1lo,s1); subum(mod,n2lo,n2hi,s2); kmulum(mod,s1,s2,mid2);
                    397:        addum(mod,mid1,mid2,mid);
                    398:        if ( mid ) {
                    399:                len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
                    400:                bzero((char *)COEF(t2),len*sizeof(int));
                    401:                bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
                    402:                addum(mod,t1,t2,nr);
                    403:        } else
                    404:                copyum(t1,nr);
                    405: }
                    406:
                    407: void ksquareummain(mod,n1,nr)
                    408: int mod;
                    409: UM n1,nr;
                    410: {
                    411:        int d1,h,len;
                    412:        UM n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
                    413:
                    414:        d1 = DEG(n1)+1; h = (d1+1)/2;
                    415:        n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
                    416:        lo = W_UMALLOC(2*d1+1); hi = W_UMALLOC(2*d1+1);
                    417:        mid1 = W_UMALLOC(2*d1+1); mid2 = W_UMALLOC(2*d1+1);
                    418:        mid = W_UMALLOC(2*d1+1);
                    419:        s1 = W_UMALLOC(2*d1+1);
                    420:        extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
                    421:        ksquareum(mod,n1hi,hi); ksquareum(mod,n1lo,lo);
                    422:        len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
                    423:        bzero((char *)COEF(t1),len*sizeof(int));
                    424:        if ( lo )
                    425:                bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
                    426:        if ( hi )
                    427:                bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
                    428:
                    429:        addum(mod,hi,lo,mid1);
                    430:        subum(mod,n1hi,n1lo,s1); ksquareum(mod,s1,mid2);
                    431:        subum(mod,mid1,mid2,mid);
                    432:        if ( mid ) {
                    433:                len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
                    434:                bzero((char *)COEF(t2),len*sizeof(int));
                    435:                bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
                    436:                addum(mod,t1,t2,nr);
                    437:        } else
                    438:                copyum(t1,nr);
                    439: }

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