xref: /netbsd-src/external/lgpl3/gmp/dist/mpn/cray/mulww.f (revision d536862b7d93d77932ef5de7eebdc48d76921b77)
1c    Helper for mpn_mul_1, mpn_addmul_1, and mpn_submul_1 for Cray PVP.
2
3c    Copyright 1996, 2000 Free Software Foundation, Inc.
4
5c    This file is part of the GNU MP Library.
6c
7c    The GNU MP Library is free software; you can redistribute it and/or modify
8c    it under the terms of either:
9c
10c      * the GNU Lesser General Public License as published by the Free
11c        Software Foundation; either version 3 of the License, or (at your
12c        option) any later version.
13c
14c    or
15c
16c      * the GNU General Public License as published by the Free Software
17c        Foundation; either version 2 of the License, or (at your option) any
18c        later version.
19c
20c    or both in parallel, as here.
21c
22c    The GNU MP Library is distributed in the hope that it will be useful, but
23c    WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24c    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
25c    for more details.
26c
27c    You should have received copies of the GNU General Public License and the
28c    GNU Lesser General Public License along with the GNU MP Library.  If not,
29c    see https://www.gnu.org/licenses/.
30
31c    p1[] = hi(a[]*s); the upper limbs of each product
32c    p0[] = low(a[]*s); the corresponding lower limbs
33c    n is number of limbs in the vectors
34
35      subroutine gmpn_mulww(p1,p0,a,n,s)
36      integer*8 p1(0:*),p0(0:*),a(0:*),s
37      integer n
38
39      integer*8 a0,a1,a2,s0,s1,s2,c
40      integer*8 ai,t0,t1,t2,t3,t4
41
42      s0 = shiftl(and(s,4194303),24)
43      s1 = shiftl(and(shiftr(s,22),4194303),24)
44      s2 = shiftl(and(shiftr(s,44),4194303),24)
45
46      do i = 0,n-1
47         ai = a(i)
48         a0 = shiftl(and(ai,4194303),24)
49         a1 = shiftl(and(shiftr(ai,22),4194303),24)
50         a2 = shiftl(and(shiftr(ai,44),4194303),24)
51
52         t0 = i24mult(a0,s0)
53         t1 = i24mult(a0,s1)+i24mult(a1,s0)
54         t2 = i24mult(a0,s2)+i24mult(a1,s1)+i24mult(a2,s0)
55         t3 = i24mult(a1,s2)+i24mult(a2,s1)
56         t4 = i24mult(a2,s2)
57
58         p0(i)=shiftl(t2,44)+shiftl(t1,22)+t0
59         c=shiftr(shiftr(t0,22)+and(t1,4398046511103)+
60     $        shiftl(and(t2,1048575),22),42)
61         p1(i)=shiftl(t4,24)+shiftl(t3,2)+shiftr(t2,20)+shiftr(t1,42)+c
62      end do
63      end
64