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