0000: 2f 2a 20 73 74 61 74 6c 69 62 2e 63 20 2d 2d 20 /* statlib.c --
0010: 53 74 61 74 69 73 74 69 63 61 6c 20 66 75 6e 63 Statistical func
0020: 74 69 6f 6e 73 20 66 6f 72 20 74 65 73 74 69 6e tions for testin
0030: 67 20 74 68 65 20 72 61 6e 64 6f 6d 6e 65 73 73 g the randomness
0040: 20 6f 66 0a 20 20 20 6e 75 6d 62 65 72 20 73 65 of. number se
0050: 71 75 65 6e 63 65 73 2e 20 2a 2f 0a 0a 2f 2a 0a quences. */../*.
0060: 43 6f 70 79 72 69 67 68 74 20 31 39 39 39 2c 20 Copyright 1999,
0070: 32 30 30 30 20 46 72 65 65 20 53 6f 66 74 77 61 2000 Free Softwa
0080: 72 65 20 46 6f 75 6e 64 61 74 69 6f 6e 2c 20 49 re Foundation, I
0090: 6e 63 2e 0a 0a 54 68 69 73 20 66 69 6c 65 20 69 nc...This file i
00a0: 73 20 70 61 72 74 20 6f 66 20 74 68 65 20 47 4e s part of the GN
00b0: 55 20 4d 50 20 4c 69 62 72 61 72 79 2e 0a 0a 54 U MP Library...T
00c0: 68 65 20 47 4e 55 20 4d 50 20 4c 69 62 72 61 72 he GNU MP Librar
00d0: 79 20 69 73 20 66 72 65 65 20 73 6f 66 74 77 61 y is free softwa
00e0: 72 65 3b 20 79 6f 75 20 63 61 6e 20 72 65 64 69 re; you can redi
00f0: 73 74 72 69 62 75 74 65 20 69 74 20 61 6e 64 2f stribute it and/
0100: 6f 72 20 6d 6f 64 69 66 79 0a 69 74 20 75 6e 64 or modify.it und
0110: 65 72 20 74 68 65 20 74 65 72 6d 73 20 6f 66 20 er the terms of
0120: 74 68 65 20 47 4e 55 20 4c 65 73 73 65 72 20 47 the GNU Lesser G
0130: 65 6e 65 72 61 6c 20 50 75 62 6c 69 63 20 4c 69 eneral Public Li
0140: 63 65 6e 73 65 20 61 73 20 70 75 62 6c 69 73 68 cense as publish
0150: 65 64 20 62 79 0a 74 68 65 20 46 72 65 65 20 53 ed by.the Free S
0160: 6f 66 74 77 61 72 65 20 46 6f 75 6e 64 61 74 69 oftware Foundati
0170: 6f 6e 3b 20 65 69 74 68 65 72 20 76 65 72 73 69 on; either versi
0180: 6f 6e 20 32 2e 31 20 6f 66 20 74 68 65 20 4c 69 on 2.1 of the Li
0190: 63 65 6e 73 65 2c 20 6f 72 20 28 61 74 20 79 6f cense, or (at yo
01a0: 75 72 0a 6f 70 74 69 6f 6e 29 20 61 6e 79 20 6c ur.option) any l
01b0: 61 74 65 72 20 76 65 72 73 69 6f 6e 2e 0a 0a 54 ater version...T
01c0: 68 65 20 47 4e 55 20 4d 50 20 4c 69 62 72 61 72 he GNU MP Librar
01d0: 79 20 69 73 20 64 69 73 74 72 69 62 75 74 65 64 y is distributed
01e0: 20 69 6e 20 74 68 65 20 68 6f 70 65 20 74 68 61 in the hope tha
01f0: 74 20 69 74 20 77 69 6c 6c 20 62 65 20 75 73 65 t it will be use
0200: 66 75 6c 2c 20 62 75 74 0a 57 49 54 48 4f 55 54 ful, but.WITHOUT
0210: 20 41 4e 59 20 57 41 52 52 41 4e 54 59 3b 20 77 ANY WARRANTY; w
0220: 69 74 68 6f 75 74 20 65 76 65 6e 20 74 68 65 20 ithout even the
0230: 69 6d 70 6c 69 65 64 20 77 61 72 72 61 6e 74 79 implied warranty
0240: 20 6f 66 20 4d 45 52 43 48 41 4e 54 41 42 49 4c of MERCHANTABIL
0250: 49 54 59 0a 6f 72 20 46 49 54 4e 45 53 53 20 46 ITY.or FITNESS F
0260: 4f 52 20 41 20 50 41 52 54 49 43 55 4c 41 52 20 OR A PARTICULAR
0270: 50 55 52 50 4f 53 45 2e 20 20 53 65 65 20 74 68 PURPOSE. See th
0280: 65 20 47 4e 55 20 4c 65 73 73 65 72 20 47 65 6e e GNU Lesser Gen
0290: 65 72 61 6c 20 50 75 62 6c 69 63 0a 4c 69 63 65 eral Public.Lice
02a0: 6e 73 65 20 66 6f 72 20 6d 6f 72 65 20 64 65 74 nse for more det
02b0: 61 69 6c 73 2e 0a 0a 59 6f 75 20 73 68 6f 75 6c ails...You shoul
02c0: 64 20 68 61 76 65 20 72 65 63 65 69 76 65 64 20 d have received
02d0: 61 20 63 6f 70 79 20 6f 66 20 74 68 65 20 47 4e a copy of the GN
02e0: 55 20 4c 65 73 73 65 72 20 47 65 6e 65 72 61 6c U Lesser General
02f0: 20 50 75 62 6c 69 63 20 4c 69 63 65 6e 73 65 0a Public License.
0300: 61 6c 6f 6e 67 20 77 69 74 68 20 74 68 65 20 47 along with the G
0310: 4e 55 20 4d 50 20 4c 69 62 72 61 72 79 3b 20 73 NU MP Library; s
0320: 65 65 20 74 68 65 20 66 69 6c 65 20 43 4f 50 59 ee the file COPY
0330: 49 4e 47 2e 4c 49 42 2e 20 20 49 66 20 6e 6f 74 ING.LIB. If not
0340: 2c 20 77 72 69 74 65 20 74 6f 0a 74 68 65 20 46 , write to.the F
0350: 72 65 65 20 53 6f 66 74 77 61 72 65 20 46 6f 75 ree Software Fou
0360: 6e 64 61 74 69 6f 6e 2c 20 49 6e 63 2e 2c 20 35 ndation, Inc., 5
0370: 39 20 54 65 6d 70 6c 65 20 50 6c 61 63 65 20 2d 9 Temple Place -
0380: 20 53 75 69 74 65 20 33 33 30 2c 20 42 6f 73 74 Suite 330, Bost
0390: 6f 6e 2c 0a 4d 41 20 30 32 31 31 31 2d 31 33 30 on,.MA 02111-130
03a0: 37 2c 20 55 53 41 2e 0a 2a 2f 0a 0a 2f 2a 20 54 7, USA..*/../* T
03b0: 68 65 20 74 68 65 6f 72 69 65 73 20 66 6f 72 20 he theories for
03c0: 74 68 65 73 65 20 66 75 6e 63 74 69 6f 6e 73 20 these functions
03d0: 61 72 65 20 74 61 6b 65 6e 20 66 72 6f 6d 20 44 are taken from D
03e0: 2e 20 4b 6e 75 74 68 27 73 20 22 54 68 65 20 41 . Knuth's "The A
03f0: 72 74 0a 6f 66 20 43 6f 6d 70 75 74 65 72 20 50 rt.of Computer P
0400: 72 6f 67 72 61 6d 6d 69 6e 67 3a 20 56 6f 6c 75 rogramming: Volu
0410: 6d 65 20 32 2c 20 53 65 6d 69 6e 75 6d 65 72 69 me 2, Seminumeri
0420: 63 61 6c 20 41 6c 67 6f 72 69 74 68 6d 73 22 2c cal Algorithms",
0430: 20 54 68 69 72 64 0a 45 64 69 74 69 6f 6e 2c 20 Third.Edition,
0440: 41 64 64 69 73 6f 6e 20 57 65 73 6c 65 79 2c 20 Addison Wesley,
0450: 31 39 39 38 2e 20 2a 2f 0a 0a 2f 2a 20 49 6d 70 1998. */../* Imp
0460: 6c 65 6d 65 6e 74 61 74 69 6f 6e 20 6e 6f 74 65 lementation note
0470: 73 2e 0a 0a 54 68 65 20 4b 6f 6c 6d 6f 67 6f 72 s...The Kolmogor
0480: 6f 76 2d 53 6d 69 72 6e 6f 76 20 74 65 73 74 2e ov-Smirnov test.
0490: 0a 0a 45 71 2e 20 28 31 33 29 20 69 6e 20 4b 6e ..Eq. (13) in Kn
04a0: 75 74 68 2c 20 70 2e 20 35 30 2c 20 73 61 79 73 uth, p. 50, says
04b0: 20 74 68 61 74 20 69 66 20 58 31 2c 20 58 32 2c that if X1, X2,
04c0: 20 2e 2e 2e 2c 20 58 6e 20 61 72 65 20 69 6e 64 ..., Xn are ind
04d0: 65 70 65 6e 64 65 6e 74 0a 6f 62 73 65 72 76 61 ependent.observa
04e0: 74 69 6f 6e 73 20 61 72 72 61 6e 67 65 64 20 69 tions arranged i
04f0: 6e 74 6f 20 61 73 63 65 6e 64 69 6e 67 20 6f 72 nto ascending or
0500: 64 65 72 0a 0a 09 4b 70 20 3d 20 73 71 72 28 6e der...Kp = sqr(n
0510: 29 20 2a 20 6d 61 78 28 6a 2f 6e 20 2d 20 46 28 ) * max(j/n - F(
0520: 58 6a 29 29 09 09 66 6f 72 20 61 6c 6c 20 31 3c Xj))..for all 1<
0530: 3d 6a 3c 3d 6e 0a 09 4b 6d 20 3d 20 73 71 72 28 =j<=n..Km = sqr(
0540: 6e 29 20 2a 20 6d 61 78 28 46 28 58 6a 29 20 2d n) * max(F(Xj) -
0550: 20 28 6a 2d 31 29 2f 6e 29 29 09 66 6f 72 20 61 (j-1)/n)).for a
0560: 6c 6c 20 31 3c 3d 6a 3c 3d 6e 0a 0a 77 68 65 72 ll 1<=j<=n..wher
0570: 65 20 46 28 78 29 20 3d 20 50 72 28 58 20 3c 3d e F(x) = Pr(X <=
0580: 20 78 29 20 3d 20 70 72 6f 62 61 62 69 6c 69 74 x) = probabilit
0590: 79 20 74 68 61 74 20 28 58 20 3c 3d 20 78 29 2c y that (X <= x),
05a0: 20 77 68 69 63 68 20 66 6f 72 20 61 0a 75 6e 69 which for a.uni
05b0: 66 6f 72 6d 6c 79 20 64 69 73 74 72 69 62 75 74 formly distribut
05c0: 65 64 20 72 61 6e 64 6f 6d 20 72 65 61 6c 20 6e ed random real n
05d0: 75 6d 62 65 72 20 62 65 74 77 65 65 6e 20 7a 65 umber between ze
05e0: 72 6f 20 61 6e 64 20 6f 6e 65 20 69 73 0a 65 78 ro and one is.ex
05f0: 61 63 74 6c 79 20 74 68 65 20 6e 75 6d 62 65 72 actly the number
0600: 20 69 74 73 65 6c 66 20 28 78 29 2e 0a 0a 0a 54 itself (x)....T
0610: 68 65 20 61 6e 73 77 65 72 20 74 6f 20 65 78 65 he answer to exe
0620: 72 63 69 73 65 20 32 33 20 67 69 76 65 73 20 74 rcise 23 gives t
0630: 68 65 20 66 6f 6c 6c 6f 77 69 6e 67 20 69 6d 70 he following imp
0640: 6c 65 6d 65 6e 74 61 74 69 6f 6e 2c 20 77 68 69 lementation, whi
0650: 63 68 0a 64 6f 65 73 6e 27 74 20 6e 65 65 64 20 ch.doesn't need
0660: 74 68 65 20 6f 62 73 65 72 76 61 74 69 6f 6e 73 the observations
0670: 20 74 6f 20 62 65 20 73 6f 72 74 65 64 20 69 6e to be sorted in
0680: 20 61 73 63 65 6e 64 69 6e 67 20 6f 72 64 65 72 ascending order
0690: 3a 0a 0a 66 6f 72 20 28 6b 20 3d 20 30 3b 20 6b :..for (k = 0; k
06a0: 20 3c 20 6d 3b 20 6b 2b 2b 29 0a 09 61 5b 6b 5d < m; k++)..a[k]
06b0: 20 3d 20 31 2e 30 0a 09 62 5b 6b 5d 20 3d 20 30 = 1.0..b[k] = 0
06c0: 2e 30 0a 09 63 5b 6b 5d 20 3d 20 30 0a 0a 66 6f .0..c[k] = 0..fo
06d0: 72 20 28 65 61 63 68 20 6f 62 73 65 72 76 61 74 r (each observat
06e0: 69 6f 6e 20 58 6a 29 0a 09 59 20 3d 20 46 28 58 ion Xj)..Y = F(X
06f0: 6a 29 0a 09 6b 20 3d 20 66 6c 6f 6f 72 20 28 6d j)..k = floor (m
0700: 20 2a 20 59 29 0a 09 61 5b 6b 5d 20 3d 20 6d 69 * Y)..a[k] = mi
0710: 6e 20 28 61 5b 6b 5d 2c 20 59 29 0a 09 62 5b 6b n (a[k], Y)..b[k
0720: 5d 20 3d 20 6d 61 78 20 28 62 5b 6b 5d 2c 20 59 ] = max (b[k], Y
0730: 29 0a 09 63 5b 6b 5d 20 2b 3d 20 31 0a 0a 09 6a )..c[k] += 1...j
0740: 20 3d 20 30 0a 09 72 70 20 3d 20 72 6d 20 3d 20 = 0..rp = rm =
0750: 30 0a 09 66 6f 72 20 28 6b 20 3d 20 30 3b 20 6b 0..for (k = 0; k
0760: 20 3c 20 6d 3b 20 6b 2b 2b 29 0a 09 09 69 66 20 < m; k++)...if
0770: 28 63 5b 6b 5d 20 3e 20 30 29 0a 09 09 09 72 6d (c[k] > 0)....rm
0780: 20 3d 20 6d 61 78 20 28 72 6d 2c 20 61 5b 6b 5d = max (rm, a[k]
0790: 20 2d 20 6a 2f 6e 29 0a 09 09 09 6a 20 2b 3d 20 - j/n)....j +=
07a0: 63 5b 6b 5d 0a 09 09 09 72 70 20 3d 20 6d 61 78 c[k]....rp = max
07b0: 20 28 72 70 2c 20 6a 2f 6e 20 2d 20 62 5b 6b 5d (rp, j/n - b[k]
07c0: 29 0a 0a 4b 70 20 3d 20 73 71 72 20 28 6e 29 20 )..Kp = sqr (n)
07d0: 2a 20 72 70 0a 4b 6d 20 3d 20 73 71 72 20 28 6e * rp.Km = sqr (n
07e0: 29 20 2a 20 72 6d 20 0a 0a 2a 2f 0a 0a 23 69 6e ) * rm ..*/..#in
07f0: 63 6c 75 64 65 20 3c 73 74 64 69 6f 2e 68 3e 0a clude <stdio.h>.
0800: 23 69 6e 63 6c 75 64 65 20 3c 73 74 64 6c 69 62 #include <stdlib
0810: 2e 68 3e 0a 23 69 6e 63 6c 75 64 65 20 3c 6d 61 .h>.#include <ma
0820: 74 68 2e 68 3e 0a 0a 23 69 6e 63 6c 75 64 65 20 th.h>..#include
0830: 22 67 6d 70 2e 68 22 0a 23 69 6e 63 6c 75 64 65 "gmp.h".#include
0840: 20 22 67 6d 70 73 74 61 74 2e 68 22 0a 0a 2f 2a "gmpstat.h"../*
0850: 20 6b 73 20 28 4b 70 2c 20 4b 6d 2c 20 58 2c 20 ks (Kp, Km, X,
0860: 50 2c 20 6e 29 20 2d 2d 20 50 65 72 66 6f 72 6d P, n) -- Perform
0870: 20 61 20 4b 6f 6c 6d 6f 67 6f 72 6f 76 2d 53 6d a Kolmogorov-Sm
0880: 69 72 6e 6f 76 20 74 65 73 74 20 6f 6e 20 74 68 irnov test on th
0890: 65 20 4e 0a 20 20 20 72 65 61 6c 20 6e 75 6d 62 e N. real numb
08a0: 65 72 73 20 62 65 74 77 65 65 6e 20 7a 65 72 6f ers between zero
08b0: 20 61 6e 64 20 6f 6e 65 20 69 6e 20 76 65 63 74 and one in vect
08c0: 6f 72 20 58 2e 20 20 50 20 69 73 20 74 68 65 0a or X. P is the.
08d0: 20 20 20 64 69 73 74 72 69 62 75 74 69 6f 6e 20 distribution
08e0: 66 75 6e 63 74 69 6f 6e 2c 20 63 61 6c 6c 65 64 function, called
08f0: 20 66 6f 72 20 65 61 63 68 20 65 6e 74 72 79 20 for each entry
0900: 69 6e 20 58 2c 20 77 68 69 63 68 20 73 68 6f 75 in X, which shou
0910: 6c 64 0a 20 20 20 63 61 6c 63 75 6c 61 74 65 20 ld. calculate
0920: 74 68 65 20 70 72 6f 62 61 62 69 6c 69 74 79 20 the probability
0930: 6f 66 20 58 20 62 65 69 6e 67 20 67 72 65 61 74 of X being great
0940: 65 72 20 74 68 61 6e 20 6f 72 20 65 71 75 61 6c er than or equal
0950: 20 74 6f 20 61 6e 79 0a 20 20 20 6e 75 6d 62 65 to any. numbe
0960: 72 20 69 6e 20 74 68 65 20 73 65 71 75 65 6e 63 r in the sequenc
0970: 65 2e 20 20 28 46 6f 72 20 61 20 75 6e 69 66 6f e. (For a unifo
0980: 72 6d 6c 79 20 64 69 73 74 72 69 62 75 74 65 64 rmly distributed
0990: 20 73 65 71 75 65 6e 63 65 20 6f 66 0a 20 20 20 sequence of.
09a0: 72 65 61 6c 20 6e 75 6d 62 65 72 73 20 62 65 74 real numbers bet
09b0: 77 65 65 6e 20 7a 65 72 6f 20 61 6e 64 20 6f 6e ween zero and on
09c0: 65 2c 20 74 68 69 73 20 69 73 20 73 69 6d 70 6c e, this is simpl
09d0: 79 20 65 71 75 61 6c 20 74 6f 20 58 2e 29 20 20 y equal to X.)
09e0: 54 68 65 0a 20 20 20 72 65 73 75 6c 74 20 69 73 The. result is
09f0: 20 70 75 74 20 69 6e 20 4b 70 20 61 6e 64 20 4b put in Kp and K
0a00: 6d 2e 20 20 2a 2f 0a 0a 76 6f 69 64 0a 6b 73 20 m. */..void.ks
0a10: 28 6d 70 66 5f 74 20 4b 70 2c 0a 20 20 20 20 6d (mpf_t Kp,. m
0a20: 70 66 5f 74 20 4b 6d 2c 0a 20 20 20 20 6d 70 66 pf_t Km,. mpf
0a30: 5f 74 20 58 5b 5d 2c 0a 20 20 20 20 76 6f 69 64 _t X[],. void
0a40: 20 28 50 29 20 28 6d 70 66 5f 74 2c 20 6d 70 66 (P) (mpf_t, mpf
0a50: 5f 74 29 2c 0a 20 20 20 20 75 6e 73 69 67 6e 65 _t),. unsigne
0a60: 64 20 6c 6f 6e 67 20 69 6e 74 20 6e 29 0a 7b 0a d long int n).{.
0a70: 20 20 6d 70 66 5f 74 20 4b 74 3b 09 09 09 2f 2a mpf_t Kt;.../*
0a80: 20 74 65 6d 70 20 2a 2f 0a 20 20 6d 70 66 5f 74 temp */. mpf_t
0a90: 20 66 5f 78 3b 0a 20 20 6d 70 66 5f 74 20 66 5f f_x;. mpf_t f_
0aa0: 6a 3b 09 09 09 2f 2a 20 6a 20 2a 2f 0a 20 20 6d j;.../* j */. m
0ab0: 70 66 5f 74 20 66 5f 6a 6e 71 3b 09 09 09 2f 2a pf_t f_jnq;.../*
0ac0: 20 6a 2f 6e 20 6f 72 20 28 6a 2d 31 29 2f 6e 20 j/n or (j-1)/n
0ad0: 2a 2f 0a 20 20 75 6e 73 69 67 6e 65 64 20 6c 6f */. unsigned lo
0ae0: 6e 67 20 69 6e 74 20 6a 3b 0a 0a 20 20 2f 2a 20 ng int j;.. /*
0af0: 53 6f 72 74 20 74 68 65 20 76 65 63 74 6f 72 20 Sort the vector
0b00: 69 6e 20 61 73 63 65 6e 64 69 6e 67 20 6f 72 64 in ascending ord
0b10: 65 72 2e 20 2a 2f 20 20 0a 20 20 71 73 6f 72 74 er. */ . qsort
0b20: 20 28 58 2c 20 6e 2c 20 73 69 7a 65 6f 66 20 28 (X, n, sizeof (
0b30: 5f 5f 6d 70 66 5f 73 74 72 75 63 74 29 2c 20 6d __mpf_struct), m
0b40: 70 66 5f 63 6d 70 29 3b 0a 0a 20 20 2f 2a 20 4b pf_cmp);.. /* K
0b50: 2d 53 20 74 65 73 74 2e 20 2a 2f 0a 20 20 2f 2a -S test. */. /*
0b60: 20 09 4b 70 20 3d 20 73 71 72 28 6e 29 20 2a 20 .Kp = sqr(n) *
0b70: 6d 61 78 28 6a 2f 6e 20 2d 20 46 28 58 6a 29 29 max(j/n - F(Xj))
0b80: 09 09 66 6f 72 20 61 6c 6c 20 31 3c 3d 6a 3c 3d ..for all 1<=j<=
0b90: 6e 0a 09 4b 6d 20 3d 20 73 71 72 28 6e 29 20 2a n..Km = sqr(n) *
0ba0: 20 6d 61 78 28 46 28 58 6a 29 20 2d 20 28 6a 2d max(F(Xj) - (j-
0bb0: 31 29 2f 6e 29 29 09 66 6f 72 20 61 6c 6c 20 31 1)/n)).for all 1
0bc0: 3c 3d 6a 3c 3d 6e 0a 20 20 2a 2f 0a 0a 20 20 6d <=j<=n. */.. m
0bd0: 70 66 5f 69 6e 69 74 20 28 4b 74 29 3b 20 6d 70 pf_init (Kt); mp
0be0: 66 5f 69 6e 69 74 20 28 66 5f 78 29 3b 20 6d 70 f_init (f_x); mp
0bf0: 66 5f 69 6e 69 74 20 28 66 5f 6a 29 3b 20 6d 70 f_init (f_j); mp
0c00: 66 5f 69 6e 69 74 20 28 66 5f 6a 6e 71 29 3b 20 f_init (f_jnq);
0c10: 0a 20 20 6d 70 66 5f 73 65 74 5f 75 69 20 28 4b . mpf_set_ui (K
0c20: 70 2c 20 30 29 3b 20 20 6d 70 66 5f 73 65 74 5f p, 0); mpf_set_
0c30: 75 69 20 28 4b 6d 2c 20 30 29 3b 0a 20 20 66 6f ui (Km, 0);. fo
0c40: 72 20 28 6a 20 3d 20 31 3b 20 6a 20 3c 3d 20 6e r (j = 1; j <= n
0c50: 3b 20 6a 2b 2b 29 0a 20 20 20 20 7b 0a 20 20 20 ; j++). {.
0c60: 20 20 20 50 20 28 66 5f 78 2c 20 58 5b 6a 2d 31 P (f_x, X[j-1
0c70: 5d 29 3b 0a 20 20 20 20 20 20 6d 70 66 5f 73 65 ]);. mpf_se
0c80: 74 5f 75 69 20 28 66 5f 6a 2c 20 6a 29 3b 0a 0a t_ui (f_j, j);..
0c90: 20 20 20 20 20 20 6d 70 66 5f 64 69 76 5f 75 69 mpf_div_ui
0ca0: 20 28 66 5f 6a 6e 71 2c 20 66 5f 6a 2c 20 6e 29 (f_jnq, f_j, n)
0cb0: 3b 0a 20 20 20 20 20 20 6d 70 66 5f 73 75 62 20 ;. mpf_sub
0cc0: 28 4b 74 2c 20 66 5f 6a 6e 71 2c 20 66 5f 78 29 (Kt, f_jnq, f_x)
0cd0: 3b 0a 20 20 20 20 20 20 69 66 20 28 6d 70 66 5f ;. if (mpf_
0ce0: 63 6d 70 20 28 4b 74 2c 20 4b 70 29 20 3e 20 30 cmp (Kt, Kp) > 0
0cf0: 29 0a 09 6d 70 66 5f 73 65 74 20 28 4b 70 2c 20 )..mpf_set (Kp,
0d00: 4b 74 29 3b 0a 20 20 20 20 20 20 69 66 20 28 67 Kt);. if (g
0d10: 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 5f 32 _debug > DEBUG_2
0d20: 29 0a 09 7b 0a 09 20 20 70 72 69 6e 74 66 20 28 )..{.. printf (
0d30: 22 6a 3d 25 6c 75 20 22 2c 20 6a 29 3b 0a 09 20 "j=%lu ", j);..
0d40: 20 70 72 69 6e 74 66 20 28 22 50 28 29 3d 22 29 printf ("P()=")
0d50: 3b 20 6d 70 66 5f 6f 75 74 5f 73 74 72 20 28 73 ; mpf_out_str (s
0d60: 74 64 6f 75 74 2c 20 31 30 2c 20 32 2c 20 66 5f tdout, 10, 2, f_
0d70: 78 29 3b 20 70 72 69 6e 74 66 20 28 22 5c 74 22 x); printf ("\t"
0d80: 29 3b 0a 0a 09 20 20 70 72 69 6e 74 66 20 28 22 );... printf ("
0d90: 6a 6e 71 3d 22 29 3b 20 6d 70 66 5f 6f 75 74 5f jnq="); mpf_out_
0da0: 73 74 72 20 28 73 74 64 6f 75 74 2c 20 31 30 2c str (stdout, 10,
0db0: 20 32 2c 20 66 5f 6a 6e 71 29 3b 20 70 72 69 6e 2, f_jnq); prin
0dc0: 74 66 20 28 22 20 22 29 3b 0a 09 20 20 70 72 69 tf (" ");.. pri
0dd0: 6e 74 66 20 28 22 64 69 66 66 3d 22 29 3b 20 6d ntf ("diff="); m
0de0: 70 66 5f 6f 75 74 5f 73 74 72 20 28 73 74 64 6f pf_out_str (stdo
0df0: 75 74 2c 20 31 30 2c 20 32 2c 20 4b 74 29 3b 20 ut, 10, 2, Kt);
0e00: 70 72 69 6e 74 66 20 28 22 20 22 29 3b 0a 09 20 printf (" ");..
0e10: 20 70 72 69 6e 74 66 20 28 22 4b 70 3d 22 29 3b printf ("Kp=");
0e20: 20 6d 70 66 5f 6f 75 74 5f 73 74 72 20 28 73 74 mpf_out_str (st
0e30: 64 6f 75 74 2c 20 31 30 2c 20 32 2c 20 4b 70 29 dout, 10, 2, Kp)
0e40: 3b 20 70 72 69 6e 74 66 20 28 22 5c 74 22 29 3b ; printf ("\t");
0e50: 0a 09 7d 0a 20 20 20 20 20 20 6d 70 66 5f 73 75 ..}. mpf_su
0e60: 62 5f 75 69 20 28 66 5f 6a 2c 20 66 5f 6a 2c 20 b_ui (f_j, f_j,
0e70: 31 29 3b 0a 20 20 20 20 20 20 6d 70 66 5f 64 69 1);. mpf_di
0e80: 76 5f 75 69 20 28 66 5f 6a 6e 71 2c 20 66 5f 6a v_ui (f_jnq, f_j
0e90: 2c 20 6e 29 3b 0a 20 20 20 20 20 20 6d 70 66 5f , n);. mpf_
0ea0: 73 75 62 20 28 4b 74 2c 20 66 5f 78 2c 20 66 5f sub (Kt, f_x, f_
0eb0: 6a 6e 71 29 3b 0a 20 20 20 20 20 20 69 66 20 28 jnq);. if (
0ec0: 6d 70 66 5f 63 6d 70 20 28 4b 74 2c 20 4b 6d 29 mpf_cmp (Kt, Km)
0ed0: 20 3e 20 30 29 0a 09 6d 70 66 5f 73 65 74 20 28 > 0)..mpf_set (
0ee0: 4b 6d 2c 20 4b 74 29 3b 0a 0a 20 20 20 20 20 20 Km, Kt);..
0ef0: 69 66 20 28 67 5f 64 65 62 75 67 20 3e 20 44 45 if (g_debug > DE
0f00: 42 55 47 5f 32 29 0a 09 7b 0a 09 20 20 70 72 69 BUG_2)..{.. pri
0f10: 6e 74 66 20 28 22 6a 6e 71 3d 22 29 3b 20 6d 70 ntf ("jnq="); mp
0f20: 66 5f 6f 75 74 5f 73 74 72 20 28 73 74 64 6f 75 f_out_str (stdou
0f30: 74 2c 20 31 30 2c 20 32 2c 20 66 5f 6a 6e 71 29 t, 10, 2, f_jnq)
0f40: 3b 20 70 72 69 6e 74 66 20 28 22 20 22 29 3b 0a ; printf (" ");.
0f50: 09 20 20 70 72 69 6e 74 66 20 28 22 64 69 66 66 . printf ("diff
0f60: 3d 22 29 3b 20 6d 70 66 5f 6f 75 74 5f 73 74 72 ="); mpf_out_str
0f70: 20 28 73 74 64 6f 75 74 2c 20 31 30 2c 20 32 2c (stdout, 10, 2,
0f80: 20 4b 74 29 3b 20 70 72 69 6e 74 66 20 28 22 20 Kt); printf ("
0f90: 22 29 3b 0a 09 20 20 70 72 69 6e 74 66 20 28 22 ");.. printf ("
0fa0: 4b 6d 3d 22 29 3b 20 6d 70 66 5f 6f 75 74 5f 73 Km="); mpf_out_s
0fb0: 74 72 20 28 73 74 64 6f 75 74 2c 20 31 30 2c 20 tr (stdout, 10,
0fc0: 32 2c 20 4b 6d 29 3b 20 70 72 69 6e 74 66 20 28 2, Km); printf (
0fd0: 22 20 22 29 3b 0a 09 20 20 70 72 69 6e 74 66 20 " ");.. printf
0fe0: 28 22 5c 6e 22 29 3b 0a 09 7d 0a 20 20 20 20 7d ("\n");..}. }
0ff0: 0a 20 20 6d 70 66 5f 73 71 72 74 5f 75 69 20 28 . mpf_sqrt_ui (
1000: 4b 74 2c 20 6e 29 3b 0a 20 20 6d 70 66 5f 6d 75 Kt, n);. mpf_mu
1010: 6c 20 28 4b 70 2c 20 4b 70 2c 20 4b 74 29 3b 0a l (Kp, Kp, Kt);.
1020: 20 20 6d 70 66 5f 6d 75 6c 20 28 4b 6d 2c 20 4b mpf_mul (Km, K
1030: 6d 2c 20 4b 74 29 3b 0a 0a 20 20 6d 70 66 5f 63 m, Kt);.. mpf_c
1040: 6c 65 61 72 20 28 4b 74 29 3b 20 6d 70 66 5f 63 lear (Kt); mpf_c
1050: 6c 65 61 72 20 28 66 5f 78 29 3b 20 6d 70 66 5f lear (f_x); mpf_
1060: 63 6c 65 61 72 20 28 66 5f 6a 29 3b 20 6d 70 66 clear (f_j); mpf
1070: 5f 63 6c 65 61 72 20 28 66 5f 6a 6e 71 29 3b 20 _clear (f_jnq);
1080: 0a 7d 0a 0a 2f 2a 20 6b 73 5f 74 61 62 6c 65 28 .}../* ks_table(
1090: 76 61 6c 2c 20 6e 29 20 2d 2d 20 63 61 6c 63 75 val, n) -- calcu
10a0: 6c 61 74 65 20 70 72 6f 62 61 62 69 6c 69 74 79 late probability
10b0: 20 66 6f 72 20 4b 70 2f 4b 6d 20 6c 65 73 73 20 for Kp/Km less
10c0: 74 68 61 6e 20 6f 72 0a 20 20 20 65 71 75 61 6c than or. equal
10d0: 20 74 6f 20 56 41 4c 20 77 69 74 68 20 4e 20 6f to VAL with N o
10e0: 62 73 65 72 76 61 74 69 6f 6e 73 2e 20 20 53 65 bservations. Se
10f0: 65 20 5b 4b 6e 75 74 68 20 73 65 63 74 69 6f 6e e [Knuth section
1100: 20 33 2e 33 2e 31 5d 20 2a 2f 0a 0a 76 6f 69 64 3.3.1] */..void
1110: 0a 6b 73 5f 74 61 62 6c 65 20 28 6d 70 66 5f 74 .ks_table (mpf_t
1120: 20 70 2c 20 6d 70 66 5f 74 20 76 61 6c 2c 20 63 p, mpf_t val, c
1130: 6f 6e 73 74 20 75 6e 73 69 67 6e 65 64 20 69 6e onst unsigned in
1140: 74 20 6e 29 0a 7b 0a 20 20 2f 2a 20 57 65 20 75 t n).{. /* We u
1150: 73 65 20 45 71 2e 20 28 32 37 29 2c 20 4b 6e 75 se Eq. (27), Knu
1160: 74 68 20 70 2e 35 38 2c 20 73 6b 69 70 70 69 6e th p.58, skippin
1170: 67 20 4f 28 31 2f 6e 29 20 66 6f 72 20 73 69 6d g O(1/n) for sim
1180: 70 6c 69 63 69 74 79 2e 0a 20 20 20 20 20 54 68 plicity.. Th
1190: 69 73 20 73 68 6f 72 74 63 75 74 20 77 69 6c 6c is shortcut will
11a0: 20 72 65 73 75 6c 74 20 69 6e 20 74 6f 6f 20 68 result in too h
11b0: 69 67 68 20 70 72 6f 62 61 62 69 6c 69 74 69 65 igh probabilitie
11c0: 73 2c 20 65 73 70 65 63 69 61 6c 6c 79 0a 20 20 s, especially.
11d0: 20 20 20 77 68 65 6e 20 6e 20 69 73 20 73 6d 61 when n is sma
11e0: 6c 6c 2e 0a 0a 20 20 20 20 20 50 72 28 4b 70 28 ll... Pr(Kp(
11f0: 6e 29 20 3c 3d 20 73 29 20 3d 20 31 20 2d 20 70 n) <= s) = 1 - p
1200: 6f 77 28 65 2c 20 2d 32 2a 73 5e 32 29 20 2a 20 ow(e, -2*s^2) *
1210: 28 31 20 2d 20 32 2f 33 2a 73 2f 73 71 72 74 28 (1 - 2/3*s/sqrt(
1220: 6e 29 20 2b 20 4f 28 31 2f 6e 29 29 20 2a 2f 0a n) + O(1/n)) */.
1230: 0a 20 20 2f 2a 20 57 65 20 68 61 76 65 20 27 73 . /* We have 's
1240: 27 20 69 6e 20 76 61 72 69 61 62 6c 65 20 56 41 ' in variable VA
1250: 4c 20 61 6e 64 20 73 74 6f 72 65 20 74 68 65 20 L and store the
1260: 72 65 73 75 6c 74 20 69 6e 20 50 2e 20 2a 2f 0a result in P. */.
1270: 0a 20 20 6d 70 66 5f 74 20 74 31 2c 20 74 32 3b . mpf_t t1, t2;
1280: 0a 0a 20 20 6d 70 66 5f 69 6e 69 74 20 28 74 31 .. mpf_init (t1
1290: 29 3b 20 6d 70 66 5f 69 6e 69 74 20 28 74 32 29 ); mpf_init (t2)
12a0: 3b 0a 0a 20 20 2f 2a 20 74 31 20 3d 20 31 20 2d ;.. /* t1 = 1 -
12b0: 20 32 2f 33 20 2a 20 73 2f 73 71 72 74 28 6e 29 2/3 * s/sqrt(n)
12c0: 20 2a 2f 0a 20 20 6d 70 66 5f 73 71 72 74 5f 75 */. mpf_sqrt_u
12d0: 69 20 28 74 31 2c 20 6e 29 3b 0a 20 20 6d 70 66 i (t1, n);. mpf
12e0: 5f 64 69 76 20 28 74 31 2c 20 76 61 6c 2c 20 74 _div (t1, val, t
12f0: 31 29 3b 0a 20 20 6d 70 66 5f 6d 75 6c 5f 75 69 1);. mpf_mul_ui
1300: 20 28 74 31 2c 20 74 31 2c 20 32 29 3b 0a 20 20 (t1, t1, 2);.
1310: 6d 70 66 5f 64 69 76 5f 75 69 20 28 74 31 2c 20 mpf_div_ui (t1,
1320: 74 31 2c 20 33 29 3b 0a 20 20 6d 70 66 5f 75 69 t1, 3);. mpf_ui
1330: 5f 73 75 62 20 28 74 31 2c 20 31 2c 20 74 31 29 _sub (t1, 1, t1)
1340: 3b 0a 0a 20 20 2f 2a 20 74 32 20 3d 20 70 6f 77 ;.. /* t2 = pow
1350: 28 65 2c 20 2d 32 2a 73 5e 32 29 20 2a 2f 0a 23 (e, -2*s^2) */.#
1360: 69 66 6e 64 65 66 20 4f 4c 44 47 4d 50 0a 20 20 ifndef OLDGMP.
1370: 6d 70 66 5f 70 6f 77 5f 75 69 20 28 74 32 2c 20 mpf_pow_ui (t2,
1380: 76 61 6c 2c 20 32 29 3b 09 2f 2a 20 74 32 20 3d val, 2);./* t2 =
1390: 20 73 5e 32 20 2a 2f 0a 20 20 6d 70 66 5f 73 65 s^2 */. mpf_se
13a0: 74 5f 64 20 28 74 32 2c 20 65 78 70 20 28 2d 28 t_d (t2, exp (-(
13b0: 32 2e 30 20 2a 20 6d 70 66 5f 67 65 74 5f 64 20 2.0 * mpf_get_d
13c0: 28 74 32 29 29 29 29 3b 0a 23 65 6c 73 65 0a 20 (t2))));.#else.
13d0: 20 2f 2a 20 68 6d 6d 6d 2c 20 67 6d 70 20 64 6f /* hmmm, gmp do
13e0: 65 73 6e 27 74 20 68 61 76 65 20 70 6f 77 28 29 esn't have pow()
13f0: 20 66 6f 72 20 66 6c 6f 61 74 73 2e 20 20 75 73 for floats. us
1400: 65 20 64 6f 75 62 6c 65 73 2e 20 2a 2f 0a 20 20 e doubles. */.
1410: 6d 70 66 5f 73 65 74 5f 64 20 28 74 32 2c 20 70 mpf_set_d (t2, p
1420: 6f 77 20 28 4d 5f 45 2c 20 2d 28 32 20 2a 20 70 ow (M_E, -(2 * p
1430: 6f 77 20 28 6d 70 66 5f 67 65 74 5f 64 20 28 76 ow (mpf_get_d (v
1440: 61 6c 29 2c 20 32 29 29 29 29 3b 0a 23 65 6e 64 al), 2))));.#end
1450: 69 66 20 20 0a 0a 20 20 2f 2a 20 70 20 3d 20 31 if .. /* p = 1
1460: 20 2d 20 74 31 20 2a 20 74 32 20 2a 2f 0a 20 20 - t1 * t2 */.
1470: 6d 70 66 5f 6d 75 6c 20 28 74 31 2c 20 74 31 2c mpf_mul (t1, t1,
1480: 20 74 32 29 3b 0a 20 20 6d 70 66 5f 75 69 5f 73 t2);. mpf_ui_s
1490: 75 62 20 28 70 2c 20 31 2c 20 74 31 29 3b 0a 0a ub (p, 1, t1);..
14a0: 20 20 6d 70 66 5f 63 6c 65 61 72 20 28 74 31 29 mpf_clear (t1)
14b0: 3b 20 6d 70 66 5f 63 6c 65 61 72 20 28 74 32 29 ; mpf_clear (t2)
14c0: 3b 0a 7d 0a 0a 73 74 61 74 69 63 20 64 6f 75 62 ;.}..static doub
14d0: 6c 65 20 78 32 5f 74 61 62 6c 65 5f 58 5b 5d 5b le x2_table_X[][
14e0: 37 5d 20 3d 20 7b 0a 20 20 7b 20 2d 32 2e 33 33 7] = {. { -2.33
14f0: 2c 20 2d 31 2e 36 34 2c 20 2d 2e 36 37 34 2c 20 , -1.64, -.674,
1500: 30 2e 30 2c 20 30 2e 36 37 34 2c 20 31 2e 36 34 0.0, 0.674, 1.64
1510: 2c 20 32 2e 33 33 20 7d 2c 20 2f 2a 20 78 20 2a , 2.33 }, /* x *
1520: 2f 0a 20 20 7b 20 35 2e 34 32 38 39 2c 20 32 2e /. { 5.4289, 2.
1530: 36 38 39 36 2c 20 2e 34 35 34 32 37 36 2c 20 30 6896, .454276, 0
1540: 2e 30 2c 20 2e 34 35 34 32 37 36 2c 20 32 2e 36 .0, .454276, 2.6
1550: 38 39 36 2c 20 35 2e 34 32 38 39 7d 20 2f 2a 20 896, 5.4289} /*
1560: 78 5e 32 20 2a 2f 0a 7d 3b 0a 0a 23 64 65 66 69 x^2 */.};..#defi
1570: 6e 65 20 5f 32 44 33 20 28 28 64 6f 75 62 6c 65 ne _2D3 ((double
1580: 29 20 2e 36 36 36 36 36 36 36 36 36 36 29 0a 0a ) .6666666666)..
1590: 2f 2a 20 78 32 5f 74 61 62 6c 65 20 28 74 2c 20 /* x2_table (t,
15a0: 76 2c 20 6e 29 20 2d 2d 20 72 65 74 75 72 6e 20 v, n) -- return
15b0: 63 68 69 2d 73 71 75 61 72 65 20 74 61 62 6c 65 chi-square table
15c0: 20 72 6f 77 20 66 6f 72 20 56 20 69 6e 20 54 5b row for V in T[
15d0: 5d 2e 20 2a 2f 0a 76 6f 69 64 0a 78 32 5f 74 61 ]. */.void.x2_ta
15e0: 62 6c 65 20 28 64 6f 75 62 6c 65 20 74 5b 5d 2c ble (double t[],
15f0: 0a 09 20 20 75 6e 73 69 67 6e 65 64 20 69 6e 74 .. unsigned int
1600: 20 76 29 0a 7b 0a 20 20 69 6e 74 20 66 3b 0a 0a v).{. int f;..
1610: 0a 20 20 2f 2a 20 46 49 58 4d 45 3a 20 44 6f 20 . /* FIXME: Do
1620: 61 20 74 61 62 6c 65 20 6c 6f 6f 6b 75 70 20 66 a table lookup f
1630: 6f 72 20 76 20 3c 3d 20 33 30 20 73 69 6e 63 65 or v <= 30 since
1640: 20 74 68 65 20 66 6f 6c 6c 6f 77 69 6e 67 20 66 the following f
1650: 6f 72 6d 75 6c 61 0a 20 20 20 20 20 5b 4b 6e 75 ormula. [Knu
1660: 74 68 2c 20 76 6f 6c 20 32 2c 20 33 2e 33 2e 31 th, vol 2, 3.3.1
1670: 5d 20 69 73 20 6f 6e 6c 79 20 67 6f 6f 64 20 66 ] is only good f
1680: 6f 72 20 76 20 3e 20 33 30 2e 20 2a 2f 0a 0a 20 or v > 30. */..
1690: 20 2f 2a 20 76 61 6c 75 65 20 3d 20 76 20 2b 20 /* value = v +
16a0: 73 71 72 74 28 32 2a 76 29 20 2a 20 58 5b 70 5d sqrt(2*v) * X[p]
16b0: 20 2b 20 28 32 2f 33 29 20 2a 20 58 5b 70 5d 5e + (2/3) * X[p]^
16c0: 32 20 2d 20 32 2f 33 20 2b 20 4f 28 31 2f 73 71 2 - 2/3 + O(1/sq
16d0: 72 74 28 74 29 20 2a 2f 0a 20 20 2f 2a 20 4e 4f rt(t) */. /* NO
16e0: 54 45 3a 20 54 68 65 20 4f 28 29 20 74 65 72 6d TE: The O() term
16f0: 20 69 73 20 69 67 6e 6f 72 65 64 20 66 6f 72 20 is ignored for
1700: 73 69 6d 70 6c 69 63 69 74 79 2e 20 2a 2f 0a 20 simplicity. */.
1710: 20 0a 20 20 66 6f 72 20 28 66 20 3d 20 30 3b 20 . for (f = 0;
1720: 66 20 3c 20 37 3b 20 66 2b 2b 29 0a 20 20 20 20 f < 7; f++).
1730: 20 20 74 5b 66 5d 20 3d 0a 09 76 20 2b 0a 09 73 t[f] =..v +..s
1740: 71 72 74 20 28 32 20 2a 20 76 29 20 2a 20 78 32 qrt (2 * v) * x2
1750: 5f 74 61 62 6c 65 5f 58 5b 30 5d 5b 66 5d 20 2b _table_X[0][f] +
1760: 0a 09 5f 32 44 33 20 2a 20 78 32 5f 74 61 62 6c .._2D3 * x2_tabl
1770: 65 5f 58 5b 31 5d 5b 66 5d 20 2d 20 5f 32 44 33 e_X[1][f] - _2D3
1780: 3b 0a 7d 0a 0a 0a 2f 2a 20 50 28 70 2c 20 78 29 ;.}.../* P(p, x)
1790: 20 2d 2d 20 44 69 73 74 72 69 62 75 74 69 6f 6e -- Distribution
17a0: 20 66 75 6e 63 74 69 6f 6e 2e 20 20 43 61 6c 63 function. Calc
17b0: 75 6c 61 74 65 20 74 68 65 20 70 72 6f 62 61 62 ulate the probab
17c0: 69 6c 69 74 79 20 6f 66 20 58 0a 62 65 69 6e 67 ility of X.being
17d0: 20 67 72 65 61 74 65 72 20 74 68 61 6e 20 6f 72 greater than or
17e0: 20 65 71 75 61 6c 20 74 6f 20 61 6e 79 20 6e 75 equal to any nu
17f0: 6d 62 65 72 20 69 6e 20 74 68 65 20 73 65 71 75 mber in the sequ
1800: 65 6e 63 65 2e 20 20 46 6f 72 20 61 0a 72 61 6e ence. For a.ran
1810: 64 6f 6d 20 72 65 61 6c 20 6e 75 6d 62 65 72 20 dom real number
1820: 62 65 74 77 65 65 6e 20 7a 65 72 6f 20 61 6e 64 between zero and
1830: 20 6f 6e 65 20 67 69 76 65 6e 20 62 79 20 61 20 one given by a
1840: 75 6e 69 66 6f 72 6d 6c 79 0a 64 69 73 74 72 69 uniformly.distri
1850: 62 75 74 65 64 20 72 61 6e 64 6f 6d 20 6e 75 6d buted random num
1860: 62 65 72 20 67 65 6e 65 72 61 74 6f 72 2c 20 74 ber generator, t
1870: 68 69 73 20 69 73 20 73 69 6d 70 6c 79 20 65 71 his is simply eq
1880: 75 61 6c 20 74 6f 20 58 2e 20 2a 2f 0a 0a 73 74 ual to X. */..st
1890: 61 74 69 63 20 76 6f 69 64 20 0a 50 20 28 6d 70 atic void .P (mp
18a0: 66 5f 74 20 70 2c 20 6d 70 66 5f 74 20 78 29 0a f_t p, mpf_t x).
18b0: 7b 0a 20 20 6d 70 66 5f 73 65 74 20 28 70 2c 20 {. mpf_set (p,
18c0: 78 29 3b 0a 7d 0a 0a 2f 2a 20 6d 70 66 5f 66 72 x);.}../* mpf_fr
18d0: 65 71 74 28 29 20 2d 2d 20 46 72 65 71 75 65 6e eqt() -- Frequen
18e0: 63 79 20 74 65 73 74 20 75 73 69 6e 67 20 4b 53 cy test using KS
18f0: 20 6f 6e 20 4e 20 72 65 61 6c 20 6e 75 6d 62 65 on N real numbe
1900: 72 73 20 62 65 74 77 65 65 6e 20 7a 65 72 6f 0a rs between zero.
1910: 20 20 20 61 6e 64 20 6f 6e 65 2e 20 20 53 65 65 and one. See
1920: 20 5b 4b 6e 75 74 68 20 76 6f 6c 20 32 2c 20 70 [Knuth vol 2, p
1930: 2e 36 31 5d 2e 20 2a 2f 0a 76 6f 69 64 0a 6d 70 .61]. */.void.mp
1940: 66 5f 66 72 65 71 74 20 28 6d 70 66 5f 74 20 4b f_freqt (mpf_t K
1950: 70 2c 0a 09 20 20 20 6d 70 66 5f 74 20 4b 6d 2c p,.. mpf_t Km,
1960: 0a 09 20 20 20 6d 70 66 5f 74 20 58 5b 5d 2c 0a .. mpf_t X[],.
1970: 09 20 20 20 63 6f 6e 73 74 20 75 6e 73 69 67 6e . const unsign
1980: 65 64 20 6c 6f 6e 67 20 69 6e 74 20 6e 29 0a 7b ed long int n).{
1990: 0a 20 20 6b 73 20 28 4b 70 2c 20 4b 6d 2c 20 58 . ks (Kp, Km, X
19a0: 2c 20 50 2c 20 6e 29 3b 0a 7d 0a 0a 0a 2f 2a 20 , P, n);.}.../*
19b0: 54 68 65 20 43 68 69 2d 73 71 75 61 72 65 20 74 The Chi-square t
19c0: 65 73 74 2e 20 20 45 71 2e 20 28 38 29 20 69 6e est. Eq. (8) in
19d0: 20 4b 6e 75 74 68 20 76 6f 6c 2e 20 32 20 73 61 Knuth vol. 2 sa
19e0: 79 73 20 74 68 61 74 20 69 66 20 59 5b 5d 0a 20 ys that if Y[].
19f0: 20 20 68 6f 6c 64 73 20 74 68 65 20 6f 62 73 65 holds the obse
1a00: 72 76 61 74 69 6f 6e 73 20 61 6e 64 20 70 5b 5d rvations and p[]
1a10: 20 69 73 20 74 68 65 20 70 72 6f 62 61 62 69 6c is the probabil
1a20: 69 74 79 20 66 6f 72 2e 2e 20 28 74 6f 20 62 65 ity for.. (to be
1a30: 0a 20 20 20 63 6f 6e 74 69 6e 75 65 64 21 29 0a . continued!).
1a40: 0a 20 20 20 56 20 3d 20 31 2f 6e 20 2a 20 73 75 . V = 1/n * su
1a50: 6d 28 28 73 3d 31 20 74 6f 20 6b 29 20 59 5b 73 m((s=1 to k) Y[s
1a60: 5d 5e 32 20 2f 20 70 5b 73 5d 29 20 2d 20 6e 20 ]^2 / p[s]) - n
1a70: 2a 2f 0a 0a 76 6f 69 64 0a 78 32 20 28 6d 70 66 */..void.x2 (mpf
1a80: 5f 74 20 56 2c 09 09 09 2f 2a 20 72 65 73 75 6c _t V,.../* resul
1a90: 74 20 2a 2f 0a 20 20 20 20 75 6e 73 69 67 6e 65 t */. unsigne
1aa0: 64 20 6c 6f 6e 67 20 69 6e 74 20 58 5b 5d 2c 09 d long int X[],.
1ab0: 2f 2a 20 64 61 74 61 20 2a 2f 0a 20 20 20 20 75 /* data */. u
1ac0: 6e 73 69 67 6e 65 64 20 69 6e 74 20 6b 2c 09 09 nsigned int k,..
1ad0: 2f 2a 20 23 6f 66 20 63 61 74 65 67 6f 72 69 65 /* #of categorie
1ae0: 73 20 2a 2f 0a 20 20 20 20 76 6f 69 64 20 28 50 s */. void (P
1af0: 29 20 28 6d 70 66 5f 74 2c 20 75 6e 73 69 67 6e ) (mpf_t, unsign
1b00: 65 64 20 6c 6f 6e 67 20 69 6e 74 2c 20 76 6f 69 ed long int, voi
1b10: 64 20 2a 29 2c 20 2f 2a 20 70 72 6f 62 61 62 69 d *), /* probabi
1b20: 6c 69 74 79 20 66 75 6e 63 20 2a 2f 0a 20 20 20 lity func */.
1b30: 20 76 6f 69 64 20 2a 78 2c 09 09 09 2f 2a 20 65 void *x,.../* e
1b40: 78 74 72 61 20 75 73 65 72 20 64 61 74 61 20 70 xtra user data p
1b50: 61 73 73 65 64 20 74 6f 20 50 28 29 20 2a 2f 0a assed to P() */.
1b60: 20 20 20 20 75 6e 73 69 67 6e 65 64 20 6c 6f 6e unsigned lon
1b70: 67 20 69 6e 74 20 6e 29 09 2f 2a 20 23 6f 66 20 g int n)./* #of
1b80: 73 61 6d 70 6c 65 73 20 2a 2f 0a 7b 0a 20 20 75 samples */.{. u
1b90: 6e 73 69 67 6e 65 64 20 69 6e 74 20 66 3b 0a 20 nsigned int f;.
1ba0: 20 6d 70 66 5f 74 20 66 5f 74 2c 20 66 5f 74 32 mpf_t f_t, f_t2
1bb0: 3b 09 09 2f 2a 20 74 65 6d 70 20 66 6c 6f 61 74 ;../* temp float
1bc0: 73 20 2a 2f 0a 0a 20 20 6d 70 66 5f 69 6e 69 74 s */.. mpf_init
1bd0: 20 28 66 5f 74 29 3b 20 6d 70 66 5f 69 6e 69 74 (f_t); mpf_init
1be0: 20 28 66 5f 74 32 29 3b 0a 0a 0a 20 20 6d 70 66 (f_t2);... mpf
1bf0: 5f 73 65 74 5f 75 69 20 28 56 2c 20 30 29 3b 0a _set_ui (V, 0);.
1c00: 20 20 66 6f 72 20 28 66 20 3d 20 30 3b 20 66 20 for (f = 0; f
1c10: 3c 20 6b 3b 20 66 2b 2b 29 0a 20 20 20 20 7b 0a < k; f++). {.
1c20: 20 20 20 20 20 20 69 66 20 28 67 5f 64 65 62 75 if (g_debu
1c30: 67 20 3e 20 44 45 42 55 47 5f 32 29 0a 09 66 70 g > DEBUG_2)..fp
1c40: 72 69 6e 74 66 20 28 73 74 64 65 72 72 2c 20 22 rintf (stderr, "
1c50: 25 75 3a 20 50 28 29 3d 22 2c 20 66 29 3b 0a 20 %u: P()=", f);.
1c60: 20 20 20 20 20 6d 70 66 5f 73 65 74 5f 75 69 20 mpf_set_ui
1c70: 28 66 5f 74 2c 20 58 5b 66 5d 29 3b 0a 20 20 20 (f_t, X[f]);.
1c80: 20 20 20 6d 70 66 5f 6d 75 6c 20 28 66 5f 74 2c mpf_mul (f_t,
1c90: 20 66 5f 74 2c 20 66 5f 74 29 3b 09 2f 2a 20 66 f_t, f_t);./* f
1ca0: 5f 74 20 3d 20 58 5b 66 5d 5e 32 20 2a 2f 0a 20 _t = X[f]^2 */.
1cb0: 20 20 20 20 20 50 20 28 66 5f 74 32 2c 20 66 2c P (f_t2, f,
1cc0: 20 78 29 3b 09 09 2f 2a 20 66 5f 74 32 20 3d 20 x);../* f_t2 =
1cd0: 50 72 28 66 29 20 2a 2f 0a 20 20 20 20 20 20 69 Pr(f) */. i
1ce0: 66 20 28 67 5f 64 65 62 75 67 20 3e 20 44 45 42 f (g_debug > DEB
1cf0: 55 47 5f 32 29 0a 09 6d 70 66 5f 6f 75 74 5f 73 UG_2)..mpf_out_s
1d00: 74 72 20 28 73 74 64 65 72 72 2c 20 31 30 2c 20 tr (stderr, 10,
1d10: 32 2c 20 66 5f 74 32 29 3b 0a 20 20 20 20 20 20 2, f_t2);.
1d20: 6d 70 66 5f 64 69 76 20 28 66 5f 74 2c 20 66 5f mpf_div (f_t, f_
1d30: 74 2c 20 66 5f 74 32 29 3b 0a 20 20 20 20 20 20 t, f_t2);.
1d40: 6d 70 66 5f 61 64 64 20 28 56 2c 20 56 2c 20 66 mpf_add (V, V, f
1d50: 5f 74 29 3b 0a 20 20 20 20 20 20 69 66 20 28 67 _t);. if (g
1d60: 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 5f 32 _debug > DEBUG_2
1d70: 29 0a 09 7b 0a 09 20 20 66 70 72 69 6e 74 66 20 )..{.. fprintf
1d80: 28 73 74 64 65 72 72 2c 20 22 5c 74 56 3d 22 29 (stderr, "\tV=")
1d90: 3b 0a 09 20 20 6d 70 66 5f 6f 75 74 5f 73 74 72 ;.. mpf_out_str
1da0: 20 28 73 74 64 65 72 72 2c 20 31 30 2c 20 32 2c (stderr, 10, 2,
1db0: 20 56 29 3b 0a 09 20 20 66 70 72 69 6e 74 66 20 V);.. fprintf
1dc0: 28 73 74 64 65 72 72 2c 20 22 5c 74 22 29 3b 0a (stderr, "\t");.
1dd0: 09 7d 0a 20 20 20 20 7d 0a 20 20 69 66 20 28 67 .}. }. if (g
1de0: 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 5f 32 _debug > DEBUG_2
1df0: 29 0a 20 20 20 20 66 70 72 69 6e 74 66 20 28 73 ). fprintf (s
1e00: 74 64 65 72 72 2c 20 22 5c 6e 22 29 3b 0a 20 20 tderr, "\n");.
1e10: 6d 70 66 5f 64 69 76 5f 75 69 20 28 56 2c 20 56 mpf_div_ui (V, V
1e20: 2c 20 6e 29 3b 0a 20 20 6d 70 66 5f 73 75 62 5f , n);. mpf_sub_
1e30: 75 69 20 28 56 2c 20 56 2c 20 6e 29 3b 0a 20 20 ui (V, V, n);.
1e40: 0a 20 20 6d 70 66 5f 63 6c 65 61 72 20 28 66 5f . mpf_clear (f_
1e50: 74 29 3b 20 6d 70 66 5f 63 6c 65 61 72 20 28 66 t); mpf_clear (f
1e60: 5f 74 32 29 3b 0a 7d 0a 0a 2f 2a 20 50 7a 66 28 _t2);.}../* Pzf(
1e70: 70 2c 20 73 2c 20 78 29 20 2d 2d 20 50 72 6f 62 p, s, x) -- Prob
1e80: 61 62 69 6c 69 74 79 20 66 6f 72 20 63 61 74 65 ability for cate
1e90: 67 6f 72 79 20 53 20 69 6e 20 6d 70 7a 5f 66 72 gory S in mpz_fr
1ea0: 65 71 74 28 29 2e 20 20 49 74 27 73 0a 20 20 20 eqt(). It's.
1eb0: 31 2f 64 20 66 6f 72 20 61 6c 6c 20 53 2e 20 20 1/d for all S.
1ec0: 58 20 69 73 20 61 20 70 6f 69 6e 74 65 72 20 74 X is a pointer t
1ed0: 6f 20 61 6e 20 75 6e 73 69 67 6e 65 64 20 69 6e o an unsigned in
1ee0: 74 20 68 6f 6c 64 69 6e 67 20 27 64 27 2e 20 2a t holding 'd'. *
1ef0: 2f 0a 73 74 61 74 69 63 20 76 6f 69 64 0a 50 7a /.static void.Pz
1f00: 66 20 28 6d 70 66 5f 74 20 70 2c 20 75 6e 73 69 f (mpf_t p, unsi
1f10: 67 6e 65 64 20 6c 6f 6e 67 20 69 6e 74 20 73 2c gned long int s,
1f20: 20 76 6f 69 64 20 2a 78 29 0a 7b 0a 20 20 6d 70 void *x).{. mp
1f30: 66 5f 73 65 74 5f 75 69 20 28 70 2c 20 31 29 3b f_set_ui (p, 1);
1f40: 0a 20 20 6d 70 66 5f 64 69 76 5f 75 69 20 28 70 . mpf_div_ui (p
1f50: 2c 20 70 2c 20 2a 28 28 75 6e 73 69 67 6e 65 64 , p, *((unsigned
1f60: 20 69 6e 74 20 2a 29 20 78 29 29 3b 0a 7d 0a 0a int *) x));.}..
1f70: 2f 2a 20 6d 70 7a 5f 66 72 65 71 74 28 56 2c 20 /* mpz_freqt(V,
1f80: 58 2c 20 69 6d 61 78 2c 20 6e 29 20 2d 2d 20 46 X, imax, n) -- F
1f90: 72 65 71 75 65 6e 63 79 20 74 65 73 74 20 6f 6e requency test on
1fa0: 20 69 6e 74 65 67 65 72 73 2e 20 20 5b 4b 6e 75 integers. [Knu
1fb0: 74 68 2c 0a 20 20 20 76 6f 6c 20 32 2c 20 33 2e th,. vol 2, 3.
1fc0: 33 2e 32 5d 2e 20 20 4b 65 65 70 20 49 4d 41 58 3.2]. Keep IMAX
1fd0: 20 6c 6f 77 20 6f 6e 20 74 68 69 73 20 6f 6e 65 low on this one
1fe0: 2c 20 73 69 6e 63 65 20 77 65 20 6c 6f 6f 70 20 , since we loop
1ff0: 66 72 6f 6d 20 30 20 74 6f 0a 20 20 20 49 4d 41 from 0 to. IMA
2000: 58 2e 20 20 31 32 38 20 6f 72 20 32 35 36 20 63 X. 128 or 256 c
2010: 6f 75 6c 64 20 62 65 20 6e 69 63 65 2e 0a 0a 20 ould be nice...
2020: 20 20 58 5b 5d 20 6d 75 73 74 20 6e 6f 74 20 63 X[] must not c
2030: 6f 6e 74 61 69 6e 20 6e 75 6d 62 65 72 73 20 6f ontain numbers o
2040: 75 74 73 69 64 65 20 74 68 65 20 72 61 6e 67 65 utside the range
2050: 20 30 20 3c 3d 20 58 20 3c 3d 20 49 4d 41 58 2e 0 <= X <= IMAX.
2060: 0a 0a 20 20 20 52 65 74 75 72 6e 20 76 61 6c 75 .. Return valu
2070: 65 20 69 73 20 6e 75 6d 62 65 72 20 6f 66 20 6f e is number of o
2080: 62 73 65 72 76 61 74 69 6f 6e 73 20 61 63 74 61 bservations acta
2090: 6c 6c 79 20 75 73 65 64 2c 20 61 66 74 65 72 0a lly used, after.
20a0: 20 20 20 64 69 73 63 61 72 64 69 6e 67 20 65 6e discarding en
20b0: 74 72 69 65 73 20 6f 75 74 20 6f 66 20 72 61 6e tries out of ran
20c0: 67 65 2e 0a 0a 20 20 20 53 69 6e 63 65 20 58 5b ge... Since X[
20d0: 5d 20 63 6f 6e 74 61 69 6e 73 20 69 6e 74 65 67 ] contains integ
20e0: 65 72 73 20 62 65 74 77 65 65 6e 20 7a 65 72 6f ers between zero
20f0: 20 61 6e 64 20 49 4d 41 58 2c 20 69 6e 63 6c 75 and IMAX, inclu
2100: 73 69 76 65 2c 20 77 65 0a 20 20 20 68 61 76 65 sive, we. have
2110: 20 49 4d 41 58 2b 31 20 63 61 74 65 67 6f 72 69 IMAX+1 categori
2120: 65 73 2e 0a 0a 20 20 20 4e 6f 74 65 20 74 68 61 es... Note tha
2130: 74 20 4e 20 73 68 6f 75 6c 64 20 62 65 20 61 74 t N should be at
2140: 20 6c 65 61 73 74 20 35 2a 49 4d 41 58 2e 20 20 least 5*IMAX.
2150: 52 65 73 75 6c 74 20 69 73 20 70 75 74 20 69 6e Result is put in
2160: 20 56 20 61 6e 64 20 63 61 6e 0a 20 20 20 62 65 V and can. be
2170: 20 63 6f 6d 70 61 72 65 64 20 74 6f 20 6f 75 74 compared to out
2180: 70 75 74 20 66 72 6f 6d 20 78 32 5f 74 61 62 6c put from x2_tabl
2190: 65 20 28 76 3d 49 4d 41 58 29 2e 20 2a 2f 0a 0a e (v=IMAX). */..
21a0: 75 6e 73 69 67 6e 65 64 20 6c 6f 6e 67 20 69 6e unsigned long in
21b0: 74 0a 6d 70 7a 5f 66 72 65 71 74 20 28 6d 70 66 t.mpz_freqt (mpf
21c0: 5f 74 20 56 2c 0a 09 20 20 20 6d 70 7a 5f 74 20 _t V,.. mpz_t
21d0: 58 5b 5d 2c 0a 09 20 20 20 75 6e 73 69 67 6e 65 X[],.. unsigne
21e0: 64 20 69 6e 74 20 69 6d 61 78 2c 0a 09 20 20 20 d int imax,..
21f0: 63 6f 6e 73 74 20 75 6e 73 69 67 6e 65 64 20 6c const unsigned l
2200: 6f 6e 67 20 69 6e 74 20 6e 29 0a 7b 0a 20 20 75 ong int n).{. u
2210: 6e 73 69 67 6e 65 64 20 6c 6f 6e 67 20 69 6e 74 nsigned long int
2220: 20 2a 76 3b 09 09 2f 2a 20 72 65 73 75 6c 74 20 *v;../* result
2230: 2a 2f 0a 20 20 75 6e 73 69 67 6e 65 64 20 69 6e */. unsigned in
2240: 74 20 66 3b 0a 20 20 75 6e 73 69 67 6e 65 64 20 t f;. unsigned
2250: 69 6e 74 20 64 3b 09 09 2f 2a 20 6e 75 6d 62 65 int d;../* numbe
2260: 72 20 6f 66 20 63 61 74 65 67 6f 72 69 65 73 20 r of categories
2270: 3d 20 69 6d 61 78 2b 31 20 2a 2f 0a 20 20 75 6e = imax+1 */. un
2280: 73 69 67 6e 65 64 20 69 6e 74 20 75 69 74 65 6d signed int uitem
2290: 70 3b 0a 20 20 75 6e 73 69 67 6e 65 64 20 6c 6f p;. unsigned lo
22a0: 6e 67 20 69 6e 74 20 75 73 65 64 6e 3b 0a 0a 0a ng int usedn;...
22b0: 20 20 64 20 3d 20 69 6d 61 78 20 2b 20 31 3b 0a d = imax + 1;.
22c0: 0a 20 20 76 20 3d 20 28 75 6e 73 69 67 6e 65 64 . v = (unsigned
22d0: 20 6c 6f 6e 67 20 69 6e 74 20 2a 29 20 63 61 6c long int *) cal
22e0: 6c 6f 63 20 28 69 6d 61 78 20 2b 20 31 2c 20 73 loc (imax + 1, s
22f0: 69 7a 65 6f 66 20 28 75 6e 73 69 67 6e 65 64 20 izeof (unsigned
2300: 6c 6f 6e 67 20 69 6e 74 29 29 3b 0a 20 20 69 66 long int));. if
2310: 20 28 4e 55 4c 4c 20 3d 3d 20 76 29 0a 20 20 20 (NULL == v).
2320: 20 7b 0a 20 20 20 20 20 20 66 70 72 69 6e 74 66 {. fprintf
2330: 20 28 73 74 64 65 72 72 2c 20 22 6d 70 7a 5f 66 (stderr, "mpz_f
2340: 72 65 71 74 28 29 3a 20 6f 75 74 20 6f 66 20 6d reqt(): out of m
2350: 65 6d 6f 72 79 5c 6e 22 29 3b 0a 20 20 20 20 20 emory\n");.
2360: 20 65 78 69 74 20 28 31 29 3b 0a 20 20 20 20 7d exit (1);. }
2370: 0a 0a 20 20 2f 2a 20 63 6f 75 6e 74 20 2a 2f 0a .. /* count */.
2380: 20 20 75 73 65 64 6e 20 3d 20 6e 3b 09 09 09 2f usedn = n;.../
2390: 2a 20 61 63 74 75 61 6c 20 6e 75 6d 62 65 72 20 * actual number
23a0: 6f 66 20 6f 62 73 65 72 76 61 74 69 6f 6e 73 20 of observations
23b0: 2a 2f 0a 20 20 66 6f 72 20 28 66 20 3d 20 30 3b */. for (f = 0;
23c0: 20 66 20 3c 20 6e 3b 20 66 2b 2b 29 0a 20 20 20 f < n; f++).
23d0: 20 7b 0a 20 20 20 20 20 20 75 69 74 65 6d 70 20 {. uitemp
23e0: 3d 20 6d 70 7a 5f 67 65 74 5f 75 69 28 58 5b 66 = mpz_get_ui(X[f
23f0: 5d 29 3b 0a 20 20 20 20 20 20 69 66 20 28 75 69 ]);. if (ui
2400: 74 65 6d 70 20 3e 20 69 6d 61 78 29 09 2f 2a 20 temp > imax)./*
2410: 73 61 6e 69 74 79 20 63 68 65 63 6b 20 2a 2f 0a sanity check */.
2420: 09 7b 0a 09 20 20 69 66 20 28 67 5f 64 65 62 75 .{.. if (g_debu
2430: 67 29 0a 09 20 20 20 20 66 70 72 69 6e 74 66 20 g).. fprintf
2440: 28 73 74 64 65 72 72 2c 20 22 6d 70 7a 5f 66 72 (stderr, "mpz_fr
2450: 65 71 74 28 29 3a 20 77 61 72 6e 69 6e 67 3a 20 eqt(): warning:
2460: 69 6e 70 75 74 20 69 6e 73 61 6e 69 74 79 3a 20 input insanity:
2470: 25 75 2c 20 22 5c 0a 09 09 20 20 20 20 20 22 69 %u, "\... "i
2480: 67 6e 6f 72 65 64 2e 5c 6e 22 2c 20 75 69 74 65 gnored.\n", uite
2490: 6d 70 29 3b 0a 09 20 20 75 73 65 64 6e 2d 2d 3b mp);.. usedn--;
24a0: 0a 09 20 20 63 6f 6e 74 69 6e 75 65 3b 0a 09 7d .. continue;..}
24b0: 0a 20 20 20 20 20 20 76 5b 75 69 74 65 6d 70 5d . v[uitemp]
24c0: 2b 2b 3b 0a 20 20 20 20 7d 0a 0a 20 20 69 66 20 ++;. }.. if
24d0: 28 67 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 (g_debug > DEBUG
24e0: 5f 32 29 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 _2). {.
24f0: 66 70 72 69 6e 74 66 20 28 73 74 64 65 72 72 2c fprintf (stderr,
2500: 20 22 63 6f 75 6e 74 73 3a 5c 6e 22 29 3b 0a 20 "counts:\n");.
2510: 20 20 20 20 20 66 6f 72 20 28 66 20 3d 20 30 3b for (f = 0;
2520: 20 66 20 3c 3d 20 69 6d 61 78 3b 20 66 2b 2b 29 f <= imax; f++)
2530: 0a 09 66 70 72 69 6e 74 66 20 28 73 74 64 65 72 ..fprintf (stder
2540: 72 2c 20 22 25 75 3a 5c 74 25 6c 75 5c 6e 22 2c r, "%u:\t%lu\n",
2550: 20 66 2c 20 76 5b 66 5d 29 3b 0a 20 20 20 20 7d f, v[f]);. }
2560: 0a 0a 20 20 2f 2a 20 63 68 69 2d 73 71 75 61 72 .. /* chi-squar
2570: 65 20 77 69 74 68 20 6b 3d 69 6d 61 78 2b 31 20 e with k=imax+1
2580: 61 6e 64 20 50 28 78 29 3d 31 2f 28 69 6d 61 78 and P(x)=1/(imax
2590: 2b 31 29 20 66 6f 72 20 61 6c 6c 20 78 2e 2a 2f +1) for all x.*/
25a0: 0a 20 20 78 32 20 28 56 2c 20 76 2c 20 64 2c 20 . x2 (V, v, d,
25b0: 50 7a 66 2c 20 28 76 6f 69 64 20 2a 29 20 26 64 Pzf, (void *) &d
25c0: 2c 20 75 73 65 64 6e 29 3b 0a 0a 20 20 66 72 65 , usedn);.. fre
25d0: 65 20 28 76 29 3b 0a 20 20 72 65 74 75 72 6e 20 e (v);. return
25e0: 28 75 73 65 64 6e 29 3b 0a 7d 0a 0a 2f 2a 20 64 (usedn);.}../* d
25f0: 65 62 75 67 20 64 75 6d 6d 79 20 74 6f 20 64 72 ebug dummy to dr
2600: 61 67 20 69 6e 20 64 75 6d 70 20 66 75 6e 63 73 ag in dump funcs
2610: 20 2a 2f 0a 76 6f 69 64 0a 66 6f 6f 5f 64 65 62 */.void.foo_deb
2620: 75 67 20 28 29 20 0a 7b 0a 20 20 69 66 20 28 30 ug () .{. if (0
2630: 29 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 6d 70 ). {. mp
2640: 66 5f 64 75 6d 70 20 28 30 29 3b 20 0a 23 69 66 f_dump (0); .#if
2650: 6e 64 65 66 20 4f 4c 44 47 4d 50 0a 20 20 20 20 ndef OLDGMP.
2660: 20 20 6d 70 7a 5f 64 75 6d 70 20 28 30 29 3b 0a mpz_dump (0);.
2670: 23 65 6e 64 69 66 20 20 20 20 20 20 0a 20 20 20 #endif .
2680: 20 7d 0a 7d 0a 0a 2f 2a 20 6d 65 72 69 74 20 28 }.}../* merit (
2690: 72 6f 70 2c 20 74 2c 20 76 2c 20 6d 29 20 2d 2d rop, t, v, m) --
26a0: 20 63 61 6c 63 75 6c 61 74 65 20 6d 65 72 69 74 calculate merit
26b0: 20 66 6f 72 20 73 70 65 63 74 72 61 6c 20 74 65 for spectral te
26c0: 73 74 20 72 65 73 75 6c 74 20 69 6e 0a 20 20 20 st result in.
26d0: 64 69 6d 65 6e 73 69 6f 6e 20 54 2c 20 73 65 65 dimension T, see
26e0: 20 4b 6e 75 74 68 20 70 2e 20 31 30 35 2e 20 20 Knuth p. 105.
26f0: 42 55 47 53 3a 20 4f 6e 6c 79 20 76 61 6c 69 64 BUGS: Only valid
2700: 20 66 6f 72 20 32 20 3c 3d 20 54 20 3c 3d 0a 20 for 2 <= T <=.
2710: 20 20 36 2e 20 2a 2f 0a 76 6f 69 64 0a 6d 65 72 6. */.void.mer
2720: 69 74 20 28 6d 70 66 5f 74 20 72 6f 70 2c 20 75 it (mpf_t rop, u
2730: 6e 73 69 67 6e 65 64 20 69 6e 74 20 74 2c 20 6d nsigned int t, m
2740: 70 66 5f 74 20 76 2c 20 6d 70 7a 5f 74 20 6d 29 pf_t v, mpz_t m)
2750: 0a 7b 0a 20 20 69 6e 74 20 66 3b 0a 20 20 6d 70 .{. int f;. mp
2760: 66 5f 74 20 66 5f 6d 2c 20 66 5f 63 6f 6e 73 74 f_t f_m, f_const
2770: 2c 20 66 5f 70 69 3b 0a 0a 20 20 6d 70 66 5f 69 , f_pi;.. mpf_i
2780: 6e 69 74 20 28 66 5f 6d 29 3b 0a 20 20 6d 70 66 nit (f_m);. mpf
2790: 5f 73 65 74 5f 7a 20 28 66 5f 6d 2c 20 6d 29 3b _set_z (f_m, m);
27a0: 0a 20 20 6d 70 66 5f 69 6e 69 74 5f 73 65 74 5f . mpf_init_set_
27b0: 64 20 28 66 5f 63 6f 6e 73 74 2c 20 4d 5f 50 49 d (f_const, M_PI
27c0: 29 3b 0a 20 20 6d 70 66 5f 69 6e 69 74 5f 73 65 );. mpf_init_se
27d0: 74 5f 64 20 28 66 5f 70 69 2c 20 4d 5f 50 49 29 t_d (f_pi, M_PI)
27e0: 3b 0a 20 20 0a 20 20 73 77 69 74 63 68 20 28 74 ;. . switch (t
27f0: 29 0a 20 20 20 20 7b 0a 20 20 20 20 63 61 73 65 ). {. case
2800: 20 32 3a 09 09 09 2f 2a 20 50 49 20 2a 2f 0a 20 2:.../* PI */.
2810: 20 20 20 20 20 62 72 65 61 6b 3b 0a 20 20 20 20 break;.
2820: 63 61 73 65 20 33 3a 09 09 09 2f 2a 20 50 49 20 case 3:.../* PI
2830: 2a 20 34 2f 33 20 2a 2f 0a 20 20 20 20 20 20 6d * 4/3 */. m
2840: 70 66 5f 6d 75 6c 5f 75 69 20 28 66 5f 63 6f 6e pf_mul_ui (f_con
2850: 73 74 2c 20 66 5f 63 6f 6e 73 74 2c 20 34 29 3b st, f_const, 4);
2860: 0a 20 20 20 20 20 20 6d 70 66 5f 64 69 76 5f 75 . mpf_div_u
2870: 69 20 28 66 5f 63 6f 6e 73 74 2c 20 66 5f 63 6f i (f_const, f_co
2880: 6e 73 74 2c 20 33 29 3b 0a 20 20 20 20 20 20 62 nst, 3);. b
2890: 72 65 61 6b 3b 0a 20 20 20 20 63 61 73 65 20 34 reak;. case 4
28a0: 3a 09 09 09 2f 2a 20 50 49 5e 32 20 2a 20 31 2f :.../* PI^2 * 1/
28b0: 32 20 2a 2f 0a 20 20 20 20 20 20 6d 70 66 5f 6d 2 */. mpf_m
28c0: 75 6c 20 28 66 5f 63 6f 6e 73 74 2c 20 66 5f 63 ul (f_const, f_c
28d0: 6f 6e 73 74 2c 20 66 5f 70 69 29 3b 0a 20 20 20 onst, f_pi);.
28e0: 20 20 20 6d 70 66 5f 64 69 76 5f 75 69 20 28 66 mpf_div_ui (f
28f0: 5f 63 6f 6e 73 74 2c 20 66 5f 63 6f 6e 73 74 2c _const, f_const,
2900: 20 32 29 3b 0a 20 20 20 20 20 20 62 72 65 61 6b 2);. break
2910: 3b 0a 20 20 20 20 63 61 73 65 20 35 3a 09 09 09 ;. case 5:...
2920: 2f 2a 20 50 49 5e 32 20 2a 20 38 2f 31 35 20 2a /* PI^2 * 8/15 *
2930: 2f 0a 20 20 20 20 20 20 6d 70 66 5f 6d 75 6c 20 /. mpf_mul
2940: 28 66 5f 63 6f 6e 73 74 2c 20 66 5f 63 6f 6e 73 (f_const, f_cons
2950: 74 2c 20 66 5f 70 69 29 3b 0a 20 20 20 20 20 20 t, f_pi);.
2960: 6d 70 66 5f 6d 75 6c 5f 75 69 20 28 66 5f 63 6f mpf_mul_ui (f_co
2970: 6e 73 74 2c 20 66 5f 63 6f 6e 73 74 2c 20 38 29 nst, f_const, 8)
2980: 3b 0a 20 20 20 20 20 20 6d 70 66 5f 64 69 76 5f ;. mpf_div_
2990: 75 69 20 28 66 5f 63 6f 6e 73 74 2c 20 66 5f 63 ui (f_const, f_c
29a0: 6f 6e 73 74 2c 20 31 35 29 3b 0a 20 20 20 20 20 onst, 15);.
29b0: 20 62 72 65 61 6b 3b 0a 20 20 20 20 63 61 73 65 break;. case
29c0: 20 36 3a 09 09 09 2f 2a 20 50 49 5e 33 20 2a 20 6:.../* PI^3 *
29d0: 31 2f 36 20 2a 2f 0a 20 20 20 20 20 20 6d 70 66 1/6 */. mpf
29e0: 5f 6d 75 6c 20 28 66 5f 63 6f 6e 73 74 2c 20 66 _mul (f_const, f
29f0: 5f 63 6f 6e 73 74 2c 20 66 5f 70 69 29 3b 0a 20 _const, f_pi);.
2a00: 20 20 20 20 20 6d 70 66 5f 6d 75 6c 20 28 66 5f mpf_mul (f_
2a10: 63 6f 6e 73 74 2c 20 66 5f 63 6f 6e 73 74 2c 20 const, f_const,
2a20: 66 5f 70 69 29 3b 0a 20 20 20 20 20 20 6d 70 66 f_pi);. mpf
2a30: 5f 64 69 76 5f 75 69 20 28 66 5f 63 6f 6e 73 74 _div_ui (f_const
2a40: 2c 20 66 5f 63 6f 6e 73 74 2c 20 36 29 3b 0a 20 , f_const, 6);.
2a50: 20 20 20 20 20 62 72 65 61 6b 3b 0a 20 20 20 20 break;.
2a60: 64 65 66 61 75 6c 74 3a 0a 20 20 20 20 20 20 66 default:. f
2a70: 70 72 69 6e 74 66 20 28 73 74 64 65 72 72 2c 0a printf (stderr,.
2a80: 09 20 20 20 20 20 20 20 22 73 70 65 63 74 20 28 . "spect (
2a90: 6d 65 72 69 74 29 3a 20 63 61 6e 27 74 20 63 61 merit): can't ca
2aa0: 6c 63 75 6c 61 74 65 20 6d 65 72 69 74 20 66 6f lculate merit fo
2ab0: 72 20 64 69 6d 65 6e 73 69 6f 6e 73 20 3e 20 36 r dimensions > 6
2ac0: 5c 6e 22 29 3b 0a 20 20 20 20 20 20 6d 70 66 5f \n");. mpf_
2ad0: 73 65 74 5f 75 69 20 28 66 5f 63 6f 6e 73 74 2c set_ui (f_const,
2ae0: 20 30 29 3b 0a 20 20 20 20 20 20 62 72 65 61 6b 0);. break
2af0: 3b 0a 20 20 20 20 7d 0a 0a 20 20 2f 2a 20 72 6f ;. }.. /* ro
2b00: 70 20 3d 20 76 5e 74 20 2a 2f 0a 20 20 6d 70 66 p = v^t */. mpf
2b10: 5f 73 65 74 20 28 72 6f 70 2c 20 76 29 3b 0a 20 _set (rop, v);.
2b20: 20 66 6f 72 20 28 66 20 3d 20 31 3b 20 66 20 3c for (f = 1; f <
2b30: 20 74 3b 20 66 2b 2b 29 0a 20 20 20 20 6d 70 66 t; f++). mpf
2b40: 5f 6d 75 6c 20 28 72 6f 70 2c 20 72 6f 70 2c 20 _mul (rop, rop,
2b50: 76 29 3b 0a 20 20 6d 70 66 5f 6d 75 6c 20 28 72 v);. mpf_mul (r
2b60: 6f 70 2c 20 72 6f 70 2c 20 66 5f 63 6f 6e 73 74 op, rop, f_const
2b70: 29 3b 0a 20 20 6d 70 66 5f 64 69 76 20 28 72 6f );. mpf_div (ro
2b80: 70 2c 20 72 6f 70 2c 20 66 5f 6d 29 3b 0a 0a 20 p, rop, f_m);..
2b90: 20 6d 70 66 5f 63 6c 65 61 72 20 28 66 5f 6d 29 mpf_clear (f_m)
2ba0: 3b 0a 20 20 6d 70 66 5f 63 6c 65 61 72 20 28 66 ;. mpf_clear (f
2bb0: 5f 63 6f 6e 73 74 29 3b 0a 20 20 6d 70 66 5f 63 _const);. mpf_c
2bc0: 6c 65 61 72 20 28 66 5f 70 69 29 3b 0a 7d 0a 0a lear (f_pi);.}..
2bd0: 64 6f 75 62 6c 65 0a 6d 65 72 69 74 5f 75 20 28 double.merit_u (
2be0: 75 6e 73 69 67 6e 65 64 20 69 6e 74 20 74 2c 20 unsigned int t,
2bf0: 6d 70 66 5f 74 20 76 2c 20 6d 70 7a 5f 74 20 6d mpf_t v, mpz_t m
2c00: 29 0a 7b 0a 20 20 6d 70 66 5f 74 20 72 6f 70 3b ).{. mpf_t rop;
2c10: 0a 20 20 64 6f 75 62 6c 65 20 72 65 73 3b 0a 20 . double res;.
2c20: 20 0a 20 20 6d 70 66 5f 69 6e 69 74 20 28 72 6f . mpf_init (ro
2c30: 70 29 3b 0a 20 20 6d 65 72 69 74 20 28 72 6f 70 p);. merit (rop
2c40: 2c 20 74 2c 20 76 2c 20 6d 29 3b 0a 20 20 72 65 , t, v, m);. re
2c50: 73 20 3d 20 6d 70 66 5f 67 65 74 5f 64 20 28 72 s = mpf_get_d (r
2c60: 6f 70 29 3b 0a 20 20 6d 70 66 5f 63 6c 65 61 72 op);. mpf_clear
2c70: 20 28 72 6f 70 29 3b 0a 20 20 72 65 74 75 72 6e (rop);. return
2c80: 20 72 65 73 3b 0a 7d 0a 0a 2f 2a 20 66 5f 66 6c res;.}../* f_fl
2c90: 6f 6f 72 20 28 72 6f 70 2c 20 6f 70 29 20 2d 2d oor (rop, op) --
2ca0: 20 53 65 74 20 72 6f 70 20 3d 20 66 6c 6f 6f 72 Set rop = floor
2cb0: 20 28 6f 70 29 2e 20 2a 2f 0a 76 6f 69 64 0a 66 (op). */.void.f
2cc0: 5f 66 6c 6f 6f 72 20 28 6d 70 66 5f 74 20 72 6f _floor (mpf_t ro
2cd0: 70 2c 20 6d 70 66 5f 74 20 6f 70 29 0a 7b 0a 20 p, mpf_t op).{.
2ce0: 20 6d 70 7a 5f 74 20 7a 3b 0a 0a 20 20 6d 70 7a mpz_t z;.. mpz
2cf0: 5f 69 6e 69 74 20 28 7a 29 3b 0a 0a 20 20 2f 2a _init (z);.. /*
2d00: 20 4e 6f 20 6d 70 66 5f 66 6c 6f 6f 72 28 29 2e No mpf_floor().
2d10: 20 20 43 6f 6e 76 65 72 74 20 74 6f 20 6d 70 7a Convert to mpz
2d20: 20 61 6e 64 20 62 61 63 6b 2e 20 2a 2f 0a 20 20 and back. */.
2d30: 6d 70 7a 5f 73 65 74 5f 66 20 28 7a 2c 20 6f 70 mpz_set_f (z, op
2d40: 29 3b 0a 20 20 6d 70 66 5f 73 65 74 5f 7a 20 28 );. mpf_set_z (
2d50: 72 6f 70 2c 20 7a 29 3b 0a 0a 20 20 6d 70 7a 5f rop, z);.. mpz_
2d60: 63 6c 65 61 72 20 28 7a 29 3b 0a 7d 0a 0a 0a 2f clear (z);.}.../
2d70: 2a 20 76 7a 5f 64 6f 74 20 28 72 6f 70 2c 20 76 * vz_dot (rop, v
2d80: 31 2c 20 76 32 2c 20 6e 65 6c 65 6d 29 20 2d 2d 1, v2, nelem) --
2d90: 20 63 6f 6d 70 75 74 65 20 64 6f 74 20 70 72 6f compute dot pro
2da0: 64 75 63 74 20 6f 66 20 7a 2d 76 65 63 74 6f 72 duct of z-vector
2db0: 73 20 56 31 2c 0a 20 20 20 56 32 2e 20 20 4e 20 s V1,. V2. N
2dc0: 69 73 20 6e 75 6d 62 65 72 20 6f 66 20 65 6c 65 is number of ele
2dd0: 6d 65 6e 74 73 20 69 6e 20 76 65 63 74 6f 72 73 ments in vectors
2de0: 20 56 31 20 61 6e 64 20 56 32 2e 20 2a 2f 0a 0a V1 and V2. */..
2df0: 76 6f 69 64 0a 76 7a 5f 64 6f 74 20 28 6d 70 7a void.vz_dot (mpz
2e00: 5f 74 20 72 6f 70 2c 20 6d 70 7a 5f 74 20 56 31 _t rop, mpz_t V1
2e10: 5b 5d 2c 20 6d 70 7a 5f 74 20 56 32 5b 5d 2c 20 [], mpz_t V2[],
2e20: 75 6e 73 69 67 6e 65 64 20 69 6e 74 20 6e 29 0a unsigned int n).
2e30: 7b 0a 20 20 6d 70 7a 5f 74 20 74 3b 0a 0a 20 20 {. mpz_t t;..
2e40: 6d 70 7a 5f 69 6e 69 74 20 28 74 29 3b 0a 20 20 mpz_init (t);.
2e50: 6d 70 7a 5f 73 65 74 5f 75 69 20 28 72 6f 70 2c mpz_set_ui (rop,
2e60: 20 30 29 3b 0a 20 20 77 68 69 6c 65 20 28 6e 2d 0);. while (n-
2e70: 2d 29 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 6d -). {. m
2e80: 70 7a 5f 6d 75 6c 20 28 74 2c 20 56 31 5b 6e 5d pz_mul (t, V1[n]
2e90: 2c 20 56 32 5b 6e 5d 29 3b 0a 20 20 20 20 20 20 , V2[n]);.
2ea0: 6d 70 7a 5f 61 64 64 20 28 72 6f 70 2c 20 72 6f mpz_add (rop, ro
2eb0: 70 2c 20 74 29 3b 0a 20 20 20 20 7d 0a 0a 20 20 p, t);. }..
2ec0: 6d 70 7a 5f 63 6c 65 61 72 20 28 74 29 3b 0a 7d mpz_clear (t);.}
2ed0: 0a 0a 76 6f 69 64 0a 73 70 65 63 74 72 61 6c 5f ..void.spectral_
2ee0: 74 65 73 74 20 28 6d 70 66 5f 74 20 72 6f 70 5b test (mpf_t rop[
2ef0: 5d 2c 20 75 6e 73 69 67 6e 65 64 20 69 6e 74 20 ], unsigned int
2f00: 54 2c 20 6d 70 7a 5f 74 20 61 2c 20 6d 70 7a 5f T, mpz_t a, mpz_
2f10: 74 20 6d 29 0a 7b 0a 20 20 2f 2a 20 4b 6e 75 74 t m).{. /* Knut
2f20: 68 20 22 53 65 6d 69 6e 75 6d 65 72 69 63 61 6c h "Seminumerical
2f30: 20 41 6c 67 6f 72 69 74 68 6d 73 2c 20 54 68 69 Algorithms, Thi
2f40: 72 64 20 45 64 69 74 69 6f 6e 22 2c 20 73 65 63 rd Edition", sec
2f50: 74 69 6f 6e 20 33 2e 33 2e 34 0a 20 20 20 20 20 tion 3.3.4.
2f60: 28 70 70 2e 20 31 30 31 2d 31 30 33 29 2e 20 2a (pp. 101-103). *
2f70: 2f 0a 0a 20 20 2f 2a 20 76 5b 74 5d 20 3d 20 6d /.. /* v[t] = m
2f80: 69 6e 20 7b 20 73 71 72 74 20 28 78 5b 31 5d 5e in { sqrt (x[1]^
2f90: 32 20 2b 20 2e 2e 2e 20 2b 20 78 5b 74 5d 5e 32 2 + ... + x[t]^2
2fa0: 29 20 7c 0a 20 20 20 20 20 78 5b 31 5d 20 2b 20 ) |. x[1] +
2fb0: 61 2a 78 5b 32 5d 20 2b 20 2e 2e 2e 20 2b 20 70 a*x[2] + ... + p
2fc0: 6f 77 20 28 61 2c 20 74 2d 31 29 20 2a 20 78 5b ow (a, t-1) * x[
2fd0: 74 5d 20 69 73 20 63 6f 6e 67 72 75 65 6e 74 20 t] is congruent
2fe0: 74 6f 20 30 20 28 6d 6f 64 20 6d 29 20 7d 20 2a to 0 (mod m) } *
2ff0: 2f 0a 0a 0a 20 20 2f 2a 20 56 61 72 69 61 62 6c /... /* Variabl
3000: 65 73 2e 20 2a 2f 0a 20 20 75 6e 73 69 67 6e 65 es. */. unsigne
3010: 64 20 69 6e 74 20 75 69 5f 74 3b 0a 20 20 75 6e d int ui_t;. un
3020: 73 69 67 6e 65 64 20 69 6e 74 20 75 69 5f 69 2c signed int ui_i,
3030: 20 75 69 5f 6a 2c 20 75 69 5f 6b 2c 20 75 69 5f ui_j, ui_k, ui_
3040: 6c 3b 0a 20 20 6d 70 66 5f 74 20 66 5f 74 6d 70 l;. mpf_t f_tmp
3050: 31 2c 20 66 5f 74 6d 70 32 3b 0a 20 20 6d 70 7a 1, f_tmp2;. mpz
3060: 5f 74 20 74 6d 70 31 2c 20 74 6d 70 32 2c 20 74 _t tmp1, tmp2, t
3070: 6d 70 33 3b 0a 20 20 6d 70 7a 5f 74 20 55 5b 47 mp3;. mpz_t U[G
3080: 4d 50 5f 53 50 45 43 54 5f 4d 41 58 54 5d 5b 47 MP_SPECT_MAXT][G
3090: 4d 50 5f 53 50 45 43 54 5f 4d 41 58 54 5d 2c 0a MP_SPECT_MAXT],.
30a0: 20 20 20 20 56 5b 47 4d 50 5f 53 50 45 43 54 5f V[GMP_SPECT_
30b0: 4d 41 58 54 5d 5b 47 4d 50 5f 53 50 45 43 54 5f MAXT][GMP_SPECT_
30c0: 4d 41 58 54 5d 2c 0a 20 20 20 20 58 5b 47 4d 50 MAXT],. X[GMP
30d0: 5f 53 50 45 43 54 5f 4d 41 58 54 5d 2c 0a 20 20 _SPECT_MAXT],.
30e0: 20 20 59 5b 47 4d 50 5f 53 50 45 43 54 5f 4d 41 Y[GMP_SPECT_MA
30f0: 58 54 5d 2c 0a 20 20 20 20 5a 5b 47 4d 50 5f 53 XT],. Z[GMP_S
3100: 50 45 43 54 5f 4d 41 58 54 5d 3b 0a 20 20 6d 70 PECT_MAXT];. mp
3110: 7a 5f 74 20 68 2c 20 68 70 2c 20 72 2c 20 73 2c z_t h, hp, r, s,
3120: 20 70 2c 20 70 70 2c 20 71 2c 20 75 2c 20 76 3b p, pp, q, u, v;
3130: 0a 0a 20 20 2f 2a 20 47 4d 50 20 69 6e 69 74 73 .. /* GMP inits
3140: 2e 20 2a 2f 0a 20 20 6d 70 66 5f 69 6e 69 74 20 . */. mpf_init
3150: 28 66 5f 74 6d 70 31 29 3b 0a 20 20 6d 70 66 5f (f_tmp1);. mpf_
3160: 69 6e 69 74 20 28 66 5f 74 6d 70 32 29 3b 0a 20 init (f_tmp2);.
3170: 20 66 6f 72 20 28 75 69 5f 69 20 3d 20 30 3b 20 for (ui_i = 0;
3180: 75 69 5f 69 20 3c 20 47 4d 50 5f 53 50 45 43 54 ui_i < GMP_SPECT
3190: 5f 4d 41 58 54 3b 20 75 69 5f 69 2b 2b 29 0a 20 _MAXT; ui_i++).
31a0: 20 20 20 7b 0a 20 20 20 20 20 20 66 6f 72 20 28 {. for (
31b0: 75 69 5f 6a 20 3d 20 30 3b 20 75 69 5f 6a 20 3c ui_j = 0; ui_j <
31c0: 20 47 4d 50 5f 53 50 45 43 54 5f 4d 41 58 54 3b GMP_SPECT_MAXT;
31d0: 20 75 69 5f 6a 2b 2b 29 0a 09 7b 0a 09 20 20 6d ui_j++)..{.. m
31e0: 70 7a 5f 69 6e 69 74 5f 73 65 74 5f 75 69 20 28 pz_init_set_ui (
31f0: 55 5b 75 69 5f 69 5d 5b 75 69 5f 6a 5d 2c 20 30 U[ui_i][ui_j], 0
3200: 29 3b 0a 09 20 20 6d 70 7a 5f 69 6e 69 74 5f 73 );.. mpz_init_s
3210: 65 74 5f 75 69 20 28 56 5b 75 69 5f 69 5d 5b 75 et_ui (V[ui_i][u
3220: 69 5f 6a 5d 2c 20 30 29 3b 0a 09 7d 0a 20 20 20 i_j], 0);..}.
3230: 20 20 20 6d 70 7a 5f 69 6e 69 74 5f 73 65 74 5f mpz_init_set_
3240: 75 69 20 28 58 5b 75 69 5f 69 5d 2c 20 30 29 3b ui (X[ui_i], 0);
3250: 0a 20 20 20 20 20 20 6d 70 7a 5f 69 6e 69 74 5f . mpz_init_
3260: 73 65 74 5f 75 69 20 28 59 5b 75 69 5f 69 5d 2c set_ui (Y[ui_i],
3270: 20 30 29 3b 0a 20 20 20 20 20 20 6d 70 7a 5f 69 0);. mpz_i
3280: 6e 69 74 20 28 5a 5b 75 69 5f 69 5d 29 3b 0a 20 nit (Z[ui_i]);.
3290: 20 20 20 7d 0a 20 20 6d 70 7a 5f 69 6e 69 74 20 }. mpz_init
32a0: 28 74 6d 70 31 29 3b 0a 20 20 6d 70 7a 5f 69 6e (tmp1);. mpz_in
32b0: 69 74 20 28 74 6d 70 32 29 3b 0a 20 20 6d 70 7a it (tmp2);. mpz
32c0: 5f 69 6e 69 74 20 28 74 6d 70 33 29 3b 0a 20 20 _init (tmp3);.
32d0: 6d 70 7a 5f 69 6e 69 74 20 28 68 29 3b 0a 20 20 mpz_init (h);.
32e0: 6d 70 7a 5f 69 6e 69 74 20 28 68 70 29 3b 0a 20 mpz_init (hp);.
32f0: 20 6d 70 7a 5f 69 6e 69 74 20 28 72 29 3b 0a 20 mpz_init (r);.
3300: 20 6d 70 7a 5f 69 6e 69 74 20 28 73 29 3b 0a 20 mpz_init (s);.
3310: 20 6d 70 7a 5f 69 6e 69 74 20 28 70 29 3b 0a 20 mpz_init (p);.
3320: 20 6d 70 7a 5f 69 6e 69 74 20 28 70 70 29 3b 0a mpz_init (pp);.
3330: 20 20 6d 70 7a 5f 69 6e 69 74 20 28 71 29 3b 0a mpz_init (q);.
3340: 20 20 6d 70 7a 5f 69 6e 69 74 20 28 75 29 3b 0a mpz_init (u);.
3350: 20 20 6d 70 7a 5f 69 6e 69 74 20 28 76 29 3b 0a mpz_init (v);.
3360: 0a 20 20 2f 2a 20 49 6d 70 6c 65 6d 65 6e 74 61 . /* Implementa
3370: 74 69 6f 6e 20 69 6e 69 74 73 2e 20 2a 2f 0a 20 tion inits. */.
3380: 20 69 66 20 28 54 20 3e 20 47 4d 50 5f 53 50 45 if (T > GMP_SPE
3390: 43 54 5f 4d 41 58 54 29 0a 20 20 20 20 54 20 3d CT_MAXT). T =
33a0: 20 47 4d 50 5f 53 50 45 43 54 5f 4d 41 58 54 3b GMP_SPECT_MAXT;
33b0: 09 09 09 2f 2a 20 46 49 58 4d 45 3a 20 4c 61 7a .../* FIXME: Laz
33c0: 79 2e 20 2a 2f 0a 0a 20 20 2f 2a 20 53 31 20 5b y. */.. /* S1 [
33d0: 49 6e 69 74 69 61 6c 69 7a 65 2e 5d 20 2a 2f 0a Initialize.] */.
33e0: 20 20 75 69 5f 74 20 3d 20 32 20 2d 20 31 3b 09 ui_t = 2 - 1;.
33f0: 09 09 2f 2a 20 4e 4f 54 45 3a 20 60 74 27 20 69 ../* NOTE: `t' i
3400: 6e 20 64 65 73 63 72 69 70 74 69 6f 6e 20 3d 3d n description ==
3410: 20 75 69 5f 74 20 2b 20 31 0a 20 20 20 20 20 20 ui_t + 1.
3420: 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20
3430: 20 20 20 20 20 20 20 20 20 20 20 20 20 66 6f 72 for
3440: 20 65 61 73 79 20 69 6e 64 65 78 69 6e 67 20 2a easy indexing *
3450: 2f 0a 20 20 6d 70 7a 5f 73 65 74 20 28 68 2c 20 /. mpz_set (h,
3460: 61 29 3b 0a 20 20 6d 70 7a 5f 73 65 74 20 28 68 a);. mpz_set (h
3470: 70 2c 20 6d 29 3b 0a 20 20 6d 70 7a 5f 73 65 74 p, m);. mpz_set
3480: 5f 75 69 20 28 70 2c 20 31 29 3b 0a 20 20 6d 70 _ui (p, 1);. mp
3490: 7a 5f 73 65 74 5f 75 69 20 28 70 70 2c 20 30 29 z_set_ui (pp, 0)
34a0: 3b 0a 20 20 6d 70 7a 5f 73 65 74 20 28 72 2c 20 ;. mpz_set (r,
34b0: 61 29 3b 0a 20 20 6d 70 7a 5f 70 6f 77 5f 75 69 a);. mpz_pow_ui
34c0: 20 28 73 2c 20 61 2c 20 32 29 3b 0a 20 20 6d 70 (s, a, 2);. mp
34d0: 7a 5f 61 64 64 5f 75 69 20 28 73 2c 20 73 2c 20 z_add_ui (s, s,
34e0: 31 29 3b 09 09 2f 2a 20 73 20 3d 20 31 20 2b 20 1);../* s = 1 +
34f0: 61 5e 32 20 2a 2f 0a 0a 20 20 2f 2a 20 53 32 20 a^2 */.. /* S2
3500: 5b 45 75 63 6c 69 64 65 61 6e 20 73 74 65 70 2e [Euclidean step.
3510: 5d 20 2a 2f 0a 20 20 77 68 69 6c 65 20 28 31 29 ] */. while (1)
3520: 0a 20 20 20 20 7b 0a 20 20 20 20 20 20 69 66 20 . {. if
3530: 28 67 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 (g_debug > DEBUG
3540: 5f 31 29 0a 09 7b 0a 09 20 20 6d 70 7a 5f 6d 75 _1)..{.. mpz_mu
3550: 6c 20 28 74 6d 70 31 2c 20 68 2c 20 70 70 29 3b l (tmp1, h, pp);
3560: 0a 09 20 20 6d 70 7a 5f 6d 75 6c 20 28 74 6d 70 .. mpz_mul (tmp
3570: 32 2c 20 68 70 2c 20 70 29 3b 0a 09 20 20 6d 70 2, hp, p);.. mp
3580: 7a 5f 73 75 62 20 28 74 6d 70 31 2c 20 74 6d 70 z_sub (tmp1, tmp
3590: 31 2c 20 74 6d 70 32 29 3b 0a 09 20 20 69 66 20 1, tmp2);.. if
35a0: 28 6d 70 7a 5f 63 6d 70 61 62 73 20 28 6d 2c 20 (mpz_cmpabs (m,
35b0: 74 6d 70 31 29 29 0a 09 20 20 20 20 7b 0a 09 20 tmp1)).. {..
35c0: 20 20 20 20 20 70 72 69 6e 74 66 20 28 22 2a 2a printf ("**
35d0: 2a 42 55 47 2a 2a 2a 3a 20 68 2a 70 70 20 2d 20 *BUG***: h*pp -
35e0: 68 70 2a 70 20 3d 20 22 29 3b 0a 09 20 20 20 20 hp*p = ");..
35f0: 20 20 6d 70 7a 5f 6f 75 74 5f 73 74 72 20 28 73 mpz_out_str (s
3600: 74 64 6f 75 74 2c 20 31 30 2c 20 74 6d 70 31 29 tdout, 10, tmp1)
3610: 3b 0a 09 20 20 20 20 20 20 70 72 69 6e 74 66 20 ;.. printf
3620: 28 22 5c 6e 22 29 3b 0a 09 20 20 20 20 7d 0a 09 ("\n");.. }..
3630: 7d 0a 20 20 20 20 20 20 69 66 20 28 67 5f 64 65 }. if (g_de
3640: 62 75 67 20 3e 20 44 45 42 55 47 5f 32 29 0a 09 bug > DEBUG_2)..
3650: 7b 0a 09 20 20 70 72 69 6e 74 66 20 28 22 68 70 {.. printf ("hp
3660: 20 3d 20 22 29 3b 0a 09 20 20 6d 70 7a 5f 6f 75 = ");.. mpz_ou
3670: 74 5f 73 74 72 20 28 73 74 64 6f 75 74 2c 20 31 t_str (stdout, 1
3680: 30 2c 20 68 70 29 3b 0a 09 20 20 70 72 69 6e 74 0, hp);.. print
3690: 66 20 28 22 5c 6e 68 20 3d 20 22 29 3b 0a 09 20 f ("\nh = ");..
36a0: 20 6d 70 7a 5f 6f 75 74 5f 73 74 72 20 28 73 74 mpz_out_str (st
36b0: 64 6f 75 74 2c 20 31 30 2c 20 68 29 3b 0a 09 20 dout, 10, h);..
36c0: 20 70 72 69 6e 74 66 20 28 22 5c 6e 22 29 3b 0a printf ("\n");.
36d0: 09 20 20 66 66 6c 75 73 68 20 28 73 74 64 6f 75 . fflush (stdou
36e0: 74 29 3b 0a 09 7d 0a 0a 20 20 20 20 20 20 69 66 t);..}.. if
36f0: 20 28 6d 70 7a 5f 73 67 6e 20 28 68 29 29 0a 09 (mpz_sgn (h))..
3700: 6d 70 7a 5f 74 64 69 76 5f 71 20 28 71 2c 20 68 mpz_tdiv_q (q, h
3710: 70 2c 20 68 29 3b 09 2f 2a 20 71 20 3d 20 66 6c p, h);./* q = fl
3720: 6f 6f 72 28 68 70 2f 68 29 20 2a 2f 0a 20 20 20 oor(hp/h) */.
3730: 20 20 20 65 6c 73 65 0a 09 6d 70 7a 5f 73 65 74 else..mpz_set
3740: 5f 75 69 20 28 71 2c 20 31 29 3b 0a 0a 20 20 20 _ui (q, 1);..
3750: 20 20 20 69 66 20 28 67 5f 64 65 62 75 67 20 3e if (g_debug >
3760: 20 44 45 42 55 47 5f 32 29 0a 09 7b 0a 09 20 20 DEBUG_2)..{..
3770: 70 72 69 6e 74 66 20 28 22 71 20 3d 20 22 29 3b printf ("q = ");
3780: 0a 09 20 20 6d 70 7a 5f 6f 75 74 5f 73 74 72 20 .. mpz_out_str
3790: 28 73 74 64 6f 75 74 2c 20 31 30 2c 20 71 29 3b (stdout, 10, q);
37a0: 0a 09 20 20 70 72 69 6e 74 66 20 28 22 5c 6e 22 .. printf ("\n"
37b0: 29 3b 0a 09 20 20 66 66 6c 75 73 68 20 28 73 74 );.. fflush (st
37c0: 64 6f 75 74 29 3b 0a 09 7d 0a 0a 20 20 20 20 20 dout);..}..
37d0: 20 6d 70 7a 5f 6d 75 6c 20 28 74 6d 70 31 2c 20 mpz_mul (tmp1,
37e0: 71 2c 20 68 29 3b 0a 20 20 20 20 20 20 6d 70 7a q, h);. mpz
37f0: 5f 73 75 62 20 28 75 2c 20 68 70 2c 20 74 6d 70 _sub (u, hp, tmp
3800: 31 29 3b 09 2f 2a 20 75 20 3d 20 68 70 20 2d 20 1);./* u = hp -
3810: 71 2a 68 20 2a 2f 0a 0a 20 20 20 20 20 20 6d 70 q*h */.. mp
3820: 7a 5f 6d 75 6c 20 28 74 6d 70 31 2c 20 71 2c 20 z_mul (tmp1, q,
3830: 70 29 3b 0a 20 20 20 20 20 20 6d 70 7a 5f 73 75 p);. mpz_su
3840: 62 20 28 76 2c 20 70 70 2c 20 74 6d 70 31 29 3b b (v, pp, tmp1);
3850: 09 2f 2a 20 76 20 3d 20 70 70 20 2d 20 71 2a 70 ./* v = pp - q*p
3860: 20 2a 2f 0a 20 20 0a 20 20 20 20 20 20 6d 70 7a */. . mpz
3870: 5f 70 6f 77 5f 75 69 20 28 74 6d 70 31 2c 20 75 _pow_ui (tmp1, u
3880: 2c 20 32 29 3b 0a 20 20 20 20 20 20 6d 70 7a 5f , 2);. mpz_
3890: 70 6f 77 5f 75 69 20 28 74 6d 70 32 2c 20 76 2c pow_ui (tmp2, v,
38a0: 20 32 29 3b 0a 20 20 20 20 20 20 6d 70 7a 5f 61 2);. mpz_a
38b0: 64 64 20 28 74 6d 70 31 2c 20 74 6d 70 31 2c 20 dd (tmp1, tmp1,
38c0: 74 6d 70 32 29 3b 0a 20 20 20 20 20 20 69 66 20 tmp2);. if
38d0: 28 6d 70 7a 5f 63 6d 70 20 28 74 6d 70 31 2c 20 (mpz_cmp (tmp1,
38e0: 73 29 20 3c 20 30 29 0a 09 7b 0a 09 20 20 6d 70 s) < 0)..{.. mp
38f0: 7a 5f 73 65 74 20 28 73 2c 20 74 6d 70 31 29 3b z_set (s, tmp1);
3900: 09 2f 2a 20 73 20 3d 20 75 5e 32 20 2b 20 76 5e ./* s = u^2 + v^
3910: 32 20 2a 2f 0a 09 20 20 6d 70 7a 5f 73 65 74 20 2 */.. mpz_set
3920: 28 68 70 2c 20 68 29 3b 09 2f 2a 20 68 70 20 3d (hp, h);./* hp =
3930: 20 68 20 2a 2f 0a 09 20 20 6d 70 7a 5f 73 65 74 h */.. mpz_set
3940: 20 28 68 2c 20 75 29 3b 09 2f 2a 20 68 20 3d 20 (h, u);./* h =
3950: 75 20 2a 2f 0a 09 20 20 6d 70 7a 5f 73 65 74 20 u */.. mpz_set
3960: 28 70 70 2c 20 70 29 3b 09 2f 2a 20 70 70 20 3d (pp, p);./* pp =
3970: 20 70 20 2a 2f 0a 09 20 20 6d 70 7a 5f 73 65 74 p */.. mpz_set
3980: 20 28 70 2c 20 76 29 3b 09 2f 2a 20 70 20 3d 20 (p, v);./* p =
3990: 76 20 2a 2f 0a 09 7d 0a 20 20 20 20 20 20 65 6c v */..}. el
39a0: 73 65 0a 09 62 72 65 61 6b 3b 0a 20 20 20 20 7d se..break;. }
39b0: 0a 0a 20 20 2f 2a 20 53 33 20 5b 43 6f 6d 70 75 .. /* S3 [Compu
39c0: 74 65 20 76 32 2e 5d 20 2a 2f 0a 20 20 6d 70 7a te v2.] */. mpz
39d0: 5f 73 75 62 20 28 75 2c 20 75 2c 20 68 29 3b 0a _sub (u, u, h);.
39e0: 20 20 6d 70 7a 5f 73 75 62 20 28 76 2c 20 76 2c mpz_sub (v, v,
39f0: 20 70 29 3b 0a 0a 20 20 6d 70 7a 5f 70 6f 77 5f p);.. mpz_pow_
3a00: 75 69 20 28 74 6d 70 31 2c 20 75 2c 20 32 29 3b ui (tmp1, u, 2);
3a10: 0a 20 20 6d 70 7a 5f 70 6f 77 5f 75 69 20 28 74 . mpz_pow_ui (t
3a20: 6d 70 32 2c 20 76 2c 20 32 29 3b 0a 20 20 6d 70 mp2, v, 2);. mp
3a30: 7a 5f 61 64 64 20 28 74 6d 70 31 2c 20 74 6d 70 z_add (tmp1, tmp
3a40: 31 2c 20 74 6d 70 32 29 3b 0a 20 20 69 66 20 28 1, tmp2);. if (
3a50: 6d 70 7a 5f 63 6d 70 20 28 74 6d 70 31 2c 20 73 mpz_cmp (tmp1, s
3a60: 29 20 3c 20 30 29 0a 20 20 20 20 7b 0a 20 20 20 ) < 0). {.
3a70: 20 20 20 6d 70 7a 5f 73 65 74 20 28 73 2c 20 74 mpz_set (s, t
3a80: 6d 70 31 29 3b 09 2f 2a 20 73 20 3d 20 75 5e 32 mp1);./* s = u^2
3a90: 20 2b 20 76 5e 32 20 2a 2f 0a 20 20 20 20 20 20 + v^2 */.
3aa0: 6d 70 7a 5f 73 65 74 20 28 68 70 2c 20 75 29 3b mpz_set (hp, u);
3ab0: 0a 20 20 20 20 20 20 6d 70 7a 5f 73 65 74 20 28 . mpz_set (
3ac0: 70 70 2c 20 76 29 3b 0a 20 20 20 20 7d 0a 20 20 pp, v);. }.
3ad0: 6d 70 66 5f 73 65 74 5f 7a 20 28 66 5f 74 6d 70 mpf_set_z (f_tmp
3ae0: 31 2c 20 73 29 3b 0a 20 20 6d 70 66 5f 73 71 72 1, s);. mpf_sqr
3af0: 74 20 28 72 6f 70 5b 75 69 5f 74 20 2d 20 31 5d t (rop[ui_t - 1]
3b00: 2c 20 66 5f 74 6d 70 31 29 3b 0a 20 20 20 20 20 , f_tmp1);.
3b10: 20 0a 20 20 2f 2a 20 53 34 20 5b 41 64 76 61 6e . /* S4 [Advan
3b20: 63 65 20 74 2e 5d 20 2a 2f 0a 20 20 6d 70 7a 5f ce t.] */. mpz_
3b30: 6e 65 67 20 28 55 5b 30 5d 5b 30 5d 2c 20 68 29 neg (U[0][0], h)
3b40: 3b 0a 20 20 6d 70 7a 5f 73 65 74 20 28 55 5b 30 ;. mpz_set (U[0
3b50: 5d 5b 31 5d 2c 20 70 29 3b 0a 20 20 6d 70 7a 5f ][1], p);. mpz_
3b60: 6e 65 67 20 28 55 5b 31 5d 5b 30 5d 2c 20 68 70 neg (U[1][0], hp
3b70: 29 3b 0a 20 20 6d 70 7a 5f 73 65 74 20 28 55 5b );. mpz_set (U[
3b80: 31 5d 5b 31 5d 2c 20 70 70 29 3b 0a 20 20 0a 20 1][1], pp);. .
3b90: 20 6d 70 7a 5f 73 65 74 20 28 56 5b 30 5d 5b 30 mpz_set (V[0][0
3ba0: 5d 2c 20 70 70 29 3b 0a 20 20 6d 70 7a 5f 73 65 ], pp);. mpz_se
3bb0: 74 20 28 56 5b 30 5d 5b 31 5d 2c 20 68 70 29 3b t (V[0][1], hp);
3bc0: 0a 20 20 6d 70 7a 5f 6e 65 67 20 28 56 5b 31 5d . mpz_neg (V[1]
3bd0: 5b 30 5d 2c 20 70 29 3b 0a 20 20 6d 70 7a 5f 6e [0], p);. mpz_n
3be0: 65 67 20 28 56 5b 31 5d 5b 31 5d 2c 20 68 29 3b eg (V[1][1], h);
3bf0: 0a 20 20 69 66 20 28 6d 70 7a 5f 63 6d 70 5f 75 . if (mpz_cmp_u
3c00: 69 20 28 70 70 2c 20 30 29 20 3e 20 30 29 0a 20 i (pp, 0) > 0).
3c10: 20 20 20 7b 0a 20 20 20 20 20 20 6d 70 7a 5f 6e {. mpz_n
3c20: 65 67 20 28 56 5b 30 5d 5b 30 5d 2c 20 56 5b 30 eg (V[0][0], V[0
3c30: 5d 5b 30 5d 29 3b 0a 20 20 20 20 20 20 6d 70 7a ][0]);. mpz
3c40: 5f 6e 65 67 20 28 56 5b 30 5d 5b 31 5d 2c 20 56 _neg (V[0][1], V
3c50: 5b 30 5d 5b 31 5d 29 3b 0a 20 20 20 20 20 20 6d [0][1]);. m
3c60: 70 7a 5f 6e 65 67 20 28 56 5b 31 5d 5b 30 5d 2c pz_neg (V[1][0],
3c70: 20 56 5b 31 5d 5b 30 5d 29 3b 0a 20 20 20 20 20 V[1][0]);.
3c80: 20 6d 70 7a 5f 6e 65 67 20 28 56 5b 31 5d 5b 31 mpz_neg (V[1][1
3c90: 5d 2c 20 56 5b 31 5d 5b 31 5d 29 3b 0a 20 20 20 ], V[1][1]);.
3ca0: 20 7d 0a 0a 20 20 77 68 69 6c 65 20 28 75 69 5f }.. while (ui_
3cb0: 74 20 2b 20 31 20 21 3d 20 54 29 09 09 2f 2a 20 t + 1 != T)../*
3cc0: 53 34 20 6c 6f 6f 70 20 2a 2f 0a 20 20 20 20 7b S4 loop */. {
3cd0: 0a 20 20 20 20 20 20 75 69 5f 74 2b 2b 3b 0a 20 . ui_t++;.
3ce0: 20 20 20 20 20 6d 70 7a 5f 6d 75 6c 20 28 72 2c mpz_mul (r,
3cf0: 20 61 2c 20 72 29 3b 0a 20 20 20 20 20 20 6d 70 a, r);. mp
3d00: 7a 5f 6d 6f 64 20 28 72 2c 20 72 2c 20 6d 29 3b z_mod (r, r, m);
3d10: 0a 0a 20 20 20 20 20 20 2f 2a 20 41 64 64 20 6e .. /* Add n
3d20: 65 77 20 72 6f 77 20 61 6e 64 20 63 6f 6c 75 6d ew row and colum
3d30: 6e 20 74 6f 20 55 20 61 6e 64 20 56 2e 20 20 54 n to U and V. T
3d40: 68 65 79 20 61 72 65 20 69 6e 69 74 69 61 6c 69 hey are initiali
3d50: 7a 65 64 20 77 69 74 68 0a 09 20 61 6c 6c 20 65 zed with.. all e
3d60: 6c 65 6d 65 6e 74 73 20 73 65 74 20 74 6f 20 7a lements set to z
3d70: 65 72 6f 2c 20 73 6f 20 63 6c 65 61 72 69 6e 67 ero, so clearing
3d80: 20 69 73 20 6e 6f 74 20 6e 65 63 65 73 73 61 72 is not necessar
3d90: 79 2e 20 2a 2f 0a 0a 20 20 20 20 20 20 6d 70 7a y. */.. mpz
3da0: 5f 6e 65 67 20 28 55 5b 75 69 5f 74 5d 5b 30 5d _neg (U[ui_t][0]
3db0: 2c 20 72 29 3b 20 2f 2a 20 55 3a 20 46 69 72 73 , r); /* U: Firs
3dc0: 74 20 63 6f 6c 20 69 6e 20 6e 65 77 20 72 6f 77 t col in new row
3dd0: 2e 20 2a 2f 0a 20 20 20 20 20 20 6d 70 7a 5f 73 . */. mpz_s
3de0: 65 74 5f 75 69 20 28 55 5b 75 69 5f 74 5d 5b 75 et_ui (U[ui_t][u
3df0: 69 5f 74 5d 2c 20 31 29 3b 20 2f 2a 20 55 3a 20 i_t], 1); /* U:
3e00: 4c 61 73 74 20 63 6f 6c 20 69 6e 20 6e 65 77 20 Last col in new
3e10: 72 6f 77 2e 20 2a 2f 0a 0a 20 20 20 20 20 20 6d row. */.. m
3e20: 70 7a 5f 73 65 74 20 28 56 5b 75 69 5f 74 5d 5b pz_set (V[ui_t][
3e30: 75 69 5f 74 5d 2c 20 6d 29 3b 20 2f 2a 20 56 3a ui_t], m); /* V:
3e40: 20 4c 61 73 74 20 63 6f 6c 20 69 6e 20 6e 65 77 Last col in new
3e50: 20 72 6f 77 2e 20 2a 2f 0a 20 20 20 20 20 20 0a row. */. .
3e60: 20 20 20 20 20 20 2f 2a 20 22 46 69 6e 61 6c 6c /* "Finall
3e70: 79 2c 20 66 6f 72 20 31 20 3c 3d 20 69 20 3c 20 y, for 1 <= i <
3e80: 74 2c 0a 09 20 20 20 73 65 74 20 71 20 3d 20 72 t,.. set q = r
3e90: 6f 75 6e 64 20 28 76 69 31 20 2a 20 72 20 2f 20 ound (vi1 * r /
3ea0: 6d 29 2c 0a 09 20 20 20 76 69 74 20 3d 20 76 69 m),.. vit = vi
3eb0: 31 2a 72 20 2d 20 71 2a 6d 2c 0a 09 20 20 20 61 1*r - q*m,.. a
3ec0: 6e 64 20 55 74 3d 55 74 2b 71 2a 55 69 20 2a 2f nd Ut=Ut+q*Ui */
3ed0: 0a 0a 20 20 20 20 20 20 66 6f 72 20 28 75 69 5f .. for (ui_
3ee0: 69 20 3d 20 30 3b 20 75 69 5f 69 20 3c 20 75 69 i = 0; ui_i < ui
3ef0: 5f 74 3b 20 75 69 5f 69 2b 2b 29 0a 09 7b 0a 09 _t; ui_i++)..{..
3f00: 20 20 6d 70 7a 5f 6d 75 6c 20 28 74 6d 70 31 2c mpz_mul (tmp1,
3f10: 20 56 5b 75 69 5f 69 5d 5b 30 5d 2c 20 72 29 3b V[ui_i][0], r);
3f20: 20 2f 2a 20 74 6d 70 31 3d 76 69 31 2a 72 20 2a /* tmp1=vi1*r *
3f30: 2f 0a 09 20 20 7a 64 69 76 5f 72 6f 75 6e 64 20 /.. zdiv_round
3f40: 28 71 2c 20 74 6d 70 31 2c 20 6d 29 3b 20 2f 2a (q, tmp1, m); /*
3f50: 20 71 3d 72 6f 75 6e 64 28 76 69 31 2a 72 2f 6d q=round(vi1*r/m
3f60: 29 20 2a 2f 0a 09 20 20 6d 70 7a 5f 6d 75 6c 20 ) */.. mpz_mul
3f70: 28 74 6d 70 32 2c 20 71 2c 20 6d 29 3b 09 2f 2a (tmp2, q, m);./*
3f80: 20 74 6d 70 32 3d 71 2a 6d 20 2a 2f 0a 09 20 20 tmp2=q*m */..
3f90: 6d 70 7a 5f 73 75 62 20 28 56 5b 75 69 5f 69 5d mpz_sub (V[ui_i]
3fa0: 5b 75 69 5f 74 5d 2c 20 74 6d 70 31 2c 20 74 6d [ui_t], tmp1, tm
3fb0: 70 32 29 3b 0a 0a 09 20 20 66 6f 72 20 28 75 69 p2);... for (ui
3fc0: 5f 6a 20 3d 20 30 3b 20 75 69 5f 6a 20 3c 3d 20 _j = 0; ui_j <=
3fd0: 75 69 5f 74 3b 20 75 69 5f 6a 2b 2b 29 20 2f 2a ui_t; ui_j++) /*
3fe0: 20 55 5b 74 5d 20 3d 20 55 5b 74 5d 20 2b 20 71 U[t] = U[t] + q
3ff0: 2a 55 5b 69 5d 20 2a 2f 0a 09 20 20 20 20 7b 0a *U[i] */.. {.
4000: 09 20 20 20 20 20 20 6d 70 7a 5f 6d 75 6c 20 28 . mpz_mul (
4010: 74 6d 70 31 2c 20 71 2c 20 55 5b 75 69 5f 69 5d tmp1, q, U[ui_i]
4020: 5b 75 69 5f 6a 5d 29 3b 09 2f 2a 20 74 6d 70 3d [ui_j]);./* tmp=
4030: 71 2a 75 69 6a 20 2a 2f 0a 09 20 20 20 20 20 20 q*uij */..
4040: 6d 70 7a 5f 61 64 64 20 28 55 5b 75 69 5f 74 5d mpz_add (U[ui_t]
4050: 5b 75 69 5f 6a 5d 2c 20 55 5b 75 69 5f 74 5d 5b [ui_j], U[ui_t][
4060: 75 69 5f 6a 5d 2c 20 74 6d 70 31 29 3b 20 2f 2a ui_j], tmp1); /*
4070: 20 75 74 6a 20 3d 20 75 74 6a 20 2b 20 71 2a 75 utj = utj + q*u
4080: 69 6a 20 2a 2f 0a 09 20 20 20 20 7d 0a 09 7d 0a ij */.. }..}.
4090: 0a 20 20 20 20 20 20 2f 2a 20 73 20 3d 20 6d 69 . /* s = mi
40a0: 6e 20 28 73 2c 20 7a 64 6f 74 20 28 55 5b 74 5d n (s, zdot (U[t]
40b0: 2c 20 55 5b 74 5d 29 20 2a 2f 0a 20 20 20 20 20 , U[t]) */.
40c0: 20 76 7a 5f 64 6f 74 20 28 74 6d 70 31 2c 20 55 vz_dot (tmp1, U
40d0: 5b 75 69 5f 74 5d 2c 20 55 5b 75 69 5f 74 5d 2c [ui_t], U[ui_t],
40e0: 20 75 69 5f 74 20 2b 20 31 29 3b 0a 20 20 20 20 ui_t + 1);.
40f0: 20 20 69 66 20 28 6d 70 7a 5f 63 6d 70 20 28 74 if (mpz_cmp (t
4100: 6d 70 31 2c 20 73 29 20 3c 20 30 29 0a 09 6d 70 mp1, s) < 0)..mp
4110: 7a 5f 73 65 74 20 28 73 2c 20 74 6d 70 31 29 3b z_set (s, tmp1);
4120: 0a 0a 20 20 20 20 20 20 75 69 5f 6b 20 3d 20 75 .. ui_k = u
4130: 69 5f 74 3b 0a 20 20 20 20 20 20 75 69 5f 6a 20 i_t;. ui_j
4140: 3d 20 30 3b 09 09 09 2f 2a 20 57 41 52 4e 49 4e = 0;.../* WARNIN
4150: 47 3a 20 75 69 5f 6a 20 6e 6f 20 6c 6f 6e 67 65 G: ui_j no longe
4160: 72 20 61 20 74 65 6d 70 2e 20 2a 2f 0a 0a 20 20 r a temp. */..
4170: 20 20 20 20 2f 2a 20 53 35 20 5b 54 72 61 6e 73 /* S5 [Trans
4180: 66 6f 72 6d 2e 5d 20 2a 2f 0a 20 20 20 20 20 20 form.] */.
4190: 69 66 20 28 67 5f 64 65 62 75 67 20 3e 20 44 45 if (g_debug > DE
41a0: 42 55 47 5f 32 29 0a 09 70 72 69 6e 74 66 20 28 BUG_2)..printf (
41b0: 22 28 74 2c 20 6b 2c 20 6a 2c 20 71 31 2c 20 71 "(t, k, j, q1, q
41c0: 32 2c 20 2e 2e 2e 29 5c 6e 22 29 3b 0a 20 20 20 2, ...)\n");.
41d0: 20 20 20 64 6f 20 0a 09 7b 0a 09 20 20 69 66 20 do ..{.. if
41e0: 28 67 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 (g_debug > DEBUG
41f0: 5f 32 29 0a 09 20 20 20 20 70 72 69 6e 74 66 20 _2).. printf
4200: 28 22 28 25 75 2c 20 25 75 2c 20 25 75 22 2c 20 ("(%u, %u, %u",
4210: 75 69 5f 74 20 2b 20 31 2c 20 75 69 5f 6b 20 2b ui_t + 1, ui_k +
4220: 20 31 2c 20 75 69 5f 6a 20 2b 20 31 29 3b 0a 0a 1, ui_j + 1);..
4230: 09 20 20 66 6f 72 20 28 75 69 5f 69 20 3d 20 30 . for (ui_i = 0
4240: 3b 20 75 69 5f 69 20 3c 3d 20 75 69 5f 74 3b 20 ; ui_i <= ui_t;
4250: 75 69 5f 69 2b 2b 29 0a 09 20 20 20 20 7b 0a 09 ui_i++).. {..
4260: 20 20 20 20 20 20 69 66 20 28 75 69 5f 69 20 21 if (ui_i !
4270: 3d 20 75 69 5f 6a 29 0a 09 09 7b 0a 09 09 20 20 = ui_j)...{...
4280: 76 7a 5f 64 6f 74 20 28 74 6d 70 31 2c 20 56 5b vz_dot (tmp1, V[
4290: 75 69 5f 69 5d 2c 20 56 5b 75 69 5f 6a 5d 2c 20 ui_i], V[ui_j],
42a0: 75 69 5f 74 20 2b 20 31 29 3b 20 2f 2a 20 74 6d ui_t + 1); /* tm
42b0: 70 31 3d 64 6f 74 28 56 69 2c 56 6a 29 2e 20 2a p1=dot(Vi,Vj). *
42c0: 2f 0a 09 09 20 20 6d 70 7a 5f 61 62 73 20 28 74 /... mpz_abs (t
42d0: 6d 70 32 2c 20 74 6d 70 31 29 3b 0a 09 09 20 20 mp2, tmp1);...
42e0: 6d 70 7a 5f 6d 75 6c 5f 75 69 20 28 74 6d 70 32 mpz_mul_ui (tmp2
42f0: 2c 20 74 6d 70 32 2c 20 32 29 3b 20 2f 2a 20 74 , tmp2, 2); /* t
4300: 6d 70 32 20 3d 20 32 2a 61 62 73 28 64 6f 74 28 mp2 = 2*abs(dot(
4310: 56 69 2c 56 6a 29 20 2a 2f 0a 09 09 20 20 76 7a Vi,Vj) */... vz
4320: 5f 64 6f 74 20 28 74 6d 70 33 2c 20 56 5b 75 69 _dot (tmp3, V[ui
4330: 5f 6a 5d 2c 20 56 5b 75 69 5f 6a 5d 2c 20 75 69 _j], V[ui_j], ui
4340: 5f 74 20 2b 20 31 29 3b 20 2f 2a 20 74 6d 70 33 _t + 1); /* tmp3
4350: 3d 64 6f 74 28 56 6a 2c 56 6a 29 2e 20 2a 2f 0a =dot(Vj,Vj). */.
4360: 0a 09 09 20 20 69 66 20 28 6d 70 7a 5f 63 6d 70 ... if (mpz_cmp
4370: 20 28 74 6d 70 32 2c 20 74 6d 70 33 29 20 3e 20 (tmp2, tmp3) >
4380: 30 29 0a 09 09 20 20 20 20 7b 0a 09 09 20 20 20 0)... {...
4390: 20 20 20 7a 64 69 76 5f 72 6f 75 6e 64 20 28 71 zdiv_round (q
43a0: 2c 20 74 6d 70 31 2c 20 74 6d 70 33 29 3b 20 2f , tmp1, tmp3); /
43b0: 2a 20 71 3d 72 6f 75 6e 64 28 56 69 2e 56 6a 2f * q=round(Vi.Vj/
43c0: 56 6a 2e 56 6a 29 20 2a 2f 0a 09 09 20 20 20 20 Vj.Vj) */...
43d0: 20 20 69 66 20 28 67 5f 64 65 62 75 67 20 3e 20 if (g_debug >
43e0: 44 45 42 55 47 5f 32 29 0a 09 09 09 7b 0a 09 09 DEBUG_2)....{...
43f0: 09 20 20 70 72 69 6e 74 66 20 28 22 2c 20 22 29 . printf (", ")
4400: 3b 0a 09 09 09 20 20 6d 70 7a 5f 6f 75 74 5f 73 ;.... mpz_out_s
4410: 74 72 20 28 73 74 64 6f 75 74 2c 20 31 30 2c 20 tr (stdout, 10,
4420: 71 29 3b 0a 09 09 09 7d 0a 0a 09 09 20 20 20 20 q);....}....
4430: 20 20 66 6f 72 20 28 75 69 5f 6c 20 3d 20 30 3b for (ui_l = 0;
4440: 20 75 69 5f 6c 20 3c 3d 20 75 69 5f 74 3b 20 75 ui_l <= ui_t; u
4450: 69 5f 6c 2b 2b 29 0a 09 09 09 7b 0a 09 09 09 20 i_l++)....{....
4460: 20 6d 70 7a 5f 6d 75 6c 20 28 74 6d 70 31 2c 20 mpz_mul (tmp1,
4470: 71 2c 20 56 5b 75 69 5f 6a 5d 5b 75 69 5f 6c 5d q, V[ui_j][ui_l]
4480: 29 3b 0a 09 09 09 20 20 6d 70 7a 5f 73 75 62 20 );.... mpz_sub
4490: 28 56 5b 75 69 5f 69 5d 5b 75 69 5f 6c 5d 2c 20 (V[ui_i][ui_l],
44a0: 56 5b 75 69 5f 69 5d 5b 75 69 5f 6c 5d 2c 20 74 V[ui_i][ui_l], t
44b0: 6d 70 31 29 3b 20 2f 2a 20 56 69 3d 56 69 2d 71 mp1); /* Vi=Vi-q
44c0: 2a 56 6a 20 2a 2f 0a 09 09 09 20 20 6d 70 7a 5f *Vj */.... mpz_
44d0: 6d 75 6c 20 28 74 6d 70 31 2c 20 71 2c 20 55 5b mul (tmp1, q, U[
44e0: 75 69 5f 69 5d 5b 75 69 5f 6c 5d 29 3b 0a 09 09 ui_i][ui_l]);...
44f0: 09 20 20 6d 70 7a 5f 61 64 64 20 28 55 5b 75 69 . mpz_add (U[ui
4500: 5f 6a 5d 5b 75 69 5f 6c 5d 2c 20 55 5b 75 69 5f _j][ui_l], U[ui_
4510: 6a 5d 5b 75 69 5f 6c 5d 2c 20 74 6d 70 31 29 3b j][ui_l], tmp1);
4520: 20 2f 2a 20 55 6a 3d 55 6a 2b 71 2a 55 69 20 2a /* Uj=Uj+q*Ui *
4530: 2f 0a 09 09 09 7d 0a 09 09 20 20 20 20 20 20 0a /....}... .
4540: 09 09 20 20 20 20 20 20 76 7a 5f 64 6f 74 20 28 .. vz_dot (
4550: 74 6d 70 31 2c 20 55 5b 75 69 5f 6a 5d 2c 20 55 tmp1, U[ui_j], U
4560: 5b 75 69 5f 6a 5d 2c 20 75 69 5f 74 20 2b 20 31 [ui_j], ui_t + 1
4570: 29 3b 20 2f 2a 20 74 6d 70 31 3d 64 6f 74 28 55 ); /* tmp1=dot(U
4580: 6a 2c 55 6a 29 20 2a 2f 0a 09 09 20 20 20 20 20 j,Uj) */...
4590: 20 69 66 20 28 6d 70 7a 5f 63 6d 70 20 28 74 6d if (mpz_cmp (tm
45a0: 70 31 2c 20 73 29 20 3c 20 30 29 20 2f 2a 20 73 p1, s) < 0) /* s
45b0: 20 3d 20 6d 69 6e 28 73 2c 64 6f 74 28 55 6a 2c = min(s,dot(Uj,
45c0: 55 6a 29 29 20 2a 2f 0a 09 09 09 6d 70 7a 5f 73 Uj)) */....mpz_s
45d0: 65 74 20 28 73 2c 20 74 6d 70 31 29 3b 0a 09 09 et (s, tmp1);...
45e0: 20 20 20 20 20 20 75 69 5f 6b 20 3d 20 75 69 5f ui_k = ui_
45f0: 6a 3b 0a 09 09 20 20 20 20 7d 0a 09 09 20 20 65 j;... }... e
4600: 6c 73 65 20 69 66 20 28 67 5f 64 65 62 75 67 20 lse if (g_debug
4610: 3e 20 44 45 42 55 47 5f 32 29 0a 09 09 20 20 20 > DEBUG_2)...
4620: 20 70 72 69 6e 74 66 20 28 22 2c 20 23 22 29 3b printf (", #");
4630: 20 2f 2a 20 32 7c 56 69 2e 56 6a 7c 20 3c 3d 20 /* 2|Vi.Vj| <=
4640: 56 6a 2e 56 6a 20 2a 2f 0a 09 09 7d 0a 09 20 20 Vj.Vj */...}..
4650: 20 20 20 20 65 6c 73 65 20 69 66 20 28 67 5f 64 else if (g_d
4660: 65 62 75 67 20 3e 20 44 45 42 55 47 5f 32 29 0a ebug > DEBUG_2).
4670: 09 09 70 72 69 6e 74 66 20 28 22 2c 20 2a 22 29 ..printf (", *")
4680: 3b 09 2f 2a 20 69 20 3d 3d 20 6a 20 2a 2f 0a 09 ;./* i == j */..
4690: 20 20 20 20 7d 0a 0a 09 20 20 69 66 20 28 67 5f }... if (g_
46a0: 64 65 62 75 67 20 3e 20 44 45 42 55 47 5f 32 29 debug > DEBUG_2)
46b0: 0a 09 20 20 20 20 70 72 69 6e 74 66 20 28 22 29 .. printf (")
46c0: 5c 6e 22 29 3b 0a 0a 09 20 20 2f 2a 20 53 36 20 \n");... /* S6
46d0: 5b 41 64 76 61 6e 63 65 20 6a 2e 5d 20 2a 2f 0a [Advance j.] */.
46e0: 09 20 20 69 66 20 28 75 69 5f 6a 20 3d 3d 20 75 . if (ui_j == u
46f0: 69 5f 74 29 0a 09 20 20 20 20 75 69 5f 6a 20 3d i_t).. ui_j =
4700: 20 30 3b 0a 09 20 20 65 6c 73 65 0a 09 20 20 20 0;.. else..
4710: 20 75 69 5f 6a 2b 2b 3b 0a 09 7d 0a 20 20 20 20 ui_j++;..}.
4720: 20 20 77 68 69 6c 65 20 28 75 69 5f 6a 20 21 3d while (ui_j !=
4730: 20 75 69 5f 6b 29 3b 09 2f 2a 20 53 35 20 2a 2f ui_k);./* S5 */
4740: 0a 0a 20 20 20 20 20 20 2f 2a 20 46 72 6f 6d 20 .. /* From
4750: 4b 6e 75 74 68 20 70 2e 20 31 30 34 3a 20 22 54 Knuth p. 104: "T
4760: 68 65 20 65 78 68 61 75 73 74 69 76 65 20 73 65 he exhaustive se
4770: 61 72 63 68 20 69 6e 20 73 74 65 70 73 20 53 38 arch in steps S8
4780: 2d 53 31 30 0a 09 20 72 65 64 75 63 65 73 20 74 -S10.. reduces t
4790: 68 65 20 76 61 6c 75 65 20 6f 66 20 73 20 6f 6e he value of s on
47a0: 6c 79 20 72 61 72 65 6c 79 2e 22 20 2a 2f 0a 23 ly rarely." */.#
47b0: 69 66 64 65 66 20 44 4f 5f 53 45 41 52 43 48 0a ifdef DO_SEARCH.
47c0: 20 20 20 20 20 20 2f 2a 20 53 37 20 5b 50 72 65 /* S7 [Pre
47d0: 70 61 72 65 20 66 6f 72 20 73 65 61 72 63 68 2e pare for search.
47e0: 5d 20 2a 2f 0a 20 20 20 20 20 20 2f 2a 20 46 69 ] */. /* Fi
47f0: 6e 64 20 6d 69 6e 69 6d 75 6d 20 69 6e 20 28 78 nd minimum in (x
4800: 5b 31 5d 2c 20 2e 2e 2e 2c 20 78 5b 74 5d 29 20 [1], ..., x[t])
4810: 73 61 74 69 73 66 79 69 6e 67 20 63 6f 6e 64 69 satisfying condi
4820: 74 69 6f 6e 0a 09 20 78 5b 6b 5d 5e 32 20 3c 3d tion.. x[k]^2 <=
4830: 20 66 28 79 5b 31 5d 2c 20 2e 2e 2e 2c 79 5b 74 f(y[1], ...,y[t
4840: 5d 29 20 2a 20 64 6f 74 28 56 5b 6b 5d 2c 56 5b ]) * dot(V[k],V[
4850: 6b 5d 29 20 2a 2f 0a 0a 20 20 20 20 20 20 75 69 k]) */.. ui
4860: 5f 6b 20 3d 20 75 69 5f 74 3b 0a 20 20 20 20 20 _k = ui_t;.
4870: 20 69 66 20 28 67 5f 64 65 62 75 67 20 3e 20 44 if (g_debug > D
4880: 45 42 55 47 5f 32 29 0a 09 7b 0a 09 20 20 70 72 EBUG_2)..{.. pr
4890: 69 6e 74 66 20 28 22 73 65 61 72 63 68 69 6e 67 intf ("searching
48a0: 2e 2e 2e 22 29 3b 0a 09 20 20 2f 2a 66 6f 72 20 ...");.. /*for
48b0: 28 66 20 3d 20 30 3b 20 66 20 3c 20 75 69 5f 74 (f = 0; f < ui_t
48c0: 2a 2f 0a 09 20 20 66 66 6c 75 73 68 20 28 73 74 */.. fflush (st
48d0: 64 6f 75 74 29 3b 0a 09 7d 0a 0a 20 20 20 20 20 dout);..}..
48e0: 20 2f 2a 20 5a 5b 69 5d 20 3d 20 66 6c 6f 6f 72 /* Z[i] = floor
48f0: 20 28 73 71 72 74 20 28 66 6c 6f 6f 72 20 28 64 (sqrt (floor (d
4900: 6f 74 28 56 5b 69 5d 2c 56 5b 69 5d 29 20 2a 20 ot(V[i],V[i]) *
4910: 73 20 2f 20 6d 5e 32 29 29 29 3b 20 2a 2f 0a 20 s / m^2))); */.
4920: 20 20 20 20 20 6d 70 7a 5f 70 6f 77 5f 75 69 20 mpz_pow_ui
4930: 28 74 6d 70 31 2c 20 6d 2c 20 32 29 3b 0a 20 20 (tmp1, m, 2);.
4940: 20 20 20 20 6d 70 66 5f 73 65 74 5f 7a 20 28 66 mpf_set_z (f
4950: 5f 74 6d 70 31 2c 20 74 6d 70 31 29 3b 0a 20 20 _tmp1, tmp1);.
4960: 20 20 20 20 6d 70 66 5f 73 65 74 5f 7a 20 28 66 mpf_set_z (f
4970: 5f 74 6d 70 32 2c 20 73 29 3b 0a 20 20 20 20 20 _tmp2, s);.
4980: 20 6d 70 66 5f 64 69 76 20 28 66 5f 74 6d 70 31 mpf_div (f_tmp1
4990: 2c 20 66 5f 74 6d 70 32 2c 20 66 5f 74 6d 70 31 , f_tmp2, f_tmp1
49a0: 29 3b 09 2f 2a 20 66 5f 74 6d 70 31 20 3d 20 73 );./* f_tmp1 = s
49b0: 2f 6d 5e 32 20 2a 2f 0a 20 20 20 20 20 20 66 6f /m^2 */. fo
49c0: 72 20 28 75 69 5f 69 20 3d 20 30 3b 20 75 69 5f r (ui_i = 0; ui_
49d0: 69 20 3c 3d 20 75 69 5f 74 3b 20 75 69 5f 69 2b i <= ui_t; ui_i+
49e0: 2b 29 0a 09 7b 0a 09 20 20 76 7a 5f 64 6f 74 20 +)..{.. vz_dot
49f0: 28 74 6d 70 31 2c 20 56 5b 75 69 5f 69 5d 2c 20 (tmp1, V[ui_i],
4a00: 56 5b 75 69 5f 69 5d 2c 20 75 69 5f 74 20 2b 20 V[ui_i], ui_t +
4a10: 31 29 3b 0a 09 20 20 6d 70 66 5f 73 65 74 5f 7a 1);.. mpf_set_z
4a20: 20 28 66 5f 74 6d 70 32 2c 20 74 6d 70 31 29 3b (f_tmp2, tmp1);
4a30: 0a 09 20 20 6d 70 66 5f 6d 75 6c 20 28 66 5f 74 .. mpf_mul (f_t
4a40: 6d 70 32 2c 20 66 5f 74 6d 70 32 2c 20 66 5f 74 mp2, f_tmp2, f_t
4a50: 6d 70 31 29 3b 0a 09 20 20 66 5f 66 6c 6f 6f 72 mp1);.. f_floor
4a60: 20 28 66 5f 74 6d 70 32 2c 20 66 5f 74 6d 70 32 (f_tmp2, f_tmp2
4a70: 29 3b 0a 09 20 20 6d 70 66 5f 73 71 72 74 20 28 );.. mpf_sqrt (
4a80: 66 5f 74 6d 70 32 2c 20 66 5f 74 6d 70 32 29 3b f_tmp2, f_tmp2);
4a90: 0a 09 20 20 6d 70 7a 5f 73 65 74 5f 66 20 28 5a .. mpz_set_f (Z
4aa0: 5b 75 69 5f 69 5d 2c 20 66 5f 74 6d 70 32 29 3b [ui_i], f_tmp2);
4ab0: 0a 09 7d 0a 0a 20 20 20 20 20 20 2f 2a 20 53 38 ..}.. /* S8
4ac0: 20 5b 41 64 76 61 6e 63 65 20 58 5b 6b 5d 2e 5d [Advance X[k].]
4ad0: 20 2a 2f 0a 20 20 20 20 20 20 64 6f 20 0a 09 7b */. do ..{
4ae0: 0a 09 20 20 69 66 20 28 67 5f 64 65 62 75 67 20 .. if (g_debug
4af0: 3e 20 44 45 42 55 47 5f 32 29 0a 09 20 20 20 20 > DEBUG_2)..
4b00: 7b 0a 09 20 20 20 20 20 20 70 72 69 6e 74 66 20 {.. printf
4b10: 28 22 58 5b 25 75 5d 20 3d 20 22 2c 20 75 69 5f ("X[%u] = ", ui_
4b20: 6b 29 3b 0a 09 20 20 20 20 20 20 6d 70 7a 5f 6f k);.. mpz_o
4b30: 75 74 5f 73 74 72 20 28 73 74 64 6f 75 74 2c 20 ut_str (stdout,
4b40: 31 30 2c 20 58 5b 75 69 5f 6b 5d 29 3b 0a 09 20 10, X[ui_k]);..
4b50: 20 20 20 20 20 70 72 69 6e 74 66 20 28 22 5c 74 printf ("\t
4b60: 5a 5b 25 75 5d 20 3d 20 22 2c 20 75 69 5f 6b 29 Z[%u] = ", ui_k)
4b70: 3b 0a 09 20 20 20 20 20 20 6d 70 7a 5f 6f 75 74 ;.. mpz_out
4b80: 5f 73 74 72 20 28 73 74 64 6f 75 74 2c 20 31 30 _str (stdout, 10
4b90: 2c 20 5a 5b 75 69 5f 6b 5d 29 3b 0a 09 20 20 20 , Z[ui_k]);..
4ba0: 20 20 20 70 72 69 6e 74 66 20 28 22 5c 6e 22 29 printf ("\n")
4bb0: 3b 0a 09 20 20 20 20 20 20 66 66 6c 75 73 68 20 ;.. fflush
4bc0: 28 73 74 64 6f 75 74 29 3b 0a 09 20 20 20 20 7d (stdout);.. }
4bd0: 0a 09 20 20 20 20 20 20 0a 09 20 20 69 66 20 28 .. .. if (
4be0: 6d 70 7a 5f 63 6d 70 20 28 58 5b 75 69 5f 6b 5d mpz_cmp (X[ui_k]
4bf0: 2c 20 5a 5b 75 69 5f 6b 5d 29 29 0a 09 20 20 20 , Z[ui_k]))..
4c00: 20 7b 0a 09 20 20 20 20 20 20 6d 70 7a 5f 61 64 {.. mpz_ad
4c10: 64 5f 75 69 20 28 58 5b 75 69 5f 6b 5d 2c 20 58 d_ui (X[ui_k], X
4c20: 5b 75 69 5f 6b 5d 2c 20 31 29 3b 0a 09 20 20 20 [ui_k], 1);..
4c30: 20 20 20 66 6f 72 20 28 75 69 5f 69 20 3d 20 30 for (ui_i = 0
4c40: 3b 20 75 69 5f 69 20 3c 3d 20 75 69 5f 74 3b 20 ; ui_i <= ui_t;
4c50: 75 69 5f 69 2b 2b 29 0a 09 09 6d 70 7a 5f 61 64 ui_i++)...mpz_ad
4c60: 64 20 28 59 5b 75 69 5f 69 5d 2c 20 59 5b 75 69 d (Y[ui_i], Y[ui
4c70: 5f 69 5d 2c 20 55 5b 75 69 5f 6b 5d 5b 75 69 5f _i], U[ui_k][ui_
4c80: 69 5d 29 3b 0a 0a 09 20 20 20 20 20 20 2f 2a 20 i]);... /*
4c90: 53 39 20 5b 41 64 76 61 6e 63 65 20 6b 2e 5d 20 S9 [Advance k.]
4ca0: 2a 2f 0a 09 20 20 20 20 20 20 77 68 69 6c 65 20 */.. while
4cb0: 28 2b 2b 75 69 5f 6b 20 3c 3d 20 75 69 5f 74 29 (++ui_k <= ui_t)
4cc0: 0a 09 09 7b 0a 09 09 20 20 6d 70 7a 5f 6e 65 67 ...{... mpz_neg
4cd0: 20 28 58 5b 75 69 5f 6b 5d 2c 20 5a 5b 75 69 5f (X[ui_k], Z[ui_
4ce0: 6b 5d 29 3b 0a 09 09 20 20 6d 70 7a 5f 6d 75 6c k]);... mpz_mul
4cf0: 5f 75 69 20 28 74 6d 70 31 2c 20 5a 5b 75 69 5f _ui (tmp1, Z[ui_
4d00: 6b 5d 2c 20 32 29 3b 0a 09 09 20 20 66 6f 72 20 k], 2);... for
4d10: 28 75 69 5f 69 20 3d 20 30 3b 20 75 69 5f 69 20 (ui_i = 0; ui_i
4d20: 3c 3d 20 75 69 5f 74 3b 20 75 69 5f 69 2b 2b 29 <= ui_t; ui_i++)
4d30: 0a 09 09 20 20 20 20 7b 0a 09 09 20 20 20 20 20 ... {...
4d40: 20 6d 70 7a 5f 6d 75 6c 20 28 74 6d 70 32 2c 20 mpz_mul (tmp2,
4d50: 74 6d 70 31 2c 20 55 5b 75 69 5f 6b 5d 5b 75 69 tmp1, U[ui_k][ui
4d60: 5f 69 5d 29 3b 0a 09 09 20 20 20 20 20 20 6d 70 _i]);... mp
4d70: 7a 5f 73 75 62 20 28 59 5b 75 69 5f 69 5d 2c 20 z_sub (Y[ui_i],
4d80: 59 5b 75 69 5f 69 5d 2c 20 74 6d 70 32 29 3b 0a Y[ui_i], tmp2);.
4d90: 09 09 20 20 20 20 7d 0a 09 09 7d 0a 09 20 20 20 .. }...}..
4da0: 20 20 20 76 7a 5f 64 6f 74 20 28 74 6d 70 31 2c vz_dot (tmp1,
4db0: 20 59 2c 20 59 2c 20 75 69 5f 74 20 2b 20 31 29 Y, Y, ui_t + 1)
4dc0: 3b 0a 09 20 20 20 20 20 20 69 66 20 28 6d 70 7a ;.. if (mpz
4dd0: 5f 63 6d 70 20 28 74 6d 70 31 2c 20 73 29 20 3c _cmp (tmp1, s) <
4de0: 20 30 29 0a 09 09 6d 70 7a 5f 73 65 74 20 28 73 0)...mpz_set (s
4df0: 2c 20 74 6d 70 31 29 3b 0a 09 20 20 20 20 7d 0a , tmp1);.. }.
4e00: 09 7d 0a 20 20 20 20 20 20 77 68 69 6c 65 20 28 .}. while (
4e10: 2d 2d 75 69 5f 6b 29 3b 0a 23 65 6e 64 69 66 20 --ui_k);.#endif
4e20: 2f 2a 20 44 4f 5f 53 45 41 52 43 48 20 2a 2f 0a /* DO_SEARCH */.
4e30: 20 20 20 20 20 20 6d 70 66 5f 73 65 74 5f 7a 20 mpf_set_z
4e40: 28 66 5f 74 6d 70 31 2c 20 73 29 3b 0a 20 20 20 (f_tmp1, s);.
4e50: 20 20 20 6d 70 66 5f 73 71 72 74 20 28 72 6f 70 mpf_sqrt (rop
4e60: 5b 75 69 5f 74 20 2d 20 31 5d 2c 20 66 5f 74 6d [ui_t - 1], f_tm
4e70: 70 31 29 3b 0a 23 69 66 64 65 66 20 44 4f 5f 53 p1);.#ifdef DO_S
4e80: 45 41 52 43 48 0a 20 20 20 20 20 20 69 66 20 28 EARCH. if (
4e90: 67 5f 64 65 62 75 67 20 3e 20 44 45 42 55 47 5f g_debug > DEBUG_
4ea0: 32 29 0a 09 70 72 69 6e 74 66 20 28 22 64 6f 6e 2)..printf ("don
4eb0: 65 2e 5c 6e 22 29 3b 0a 23 65 6e 64 69 66 20 2f e.\n");.#endif /
4ec0: 2a 20 44 4f 5f 53 45 41 52 43 48 20 2a 2f 0a 20 * DO_SEARCH */.
4ed0: 20 20 20 7d 20 2f 2a 20 53 34 20 6c 6f 6f 70 20 } /* S4 loop
4ee0: 2a 2f 0a 0a 20 20 2f 2a 20 43 6c 65 61 72 20 47 */.. /* Clear G
4ef0: 4d 50 20 76 61 72 69 61 62 6c 65 73 2e 20 2a 2f MP variables. */
4f00: 0a 0a 20 20 6d 70 66 5f 63 6c 65 61 72 20 28 66 .. mpf_clear (f
4f10: 5f 74 6d 70 31 29 3b 0a 20 20 6d 70 66 5f 63 6c _tmp1);. mpf_cl
4f20: 65 61 72 20 28 66 5f 74 6d 70 32 29 3b 0a 20 20 ear (f_tmp2);.
4f30: 66 6f 72 20 28 75 69 5f 69 20 3d 20 30 3b 20 75 for (ui_i = 0; u
4f40: 69 5f 69 20 3c 20 47 4d 50 5f 53 50 45 43 54 5f i_i < GMP_SPECT_
4f50: 4d 41 58 54 3b 20 75 69 5f 69 2b 2b 29 0a 20 20 MAXT; ui_i++).
4f60: 20 20 7b 0a 20 20 20 20 20 20 66 6f 72 20 28 75 {. for (u
4f70: 69 5f 6a 20 3d 20 30 3b 20 75 69 5f 6a 20 3c 20 i_j = 0; ui_j <
4f80: 47 4d 50 5f 53 50 45 43 54 5f 4d 41 58 54 3b 20 GMP_SPECT_MAXT;
4f90: 75 69 5f 6a 2b 2b 29 0a 09 7b 0a 09 20 20 6d 70 ui_j++)..{.. mp
4fa0: 7a 5f 63 6c 65 61 72 20 28 55 5b 75 69 5f 69 5d z_clear (U[ui_i]
4fb0: 5b 75 69 5f 6a 5d 29 3b 0a 09 20 20 6d 70 7a 5f [ui_j]);.. mpz_
4fc0: 63 6c 65 61 72 20 28 56 5b 75 69 5f 69 5d 5b 75 clear (V[ui_i][u
4fd0: 69 5f 6a 5d 29 3b 0a 09 7d 0a 20 20 20 20 20 20 i_j]);..}.
4fe0: 6d 70 7a 5f 63 6c 65 61 72 20 28 58 5b 75 69 5f mpz_clear (X[ui_
4ff0: 69 5d 29 3b 0a 20 20 20 20 20 20 6d 70 7a 5f 63 i]);. mpz_c
5000: 6c 65 61 72 20 28 59 5b 75 69 5f 69 5d 29 3b 0a lear (Y[ui_i]);.
5010: 20 20 20 20 20 20 6d 70 7a 5f 63 6c 65 61 72 20 mpz_clear
5020: 28 5a 5b 75 69 5f 69 5d 29 3b 0a 20 20 20 20 7d (Z[ui_i]);. }
5030: 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 20 28 74 6d . mpz_clear (tm
5040: 70 31 29 3b 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 p1);. mpz_clear
5050: 20 28 74 6d 70 32 29 3b 0a 20 20 6d 70 7a 5f 63 (tmp2);. mpz_c
5060: 6c 65 61 72 20 28 74 6d 70 33 29 3b 0a 20 20 6d lear (tmp3);. m
5070: 70 7a 5f 63 6c 65 61 72 20 28 68 29 3b 0a 20 20 pz_clear (h);.
5080: 6d 70 7a 5f 63 6c 65 61 72 20 28 68 70 29 3b 0a mpz_clear (hp);.
5090: 20 20 6d 70 7a 5f 63 6c 65 61 72 20 28 72 29 3b mpz_clear (r);
50a0: 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 20 28 73 29 . mpz_clear (s)
50b0: 3b 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 20 28 70 ;. mpz_clear (p
50c0: 29 3b 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 20 28 );. mpz_clear (
50d0: 70 70 29 3b 0a 20 20 6d 70 7a 5f 63 6c 65 61 72 pp);. mpz_clear
50e0: 20 28 71 29 3b 0a 20 20 6d 70 7a 5f 63 6c 65 61 (q);. mpz_clea
50f0: 72 20 28 75 29 3b 0a 20 20 6d 70 7a 5f 63 6c 65 r (u);. mpz_cle
5100: 61 72 20 28 76 29 3b 0a 0a 20 20 72 65 74 75 72 ar (v);.. retur
5110: 6e 3b 0a 7d 0a 0a n;.}..