@= f:=1; repeat be_careful:=p-q; p:=be_careful+p; if p>=0 then f:=f+f+1 else begin double(f); p:=p+q; end; until f>=mpfract_one; be_careful:=p-q; if be_careful+p>=0 then incr(f) @ The dual of |make_mpfract| is |take_mpfract|, which multiplies a given integer~|q| by a fraction~|f|. When the operands are positive, it computes $p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function of |q| and~|f|. @ = function take_mpfract(@!q:integer;@!f:mpfract):integer; var @!p:integer; {the fraction so far} @!negative:boolean; {should the result be negated?} @!n:integer; {additional multiple of $q$} @!be_careful:integer; {disables certain compiler optimizations} begin @ =0| and |q>0|@>; if f ; be_careful:=n-el_gordo; if be_careful+p>0 then begin arith_error:=true; n:=el_gordo-p; end; if negative then take_mpfract:=-(n+p) else take_mpfract:=n+p; end; @ @ =0| and |q>0|@>= if f>=0 then negative:=false else begin negate(f); negative:=true; end; if q<0 then begin negate(q); negative:=not negative; end; @ The invariant relations in this case are (i)~$\lfloor(qf+p)/2^k\rfloor =\lfloor qf_0/2^{28}+{1\over2}\rfloor$, where $k$ is an integer and $f_0$ is the original value of~$f$; (ii)~$2^k\L f<2^{k+1}$. @^inner loop@> @ = p:=mpfract_half; {that's $2^{27}$; the invariants hold now with $k=28$} if q = @!randoms:array[0..54] of mpfract; {the last 55 random values generated} @!j_random:0..54; {the number of unused |randoms|} @ This array of pseudo-random numbers is set starting from a seed value, that is kept in the global integer |random_seed|. @ = @!random_seed:integer; {seed for pseudo-random number generation} @ @ = primitive("randomseed",last_item,random_seed_code);@/ @!@:randomseed_}{\.{\\randomseed} primitive@> @ @ = random_seed_code: print_esc("randomseed"); @ @ = random_seed_code: cur_val:=random_seed; @ We set the initial value from the system time. System integrators could provide a better source of pseudo-randomness. Every time a new seed value is assigned, the array has to be regenerated for consumption by routines explained a little later. @<\Prote\ initializations@>= random_seed:=sys_time; init_randoms; @ Since changing the value must trigger the redefinition of the array, a dedicated primitive is defined to take the new seed and call |init_randoms|. @ = primitive("setrandomseed",convert,set_random_seed_code);@/ @!@:setrandomseed_}{\.{\\setrandomseed} primitive@> @ @ = set_random_seed_code: print_esc("setrandomseed"); @ Once we have retrieved and redefined |random_seed|, we must regenerate the |randoms| array. @ = set_random_seed_code: begin scan_int; random_seed:=cur_val; init_randoms; end; @ @ = set_random_seed_code: print_int(random_seed); @ To consume a random fraction, the program below will say `|next_random|' and then it will fetch |randoms[j_random]|. @d next_random==if j_random=0 then new_randoms else decr(j_random) @ = procedure new_randoms; var @!k:0..54; {index into |randoms|} @!x:integer; {accumulator} begin for k:=0 to 23 do begin x:=randoms[k]-randoms[k+31]; if x<0 then x:=x+mpfract_one; randoms[k]:=x; end; for k:=24 to 54 do begin x:=randoms[k]-randoms[k-24]; if x<0 then x:=x+mpfract_one; randoms[k]:=x; end; j_random:=54; end; @ To initialize the |randoms| table, we call the following routine. @ = procedure init_randoms; var @!j,@!jj,@!k:mpfract; {more or less random integers} @!i:0..54; {index into |randoms|} begin j:=abs(random_seed); while j>=mpfract_one do j:=halfp(j); k:=1; for i:=0 to 54 do begin jj:=k; k:=j-k; j:=jj; if k<0 then k:=k+mpfract_one; randoms[(i*21)mod 55]:=j; end; new_randoms; new_randoms; new_randoms; {``warm up'' the array} end; @ To produce a uniform random number in the range |0<=u =u>x| or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here. Note that the call of |mult_integers| will produce the values 0 and~|x| with about half the probability that it will produce any other particular values between 0 and~|x|, because it rounds its answers. @ = function unif_rand(@!x:scaled):scaled; var @!y:scaled; {trial value} begin next_random; y:=take_mpfract(abs(x),randoms[j_random]); if y=abs(x) then unif_rand:=0 else if x>0 then unif_rand:=y else unif_rand:=-y; end; @ This can be used by calling the following primitive. @ = primitive("uniformdeviate",convert,uniform_deviate_code);@/ @!@:uniformdeviate_}{\.{\\uniformdeviate} primitive@> @ @ = uniform_deviate_code: print_esc("uniformdeviate"); @ It takes one integer argument obviously that will be the argument to the function. @ = uniform_deviate_code: begin scan_int; cur_val:=unif_rand(cur_val); end; @ @ = uniform_deviate_code: print_int(cur_val); @ The following somewhat different subroutine tests rigorously if $ab$ is greater than, equal to, or less than~$cd$, given integers $(a,b,c,d)$. In most cases a quick decision is reached. The result is $+1$, 0, or~$-1$ in the three respective cases. @d return_sign(#)==begin ab_vs_cd:=#; return; end @ = function ab_vs_cd(@!a,b,c,d:integer):integer; label exit; var @!q,@!r:integer; {temporary registers} begin @ =0|, |b,d>0|@>; loop@+ begin q := a div d; r := c div b; if q<>r then if q>r then return_sign(1)@+else return_sign(-1); q := a mod d; r := c mod b; if r=0 then if q=0 then return_sign(0)@+else return_sign(1); if q=0 then return_sign(-1); a:=b; b:=q; c:=d; d:=r; end; {now |a>d>0| and |c>b>0|} exit:end; @ @ = if a<0 then begin negate(a); negate(b); end; if c<0 then begin negate(c); negate(d); end; if d<=0 then begin if b>=0 then if ((a=0)or(b=0))and((c=0)or(d=0)) then return_sign(0) else return_sign(1); if d=0 then if a=0 then return_sign(0)@+else return_sign(-1); q:=a; a:=c; c:=q; q:=-b; b:=-d; d:=q; end else if b<=0 then begin if b<0 then if a>0 then return_sign(-1); if c=0 then return_sign(0) else return_sign(-1); end @ Finally, a normal deviate with mean zero and unit standard deviation can readily be obtained with the ratio method (Algorithm 3.4.1R in {\sl The Art of Computer Programming\/}). @ = function norm_rand:scaled; var @!x,@!u,@!l:integer; {what the book would call $2^{16}X$, $2^{28}U$, and $-2^{24}\ln U$} begin repeat repeat next_random; x:=take_mpfract(112429,randoms[j_random]-mpfract_half); {$2^{16}\sqrt{8/e}\approx 112428.82793$} next_random; u:=randoms[j_random]; until abs(x) =0; norm_rand:=x; end; @ This can be used by calling the following primitive. @= primitive("normaldeviate",convert,normal_deviate_code);@/ @!@:normaldeviate_}{\.{\\normaldeviate} primitive@> @ @ = normal_deviate_code: print_esc("normaldeviate"); @ @ = normal_deviate_code: cur_val:=norm_rand; @ @ = normal_deviate_code: print_int(cur_val); @*1 DVI related primitives. These primitives are related to positions in the DVI output. The \TeX\ and DVI system coordinates relate to an origin that is at the upper left corner. The \TeX\ coordinates are computed relative to an origin that has $(0,0)$ coordinates. Coordinates grow then rightward and downward. This is the {\sl page} coordinates relative to what is typeset (what \TeX\ is dealing with). But this typesetting material has to be put on what we will call {\sl paper}. The material put into shape by \TeX\ is put on the paper. On this paper, where will be put the \TeX\ origin? It is considered to be $1in$ at the right and $1in$ down from the upper left corner of the paper (see m.590, alinea 2). @d DVI_std_x_offset=4736286 {1 inch in sp} @d DVI_std_y_offset=4736286 {1 inch in sp} @ But the paper size is not specified in the DVI file and is not being dealt with by \TeX. In order to have a common reference point, and since the \.{\\lastxpos} and \.{\\lastypos} primitives originated in pdf\TeX, these two primitives give positions, in scaled points, relative to the lower left corner of the paper. Hence the need, for these primitive, to define the paper size, with the (misnamed) \.{\\pagewidth} and \.{\\pageheight}. \.{\\pagewidth} and \.{\\pageheight} are dimension parameters, initialized to $0$ by the generic \TeX\ code. @ = primitive("pagewidth",assign_dimen,dimen_base+page_width_code);@/ @!@:pagewidth_}{\.{\\pagewidth} primitive@> primitive("pageheight",assign_dimen,dimen_base+page_height_code);@/ @!@:pageheight_}{\.{\\pageheight} primitive@> @ When instructed to, the \.{h} and \.{v} last values are transformed, in the coordinates system defined above and saved in the global variables |last_saved_xpos| and |last_saved_ypos|. They are initialized to $0$ and we do not make any verification that a call to the \.{\\savepos} primitive---to come---has been made before retrieving their values. @ = @!last_saved_xpos,last_saved_ypos: scaled; {last (x,y) DVI pos saved} @ @<\Prote\ initializations@>= last_saved_xpos:=0; last_saved_ypos:=0; @ @ = last_saved_xpos:=cur_h+DVI_std_x_offset; last_saved_ypos:=page_height-(cur_v+DVI_std_y_offset); @ @ = primitive("lastxpos",last_item,last_xpos_code);@/ @!@:lastxpos_}{\.{\\lastxpos} primitive@> primitive("lastypos",last_item,last_ypos_code);@/ @!@:lastypos_}{\.{\\lastypos} primitive@> @ @ = last_xpos_code: print_esc("lastxpos"); last_ypos_code: print_esc("lastypos"); @ @ = last_xpos_code: cur_val:=last_saved_xpos; last_ypos_code: cur_val:=last_saved_ypos; @ |last_saved_xpos| and |last_saved_ypos| are only defined when instructed to by the call the the \.{\\savepos} primitive. Since the real work has to be done at \.{shipout} time, it is a case to be treated like the \.{\\special} primitive, that is it belongs to the \.{extension} class. We will add something more in the handling of the primitive: it will insert a \.{whatsit} in the DVI file so that one, using the program |dvitype|, could retrieve more than one {\it hic}. So there is a counter incremented whenever the primitive is called. @ = last_save_pos_number:integer; {identifying the order of the call} @ @<\Prote\ initializations@>= last_save_pos_number:=0; {i.e. none} @ @ = primitive("savepos",extension,save_pos_code);@/ @!@:savepos_}{\.{\\savepos} primitive@> @ @ = save_pos_code: print_esc("savepos"); @ @ = save_pos_code:@ ; @ We need the basic two words node, since we don't pass any parameter and it is just an instruction to do something. So the \.{whatsit} node is just the call. @ = begin new_whatsit(save_pos_code,small_node_size); write_stream(tail):=null; write_tokens(tail):=null; end @ @ = save_pos_code: print_esc("savepos"); @ @ = save_pos_code: begin r:=get_node(small_node_size); words:=small_node_size; end; @ @ = save_pos_code: free_node(p,small_node_size); @ So, after these trivial initializations, what will we effectively do? When the following procedure will be called, we define |last_saved_xpos|, |last_saved_ypos|, increment |last_save_pos_number|, and a |warning| followed by three |key=value| space separated definitions as a \.{\\special}, the first being prefixed by the string |__PROTE_| (shall be considered a reserved prefix) and the string |SAVEPOS_|, equal to the index of the call, and the |XPOS| and |YPOS| definitions. This is obviously, from the previous description, a variation around |special_out|. @ = procedure save_pos_out(@!p:pointer); var old_setting:0..max_selector; {holds print |selector|} @!k:pool_pointer; {index into |str_pool|} begin synch_h; synch_v; incr(last_save_pos_number); @ @; old_setting:=selector; selector:=new_string; print("warning __PROTE_");print("SAVEPOS");print_char("="); print_int(last_save_pos_number);print_char(" "); print("XPOS");print("=");print_int(last_saved_xpos);print_char(" "); print("YPOS");print("=");print_int(last_saved_ypos); selector:=old_setting; str_room(1); {abort if probably overflowed and truncated} dvi_out(xxx1); dvi_out(cur_length); {it's less than 256} for k:=str_start[str_ptr] to pool_ptr-1 do dvi_out(so(str_pool[k])); pool_ptr:=str_start[str_ptr]; {forget the not commited tentative string} end; @ @ = save_pos_code: save_pos_out(p); @* \[54] System-dependent changes. @z