1)abbrev domain SOREXPV SortedExponentVector
2++ Domain for storing information about structure of polynomials
3++ as vectors of exponents
4SortedExponentVector == U32Vector
5
6
7)abbrev domain VECREC1 VectorModularReconstructor
8++ Description: This domain supports modular methods based on
9++  evaluation and rational reconstruction.  All computation
10++  are done on polynomials modulo machine sized prime p -- p must
11++  be chosen small enough to avoid overflow in intermediate
12++  calculations. Each evaluation is supposed to produce vector of
13++  values. Once enough evaluations are known rational reconstruction
14++  produces vector of rational functions or multivariate polynomials.
15VectorModularReconstructor() : Export == Implementation where
16    PA ==> U32Vector
17    RR ==> Record(numers : PrimitiveArray PA, denoms : PrimitiveArray PA)
18    RatRec ==> Record(numer : PA, denom : PA)
19    VI ==> Vector Integer
20    PPA ==> PrimitiveArray PA
21    PDR ==> Record(nvars : Integer, offsetdata : VI, _
22                   expdata : SortedExponentVector, _
23                   coeffdata : PA)
24    Export ==> with
25        empty : (Integer, Integer) -> %
26          ++ empty(n, p) initializes reconstructor with n slots
27          ++ working modulo p
28        add_slots : (List Integer, %) -> Void
29          ++ extend reconstructor adding zeros at specified
30          ++ positions.
31        chinese_update : (PA, Integer, %) -> Void
32          ++ chinese_update(v, pt, r) informs r that
33          ++ evaluation at pt gave vector of values v
34        rational_reconstruction : % -> Union(RR, "failed")
35          ++ reconstruct vector of rational functions based on
36          ++ information stored in reconstructor.
37        rational_reconstruction : (PA, PA, Integer, Integer) _
38                            -> Union(RatRec, "failed")
39          ++ reconstruct rational function
40        repack_polys : (Integer, VI, SortedExponentVector, PPA) -> PDR
41          ++ convert polynomials represented as parallel vector
42          ++ of exponents in k variables and vector of univariate
43          ++ polynomials to parallel vector of exponents in k+1
44          ++ variables and coefficients
45        remove_denoms : (VI, PPA, PPA, Integer) -> PPA
46          ++ remove_denoms(offsets, nums, denoms, p) removes common
47          ++ denominator from vectors of rational functions.  Several
48          ++ vectors of rational functions are packed into nums
49          ++ (storing numerators) and dens (storing denominators)
50          ++ Vector i starts at position offsets(i).  Computations
51          ++ are done modulo p.
52        reconstruct : (Integer, VI, VI, SortedExponentVector, _
53                       PPA, PPA, Integer) -> PDR
54          ++ reconstruct polynomials in n + 1 variables from
55          ++ result of rational reconstruction
56        reconstruct : (%, Integer, VI, VI, SortedExponentVector) _
57                       -> Union(PDR, "failed")
58          ++ reconstruct polynomials in n + 1 variables
59
60    Implementation ==> add
61        Rep := Record(prime : Integer, lpol : PA, curj : Integer, _
62                      npoints : Integer, npolys : Integer, _
63                      palloc : Integer, polys : PrimitiveArray PA, _
64                      next_rec : Integer, rec_step : Integer, _
65                      numers : PrimitiveArray PA, _
66                      denoms : PrimitiveArray PA)
67
68        modInverse ==> invmod
69
70        import from U32VectorPolynomialOperations
71
72        empty(npoly, np) ==
73            polyvec := new(npoly::NonNegativeInteger, _
74                           empty()$U32Vector)$PrimitiveArray(U32Vector)
75            for i in 0..(npoly - 1) repeat
76                polyvec(i) := new(5, 0)$U32Vector
77            state := [np, new(5, 0)$U32Vector, 0, 0, npoly, _
78                      5, polyvec, 3, 1, empty()$PrimitiveArray(PA), _
79                      empty()$PrimitiveArray(PA)]$Rep
80            setelt!(state.lpol, 0, 1)
81            state
82
83        add_slots(ndl : List Integer, statearg : %) : Void ==
84            state := statearg::Rep
85            polyvec := state.polys
86            m := state.palloc
87            n0 := #polyvec
88            n1 := #ndl
89            npoly := n0 + n1
90            nvec := new(npoly::NonNegativeInteger, _
91                           empty()$U32Vector)$PrimitiveArray(U32Vector)
92            li := first(ndl)
93            j : Integer := 0
94            for i in 0..(npoly - 1) repeat
95                i = li =>
96                    nvec(i) := new(m ::NonNegativeInteger, 0)$PA
97                    ndl := rest(ndl)
98                    li :=
99                        empty?(ndl) => npoly
100                        first(ndl)
101                nvec(i) := polyvec(j)
102                j := j + 1
103            if not(empty?(state.numers)) then
104                state.numers := new(npoly::NonNegativeInteger, _
105                                    empty()$PA)$PrimitiveArray(U32Vector)
106                state.denoms := new(npoly::NonNegativeInteger, _
107                                    empty()$PA)$PrimitiveArray(U32Vector)
108            state.polys := polyvec
109            state.npolys := npoly
110
111        double_poly_space(statearg : %) : Void ==
112            state := statearg::Rep
113            polyvec := state.polys
114            m := state.palloc
115            n := 2*m
116            for i in 0..(state.npolys - 1) repeat
117                np := new(n::NonNegativeInteger, 0)$U32Vector
118                op := polyvec(i)
119                copy_first(np, op, m)
120                polyvec(i) := np
121            state.palloc := n
122
123        chinese_update(vec, pt, statearg) ==
124            state := statearg::Rep
125            mtvec := state.lpol
126            npt := state.npoints
127            npt1 := npt + 1
128            p := state.prime
129            mtval := eval_at(mtvec, npt, pt, p)
130            mtval = 0 => error "Duplicate point in update"
131            mtcor := modInverse(mtval, p)
132            state.npoints := npt1
133            if npt1 > state.palloc then
134                double_poly_space(statearg)
135            polyvec := state.polys
136            nn := npt - 1
137            for i in 0..(state.npolys - 1) repeat
138               pol := polyvec(i)
139               cor := vec(i) - eval_at(pol, nn, pt, p)
140               cor :=
141                   cor < 0 => cor + p
142                   cor
143               cor := positiveRemainder(cor*mtcor, p)
144               vector_add_mul(pol, mtvec, 0, npt, cor, p)
145            if #mtvec < npt1 + 1 then
146                nmt := new(2*(npt1::NonNegativeInteger), 0)$U32Vector
147                copy_first(nmt, mtvec, npt1)
148                mtvec := nmt
149            mul_by_binomial(mtvec, npt1 + 1, p - pt, p)
150            state.lpol := mtvec
151
152        rational_reconstruction(x : PA, y : PA, i : Integer, p : Integer) : _
153                            Union(RatRec, "failed") ==
154            -- invariant: r0 = t0*y + s0*x, r1 = t1*y + s1*x
155            -- we do not need t0 and t1, so we do not compute them
156            j := degree(y)
157            r0 := new(qcoerce(j+1)@NonNegativeInteger, 0)$PA
158            copy_first(r0, y, j + 1)
159            dr0 := j
160            -- s0 is 0
161            s0 := new(qcoerce(j+1)@NonNegativeInteger, 0)$PA
162            ds0 := 0$Integer
163            r1 := new(qcoerce(j+1)@NonNegativeInteger, 0)$PA
164            dr1 := degree(x)
165            copy_first(r1, x, dr1 + 1)
166            -- s1 is 1
167            s1 := new(qcoerce(j+1)@NonNegativeInteger, 0)$PA
168            s1(0) := 1
169            ds1 := 0$Integer
170            while dr1 > i repeat
171                while dr0 >= dr1 repeat
172                    delta := dr0 - dr1
173                    c1 := p - r0(dr0)
174                    c0 := r1(dr1)
175                    r0(dr0) := 0
176                    dr0 := dr0 - 1
177                    vector_combination(r0, c0, r1, c1, dr0, delta, p)
178                    while r0(dr0) = 0 repeat
179                        dr0 := dr0 - 1
180                        if dr0 < 0 then break
181                    ds0a := ds1 + delta
182                    ds0a :=
183                        ds0 > ds0a => ds0
184                        ds0a
185                    vector_combination(s0, c0, s1, c1, ds0a, delta, p)
186                    ds0 :=
187                        ds0 > ds0a => ds0
188                        ds0 < ds0a => ds0a
189                        while s0(ds0a) = 0 repeat ds0a := ds0a - 1
190                        ds0a
191                tmpp := r0
192                tmp := dr0
193                r0 := r1
194                dr0 := dr1
195                r1 := tmpp
196                dr1 := tmp
197                tmpp := s0
198                tmp := ds0
199                s0 := s1
200                ds0 := ds1
201                s1 := tmpp
202                ds1 := tmp
203            ds1 > j - i - 2 => "failed"
204            degree(gcd(s1, y, p)) ~= 0 => "failed"
205            c := s1(ds1)
206            c := modInverse(c, p)
207            mul_by_scalar(r1, dr1, c, p)
208            mul_by_scalar(s1, ds1, c, p)
209            [r1, s1]
210
211        rational_reconstruction(statearg : %) : Union(RR, "failed") ==
212            state := statearg::Rep
213            modulus := state.lpol
214            polyvec := state.polys
215            p := state.prime
216            j0 := state.curj
217            -- Compiler bug
218            -- j0 <= 3 => "failed"
219            m := state.npoints
220            m <= state.next_rec => return "failed"
221            state.next_rec := state.next_rec + state.rec_step
222            if m > 30 then
223                state.rec_step :=
224                    state.rec_step +
225                           m > 200 => 4
226                           1
227            bound := m quo 2
228            ok := true
229            pp := rational_reconstruction(polyvec(j0), modulus, bound, p)
230            pp case "failed" => "failed"
231            n := state.npolys
232            if empty?(state.numers) then
233               state.numers := new(n::NonNegativeInteger, _
234                                   empty()$PA)$PrimitiveArray(PA)
235               state.denoms := new(n::NonNegativeInteger, _
236                                   empty()$PA)$PrimitiveArray(PA)
237            nums := state.numers
238            dens := state.denoms
239            ppr := pp :: RatRec
240            nums(j0) := ppr.numer
241            dens(j0) := ppr.denom
242            cden := ppr.denom
243            j := j0
244            repeat
245               j := j + 1
246               if j >= n then j := j - n
247               j = j0 => return [nums, dens]
248               r1 := polyvec(j)
249               r1 := mul(r1, cden, p)
250               remainder!(r1, modulus, p)
251               (deg_r1 := degree(r1)) < bound =>
252                   rp := new(qcoerce(deg_r1 + 1)@NonNegativeInteger, 0)$PA
253                   copy_first(rp, r1, deg_r1 + 1)
254                   nums(j) := rp
255                   dens(j) := cden
256               pp := rational_reconstruction(r1, modulus, bound, p)
257               pp case "failed" =>
258                   state.curj := j
259                   return "failed"
260               ppr := pp :: RatRec
261               cden := mul(cden, ppr.denom, p)
262               degree(cden) > bound =>
263                   state.curj := j
264                   return "failed"
265               nums(j) := ppr.numer
266               dens(j) := cden
267
268
269        repack_polys(var_cnt : Integer, poly_offsets : VI, _
270                      exps : SortedExponentVector, _
271                      coeffs : PrimitiveArray U32Vector) : PDR ==
272            m : Integer := 0
273            n := #coeffs
274            -- count nonzero coefficients
275            for i in 0..(n - 1) repeat
276                ci := coeffs(i)
277                k := #ci
278                for j in 0..(k - 1) repeat
279                    if ci(j) ~= 0 then m := m + 1
280            nnvars := var_cnt + 1
281            nexps := new(qcoerce(m*nnvars)@NonNegativeInteger, _
282                         0)$SortedExponentVector
283            ncoeffs := new(qcoerce(m)@NonNegativeInteger, 0)$U32Vector
284            pi_cnt := #poly_offsets
285            npo := new(pi_cnt, 0)$VI
286            pi : Integer := 1
287            opi := poly_offsets(pi)
288            nm : SingleInteger := 0
289            oei : SingleInteger := 0
290            nei : SingleInteger := 0
291            for i in 0..(n - 1) repeat
292                while opi = i repeat
293                    npo(pi) := nm
294                    pi := pi + 1
295                    opi :=
296                        pi <= pi_cnt =>
297                            opi := poly_offsets(pi)
298                        -1
299                ci := coeffs(i)
300                k := #ci
301                for j in 0..(k - 1) repeat
302                    cij := ci(j)
303                    if cij ~= 0 then
304                        ncoeffs(nm) := cij
305                        nm := nm + 1
306                        oei0 := oei
307                        for i1 in 1..var_cnt repeat
308                            nexps(nei) := exps(oei0)
309                            nei := nei + 1
310                            oei0 := oei0 + 1
311                        nexps(nei) := j
312                        nei := nei + 1
313                oei := oei + qconvert(var_cnt)@SingleInteger
314            for i in pi..pi_cnt repeat npo(i) := nm
315            [nnvars, npo, nexps, ncoeffs]
316
317        remove_denoms(block_offsets : VI, nums : PPA, _
318                       dens : PPA, p : Integer) : PPA ==
319            nb := #block_offsets
320            np := #nums
321            res := new(np, empty()$PA)$PPA
322            tmpp1 : PA
323            tmpp2 : PA
324            for ib in 1..nb repeat
325                li := block_offsets(ib)
326                hi : Integer :=
327                    ib = nb => np
328                    block_offsets(ib + 1)
329                hi := hi - 1
330                cden := lcm(dens, li, hi, p)
331                dcden := degree(cden)
332                tmpp1 := new(qcoerce(dcden + 1)@NonNegativeInteger, 0)$PA
333                tmpp2 := new(qcoerce(dcden + 1)@NonNegativeInteger, 0)$PA
334                for i in li..hi repeat
335                    copy_first(tmpp1, cden, dcden + 1)
336                    for j in 0..dcden repeat tmpp2(j) := 0
337                    divide!(tmpp1, dens(i), tmpp2, p)
338                    dt := degree(tmpp2)
339                    res(i) :=
340                        dt > 0 => mul(tmpp2, nums(i), p)
341                        dt = 0 and tmpp2(0) ~= 1 =>
342                            error "remove_denoms expect quotient to be 1"
343                        nums(i)
344                cfactor := gcd(res, li, hi, p)
345                dcf := degree(cfactor)
346                dtmp := dcden
347                if dcf >= 1 then
348                    for i in li..hi repeat
349                        resi := res(i)
350                        dresi := degree(resi)
351                        if dresi > dtmp then
352                            dtmp := dresi
353                            tmpp1 := new(qcoerce(dtmp + 1)@ _
354                                          NonNegativeInteger, 0)$PA
355                            tmpp2 := new(qcoerce(dtmp + 1)@ _
356                                          NonNegativeInteger, 0)$PA
357                        else
358                            for j in 0..dtmp repeat tmpp2(j) := 0
359                        copy_first(tmpp1, resi, dresi + 1)
360                        divide!(tmpp1, cfactor, tmpp2, p)
361                        dt := degree(tmpp2)
362                        res(i) := new(qcoerce(dt + 1)@NonNegativeInteger, 0)$PA
363                        copy_first(res(i), tmpp2, dt + 1)
364            res
365
366        reconstruct(var_cnt : Integer, block_offsets : VI,
367                    poly_offsets : VI, exps : SortedExponentVector, _
368                    nums : PPA, dens : PPA, p : Integer) : PDR ==
369            ppa := remove_denoms(block_offsets, nums, dens, p)
370            repack_polys(var_cnt, poly_offsets, exps, ppa)
371
372        reconstruct(statearg : %, var_cnt : Integer, _
373                    block_offsets : VI, poly_offsets : VI, _
374                    exps : SortedExponentVector) _
375                     : Union(PDR, "failed") ==
376            pp := rational_reconstruction(statearg)
377            pp case "failed" => return "failed"
378            state := statearg::Rep
379            ppr := pp::RR
380            reconstruct(var_cnt, block_offsets, poly_offsets, _
381                        exps, ppr.numers, ppr.denoms, state.prime)
382
383
384)abbrev domain VECREC2 VectorIntegerReconstructor
385++ Description: This domain supports modular methods based on
386++  evaluation and rational reconstruction.  Each evaluation
387++  is done modulo machine sized prime p.  Both Chinese
388++  remaindering and (linear) Hensel lift are supported.
389++  Once enough evaluations are known rational reconstruction
390++  produces vector of rational numbers or integers.
391VectorIntegerReconstructor() : Export == Implementation where
392   RatRec ==> Record(num : Integer, den : Integer)
393   PAI ==> PrimitiveArray Integer
394   RR ==> Record(numers : PAI, denoms : PAI)
395   VI ==> Vector Integer
396   Export ==> with
397       empty : (Integer) -> %
398         ++ empty(n) produces reconstructor with n slots
399       chinese_update : (U32Vector, Integer, %) -> Void
400         ++ chinese_update(v, p, r) informs r about evaluation at p
401       hensel_update : (U32Vector, Integer, %) -> Void
402         ++ one step of Hensel lifting
403       rational_reconstruction : (Integer, Integer, Integer, Integer) -> _
404           Union(RatRec, "failed")
405       rational_reconstruction : % -> Union(RR, "failed")
406       remove_denoms : (VI, PAI, PAI) -> PAI
407         ++ remove common denominators in blocks
408       reconstruct : (%, VI) -> Union(PAI, "failed")
409         ++ combines rational reconstruction with removal of
410         ++ common denominators in blocks.
411
412   Implementation ==> add
413       Rep := Record(cmod : Integer, curj : Integer,  _
414                     nmods : Integer, nints : Integer, _
415                     ints : PrimitiveArray Integer,
416                     bints : PrimitiveArray Integer,
417                     bcmod : Integer, bnmods : Integer,
418                     nrecs : Integer,
419                     numers : PrimitiveArray Integer,
420                     denoms : PrimitiveArray Integer)
421
422
423       modInverse(c : Integer, p : Integer) : Integer ==
424          (extendedEuclidean(c, p, 1)::Record(coef1 : Integer, _
425               coef2 : Integer)).coef1
426
427       empty(nint) ==
428           intvec := new(nint::NonNegativeInteger, _
429                          0)$PrimitiveArray(Integer)
430           [1, 0, 0, nint, intvec, empty()$PrimitiveArray(Integer), _
431            1, 0, 0, empty()$PrimitiveArray(Integer), _
432            empty()$PrimitiveArray(Integer)]$Rep
433
434       chinese_update(vec, p, statearg) ==
435           state := statearg::Rep
436           mp := state.cmod
437           mpval := positiveRemainder(mp, p)
438           mpval = 0 => error "Duplicate modulus in update"
439           mpcor := modInverse(mpval, p)
440           mpfact := mpcor*mp
441           nmp := mp*p
442           mpfact := positiveRemainder(mpfact, nmp)
443           state.nmods := state.nmods + 1
444           intvec := state.ints
445           for i in 0..(state.nints - 1) repeat
446               ii := intvec(i)
447               cor := positiveRemainder(vec(i) - ii, p)
448               intvec(i) := positiveRemainder(ii + mpfact*cor, nmp)
449           if state.nmods >= 200 and
450             positiveRemainder(state.nmods, 100) = 0 then
451               state.bnmods := state.nmods
452               empty?(state.bints) =>
453                   state.bints := new(state.nints::NonNegativeInteger, _
454                          0)$PrimitiveArray(Integer)
455                   bintvec := state.bints
456                   for i in 0..(state.nints - 1) repeat
457                       bintvec(i) := intvec(i)
458                       intvec(i) := 0
459                   state.bcmod := nmp
460                   nmp := 1
461               bintvec := state.bints
462               bmp := state.bcmod
463               mpval := positiveRemainder(bmp, nmp)
464               mpcor := modInverse(mpval, nmp)
465               mpfact := mpcor*bmp
466               nbmp := bmp*nmp
467               mpfact := positiveRemainder(mpfact, nbmp)
468               for i in 0..(state.nints - 1) repeat
469                   ii := bintvec(i)
470                   cor := positiveRemainder(intvec(i) - ii, nmp)
471                   bintvec(i) := positiveRemainder(ii + mpfact*cor, nbmp)
472                   intvec(i) := 0
473               state.bcmod := nbmp
474               nmp := 1
475           state.cmod := nmp
476
477       hensel_update(vec, p, statearg) ==
478           state := statearg::Rep
479           mp := state.cmod
480           intvec := state.ints
481           for i in 0..(state.nints - 1) repeat
482               intvec(i) := intvec(i) + vec(i)*mp
483           state.cmod := p*mp
484
485       rational_reconstruction(x : Integer, y : Integer, i : Integer, _
486                 j : Integer) : Union(RatRec, "failed") ==
487           r0 : Integer := y
488           s0 : Integer := 0
489           r1 : Integer := positiveRemainder(x, y)
490           s1 : Integer := 1
491           while r1 > i repeat
492               qr := divide(r0, r1)
493               r0 := r1
494               r1 := qr.remainder
495               tmp := s0 - qr.quotient*s1
496               s0 := s1
497               s1 := tmp
498           if s1 < 0 then
499               s1 := -s1
500               r1 := -r1
501           s1 > j => "failed"
502           gcd(s1, y) ~= 1 => "failed"
503           [r1, s1]
504
505
506       rational_reconstruction2(statearg : %, block_offsets : VI) _
507               : Union(RR, "failed") ==
508           state := statearg::Rep
509           modulus := state.cmod
510           intvec := state.ints
511           if state.nmods >= 200 then
512               if state.nmods - state.nrecs < 150 then
513                   state.cmod ~= 1 =>
514                       positiveRemainder(state.nmods, 100) = 0 =>
515                           error "impossible"
516                       return "failed"
517               state.nrecs := state.bnmods
518               modulus := state.bcmod
519               intvec := state.bints
520           j0 := state.curj
521           bound := approxSqrt(modulus)$IntegerRoots(Integer) quo 10
522           bound2 := modulus - bound
523           ok := true
524           pp := rational_reconstruction(intvec(j0), modulus, bound, bound)
525           pp case "failed" => "failed"
526           n := state.nints
527           mm := #block_offsets
528           if empty?(state.numers) then
529               state.numers := new(n::NonNegativeInteger, _
530                                   0)$PrimitiveArray(Integer)
531               state.denoms := new(n::NonNegativeInteger, _
532                                   0)$PrimitiveArray(Integer)
533           nums := state.numers
534           dens := state.denoms
535           ppr := pp :: RatRec
536           nums(j0) := ppr.num
537           dens(j0) := ppr.den
538           cden := ppr.den
539           co : Integer := 0
540           cb : Integer := 1
541           for jj in 1..mm repeat
542              ctmp := block_offsets(jj)
543              if ctmp > j0 then
544                  cb := jj
545                  co := ctmp
546                  break
547           j := j0
548           repeat
549               j := j + 1
550               if j >= n then j := j - n
551               if j = co and mm > 1 then
552                   cden := 1
553                   cb :=
554                       cb = mm => 1
555                       cb + 1
556                   co := block_offsets(cb)
557               j = j0 => return [nums, dens]
558               r1 := positiveRemainder(cden*intvec(j), modulus)
559               r1 < bound =>
560                   nums(j) := r1
561                   dens(j) := cden
562               r1 > bound2 =>
563                   nums(j) := r1 - modulus
564                   dens(j) := cden
565               pp := rational_reconstruction(r1, modulus, bound, bound)
566               pp case "failed" =>
567                   state.curj := j
568                   return "failed"
569               ppr := pp :: RatRec
570               cden := cden*ppr.den
571               cden > bound =>
572                   state.curj := j
573                   return "failed"
574               nums(j) := ppr.num
575               dens(j) := cden
576
577       rational_reconstruction(statearg : %) : Union(RR, "failed") ==
578           rational_reconstruction2(statearg, new(1, 0)$VI)
579
580       lcm(nums : PAI, lo : Integer, hi : Integer) : Integer ==
581           res := nums(lo)
582           for i in (lo + 1)..hi repeat
583               res := lcm(res, nums(i))
584           res
585
586       gcd(nums : PAI, lo : Integer, hi : Integer) : Integer ==
587           res := nums(lo)
588           for i in (lo + 1)..hi repeat
589               res := gcd(res, nums(i))
590           res
591
592       remove_denoms(block_offsets : VI, nums : PAI, _
593                      dens : PAI) : PAI ==
594            nb := #block_offsets
595            np := #nums
596            res := new(np, 0)$PAI
597            for ib in 1..nb repeat
598                li := block_offsets(ib)
599                hi : Integer :=
600                    ib = nb => np
601                    block_offsets(ib + 1)
602                hi := hi - 1
603                cden := lcm(dens, li, hi)
604                for i in li..hi repeat
605                    tmpp2  := (cden exquo dens(i))::Integer
606                    res(i) := nums(i)*tmpp2
607                cfactor := gcd(res, li, hi)
608                if cfactor ~= 1 and cfactor ~= 0 then
609                    for i in li..hi repeat
610                        res(i) := (res(i) exquo cfactor)::Integer
611            res
612
613       reconstruct(statearg : %, block_offsets : VI) : _
614               Union(PAI, "failed") ==
615           pp := rational_reconstruction2(statearg, block_offsets)
616           pp case "failed" => return "failed"
617           ppr := pp::RR
618           remove_denoms(block_offsets, ppr.numers, ppr.denoms)
619