1/*============================================================================
2
3    This file is part of Antic.
4
5    Antic is free software: you can redistribute it and/or modify it under
6    the terms of the GNU Lesser General Public License (LGPL) as published
7    by the Free Software Foundation; either version 2.1 of the License, or
8    (at your option) any later version. See <http://www.gnu.org/licenses/>.
9
10===============================================================================*/
11
12/*
13Profiling MAGMA polynomial multiplication in Z[x].
14
15Usage: run magma with the -b flag to prevent the start up banner, i.e.
16
17   magma -b magma-profile.m > output.prof
18
19(C) 2007 David Harvey + Bill Hart, GPL
20
21*/
22
23target_name := "NFElemTrace";
24target_description := "MAGMA number field trace over various degree";
25
26BITS := 10;        // bits per coefficient
27ratio := 1.1;      // ratio between consecutive lengths
28monic := false;    // monic integral case or generic rational case
29
30R<x>:=PolynomialRing(Rationals());
31S<y>:=PolynomialRing(Integers());
32
33// Timing runs need to last at least this many microseconds to be counted:
34DURATION_THRESHOLD := 200000;
35// Microseconds per timing run that the prof2d_sample function aims for:
36DURATION_TARGET := 300000;
37
38
39forward prof2d_sample;
40
41/*
42This function should run count iterations at position (x, y),
43and return the total time in seconds, using the Cputime() function.
44*/
45function sampler(length, bits, count)
46
47    // first time random element generation + multiplication
48
49    countmod := 100;
50    if length ge 50 then
51       countmod := 10;
52    end if;
53    if length ge 500 then
54       countmod := 4;
55    end if;
56    time1 := Cputime();
57    for i := 1 to count do
58      if (i-1) mod countmod eq 0 then
59         p:=Polynomial([(-1)^Random(1)*RandomBits(bits): x in [1..length-1]]);
60         if monic then
61            p1:=S!p + x^(length-1);
62         else
63            den:=RandomBits(bits);
64            if den eq 0 then
65               den:=1;
66            end if;
67            p1:=((R!p)+(-1)^Random(1)*RandomBits(bits)*x^(length-1))/den;
68         end if;
69         if Degree(p1) eq 0 then
70            if monic then
71               p1:=x+1;
72            else
73               p1:=y+1;
74            end if;
75         end if;
76         K:=NumberField(p1:Check:=false);
77         a:=Random(K,2^bits);
78         if monic then
79            a:=Numerator(a);
80         end if;
81      end if;
82      c:=Trace(a);
83    end for;
84
85    time2 := Cputime();
86
87    // now time just the random element generation
88
89    for i := 1 to count do
90      if (i-1) mod countmod eq 0 then
91         p:=Polynomial([(-1)^Random(1)*RandomBits(bits): x in [1..length-1]]);
92         if monic then
93            p1:=S!p + x^(length-1);
94         else
95            den:=RandomBits(bits);
96            if den eq 0 then
97               den:=1;
98            end if;
99            p1:=((R!p)+(-1)^Random(1)*RandomBits(bits)*x^(length-1))/den;
100         end if;
101         if Degree(p1) eq 0 then
102            if monic then
103               p1:=x+1;
104            else
105               p1:=y+1;
106            end if;
107         end if;
108         K:=NumberField(p1:Check:=false);
109         a:=Random(K,2^bits);
110         if monic then
111            a:=Numerator(a);
112         end if;
113      end if;
114      c:=a;
115    end for;
116
117    time3 := Cputime();
118    return (time2 - time1) - (time3 - time2);
119end function;
120
121
122/*
123This function should loop over appropriate combinations of (x, y),
124and call prof2d_sample(x, y) for each one.
125*/
126procedure driver()
127   length:=2;
128   while length lt 10000 do
129      monic:=0;
130      prof2d_sample(length, BITS);
131      monic:=1;
132      prof2d_sample(length, BITS);
133      length := Ceiling(length*ratio);
134   end while;
135end procedure;
136
137
138
139/************************************************************************
140
141 This last section is the generic profiling code. Just leave this
142 stuff alone.
143
144************************************************************************/
145
146/*
147Formats in scientific notation with 3 decimal places
148*/
149function format_sci(x)
150    L := Floor(Log(10, x));
151    x := x / 10^L;
152    s := Sprintf("%.3oe", x);
153    if L lt 0 then
154      s := s cat "-";
155    else
156      s := s cat "+";
157    end if;
158
159    s := s cat Sprintf("%o", Floor(Abs(L / 10)));
160    s := s cat Sprintf("%o", (Abs(L) mod 10));
161
162    return s;
163end function;
164
165
166procedure prof2d_sample(x, y)
167    // number of timings that were at least DURATION_THRESHOLD microseconds:
168    good_count := 0;
169
170    // first try just a few loops
171    num_trials := 4;
172    last_time := sampler(x, y, num_trials) * 1000000.0;
173
174    max_time := 0;
175    min_time := 0;
176
177    // loop until we have enough good times
178    while true do
179      per_trial := last_time / num_trials;
180
181      // if the last recorded time was long enough, record it
182      if last_time gt DURATION_THRESHOLD then
183          if good_count gt 0 then
184            max_time := Max(max_time, per_trial);
185            min_time := Min(min_time, per_trial);
186          else
187            max_time := per_trial;
188            min_time := per_trial;
189          end if;
190
191          good_count := good_count + 1;
192          if good_count eq 5 then
193            // we've got enough data
194            // print it out and return
195            print Sprintf("%o\t%o\t%o\t%o", x, y, format_sci(min_time), format_sci(max_time));
196            return;
197          end if;
198      end if;
199
200      // adjust num_trials so that the elapsed time gravitates towards
201      // DURATION_TARGET; num_trials can be changed by a factor of
202      // at most 25%, and must be at least 1
203      if last_time lt 0.0001 then
204          last_time := 0.0001;
205      end if;
206      adjust_ratio := 1.0 * DURATION_TARGET / last_time;
207      if adjust_ratio gt 1.25 then
208          adjust_ratio := 1.25;
209      end if;
210      if adjust_ratio lt 0.75 then
211          adjust_ratio := 0.75;
212      end if;
213      num_trials := Ceiling(adjust_ratio * num_trials);
214      // just to be safe:
215      if num_trials eq 0 then
216          num_trials := 1;
217      end if;
218
219      // run another trial
220      last_time := sampler(x, y, num_trials) * 1000000.0;
221    end while;
222end procedure;
223
224
225procedure print_header()
226    print "FLINT profile output";
227    print "";
228    print "TIMESTAMP: (todo: write code to generate timestamp)";
229    print "MACHINE: (todo: write code to get machine from environment var)";
230
231    print "";
232    print "MODULE: magma";
233    print "TARGET:", target_name;
234    print "";
235    print "DESCRIPTION:";
236    print target_description;
237    print "";
238    print "============================================== begin data";
239
240end procedure;
241
242
243print_header();
244driver();
245
246quit;
247
248// ------------- end of file ------------------------------------
249