This file is indexed.

/usr/share/axiom-20170501/src/algebra/RADFF.spad is in axiom-source 20170501-3.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
)abbrev domain RADFF RadicalFunctionField
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 27 July 1993
++ Description:
++ Function field defined by y**n = f(x);

RadicalFunctionField(F, UP, UPUP, radicnd, n) : SIG == CODE where
  F : UniqueFactorizationDomain
  UP : UnivariatePolynomialCategory F
  UPUP : UnivariatePolynomialCategory Fraction UP
  radicnd : Fraction UP
  n : NonNegativeInteger

  N   ==> NonNegativeInteger
  Z   ==> Integer
  RF  ==> Fraction UP
  QF  ==> Fraction UPUP
  UP2 ==> SparseUnivariatePolynomial UP
  REC ==> Record(factor:UP, exponent:Z)
  MOD ==> monomial(1, n)$UPUP - radicnd::UPUP
  INIT ==> if (deref brandNew?) then startUp false

  SIG ==> FunctionFieldCategory(F, UP, UPUP)

  CODE ==> SimpleAlgebraicExtension(RF, UPUP, MOD) add

    import ChangeOfVariable(F, UP, UPUP)
    import InnerCommonDenominator(UP, RF, Vector UP, Vector RF)
    import UnivariatePolynomialCategoryFunctions2(RF, UPUP, UP, UP2)

    diag        : Vector RF -> Vector $
    startUp     : Boolean -> Void
    fullVector  : (Factored UP, N) -> PrimitiveArray UP
    iBasis      : (UP, N) -> Vector UP
    inftyBasis  : (RF, N) -> Vector RF
    basisvec    : () -> Vector RF
    char0StartUp: () -> Void
    charPStartUp: () -> Void
    getInfBasis : () -> Void
    radcand     : () -> UP
    charPintbas : (UPUP, RF, Vector RF, Vector RF) -> Void

    brandNew?:Reference(Boolean) := ref true
    discPoly:Reference(RF) := ref(0$RF)
    newrad:Reference(UP) := ref(0$UP)
    n1 := (n - 1)::N
    modulus := MOD
    ibasis:Vector(RF)     := new(n, 0)
    invibasis:Vector(RF)  := new(n, 0)
    infbasis:Vector(RF)   := new(n, 0)
    invinfbasis:Vector(RF):= new(n, 0)
    mini := minIndex ibasis

    discriminant()                   == (INIT; discPoly())

    radcand()                        == (INIT; newrad())

    integralBasis()                  == (INIT; diag ibasis)

    integralBasisAtInfinity()        == (INIT; diag infbasis)

    basisvec()                       == (INIT; ibasis)

    integralMatrix()                 == diagonalMatrix basisvec()

    integralMatrixAtInfinity()       == (INIT; diagonalMatrix infbasis)

    inverseIntegralMatrix()          == (INIT; diagonalMatrix invibasis)

    inverseIntegralMatrixAtInfinity()==(INIT;diagonalMatrix invinfbasis)

    definingPolynomial()             == modulus

    ramified?(point:F)               == zero?(radcand() point)

    branchPointAtInfinity?()  == (degree(radcand()) exquo n) case "failed"

    elliptic()     == (n = 2 and degree(radcand()) = 3 => radcand(); "failed")

    hyperelliptic() == (n=2 and odd? degree(radcand()) => radcand(); "failed")

    diag v == [reduce monomial(qelt(v,i+mini), i) for i in 0..n1]

    integralRepresents(v, d) ==
      ib := basisvec()
      represents
        [qelt(ib, i) * (qelt(v, i) /$RF d) for i in mini .. maxIndex ib]

    integralCoordinates f ==
      v  := coordinates f
      ib := basisvec()
      splitDenominator
        [qelt(v,i) / qelt(ib,i) for i in mini .. maxIndex ib]$Vector(RF)

    integralDerivationMatrix d ==
      dlogp := differentiate(radicnd, d) / (n * radicnd)
      v := basisvec()
      cd := splitDenominator(
                [(i - mini) * dlogp + differentiate(qelt(v, i), d) / qelt(v, i)
                                         for i in mini..maxIndex v]$Vector(RF))
      [diagonalMatrix(cd.num), cd.den]

    -- return (d0,...,d(n-1)) s.t. (1/d0, y/d1,...,y**(n-1)/d(n-1))
    -- is an integral basis for the curve y**d = p
    -- requires that p has no factor of multiplicity >= d
    iBasis(p, d) ==
      pl := fullVector(squareFree p, d)
      d1 := (d - 1)::N
      [*/[pl.j ** ((i * j) quo d) for j in 0..d1] for i in 0..d1]

    -- returns a vector [a0,a1,...,a_{m-1}] of length m such that
    -- p = a0^0 a1^1 ... a_{m-1}^{m-1}
    fullVector(p, m) ==
      ans:PrimitiveArray(UP) := new(m, 0)
      ans.0 := unit p
      l := factors p
      for i in 1..maxIndex ans repeat
        ans.i :=
          (u := find(s+->s.exponent = i, l)) case "failed" => 1
          (u::REC).factor
      ans

    -- return (f0,...,f(n-1)) s.t. (f0, y f1,..., y**(n-1) f(n-1))
    -- is a local integral basis at infinity for the curve y**d = p
    inftyBasis(p, m) ==
      rt := rootPoly(p(x := inv(monomial(1, 1)$UP :: RF)), m)
      m ^= rt.exponent =>
        error "Curve not irreducible after change of variable 0 -> infinity"
      a    := (rt.coef) x
      b:RF := 1
      v    := iBasis(rt.radicand, m)
      w:Vector(RF) := new(m, 0)
      for i in mini..maxIndex v repeat
        qsetelt_!(w, i, b / ((qelt(v, i)::RF) x))
        b := b * a
      w

    charPintbas(p, c, v, w) ==
      degree(p) ^= n => error "charPintbas: should not happen"
      q:UP2 := map(s+->retract(s)@UP, p)
      ib := integralBasis()$FunctionFieldIntegralBasis(UP, UP2,
                                          SimpleAlgebraicExtension(UP, UP2, q))
      not diagonal?(ib.basis)=> 
         error "charPintbas: integral basis not diagonal"
      a:RF := 1
      for i in minRowIndex(ib.basis) .. maxRowIndex(ib.basis)
        for j in minColIndex(ib.basis) .. maxColIndex(ib.basis)
          for k in mini .. maxIndex v repeat
            qsetelt_!(v, k, (qelt(ib.basis, i, j) / ib.basisDen) * a)
            qsetelt_!(w, k, qelt(ib.basisInv, i, j) * inv a)
            a := a * c
      void

    charPStartUp() ==
      r      := mkIntegral modulus
      charPintbas(r.poly, r.coef, ibasis, invibasis)
      x      := inv(monomial(1, 1)$UP :: RF)
      invmod := monomial(1, n)$UPUP - (radicnd x)::UPUP
      r      := mkIntegral invmod
      charPintbas(r.poly, (r.coef) x, infbasis, invinfbasis)

    startUp b ==
      brandNew?() := b
      if zero?(p := characteristic()$F) or p > n then char0StartUp()
                                                 else charPStartUp()
      dsc:RF := ((-1)$Z ** ((n *$N n1) quo 2::N) * (n::Z)**n)$Z *
               radicnd ** n1 *
                  */[qelt(ibasis, i) ** 2 for i in mini..maxIndex ibasis]
      discPoly() := primitivePart(numer dsc) / denom(dsc)
      void

    char0StartUp() ==
      rp          := rootPoly(radicnd, n)
      rp.exponent ^= n => 
         error "RadicalFunctionField: curve is not irreducible"
      newrad()    := rp.radicand
      ib          := iBasis(newrad(), n)
      infb        := inftyBasis(radicnd, n)
      invden:RF   := 1
      for i in mini..maxIndex ib repeat
        qsetelt_!(invibasis, i, a := qelt(ib, i) * invden)
        qsetelt_!(ibasis, i, inv a)
        invden := invden / rp.coef        -- always equals 1/rp.coef**(i-mini)
        qsetelt_!(infbasis, i, a := qelt(infb, i))
        qsetelt_!(invinfbasis, i, inv a)
      void

    ramified?(p:UP) ==
      (r := retractIfCan(p)@Union(F, "failed")) case F =>
        singular?(r::F)
      (radcand() exquo p) case UP

    singular?(p:UP) ==
      (r := retractIfCan(p)@Union(F, "failed")) case F =>
        singular?(r::F)
      (radcand() exquo(p**2)) case UP

    branchPoint?(p:UP) ==
      (r := retractIfCan(p)@Union(F, "failed")) case F =>
        branchPoint?(r::F)
      ((q := (radcand() exquo p)) case UP) and
        ((q::UP exquo p) case "failed")

    singular?(point:F) ==
      zero?(radcand()  point) and
        zero?(((radcand() exquo (monomial(1,1)$UP-point::UP))::UP) point)

    branchPoint?(point:F) ==
      zero?(radcand()  point) and not
        zero?(((radcand() exquo (monomial(1,1)$UP-point::UP))::UP) point)