Residual-Contraction log1p Kernel

Integrate this at the scalar-elementwise primitive level, initially replacing log1p inside softplus and stable log-sum-exp rather than modifying transformer weights. Inputs are floating-point tensors u; choose a safe range such as [-0.5, 1], and use ordinary range reduction for values outside that interval. For each element, run the following fixed-depth procedure: (1) set a[0]=u; (2) for k=0,...,K-1, compute a[k+1]=sqrt(1+a[k])-1; (3) initialize r[K]=1/3; (4) for k=K-1,...,0, compute r[k]=(2*r[k+1]+2+0.5*a[k+1])/(2+a[k+1])^3; (5) return u-0.5*u*u+u^3*r[0]. In code, fuse the square-root, residual-backward loop, polynomial reconstruction, and surrounding softplus or log-sum-exp operations into one CUDA or Triton kernel to avoid intermediate memory traffic. The formulas determine the recurrence and base residual; empirically estimate the smallest K meeting the target error by comparing against high-precision log1p over a dense grid and random samples, and add a fallback for inputs outside the safe interval or for subnormal and exceptional values. First benchmark a fused softplus and log-sum-exp kernel on an RTX-class GPU using tensors of sizes 1M, 16M, and transformer vocabulary logits, against PyTorch and CUDA implementations. Measure maximum and percentile absolute error, gradient error, kernel latency, achieved bandwidth, and end-to-end training throughput. A successful result is at least 1.2x primitive or fused-kernel speedup at comparable FP16/BF16 accuracy, with no overflow or visible loss increase in a small Transformer language-model run. Test stability especially for logits near zero, where cancellation affects direct polynomial approximations.