Hex Artifact Content
Not logged in

Artifact ce38d12c82e26afe630909457bfa6af8a5fb3ca2:


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;.}..